Regresión lineal múltiple
Luis Fernando Apáez álvarez
17 de Septiembre de 2022
Ejemplo Ancova: Parte II
Recordemos rápidamente que:
# Cargamos los datos
CADdata <- read.table("cadmium.txt", header=TRUE, sep=" ", dec=".")
# Convertimos la variable group a tipo factor
CADdata$group <- factor(CADdata$group, levels=c(1, 2, 3), labels=c("High", "Low", "No") )
# Filtramos los datos por edad. Consideraremos a la
# persona mas joven del grupo de exposicion alta
CADdata <- CADdata[CADdata$age >= 39, ]vemos un gráfico con los datos ya filtrados
par(mfrow=c(1,2)); par(mar=c(4,4,2,1.5))
plot(CADdata$age, CADdata$vitcap, col=c("green","red","blue")[CADdata$group], xlab = "Age",
ylab = "Vital Capacity (L)", xlim=c(35,68), ylim=c(2.5,5.6),
pch=c(16,16,16), main = "Exposure to cadmium")
legend(55,5.8, levels(CADdata$group),
col=c("green","red","blue"), pch = c(16,16,16) )
boxplot(vitcap ~ group, data = CADdata, col = "white", outline=FALSE)
stripchart(vitcap ~ group, data = CADdata,
method = "jitter",
pch = 19,
col = c("green","red","blue"),
vertical = TRUE,
add = TRUE)
Depués, buscaremos ajustar el modelo saturado (el que considera todas las variables), el cual nos brindará información para dirigirnos a obtener un mejor modelo. Modelo completo:
\[ E(y;x)= \beta_0 + \beta_1\cdot age + \beta_2 \cdot Low + \beta_3 \cdot No + \beta_4(age\cdot Low) + \beta_5(age\cdot No) \] Ajustamos
fit <- lm(vitcap ~ age * group, data = CADdata)
summary(fit)##
## Call:
## lm(formula = vitcap ~ age * group, data = CADdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1017 -0.3456 -0.0324 0.5152 1.0770
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.1834 1.0641 7.69 1.5e-09 ***
## age -0.0851 0.0211 -4.04 0.00022 ***
## groupLow 1.2369 2.0332 0.61 0.54622
## groupNo -3.3940 1.3446 -2.52 0.01547 *
## age:groupLow -0.0273 0.0426 -0.64 0.52586
## age:groupNo 0.0722 0.0269 2.69 0.01031 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.64 on 42 degrees of freedom
## Multiple R-squared: 0.394, Adjusted R-squared: 0.321
## F-statistic: 5.45 on 5 and 42 DF, p-value: 0.000592donde tenemos las esperanzas
\[\begin{align*} &E(Y;group=High, age)= \beta_0 + \beta_1\cdot age\\ &E(Y;group= Low, age)= \beta_0 + \beta_1 age +\beta_2 + \beta_4 age = (\beta_0 + \beta_2) + (\beta_1 + \beta_4) age\\ &E(Y;group= No, age)= \beta_0 + \beta_1 age + \beta_3 + \beta_5 age = (\beta_0 + \beta_3) + (\beta_1 + \beta_5) age \end{align*}\]
Este modelo considera tres diferentes rectas, una para cada grupo Se rechaza \(H_0\) en la prueba asociada a la tabla ANOVA (p-value: 0.000592), entonces el modelo tiene sentido. Es decir, la edad o el nivel de exposición ayudan a modelar \(E(Y)\).
Lo que debería seguir inmediatamente después de lo anterior es la verificación de supuestos, lo cual haremos más adelante.
Continuación
Seguimos con prueba de igualdad de pendientes (coeficientes asociados a las interacciones). Si no se rechaza, podríamos optar por un modelo con igualdad de pendientes o rectas paralelas (facilita la interpretación), donde la prueba a realizar contrastar es:
\[ H_0: \beta_4=0\ \ y \ \ \beta_5=0\ \ vs \ \ H_a: \beta_4\neq0 \ \ o\ \ \beta_5\neq0 \] para lo cual utilizaremos la librería multcomp
library(multcomp)
K <- matrix(c(0,0,0,0,1,0,
0,0,0,0,0,1), ncol=6, nrow=2, byrow=TRUE)
m <- c(0,0)
# via la prueba F
summary(glht(fit, linfct=K, rhs=m), test=Ftest())##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 1 == 0 -0.0273
## 2 == 0 0.0722
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 5.27 2 42 0.0091donde se rechaza \(H_{0}\) (p-value: 0.0091). No se puede considerar el modelo con rectas paralelas. Procedemos a ver si las tres rectas tienen pendiente diferente (intentar reducir el modelo), es decir, es la prueba simultánea que sigue a la prueba lineal general.
Notar que aquí se incluyen las tres hipótesis individuales asociadas a
- H0_1: pendientes de High y Low iguales (b4=0)
- H0_2: pendientes de High y No iguales (b5=0)
- H0_3: pendientes de Low y No iguales (b4-b5=0)
Esta prueba sólo detecta las diferencias más evidentes, la lectura se debe hacer con cuidado considerando que esta prueba tiene más chance de cometer el error tipo II
# ¿Todas las pendientes difieren?
K <- matrix(c(0,0,0,0,1,0,
0,0,0,0,0,1,
0,0,0,0,1,-1), ncol=6, nrow=3, byrow=TRUE)
m <- c(0,0,0)
# via pruebas simultaneas
summary(glht(fit, linfct=K, rhs=m))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = vitcap ~ age * group, data = CADdata)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 -0.0273 0.0426 -0.64 0.795
## 2 == 0 0.0722 0.0269 2.69 0.026 *
## 3 == 0 -0.0995 0.0406 -2.45 0.046 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)Recordemos que \(H_{0}\) será rechazado si todos los p-values apuntan a rechazar la \(H_{0}\) correspondiente a cada prueba de la que hablamos antes, pues estamos realizando la “intersección” de pruebas de hipótesis.
De forma simultánea podemos identificar que se rechaza \(H_0\) por dos diferencias:
Se rechaza \(\beta_5=0\), pues \(p-value < 0.026\). Hay evidencia sobre pendientes de High y No diferentes.
Se rechaza \(\beta_4-\beta:5=0\), pues \(p-value < 0.046\). Hay evidencia sobre pendientes de Low y No diferentes
Con esta prueba no se rechaza \(\beta_4=0\), por lo que podríamos optar por considerar un modelo con \(\beta_4=0\). Ajustemos el modelo reducido que sólo considera a (age*No), este nuevo modelo corresponde a
\[ E(y;x)= \beta_0 + \beta_1 age + \beta_2 Low + \beta_3 No + \beta_4(age\cdot No) \] por lo cual realizamos el ajuste del modelo anterior como sigue
fitred <- lm(vitcap ~ age + group + I(age*(group=="No")), data = CADdata)
summary(fitred)##
## Call:
## lm(formula = vitcap ~ age + group + I(age * (group == "No")),
## data = CADdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.1738 -0.3758 -0.0677 0.5384 1.0720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.5149 0.9230 9.22 9.4e-12 ***
## age -0.0918 0.0182 -5.05 8.7e-06 ***
## groupLow -0.0524 0.2647 -0.20 0.8439
## groupNo -3.7255 1.2323 -3.02 0.0042 **
## I(age * (group == "No")) 0.0789 0.0246 3.20 0.0026 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.63 on 43 degrees of freedom
## Multiple R-squared: 0.388, Adjusted R-squared: 0.331
## F-statistic: 6.81 on 4 and 43 DF, p-value: 0.000245donde tenemos las esperanzas:
\[\begin{align*} &E(Y;group=High, age)= \beta_0 + \beta_1\cdot age\\ &E(Y;group= Low, age)= \beta_0 + \beta_1 age +\beta_2 = (\beta_0 + \beta_2) + \beta_1+ age\\ &E(Y;group= No, age)= \beta_0 + \beta_1 age + \beta_3 + \beta_4 age = (\beta_0 + \beta_3) + (\beta_1 + \beta_4) age \end{align*}\]
y donde se rechaza \(H_0\) en la prueba asociada a la tabla ANOVA (p-value: 0.000245). Ahora podríamos proceder a reducir un poco más el modelo para tratar de facilitar la interpretación o bien a contestar las preguntas sobre los investigadores.
De la salida de fitred, se puede observar que una opción es considerar \(\beta_2=0\).
El modelo reducido quedaría como
\[ E(y;x)= \beta_0 + \beta_1 age + \beta_2 No + \beta_3(age\cdot No) \]
donde:
\[\begin{align*} &E(Y;group=High, age)= \beta_0 + \beta_1 age\\ &E(Y;group= Low, age)= \beta_0 + \beta_1 age \\ &E(Y;group= No, age)= \beta_0 + \beta_1 age + \beta_2 + \beta_3 age = (\beta_0 + \beta_2) + (\beta_1 + \beta_3) age \end{align*}\]
Con esto se tendrían rectas iguales para High y Low por otro lado, se tiene otra recta para No. Ajustamos otro modelo para el nuevo modelo reducido:
fitred2 <- lm(vitcap ~ age + I(group=="No") + I(age*(group=="No")), data = CADdata)
summary(fitred2)##
## Call:
## lm(formula = vitcap ~ age + I(group == "No") + I(age * (group ==
## "No")), data = CADdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.2015 -0.3863 -0.0677 0.5428 1.0975
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.4499 0.8534 9.90 9.1e-13 ***
## age -0.0910 0.0175 -5.19 5.2e-06 ***
## I(group == "No")TRUE -3.6605 1.1748 -3.12 0.0032 **
## I(age * (group == "No")) 0.0781 0.0240 3.25 0.0022 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.62 on 44 degrees of freedom
## Multiple R-squared: 0.387, Adjusted R-squared: 0.345
## F-statistic: 9.27 on 3 and 44 DF, p-value: 7.26e-05Se rechaza \(H_0\) en la prueba asociada a la tabla ANOVA (\(p-value: 7.26e-05\)). También se puede observar que no se podría reducir más el modelo.
Las rectas ajustadas son:
\[\begin{align*} &\hat{E}(Y;group=High, age)=\hat{E}(Y;group=Low, age)= 8.4499 -0.0910 age\\ &\hat{E}(Y;group=No, age)= (8.4499-3.6605)+(-0.0910+0.0781) age =4.789-0.0129 age \end{align*}\]
Se puede observar que la pendiente de los expuestos a cadmio es mayor en valor absoluto y negativa con lo que podemos interpretar que la capacidad vital decrece más rápido en ese grupo comparado con los no expuesto.
Esto se puede observar más fácilmente en una gráfica
# Modelo fitred2
# Definimos dos funciones apra graficar las rectas:
# * Recta para los que se expusieron altamente y bajamente
# al cadmio
fitredHyL <- function(X2) {fitred2$coefficients[1]+ fitred2$coefficients[2] * X2}
# * Funcion para graficar la recta para los que no se
# expusieron al cadmio
fitredN <- function(X2) {fitred2$coefficients[1]+fitred2$coefficients[3]+ (fitred2$coefficients[2]+fitred2$coefficients[4]) * X2}
# Scatterplot age vs vitcap
with(CADdata, plot(age, vitcap, col=c("red", "green", "blue")[CADdata[ , 1]] , pch = c(16,16,16)))
legend("topright",levels(CADdata[,1]), col=c("red", "green", "blue"), pch = c(16,16,16), pt.cex=1.5,cex = .9, y.intersp = 1.4 , bty="n" )
# agregamos las rectas al grafico anterior
curve(fitredHyL, from = min(CADdata$age), to = max(CADdata$age),
col = "red", add = T)
curve(fitredN, from = min(CADdata$age), to = max(CADdata$age),
col = "blue", add = T)
Lo cual nos da la información de interpretación necesaria para poder decir que la exposición al cadmio si determina el nivel de disminución en la capacidad respiratoria