Practica 3: Regresión múltiple, control estadístico y parcialización
Estadistica Multivariada - Sociología FACSO Universidad de Chile
Juan Carlos Castillo & Alejandro Plaza
1 Introducción
En esta práctica nos enfocaremos en la comprensión del control estadístico a partir de la visualización de datos y la regresión parcial. Vamos a realizar dos ejercicios, el primero con los datos de la sesión correspondiene, y el segundo replicando el ejemplo de Darlington & Hayes
2 Ejemplo 1: ingreso, educación y experiencia.
Este ejericio corresponde al que se detalla en la sesión 6 ver aquí
2.1 Librerías y datos
pacman::p_load(dplyr,
corrplot,
ggplot2,
texreg,
stargazer)
datos=read.csv("https://juancarloscastillo.github.io/metsoc-facsouchile/documents/presentaciones/6regmul2/ingedexp.csv", sep="")Ver datos
##
## =================================================
## Statistic N Mean St. Dev. Min Max
## -------------------------------------------------
## ID 10 6 3 1 10
## ingreso 10 495,000 323,565 100,000 1,000,000
## educacion 10 4 2 1 6
## experiencia 10 2 1 1 3
## -------------------------------------------------2.2 Correlaciones y scatters
## ingreso educacion experiencia
## ingreso 1.00 0.77 0.44
## educacion 0.77 1.00 0.27
## experiencia 0.44 0.27 1.00
2.3 Regresiones
reg_y_x1=lm(ingreso ~ educacion, data=datos)
reg_y_x2=lm(ingreso ~ experiencia, data=datos)
reg_y_x1_x2=lm(ingreso ~ educacion + experiencia , data=datos)
screenreg(list(reg_y_x1,reg_y_x2,reg_y_x1_x2))##
## =====================================================
## Model 1 Model 2 Model 3
## -----------------------------------------------------
## (Intercept) -91566.27 93442.62 -270638.30
## (183509.80) (302389.31) (241882.27)
## educacion 150401.61 ** 137092.20 *
## (43618.69) (44602.35)
## experiencia 174590.16 100425.53
## (124491.71) (90114.05)
## -----------------------------------------------------
## R^2 0.60 0.20 0.66
## Adj. R^2 0.55 0.10 0.56
## Num. obs. 10 10 10
## RMSE 217656.84 307471.84 214438.31
## =====================================================
## *** p < 0.001, ** p < 0.01, * p < 0.052.4 Parciales
reg_x1_x2=lm(educacion ~ experiencia , data=datos)
x1_fit_x2=fitted.values(reg_x1_x2)
resx1_2=residuals(reg_x1_x2)
datos=cbind(datos, x1_fit_x2,resx1_2)
datos## ID ingreso educacion experiencia x1_fit_x2 resx1_2
## 1 1 100000 2 1 3.196721 -1.1967213
## 2 2 200000 1 3 4.278689 -3.2786885
## 3 3 200000 3 3 4.278689 -1.2786885
## 4 4 300000 3 1 3.196721 -0.1967213
## 5 5 300000 4 2 3.737705 0.2622951
## 6 6 500000 6 2 3.737705 2.2622951
## 7 7 650000 5 3 4.278689 0.7213115
## 8 8 800000 4 2 3.737705 0.2622951
## 9 9 900000 5 3 4.278689 0.7213115
## 10 10 1000000 6 3 4.278689 1.7213115regy_resx1_2=lm(datos$ingreso ~ resx1_2)
screenreg(list(reg_y_x1,reg_y_x2,reg_y_x1_x2,regy_resx1_2), booktabs = TRUE, dcolumn = TRUE, doctype = FALSE, caption=" ")##
## ====================================================================
## Model 1 Model 2 Model 3 Model 4
## --------------------------------------------------------------------
## (Intercept) -91566.27 93442.62 -270638.30 495000.00 ***
## (183509.80) (302389.31) (241882.27) (79673.16)
## educacion 150401.61 ** 137092.20 *
## (43618.69) (44602.35)
## experiencia 174590.16 100425.53
## (124491.71) (90114.05)
## resx1_2 137092.20 *
## (52404.36)
## --------------------------------------------------------------------
## R^2 0.60 0.20 0.66 0.46
## Adj. R^2 0.55 0.10 0.56 0.39
## Num. obs. 10 10 10 10
## RMSE 217656.84 307471.84 214438.31 251948.67
## ====================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05¿Qué nos indica el coeficiente de regresión parcial? ¿Cómo se podría obtener el segundo coeficiente?
3 Ejercicio 2: Ejemplo de Darlington & Hayes
3.1 Datos
Los datos corresponden a la base de dato exercise que muestran Darlington y Hayes (2017) en el capitulo 3 de regresión multiple.
wtloss: variable continua de perdida de peso en gramos. food : variable continua de ingesta de comida medida en calorias exercise: promedio de horas de ejercicio a la semana.
3.2 Cargar bases de datos y librerias.
#cargar base de datos
library(readr)
library(dplyr)
library(ggplot2)
library(car)
library(stargazer)
library(texreg)
# Leer datos
exercise <- read_table2("https://juancarloscastillo.github.io/metsoc-facsouchile/documents/practicas/3regmul/exercise.txt")
exercise %>% select(id, exercise, food, wtloss)->data #filtrar casos con variables a examinar3.3 Sobre relaciones contraintuitivas
cor(data[2:4])Como se puede ver en los datos la correlación entre ejercicio y perdida de peso es de 0.86. Aquellas personas que se ejercitan más tienen una mayor reducción de peso. Por otro lado la correlación entre el consumo de de comida y la perdida de peso es de 0.047. Ignorando el hecho de que es una correlación bastante pequeña, se puede observar una relación positiva, que al menos exploratoriamente es contraintutiva. Es decir a mayor consumo de comida existiría una mayor reducción de peso.
Para entender este tipo de relaciones contraintuitivas una buena aproximación podría ser visualizar qué estaría ocurriendo con los datos. En el siguiente gráfico de dispersión se observa una relación levemente positiva entre la ingesta de comida y la perdida de peso.
#Relacion lineal entre comida y perdida de peso
ggplot(data, aes(x=food, y=wtloss)) +
geom_point() +
geom_smooth(method = "lm", se=F)+
xlim(0, 10)+
ylim(0, 14)+
theme_bw()
Ahora bien, puede que entre la relación de ingesta de comida con la perdida de peso exista la presencia de otra variable que pueda estar confundiendo, y que no esta considerada.
En este caso, la variable que puede estar confundiendo refiere al ejercicio. En el siguiente gráfico se puede observar nítidamente esta relación.
data$exercise1<- as.factor(data$exercise) #recodificamos a factor para plotear
#Relacion lineal condicionada por cantidad de ejercicio.
ggplot(data, aes(x=food, y=wtloss,shape=exercise1,color=exercise1)) +
geom_point() +
geom_smooth(method = "lm", se=F)+
xlim(0, 10)+
ylim(0, 14)+
theme_bw()
Al observar el anterior gráfico se puede constatar que si se toma en consideración las horas promedio semanal de ejercicio (en este caso, 0,2 y 4) la relación es negativa entre comida y perdida de peso. Es decir si estudiamos la relación entre comida y perdida de peso entre sujetos “comparables” en relación a su ejercicio, la relación se vuelve negativa.
3.4 Estimación de modelos
Este aspecto claramente tiene repercusiones al estimar un modelo de regresión lineal como se puede observar en el siguiente cuadro.
#Analisis de regresión.
m1<- lm(wtloss~food, data=data)
m2<- lm(wtloss~exercise, data=data)
m3<- lm(wtloss~exercise+food, data=data)
stargazer(m1,m2,m3, type = "text")##
## =================================================================================
## Dependent variable:
## -------------------------------------------------------------
## wtloss
## (1) (2) (3)
## ---------------------------------------------------------------------------------
## food 0.071 -0.500*
## (0.541) (0.252)
##
## exercise 1.750*** 2.000***
## (0.361) (0.333)
##
## Constant 7.143** 4.000*** 6.000***
## (2.923) (0.913) (1.275)
##
## ---------------------------------------------------------------------------------
## Observations 10 10 10
## R2 0.002 0.746 0.838
## Adjusted R2 -0.123 0.714 0.791
## Residual Std. Error 3.505 (df = 8) 1.768 (df = 8) 1.512 (df = 7)
## F Statistic 0.017 (df = 1; 8) 23.520*** (df = 1; 8) 18.047*** (df = 2; 7)
## =================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Cómo se puede observar en la tabla anterior la comida (food), tiene un beta estimado de 0.071, lo que expresa una relación positiva y relativamente baja. No obstante cuando al modelo se le introduce la variable ejercicio (exercise), en el modelo 3, no solamente observamos un notorio aumento de la magnitud, sino que además este coeficiente cambia su dirección. Ahora es negativo.
El anterior modelo se puede formalizar de la siguiente manera:
\[\hat{Y}= b_{0}+b_{1}X_{1} + b_{2}X_{2}\]
Lo que para consideraciones de este analísis corresponde a:
\[\hat{Peso}= b_{0}+b_{1}Ejercicio + b_{2}Comida\] y A partir de lo estimado obtenemos que:
\[\hat{Y}= 6+2X_{1} - 0.5X_{2}\] A partir de la anterior ecuación podemos estimar una predicción de perdida de peso para cada una de las siguientes personas “téoricas”
## [1] 6## [1] 4## [1] 6A modo de resumen, el modelo que representa a \(\hat{Y}\) expresa una función lineal entre las variables \(X_{1}\) y \(X_{2}\)
3.5 Coeficientes de Regresión parcial.
Los coeficientes de regresión parcial \(b_{1}\) y \(b_{2}\) son conocidos también como o pendientes de regresión parcial. Estos coeficientes cuantifican la relación entre Y, y cada variable independiente en términos constantes.
A continuación demostraremos qué significa mantener una variable constante matemáticamente a partir del proceso de parcializar una variable en otra.
Si X1 y X2 no estan corelacionadas, ninguna parcialización es requerido para estimar b1 y b2. Comunmente los regresores en el modelo están correlacionados en algún grado, sin embargo algunos regresores pueden estar más correlacionados que otros en modelos con varias regresores. Dos coeficientes de regresión para X1 y X2 pueden seguir estimandose al regresar Y en X1 y en X2 separadamente pero solo si X1 y X2 han sido parcializadas en sus relaciones mutuas. Esto implica construir una nueva medición de X1 y X2 que son independientes de otro.
Para hacer esto consideramos un modelo en el cual X2 es predicho por X1, o en este ejemplo, la ingesta de comida es predicha por el ejercicio.
En primer lugar se presenta el gráfico donde se muestra la relación entre comida y ejercicio.
ggplot(data, aes(x=exercise, y=food)) +
geom_point() +
geom_smooth(method = "lm", se=F)+
xlim(0, 5)+
ylim(0, 10)+
theme_bw()
A continuación se presenta la estimación del modelo de regresión lineal:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.0 1.0954451 3.651484 0.006480453
## exercise 0.5 0.4330127 1.154701 0.281536920Al realizar la estimación de un regresor sobre el otro tenemos la siguiente ecuación
\[\hat{X_{2}}=4+0.5X_{1}\]
Este modelo de X2 genera una estimación de X2 dada la información que provee X1 para cada caso. A partir del comando fitted.values() podemos obtener los valores estimados(predichos) de X2. Además con el comando residuals() obtenemos una variable con la información de los residuos de esta regresión.
#Estimamos la regresion exercise~food
reg_x1_x2 <- lm(exercise~food, data=data)
data$x1_2 <- fitted.values(reg_x1_x2) #variable de exercise estimada
data$resx1_2<- residuals(reg_x1_x2) #residuos de exercise estimada
#Estimamos la regresion food~exercise
reg_x2_x1 <- lm(food~exercise, data=data)
data$x2_1 <- fitted.values(reg_x2_x1) #variable de food estimada
data$resx2_1<- residuals(reg_x2_x1) #residuos de food estimada
head(data)## # A tibble: 6 x 9
## id exercise food wtloss exercise1 x1_2 resx1_2 x2_1 resx2_1
## <int> <int> <int> <int> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 1 0 2 6 0 1.14 -1.14 4.00 -2.00e+ 0
## 2 2 0 4 2 0 1.71 -1.71 4 5.55e-17
## 3 3 0 6 4 0 2.29 -2.29 4 2.00e+ 0
## 4 4 2 2 8 2 1.14 0.857 5 -3.00e+ 0
## 5 5 2 4 9 2 1.71 0.286 5 -1.00e+ 0
## 6 6 2 6 8 2 2.29 -0.286 5 1.00e+ 0Los residuos del modelo de regresión lineal no estan correlacionados con todos los regresores del modelo. A continuación lo podemos ver, donde la correlación entre ejercicio no esta correlacionada con el residuo de la estimación de comida predicha por ejercicio.
## [1] 3.022382e-17De esta manera se puede entender que este residuo es el componente de la comida que es independiente del ejercicio. Es decir el residuo de esta regresión cuantifica lo unico de la variable comida, que no puede ser explicado por la variable ejercicio.
Ahora consideramos un modelo de regresión lineal en donde se estima la perdida de peso en base al residuo de la estimación de comida (en base a ejercicio). A continuación se presenta el grafíco de dispersión.
ggplot(data, aes(x=resx2_1, y=wtloss)) +
geom_point() +
geom_smooth(method = "lm", se=F)+
xlim(-4, 4)+
ylim(0, 14)+
theme_bw()
##
## Call:
## lm(formula = wtloss ~ resx2_1, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.500 -2.125 0.000 2.125 5.500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.5000 1.0577 7.091 0.000103 ***
## resx2_1 -0.5000 0.5575 -0.897 0.395949
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.345 on 8 degrees of freedom
## Multiple R-squared: 0.09137, Adjusted R-squared: -0.02221
## F-statistic: 0.8045 on 1 and 8 DF, p-value: 0.3959Como se puede apreciar la ecuación de mínimos cuadrados estimada para este modelo corresponde a: \[\hat{Y}=7.5-0.5ResComida\]
Como se puede notar, este coeficiente de regresión es el mismo que el coeficiente de comida cuando se controla por ejercicio para predecir perdida de perso.
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.14285714 2.9225839 2.4440213 0.04031087
## food 0.07142857 0.5408484 0.1320676 0.89819199A partir de esto, se puede establecer que controlando por la frecuencia del ejercicio, o dando cuenta de las diferencias entre la frecuencias por ejercicio, o manteniendo constante el ejercicio, dos personas que difieren en 1 unidad de ingesta de comida, es estimado que se diferenciaran por 0.5 (gramos) en perdida de peso.