Regresión Lineal Múltiple
Nazly R. Hincapie Monsalve
2024-08-20
Paquetes
library(dplyr)##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(readr) #Para leer archivos .csv
library(ggplot2)
library(ggcorrplot)
library(MASS) # Criterio de Información de Akaike (AIC)##
## Adjuntando el paquete: 'MASS'## The following object is masked from 'package:dplyr':
##
## selectlibrary(nortest) # prueba de Anderson-Darling
library(lmtest) # prueba Breusch-Pagan## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericDataset
# https://www.kaggle.com/datasets/royjafari/customer-churn/data
churn_data = read_csv('Customer Churn.csv', show_col_types = FALSE)
head(churn_data)## # A tibble: 6 × 16
## `Call Failure` Complains `Subscription Length` `Charge Amount`
## <dbl> <dbl> <dbl> <dbl>
## 1 8 0 38 0
## 2 0 0 39 0
## 3 10 0 37 0
## 4 10 0 38 0
## 5 3 0 38 0
## 6 11 0 38 1
## # ℹ 12 more variables: `Seconds of Use` <dbl>, `Frequency of use` <dbl>,
## # `Frequency of SMS` <dbl>, `Distinct Called Numbers` <dbl>,
## # `Age Group` <dbl>, `Tariff Plan` <dbl>, Status <dbl>, Age <dbl>,
## # `Customer Value` <dbl>, FN <dbl>, FP <dbl>, Churn <dbl>#dim(churn_data)Descripción
Call Failure: number of call failures Complains: binary (0: No complaint, 1: complaint) Subscription Length: total months of subscription Charge Amount: ordinal attribute (0: lowest amount, 9: highest amount Seconds of Use: total seconds of calls Frequency of use: total number of calls Frequency of SMS: total number of text messages Distinct Called Numbers: total number of distinct phone calls Age Group: ordinal attribute (1: younger age, 5: older age) Tariff Plan: binary (1: Pay as you go, 2: contractual) Status: binary (1: active, 2: non-active) Age: age of customer Customer Value: the calculated value of customer Churn: class label (1: churn, 0: non-churn)
Utilizaremos el conjunto de datos de cancelación de clientes de Kaggle para calcular el valor del cliente (customer value).
¿Qué entendemos por valor del cliente? Básicamente, determina lo que vale un producto o servicio para un cliente, y podemos calcularlo de la siguiente manera:
Customer Value = Benefit-Cost.
Donde Beneficio y Coste son, respectivamente, el beneficio y el coste de un producto o un servicio.
Este valor es mayor si la empresa puede ofrecer a los consumidores mayores beneficios y menores costes o, idealmente, una combinación de ambos.
Este análisis podrÃa ayudar a la empresa a identificar la oportunidad más prometedora o la siguiente mejor acción en función del valor de un cliente determinado.
# Change the column names
names(churn_data) = gsub(" ", "_", names(churn_data))
names(churn_data) = gsub(" ", "_", names(churn_data))
head(churn_data)## # A tibble: 6 × 16
## Call_Failure Complains Subscription_Length Charge_Amount Seconds_of_Use
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 8 0 38 0 4370
## 2 0 0 39 0 318
## 3 10 0 37 0 2453
## 4 10 0 38 0 4198
## 5 3 0 38 0 2393
## 6 11 0 38 1 3775
## # ℹ 11 more variables: Frequency_of_use <dbl>, Frequency_of_SMS <dbl>,
## # Distinct_Called_Numbers <dbl>, Age_Group <dbl>, Tariff_Plan <dbl>,
## # Status <dbl>, Age <dbl>, Customer_Value <dbl>, FN <dbl>, FP <dbl>,
## # Churn <dbl>Modelo de regresión múltiple
cust_value_model = lm(formula = Customer_Value ~ Call_Failure + Complains + Subscription_Length + Charge_Amount + Seconds_of_Use + Frequency_of_use + Frequency_of_SMS + Distinct_Called_Numbers + Age_Group + Tariff_Plan + Status + Age,data = churn_data)Distribución de los residuales
# Obtener los residuales del modelo
model_residuals <- cust_value_model$residuals
# Crear un data frame para los residuales
residuals_df <- data.frame(Residuals = model_residuals)
# Crear el gráfico de histograma con lÃnea de densidad
ggplot(residuals_df, aes(x = Residuals)) +
geom_histogram(aes(y = ..density..), bins = 20, color = "black", fill = "maroon") +
geom_density(color = "lightblue", size = 1) +
ggtitle("Distribución de los Residuales con LÃnea de Densidad") +
xlab("Residuales") +
ylab("Densidad") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
El histograma parece sesgado a la izquierda, por lo que no podemos concluir la normalidad con suficiente confianza. En lugar del histograma, observemos los residuos a lo largo del gráfico Q-Q normal. Si hay normalidad, los valores deberÃan seguir una lÃnea recta.
# Crear el gráfico Q-Q con la lÃnea Q-Q
ggplot(residuals_df, aes(sample = Residuals)) +
stat_qq() +
stat_qq_line(color = "red") +
ggtitle("Q-Q Plot de los Residuales") +
xlab("Cuantiles Teóricos") +
ylab("Cuantiles de los Residuales") +
theme_minimal()
A partir del gráfico, podemos observar que algunas partes de los residuos se sitúan en una lÃnea recta. Entonces podemos suponer que los residuos del modelo no siguen una distribución normal.
Matriz de correlación
# Remove the Customer Value column
reduced_data <- subset(churn_data, select = - c(Customer_Value,FN,FP))
# Compute correlation at 2 decimal places
corr_matrix = round(cor(reduced_data), 2)
# Compute and show the result
ggcorrplot(corr_matrix, hc.order = TRUE, type = "lower",
lab = TRUE,lab_size = 2)
Podemos observar dos correlaciones fuertes porque su valor es superior a 0,8.
Age y Age_Group: 0.96 Frequency_of_Use y Second_of_Use: 0.95
Este resultado tiene sentido porque Age_Group se calcula a partir de Age. Además, el número total de segundos (Second_of_Use) se obtiene a partir del número total de llamadas (Frequency_of_Use).
En este caso, podemos deshacernos de Age_Group y Second_of_Use en el conjunto de datos.
second_model = lm(formula = Customer_Value ~ Call_Failure + Complains + Subscription_Length + Charge_Amount + Frequency_of_use + Frequency_of_SMS + Distinct_Called_Numbers + Tariff_Plan + Status + Age, data = churn_data)QQPlot
second_model_residuals <- second_model$residuals
# Crear un data frame para los residuales
residuals_df2 <- data.frame(Residuals = second_model_residuals)
ggplot(residuals_df2, aes(sample = Residuals)) +
stat_qq() +
stat_qq_line(color = "red") +
ggtitle("Q-Q Plot de los Residuales") +
xlab("Cuantiles Teóricos") +
ylab("Cuantiles de los Residuales") +
theme_minimal()
Podemos observar que más valores residuales están en la lÃnea recta en comparación con el primer modelo.
¿Cómo saber cuál de los dos modelos es mejor?
Una forma de responder a esta pregunta es realizar una prueba de análisis de la varianza (ANOVA) de los dos modelos. Esta prueba compara la hipótesis nula (Ho), según la cual las variables que eliminamos anteriormente no son significativas, con la hipótesis alternativa (Ha), según la cual esas variables son significativas. Si el nuevo modelo es una mejora del modelo original, entonces no rechazamos Ho. Si no es asÃ, significa que esas variables eran significativas; por lo tanto, rechazamos Ho.
# Anova test
anova(cust_value_model, second_model)## Analysis of Variance Table
##
## Model 1: Customer_Value ~ Call_Failure + Complains + Subscription_Length +
## Charge_Amount + Seconds_of_Use + Frequency_of_use + Frequency_of_SMS +
## Distinct_Called_Numbers + Age_Group + Tariff_Plan + Status +
## Age
## Model 2: Customer_Value ~ Call_Failure + Complains + Subscription_Length +
## Charge_Amount + Frequency_of_use + Frequency_of_SMS + Distinct_Called_Numbers +
## Tariff_Plan + Status + Age
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 3137 15092825
## 2 3139 23969449 -2 -8876625 922.49 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A partir del resultado del ANOVA, observamos que el valor p es muy pequeño (inferior a 0,05), por lo que rechazamos la hipótesis nula, lo que significa que el segundo modelo no es una mejora del primero. Otra forma de ver las variables importantes del modelo es mediante una prueba de significación.Una variable será significativa si su valor p es inferior a 0,05. Este resultado puede generarse mediante la función summary(). Además de proporcionar esa información sobre el modelo, también muestra la R-cuadrado ajustada, que evalúa el rendimiento de los modelos entre sÃ.
# Print the result of the model
summary(cust_value_model)##
## Call:
## lm(formula = Customer_Value ~ Call_Failure + Complains + Subscription_Length +
## Charge_Amount + Seconds_of_Use + Frequency_of_use + Frequency_of_SMS +
## Distinct_Called_Numbers + Age_Group + Tariff_Plan + Status +
## Age, data = churn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -709.81 -26.48 -2.63 24.24 191.43
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 160.201207 10.656561 15.033 < 2e-16 ***
## Call_Failure -0.489519 0.290861 -1.683 0.092474 .
## Complains 7.227189 5.000052 1.445 0.148439
## Subscription_Length 0.741287 0.156189 4.746 2.17e-06 ***
## Charge_Amount -14.298831 1.428753 -10.008 < 2e-16 ***
## Seconds_of_Use 0.047845 0.001116 42.875 < 2e-16 ***
## Frequency_of_use -0.540230 0.093055 -5.805 7.06e-09 ***
## Frequency_of_SMS 4.010644 0.012108 331.234 < 2e-16 ***
## Distinct_Called_Numbers 0.363675 0.112751 3.225 0.001271 **
## Age_Group -7.254712 5.211649 -1.392 0.164015
## Tariff_Plan 75.507316 5.669970 13.317 < 2e-16 ***
## Status -12.824538 3.884456 -3.302 0.000972 ***
## Age -7.098635 0.524340 -13.538 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 69.36 on 3137 degrees of freedom
## Multiple R-squared: 0.9821, Adjusted R-squared: 0.982
## F-statistic: 1.432e+04 on 12 and 3137 DF, p-value: < 2.2e-16En la sección de coeficientes podemos observar que Call_Failure, Complaints y Age_Group no son considerados significativos por el modelo porque su valor p es superior a 0,05. Mantenerlos no aporta ningún valor adicional al modelo.
Aplicando el mismo análisis al segundo modelo, obtenemos este resultado:
summary(second_model)##
## Call:
## lm(formula = Customer_Value ~ Call_Failure + Complains + Subscription_Length +
## Charge_Amount + Frequency_of_use + Frequency_of_SMS + Distinct_Called_Numbers +
## Tariff_Plan + Status + Age, data = churn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -709.44 -40.39 -3.30 46.36 255.80
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 221.57401 13.30390 16.655 < 2e-16 ***
## Call_Failure -5.48129 0.33591 -16.318 < 2e-16 ***
## Complains 18.13506 6.27432 2.890 0.00387 **
## Subscription_Length 0.97222 0.19653 4.947 7.94e-07 ***
## Charge_Amount 16.40113 1.54597 10.609 < 2e-16 ***
## Frequency_of_use 3.11880 0.04622 67.484 < 2e-16 ***
## Frequency_of_SMS 4.01953 0.01505 267.021 < 2e-16 ***
## Distinct_Called_Numbers -0.42846 0.13827 -3.099 0.00196 **
## Tariff_Plan -4.67264 6.61666 -0.706 0.48012
## Status -0.01064 4.87862 -0.002 0.99826
## Age -8.15912 0.19387 -42.085 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 87.38 on 3139 degrees of freedom
## Multiple R-squared: 0.9715, Adjusted R-squared: 0.9714
## F-statistic: 1.071e+04 on 10 and 3139 DF, p-value: < 2.2e-16El modelo original tiene una R-cuadrado ajustada de 0,98, que es superior a la R-cuadrado ajustada del segundo modelo (0,97). Esto significa que el modelo original con todos los predictores es mejor que el segundo modelo.
Otra estrategia para elegir eficientemente los predictores relevantes es a través del Criterio de Información de Akaike (AIC).
# Realizar la selección de predictores basada en AIC
Selected_model <- stepAIC(cust_value_model, direction = "both")## Start: AIC=26720.9
## Customer_Value ~ Call_Failure + Complains + Subscription_Length +
## Charge_Amount + Seconds_of_Use + Frequency_of_use + Frequency_of_SMS +
## Distinct_Called_Numbers + Age_Group + Tariff_Plan + Status +
## Age
##
## Df Sum of Sq RSS AIC
## - Age_Group 1 9323 15102148 26721
## <none> 15092825 26721
## - Complains 1 10052 15102877 26721
## - Call_Failure 1 13628 15106453 26722
## - Distinct_Called_Numbers 1 50054 15142879 26729
## - Status 1 52442 15145267 26730
## - Subscription_Length 1 108374 15201199 26741
## - Frequency_of_use 1 162156 15254981 26753
## - Charge_Amount 1 481885 15574709 26818
## - Tariff_Plan 1 853243 15946067 26892
## - Age 1 881819 15974644 26898
## - Seconds_of_Use 1 8844161 23936986 28172
## - Frequency_of_SMS 1 527868092 542960917 38005
##
## Step: AIC=26720.85
## Customer_Value ~ Call_Failure + Complains + Subscription_Length +
## Charge_Amount + Seconds_of_Use + Frequency_of_use + Frequency_of_SMS +
## Distinct_Called_Numbers + Tariff_Plan + Status + Age
##
## Df Sum of Sq RSS AIC
## - Complains 1 8769 15110916 26721
## <none> 15102148 26721
## + Age_Group 1 9323 15092825 26721
## - Call_Failure 1 13908 15116056 26722
## - Status 1 53036 15155184 26730
## - Distinct_Called_Numbers 1 58358 15160506 26731
## - Subscription_Length 1 106356 15208504 26741
## - Frequency_of_use 1 166606 15268753 26753
## - Charge_Amount 1 498865 15601013 26821
## - Tariff_Plan 1 911397 16013544 26903
## - Seconds_of_Use 1 8867302 23969449 28174
## - Age 1 12311082 27413230 28597
## - Frequency_of_SMS 1 541037349 556139497 38078
##
## Step: AIC=26720.68
## Customer_Value ~ Call_Failure + Subscription_Length + Charge_Amount +
## Seconds_of_Use + Frequency_of_use + Frequency_of_SMS + Distinct_Called_Numbers +
## Tariff_Plan + Status + Age
##
## Df Sum of Sq RSS AIC
## <none> 15110916 26721
## - Call_Failure 1 9991 15120908 26721
## + Complains 1 8769 15102148 26721
## + Age_Group 1 8040 15102877 26721
## - Status 1 47145 15158061 26728
## - Distinct_Called_Numbers 1 59788 15170705 26731
## - Subscription_Length 1 102118 15213034 26740
## - Frequency_of_use 1 174441 15285357 26755
## - Charge_Amount 1 514452 15625368 26824
## - Tariff_Plan 1 920559 16031475 26905
## - Seconds_of_Use 1 8922326 24033242 28180
## - Age 1 12303516 27414432 28595
## - Frequency_of_SMS 1 541082605 556193521 38077summary(Selected_model)##
## Call:
## lm(formula = Customer_Value ~ Call_Failure + Subscription_Length +
## Charge_Amount + Seconds_of_Use + Frequency_of_use + Frequency_of_SMS +
## Distinct_Called_Numbers + Tariff_Plan + Status + Age, data = churn_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -706.84 -26.43 -4.15 23.63 192.64
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 159.592566 10.645123 14.992 < 2e-16 ***
## Call_Failure -0.409082 0.283955 -1.441 0.149780
## Subscription_Length 0.716772 0.155625 4.606 4.27e-06 ***
## Charge_Amount -14.652913 1.417429 -10.338 < 2e-16 ***
## Seconds_of_Use 0.047969 0.001114 43.052 < 2e-16 ***
## Frequency_of_use -0.557552 0.092621 -6.020 1.95e-09 ***
## Frequency_of_SMS 4.007718 0.011954 335.260 < 2e-16 ***
## Distinct_Called_Numbers 0.392502 0.111374 3.524 0.000431 ***
## Tariff_Plan 77.185663 5.581625 13.829 < 2e-16 ***
## Status -11.965735 3.823597 -3.129 0.001767 **
## Age -7.792786 0.154144 -50.555 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 69.38 on 3139 degrees of freedom
## Multiple R-squared: 0.982, Adjusted R-squared: 0.982
## F-statistic: 1.717e+04 on 10 and 3139 DF, p-value: < 2.2e-16El procedimiento de selección basado en AIC ha identificado un conjunto de variables que optimizan el ajuste del modelo, descartando aquellas que no contribuyen significativamente a la explicación de la variabilidad en la variable dependiente Customer_Value.
Validación de los supuestos
La validación de los supuestos en un modelo de regresión múltiple se realiza de manera similar a la validación en un modelo de regresión simple, pero con algunas consideraciones adicionales debido a la mayor cantidad de variables independientes.
Prueba de la Media de los Errores (Prueba t de Student)
t.test(Selected_model$residuals)##
## One Sample t-test
##
## data: Selected_model$residuals
## t = -7.2573e-16, df = 3149, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -2.420019 2.420019
## sample estimates:
## mean of x
## -8.957396e-16El p-valor de la prueba es mayor que el nivel de significancia de 0.05, no se rechaza la hipótesis nula, lo que sugiere que la media de los errores es cero.
Prueba de Normalidad de los Errores (Prueba de Anderson Darling n>=30)
ad.test(Selected_model$residuals)##
## Anderson-Darling normality test
##
## data: Selected_model$residuals
## A = 140.47, p-value < 2.2e-16El p-valor es menor que 0.05, se rechaza la hipótesis nula, lo que indica que los errores no siguen una distribución normal.
# Crear el gráfico de histograma con lÃnea de densidad
ggplot(data.frame(Residuals = Selected_model$residuals), aes(x = Residuals)) +
geom_histogram(aes(y = ..density..), bins = 20, color = "black", fill = "maroon") +
geom_density(color = "lightblue", size = 1) +
ggtitle("Distribución de los Residuales con LÃnea de Densidad") +
xlab("Residuales") +
ylab("Densidad") +
theme_minimal()
Otra prueba de normalidad (Kolmogorov-Smirnov Test) que compara la distribución empÃrica de los residuos con una distribución normal.
ks.test(Selected_model$residuals, "pnorm", mean = mean(Selected_model$residuals), sd = sd(Selected_model$residuals))## Warning in ks.test.default(Selected_model$residuals, "pnorm", mean =
## mean(Selected_model$residuals), : ties should not be present for the one-sample
## Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Selected_model$residuals
## D = 0.17628, p-value < 2.2e-16
## alternative hypothesis: two-sidedPrueba de Homocedasticidad (Prueba de Breusch-Pagan)
bptest(Selected_model)##
## studentized Breusch-Pagan test
##
## data: Selected_model
## BP = 392.93, df = 10, p-value < 2.2e-16El p-valor es menor que 0.05, se rechaza la hipótesis nula, sugiriendo que no hay homocedasticidad, es decir, la varianza de los errores no es constante.
Prueba de Independencia de los Errores (Prueba de Durbin-Watson)
dwtest(Selected_model,alternative = "two.sided")##
## Durbin-Watson test
##
## data: Selected_model
## DW = 2.1641, p-value = 4.845e-06
## alternative hypothesis: true autocorrelation is not 0El valor de Durbin-Watson es un valor cercano a 2 lo que sugiere que no hay autocorrelación entre los residuos. El p-valor asociado es menor que 0.05, se rechaza la hipótesis nula, sugiriendo que los errores no son independientes.
Multicolinealidad
Evalúa la presencia de multicolinealidad entre las variables independientes. Valores de VIF superiores a 10 (algunos usan un umbral de 5) pueden indicar multicolinealidad preocupante.
library(car)## Cargando paquete requerido: carData##
## Adjuntando el paquete: 'car'## The following object is masked from 'package:dplyr':
##
## recodevif(Selected_model)## Call_Failure Subscription_Length Charge_Amount
## 2.782975 1.164519 3.040712
## Seconds_of_Use Frequency_of_use Frequency_of_SMS
## 14.311248 18.497810 1.177552
## Distinct_Called_Numbers Tariff_Plan Status
## 2.405317 1.462256 1.785347
## Age
## 1.212155Se debe revisar la relación entre Frequency_of_use y Seconds_of_Use con otras variables para ver si alguna de ellas puede ser eliminada o combinada.
Verificación de la estructura funcional
RESET Test (Regression Equation Specification Error Test) se utiliza para evaluar si un modelo de regresión lineal está bien especificado, es decir, si se han omitido variables importantes, si se ha utilizado una forma funcional incorrecta, o si hay interacción entre las variables que no se ha modelado adecuadamente.
resettest(Selected_model)##
## RESET test
##
## data: Selected_model
## RESET = 602.91, df1 = 2, df2 = 3137, p-value < 2.2e-16Un valor de RESET elevado generalmente sugiere que el modelo podrÃa estar mal especificado y el valor p extremadamente bajo indica que hay suficiente evidencia para rechazar la hipótesis nula de que el modelo está correctamente especificado.
Linealidad de la relación entre las variables independientes y la dependiente
Estos gráficos ayudan a visualizar si la relación entre cada predictor y la respuesta es lineal.
crPlots(Selected_model)
