Cigarette Consumption: linear regression Analysis
Luiselli Lucas
2023-10-18
Se requiere construir una ecuación de regresión que relacione el consumo de cigarrillos per-cápita en todo el estado (variable Sales), con diversas variables socioeconómicas y demográficas, y determinar si estas variables son útiles para predecir el consumo de los cigarrillos
library(haven) # importar sas,pdf,spss
library(goftest) # Pruebas de bondad de ajuste de Cramer-Von Mises y Anderson-Darling para distribuciones univariadas continuas, utilizando algoritmos eficientes.
library(nortest) # Prueba de normalidad##
## Attaching package: 'nortest'## The following objects are masked from 'package:goftest':
##
## ad.test, cvm.testlibrary(oddsratio)
library(olsrr, verbose = FALSE) # Modelos de regresion##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riverslibrary(ggplot2)library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ lubridate 1.9.2 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(ggpubr)
library(rstatix,verbose = FALSE)##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filterImportar data
base_fumadores <- read_sav("p081.sav")
base_fumadores## # A tibble: 51 × 8
## STATE AGE HS INCOME BLACK FEMALE PRICE SALES
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AL 27 41.3 2948 26.2 51.7 42.7 89.8
## 2 AK 22.9 66.7 4644 3 45.7 41.8 121.
## 3 AZ 26.3 58.1 3665 3 50.8 38.5 115.
## 4 AR 29.1 39.9 2878 18.3 51.5 38.8 100.
## 5 CA 28.1 62.6 4493 7 50.8 39.7 123
## 6 CO 26.2 63.9 3855 3 50.7 31.1 125.
## 7 CT 29.1 56 4917 6 51.5 45.5 120
## 8 DE 26.8 54.6 4524 14.3 51.3 41.3 155
## 9 DC 28.4 55.2 5079 71.1 53.5 32.6 200.
## 10 FL 32.3 52.6 3738 15.3 51.8 43.8 124.
## # ℹ 41 more rowsCuál es el nivel de asociación lineal de las variables predictoras con la variable Sales? Comentar
# Gráfico de dispersión para la variable "AGE" vs. "SALES"
p1<-ggplot(base_fumadores, aes(x = AGE, y = SALES)) +
geom_point() +
geom_smooth(method = "lm")+
labs(x = "AGE", y = "SALES") +
ggtitle("Diagrama de Dispersión: AGE vs. SALES")
# Gráfico de dispersión para la variable "PRICE" vs. "SALES"
p2<-ggplot(base_fumadores, aes(x = PRICE, y = SALES)) +
geom_point() +
geom_smooth(method = "lm")+
labs(x = "PRICE", y = "SALES") +
ggtitle("Diagrama de Dispersión: PRICE vs. SALES")
# Gráfico de dispersión para la variable "INCOME" vs. "SALES"
p3<-ggplot(base_fumadores, aes(x = INCOME, y = SALES)) +
geom_point() +
geom_smooth(method = "lm")+
labs(x = "INCOME", y = "SALES") +
ggtitle("Diagrama de Dispersión: INCOME vs. SALES")
# Gráfico de dispersión para la variable "BLACK" vs. "SALES"
p4<-ggplot(base_fumadores, aes(x = BLACK, y = SALES)) +
geom_point() +
geom_smooth(method = "lm")+
labs(x = "BLACK", y = "SALES") +
ggtitle("Diagrama de Dispersión: BLACK vs. SALES")
plot_total<-ggarrange(p1,p2,p3,p4, ncol = 2, nrow = 2)## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'plot_total
Establecida la Correlación tomando a Sales como variable dependiente y a Age,Income,Black,Female y Price
corr_1 <- cor(base_fumadores[, c("SALES", "AGE", "INCOME", "BLACK", "FEMALE", "PRICE")])
cor_values <- corr_1["SALES", ]
cor_values## SALES AGE INCOME BLACK FEMALE PRICE
## 1.0000000 0.2265549 0.3260679 0.1895903 0.1462211 -0.3006227El diagrama de dispersión nos muestra cómo es el tipo de asociación (lineal, curvilínea, positiva o negativa) e incluso si no hay asociación entre las variables. El coeficiente de correlación de Pearson (r) nos indica que:
“si es -1, indica una correlación lineal perfectamente negativa entre dos variables” “si es 0, indica que no hay correlación lineal entre dos variables” “si es 1, indica una correlación lineal perfectamente positiva entre dos variables”
Age: 0.226. Correlación positiva debil
Income: 0.326. Correlación positiva moderada
Black: 0.189. Correlación positiva debil
Female: 0.146. Correlación positiva debil
Price: -0.3006. Correlación negativa moderada
Se Realiza una regresión lineal múltiple, seleccionando los mejores predictores entre las variables independientes disponibles, utilizando un método de selección automática.
model1 <- lm(SALES ~ AGE + INCOME + BLACK + FEMALE + PRICE, data = base_fumadores)
summary(model1)##
## Call:
## lm(formula = SALES ~ AGE + INCOME + BLACK + FEMALE + PRICE, data = base_fumadores)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.136 -12.508 -5.277 6.351 133.141
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 100.133448 239.217526 0.419 0.67751
## AGE 4.593832 3.035857 1.513 0.13722
## INCOME 0.018418 0.007371 2.499 0.01618 *
## BLACK 0.379049 0.391066 0.969 0.33759
## FEMALE -1.067119 5.496066 -0.194 0.84692
## PRICE -3.243835 1.009591 -3.213 0.00243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.86 on 45 degrees of freedom
## Multiple R-squared: 0.3208, Adjusted R-squared: 0.2453
## F-statistic: 4.25 on 5 and 45 DF, p-value: 0.002998Otro método para analizar la correlación de una manera más ágil sobre el conjunto de la base de datos en relación a la variable dependiente, es el siguiente: los valores de Zero Order nos muestran el p-value de las correlaciones que vimos anteriormente.
model2 <- ols_regress(SALES ~ AGE + INCOME + BLACK + FEMALE + PRICE, data = base_fumadores)
model2## Model Summary
## ---------------------------------------------------------------
## R 0.566 RMSE 27.861
## R-Squared 0.321 Coef. Var 22.923
## Adj. R-Squared 0.245 MSE 776.233
## Pred R-Squared 0.130 MAE 16.123
## ---------------------------------------------------------------
## RMSE: Root Mean Square Error
## MSE: Mean Square Error
## MAE: Mean Absolute Error
##
## ANOVA
## --------------------------------------------------------------------
## Sum of
## Squares DF Mean Square F Sig.
## --------------------------------------------------------------------
## Regression 16494.939 5 3298.988 4.25 0.0030
## Residual 34930.507 45 776.233
## Total 51425.445 50
## --------------------------------------------------------------------
##
## Parameter Estimates
## --------------------------------------------------------------------------------------------
## model Beta Std. Error Std. Beta t Sig lower upper
## --------------------------------------------------------------------------------------------
## (Intercept) 100.133 239.218 0.419 0.678 -381.675 581.942
## AGE 4.594 3.036 0.269 1.513 0.137 -1.521 10.708
## INCOME 0.018 0.007 0.342 2.499 0.016 0.004 0.033
## BLACK 0.379 0.391 0.149 0.969 0.338 -0.409 1.167
## FEMALE -1.067 5.496 -0.037 -0.194 0.847 -12.137 10.003
## PRICE -3.244 1.010 -0.418 -3.213 0.002 -5.277 -1.210
## --------------------------------------------------------------------------------------------El método de selección automática de predictores Stepwise, es una combinación del método Forward y Backward. Revisa si agregando o sacando una variable puede afectar al modelo puede afectar el R2, debido a que el ingreso o exclusión de una variable puede afectar la importancia de las variables que ya estaban en el modelo. Para aplicar esta función necesitamos calcular la regresión con la función del model1
ols_step_best_subset(model1)## Best Subsets Regression
## --------------------------------------------
## Model Index Predictors
## --------------------------------------------
## 1 INCOME
## 2 INCOME PRICE
## 3 AGE INCOME PRICE
## 4 AGE INCOME BLACK PRICE
## 5 AGE INCOME BLACK FEMALE PRICE
## --------------------------------------------
##
## Subsets Regression Summary
## ---------------------------------------------------------------------------------------------------------------------------------------
## Adj. Pred
## Model R-Square R-Square R-Square C(p) AIC SBIC SBC MSEP FPE HSP APC
## ---------------------------------------------------------------------------------------------------------------------------------------
## 1 0.1063 0.0881 0.0325 12.2063 497.7181 352.3936 503.5136 47835.2430 974.6969 19.5399 0.9666
## 2 0.2503 0.2190 0.153 4.6690 490.7602 346.1878 498.4875 40983.5056 850.4721 17.0899 0.8434
## 3 0.3032 0.2588 0.1888 3.1601 489.0237 345.1087 498.6828 38916.2141 822.1552 16.5731 0.8154
## 4 0.3202 0.2611 0.1807 4.0377 489.7682 346.3253 501.3592 38813.7221 834.5043 16.8888 0.8276
## 5 0.3208 0.2453 0.1299 6.0000 491.7255 348.5582 505.2483 39662.6252 867.5551 17.6417 0.8604
## ---------------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria
## SBIC: Sawa's Bayesian Information Criteria
## SBC: Schwarz Bayesian Criteria
## MSEP: Estimated error of prediction, assuming multivariate normality
## FPE: Final Prediction Error
## HSP: Hocking's Sp
## APC: Amemiya Prediction CriteriaPara clarificar los elementos del informe: el modelo seleccionado es el modelo 2 con dos predictores, siguiendo el criterio de que con una base de pocas observaciones y muchas variables, es recomendable tener un modelo con pocas predictores. Además se observa en los resultados que el incremento de R2 ajustado es mayor entre la primera y segunda predictora y disminuye en cuanto a la segunda y tercara predictora. Pareciera suficiente trabajar con dos predictores.
modelo_ajustado2 <- lm(SALES ~ INCOME + PRICE, data = base_fumadores)
summary(modelo_ajustado2)##
## Call:
## lm(formula = SALES ~ INCOME + PRICE, data = base_fumadores)
##
## Residuals:
## Min 1Q Median 3Q Max
## -62.560 -11.863 -4.597 3.796 132.755
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 153.33842 41.23890 3.718 0.000525 ***
## INCOME 0.02208 0.00690 3.200 0.002440 **
## PRICE -3.01756 0.99396 -3.036 0.003868 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.34 on 48 degrees of freedom
## Multiple R-squared: 0.2503, Adjusted R-squared: 0.219
## F-statistic: 8.012 on 2 and 48 DF, p-value: 0.0009945Con un R-cuadrado de 0.2503, es decir, alrededor del 25.03% de la variabilidad en las ventas (SALES) se explicaría por las pridictores INCOME y PRICE en este modelo ajustado. En paralelo y según las pruebas de significancia, ambos predictores son estadísticamente significativas en el modelo (p < 0.05). Esto significa que tanto INCOME como PRICE tienen un impacto estadísticamente significativo en las ventas. A partir de lo dicho podemos concluir que un aumento de INCOME se asocia con un aumento en las ventas, mientras que PRICE tiene un coeficiente negativo, lo que sugiere que un aumento en el precio se asocia con una disminución en las ventas.
-Qué información da el coeficiente de determinación?
El R-cuadrado nos brinda información sobre la proporción de variabilidad en la variable dependiente (en este caso, las ventas, SALES) puede ser explicada por las variables independientes incluidas en el modelo. Un R-cuadrado del 0.2503 sugiere que las variables INCOME y PRICE tienen cierta capacidad para explicar las fluctuaciones en las ventas, pero aún hay una cantidad significativa de variabilidad que no se captura con estas dos variables.
Un R-cuadrado ajustado de 0.219 (21.9%) no es muy alto, lo que sugiere que estas dos variables por sí solas no explican la mayoría de la variabilidad en las ventas. Esto podría indicar que se necesitan otras variables o características adicionales para mejorar la capacidad del modelo de predecir las ventas con precisión.
-Cuáles son los supuestos necesarios para definir la prueba inferencial de los estimadores de los parámeros?
Los supuestos que deben cumplirse para aplicar una regresión lineal son: que la relación sea lineal (INCOME Y PRICE) con la variable dependiente (SALE), errores con distribución normal, que varianza sea constante y la independencia de la muestra.
Normalidad: se asume que los errores de predicción (residuales) siguen una distribución normal. Como se puede ver, no hay evidencia para rechazar la hipótesis de normalidad. Se observa en la bisectriz, que los cuantiles teóricos se solapan con los muestrales
test_normalidad_model <- ols_test_normality(modelo_ajustado2)
test_normalidad_model## -----------------------------------------------
## Test Statistic pvalue
## -----------------------------------------------
## Shapiro-Wilk 0.7864 0.0000
## Kolmogorov-Smirnov 0.2073 0.0212
## Cramer-von Mises 4.8041 0.0000
## Anderson-Darling 2.9946 0.0000
## -----------------------------------------------ols_plot_resid_qq(modelo_ajustado2)
Si se considera que en el siguiente gráfico se muestran como datos normales los que estan en la banda de 2 y -2 (que representan un 95% de la distribución normal), es posible que los casos 29, 30 y 2 sean outliers en las colas, correspondientes al 2.5% en cada cola. Pero como se nombró previamente, no hay evidencia para rechazar la hipótesis de normalidad
ols_plot_resid_stud_fit(modelo_ajustado2, print_plot = TRUE)
-Analizar la bondad del ajuste del modelo obtenido, comentando los indicadores y/o test que considera.
La bondad de ajuste de un modelo de regresión se puede evaluar a través de elementos estadísticos y gráficos.
Residuos: Los residuos son las diferencias entre los valores observados y los valores predichos por el modelo. En este caso, los residuos tienen una distribución que varía desde -62.560 hasta 132.755. El primer cuartil (1Q) es -11.863, lo que significa que el 25% de los residuos son menores o iguales a este valor. La mediana de los residuos es -4.597, lo que indica que el 50% de los residuos son menores o iguales a este valor. El tercer cuartil (3Q) es 3.796, lo que indica que el 75% de los residuos son menores o iguales a este valor. Estos valores proporcionan información sobre la dispersión de los residuos y su distribución.
Coeficientes: Los coeficientes del modelo indican cómo las variables predictoras (INCOME y PRICE) están relacionadas con la variable de respuesta (SALES). Los coeficientes tienen valores estimados, errores estándar, valores t y valores p. Los valores t y los valores p se utilizan para evaluar la significancia estadística de los coeficientes. En este caso, el Intercept tiene un valor estimado significativamente diferente de cero, al igual que las variables predictoras INCOME y PRICE, según los valores p que son menores que 0.05. Esto sugiere que estas variables son significativas para el modelo.
Error Estándar Residual: El error estándar residual es una medida que indica cuánto varían los valores observados alrededor de las predicciones del modelo. En este caso, el error estándar residual es de aproximadamente 28.34.
R-cuadrado: El R-cuadrado, como ya se señaló, es 0.2503, lo que significa que aproximadamente el 25.03% de la variabilidad en las ventas (SALES) se explica mediante el modelo de regresión que incluye las variables predictoras INCOME y PRICE. Esto indica que el modelo tiene un poder explicativo limitado sobre las ventas. El R-cuadrado ajustado, es 0.219, que tiene en cuenta la complejidad del modelo y proporciona una medida más realista de la calidad del ajuste. Indica que alrededor del 21.9% de la variabilidad en SALES se explica por el modelo, teniendo en cuenta la complejidad del mismo.
Estadística F y Valor p: La estadística F es 8.012, y el valor p asociado es 0.0009945. Esta estadística y valor p se utilizan para evaluar si el modelo en su conjunto es significativo. En este caso, el valor p es muy bajo, por lo tanto podríamos decir que el modelo en su conjunto es estadísticamente significativo.
Respecto a la bondad del ajuste en base a lo desarrollado, podemos decir que el modelo de regresión muestra que las variables predictoras INCOME y PRICE tienen efectos significativos en las ventas (SALES). Sin embargo, el modelo explica solo una parte modesta de la variabilidad en las ventas, como lo indica el R-cuadrado. Es importante considerar otros factores y evaluar la adecuación del modelo en función de los objetivos de análisis y las necesidades de predicción
-Realizar un análisis de los residuos del modelo para evaluar el cumplimiento de los supuestos. Para esto, realizar gráficos de los residuos con el valor predicho.
modelo_ajustado2 <- lm(SALES ~ INCOME + PRICE, data = base_fumadores)
summary(modelo_ajustado2)##
## Call:
## lm(formula = SALES ~ INCOME + PRICE, data = base_fumadores)
##
## Residuals:
## Min 1Q Median 3Q Max
## -62.560 -11.863 -4.597 3.796 132.755
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 153.33842 41.23890 3.718 0.000525 ***
## INCOME 0.02208 0.00690 3.200 0.002440 **
## PRICE -3.01756 0.99396 -3.036 0.003868 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.34 on 48 degrees of freedom
## Multiple R-squared: 0.2503, Adjusted R-squared: 0.219
## F-statistic: 8.012 on 2 and 48 DF, p-value: 0.0009945Los residuos son las diferencias entre los valores observados y los valores predichos por el modelo. Aquí es posible ver que los residuos varían desde -62.560 hasta 132.755, sugiriendo que los residuos tienen una amplia dispersión y no parecen seguir un patrón claro. Sin embargo, como se vio en el gráfico de residuos, para verificar estos cumplen con supuesto de normalidad
El R-cuadrado es 0.2503, lo que significa que el modelo explica aproximadamente el 25.03% de la variabilidad en las ventas.
El R-cuadrado ajustado es 0.219, lo que indica que, después de ajustar por el número de variables en el modelo, el poder explicativo disminuye un poco.
El estadístico F y su p-valor se utilizan para evaluar si el modelo (en su conjunto) es estadísticamente significativo. En este caso se oberva un p-valor es muy bajo (0.0009945), lo que sugiere que el modelo es estadísticamente significativo en la explicación de las ventas.
-Analizar la colinealidad de las variables predictoras presentes en la ecuación.
colinealidad <- ols_vif_tol(modelo_ajustado2)
colinealidad## Variables Tolerance VIF
## 1 INCOME 0.9539652 1.048256
## 2 PRICE 0.9539652 1.048256TOL (Tolerancia) es una medida de cuánto de la variabilidad de una variable independiente no está correlacionada con las otras variables independientes. VIF (Variance inflation factor) es simplemente el inverso de la tolerancia. Se utiliza para cuantificar cuánto aumenta la varianza de los coeficientes de regresión debido a la multicolinealidad. Un VIF alto (generalmente mayor que 10) indica multicolinealidad significativa.
Estos valores (1.04) sugieren que no hay una fuerte multicolinealidad entre estas dos variables ya que tanto la tolerancia como el VIF están cerca de 1, lo que indica que no hay una alta correlación entre INCOME y PRICE.
-Analizar la presencia de observaciones atípicas y/o influyentes. Comentar y resolver según el caso.
Se utiliza la distancia de Cook para analizar si hay valores atípicos o puntos de influencia
ols_plot_cooksd_bar(modelo_ajustado2, print_plot = TRUE)
Suponiendo que por la lejanía observadas en el gráfico de residudos anterior, donde las posiciones 29 y 30 resultan outliers que se alejan considerablemente de la banda -2 y 2 y se retiran de la base:
base_fumadores_ajustada <- base_fumadores[c(-29,-30),]
base_fumadores_ajustada## # A tibble: 49 × 8
## STATE AGE HS INCOME BLACK FEMALE PRICE SALES
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AL 27 41.3 2948 26.2 51.7 42.7 89.8
## 2 AK 22.9 66.7 4644 3 45.7 41.8 121.
## 3 AZ 26.3 58.1 3665 3 50.8 38.5 115.
## 4 AR 29.1 39.9 2878 18.3 51.5 38.8 100.
## 5 CA 28.1 62.6 4493 7 50.8 39.7 123
## 6 CO 26.2 63.9 3855 3 50.7 31.1 125.
## 7 CT 29.1 56 4917 6 51.5 45.5 120
## 8 DE 26.8 54.6 4524 14.3 51.3 41.3 155
## 9 DC 28.4 55.2 5079 71.1 53.5 32.6 200.
## 10 FL 32.3 52.6 3738 15.3 51.8 43.8 124.
## # ℹ 39 more rowsmodel_ajustada_3 <- lm(SALES ~ INCOME + PRICE, data = base_fumadores_ajustada)
summary(model_ajustada_3)##
## Call:
## lm(formula = SALES ~ INCOME + PRICE, data = base_fumadores_ajustada)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.100 -8.694 0.040 4.469 43.319
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 153.240894 27.474847 5.577 1.24e-06 ***
## INCOME 0.018601 0.004487 4.146 0.000144 ***
## PRICE -2.780184 0.654177 -4.250 0.000103 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.2 on 46 degrees of freedom
## Multiple R-squared: 0.393, Adjusted R-squared: 0.3666
## F-statistic: 14.89 on 2 and 46 DF, p-value: 1.033e-05A partir del ajuste se puede destacar lo siguiente: - El R-cuadrado ajustado en el Modelo 2 (0.3666) es más alto que en el Modelo 1 (0.219). - El Modelo 2 tiene un error estándar residual más bajo (18.2) en comparación con el Modelo 1 (28.34). Un error estándar residual más bajo indica un mejor ajuste del modelo a los datos. - En ambos modelos, los coeficientes de INCOME y PRICE son significativos (valores p < 0.05). Sin embargo, en el Modelo 2, los coeficientes tienen un valor p aún más bajo, lo que indica una mayor significancia.