Ejemplo Ancova: Parte III

Cargamos los datos y realizamos el filtro de los datos

# 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, ]

Realizamos el ajuste de tres modelos

# Modelo saturado
fit <- lm(vitcap ~ age * group, data = CADdata) 

# Primer modelo reducido
fitred <- lm(vitcap ~ age + group + I(age*(group=="No")), data = CADdata) 

# Segundo modelo reducido
fitred2 <- lm(vitcap ~ age + I(group=="No") + I(age*(group=="No")), data = CADdata)

Análisis de los supuestos

Recordemos que, en la práctica, no debemos ser tan estrictos con los supuestos que debe cumplir nuestro modelo.

Ahora bien, realizaremos un análisis rápido de los supuestos que debe cumplir nuestro modelo, para lo cual nos apoyaremos de algunas gráficas que R base nos brinda:

par(mfrow = c(2,2), mgp = c(2,0.7,0), mar = c(3,3,1.5,1))
# linealidad
plot(fitred2, 1)

# homocedasticidad
plot(fitred2, 3)   

# normalidad
plot(fitred2, 2)   

# outliers
plot(fitred2, 5, cook.levels = c(4/(dim(CADdata)[1]-2), 0.5, 1.0))

En la primer gráfica donde trabajamos sobre los residuales, buscamos que la curva marcada de color rojo sea lo más parecida a la recta horizontal \(y=0\), aunque debemos ser cautelosos pues R exagera los comportamientos de los datos en dichas curvas rojas para que sea más sencillo la detección de patrones. Notamos que a partir del 3.5 (en la misma gráfica, en la de la izquierda de la primer fila) empiezan a haber más datos, siendo lo anterior la posible causa del porqué la curva roja tiene el comportamiento diferente alrededor del 3.5; la segunda gráfica de la primer fila es sobre la homocedasticidad, donde de nuevo vemos que a partir de un cierto valor comenzamos a tener más datos. Dentro de esta gráfica buscamos no hallar algún patrón obvio, el cual es representado vía la curva de rojo, esto es, de nuevo, nos gustaría tener lo más parecido a una recta horizontal.

Para la primer gráfica de la segunda fila, nos gustaría que los puntos siguieran el comportamiento de la recta. Finalmente para el último gráfico se analiza si tenemos algún valor atípico influyente.

# outliers
plot(fitred2, 5, cook.levels = c(4/(dim(CADdata)[1]-2), 0.5, 1.0))

Para el eje \(y\), esperamos no encontrar valores fuera del intervalo \((-2,2)\), pues tendríamos en dicho caso posibles valores influyentes; asimismo, para el eje \(x\) tenemos la distancia de los datos respecto a su valor promedio, de donde buscamos no hallar distancias grandes respecto a dicho valor promedio. Notamos que hay dos valores para los cuales deberíamos tener cuidado pues están muy alejados del resto de puntos.

De arriba a abajo, en esta gráfica, la línea curva punteada es más estrcita y nos dice que los valores que estén por encima de dicha curva son los valores más influyentes y que sí repercuten en nuestro modelo, de modo que el único valor a considerarse es el 82, para el cual se debería hacer un análisis detallado para ver si es posible quitarlo o no. Arriba de dicha curva punteada (la cual no se alcanza a ver en el gráfico) se encuentra una curva menos estricta, por lo cual el punto 82, siendo menos estrictos, no sería considerado como un valor tan influyente. En sentido general, parece ser que no tenemos problemas sobre los valores influyentes (recordemos que en la práctica debemos ser menos estrictos).

Completaremos los análisis de las gráficas anteriores mediante pruebas de hipótesis.

# Homocedasticidad
library(car)

#H0: la varianza es constante 
car::ncvTest(fitred2)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.041, Df = 1, p = 0.8

No se rechaza \(H_0\), no hay evidencia en contra de la homocedasticidad.

# Normalidad 
# Se basa en los residuales estandarizados o estudentizados
# H0: los datos provienen de una distribución normal
library(broom)
Datosfitred2=augment(fitred2)
shapiro.test(Datosfitred2$.std.resid)
## 
##  Shapiro-Wilk normality test
## 
## data:  Datosfitred2$.std.resid
## W = 1, p-value = 0.06
library(normtest)
normtest::jb.norm.test(Datosfitred2$.std.resid)
## 
##  Jarque-Bera test for normality
## 
## data:  Datosfitred2$.std.resid
## JB = 2, p-value = 0.2

No se rechaza \(H_0\) en ambas pruebas, no hay evidencia en contra de la normalidad.

# Linealidad
library(car)
residualPlots(fit)

##            Test stat Pr(>|Test stat|)
## age             0.80             0.43
## group                                
## Tukey test      0.16             0.88

Analizamos los residuales de pearson para cada variable en cuestión, de modo que no hay fuerte evidencia en contra de la linealidad. Además, en las pruebas no se está rechazando \(H_0\) (buscamos justamente no rechazar), con \(H_{0}\) indica que los residuales no dependen de la variable.

Uso del modelo

Con este modelo se pueden contestar preguntas puntuales, por ejemplo, si a los 50 años ya se nota en promedio una diferencia entre los que están expuestos a cadmio y los que no

  • \(H_0: E(Y;group= High(Low), age=50) = E(Y;group= No, age=50)\) vs
  • \(H_a: E(Y;group= High(Low), age=50)\neq E(Y;group= No, age=50)\)

Notar que estas hipótesis se pueden escribir en términos de los parámetros:

\[ H_0: \beta_0 + \beta_1 50 = (\beta_0 + \beta_2) + (\beta_1 + \beta_3) 50 \]

Se pasan los parámetros de un lado de la igualdad \(H_0: \beta_2+\beta_3 (50) = 0\).

library(multcomp)
K <- matrix(c(0,0,1,50), ncol=4, nrow=1, byrow=TRUE)
m <- c(0)
summary(glht(fitred2, linfct=K, rhs=m))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = vitcap ~ age + I(group == "No") + I(age * (group == 
##     "No")), data = CADdata)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)
## 1 == 0    0.245      0.185    1.32     0.19
## (Adjusted p values reported -- single-step method)

No se rechaza \(H_0\), es decir, a los 50 años no se encuentra evidencia que indique una diferencia promedio en la capacidad vital entre los que fueron expuestos y los que no.

A veces es importante incluir la dirección, es decir, los investigadores sospechan que a los 55 años la capacidad vital es en promedio mayor en los no expuestos comparados con los expuestos. Aquí tomamos información del contexto del problema con base en lo que sospechan los investigadores esto se expresa con una dirección en la forma de plantear las hipótesis

  • \(H_0: E(Y;group= High(Low), age=55) \geq E(Y;group= No, age=55)\) vs
  • \(H_a: E(Y;group= High(Low), age=55) < E(Y;group= No, age=55)\)

Notar que en la \(H_a\) sólo se permiten los casos \(\neq, <\) o \(>\). Esto se puede escribir en términos de los parámetros

\[ H_0: \beta_0 + \beta_1 55 \geq (\beta_0 + \beta_2) + (\beta_1 + \beta_3) 55 \]

Es decir, pasando los parámetros a un sólo lado de la desigualdad

\[ H_0: 0 \geq \beta_2 + \beta_3 55\ \ \ vs \ \ \ H_a: 0 < \beta_2 + \beta_3 55 \]

La alternativa se usa para definir la dirección \(< f(\beta_0,\beta_1,...,\beta_p)\) es greater, \(> f(\beta_0,\beta_1,...,\beta_p)\) es less.

K <- matrix(c(0,0,1,55), ncol=4, nrow=1, byrow=TRUE)
m <- c(0)
summary(glht(fitred2, linfct=K, rhs=m, alternative="greater"))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = vitcap ~ age + I(group == "No") + I(age * (group == 
##     "No")), data = CADdata)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>t)   
## 1 <= 0    0.635      0.242    2.63 0.0059 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Se rechaza \(H_0\), es decir, hay evidencia para concluir que el promedio de capacidad vital es menor en los expuestos comparados con los no expuestos a la edad de 55 años. La diferencia promedio estimada es 0.635.

Una forma de resumir los resultados es presentar una gráfica con los intervalos de confianza simultáneos. Con esto es posible identificar en que punto ya se podría observar una diferencia entre los promedios en la capacidad vital.

Hay dos opciones:

  • I. La más fácil de interpretar corresponde a incluir los datos las rectas ajustadas y los intervalos de confianza de \(E(Y)\). Dado que no considera ninguna información adicional (dirección) podría evitar distinguir claramente algo.

Calculares los intervalos de confianza simultáneos para la esperanza de la capacidad vital, para expuestos y no expuestos. Bajaremos la confianza a 90%, pues serán intervalos simultáneos de la esperanza y estos suelen ser muy conservadores

# Grafico base
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")

# Creamos una malla de valores asociados a la edad
# 39, 39.5, ... , 64.5, 65
age <- seq(from = 39, to = 65, by = .5)

# Calcularemos 53 intervalos de confianza
# length(age)

# Para la banda del grupo expuesto (high o low)
# E(Y;group=High (low), age)= b0 + b1 age
# Creamos una matriz de la forma
# 1 valor_edad 0 0
# 1 valor_edad 0 0
# etc
# correspondiente a los valores asociados a las betas
# 1*b0 b1*age  0*b2  0*b3
# 1 valor_edad   0   0
KH <- cbind(1, age, 0, 0)


# Para del grupo no expuesto
# E(Y;group=No, age)= (b0 + b2) + (b1 + b3) age
KN <- cbind(1, age, 1, age)

# Juntamos ambas matrices en la matriz K para utilizar la libreria multcomp
K <- rbind(KH, KN)

# Activamos glht sobre el modelo filtred2 con la matriz K sobre la cual nos interesa hacer las evaluaciones
fitE <- glht(fitred2, linfct = K)

# Creamos los intervalos de confianza con un nivel de confianza del 90% (se guardas los limites superior e inferior)
fitci <- confint(fitE, level = 0.90)

# Graficamos las bandas de confianza:

# Accedemos a los coeficientes del modelo pero solo para los expuestos, hasta la fila 53 y graficamos la recta
lines(age, coef(fitE)[1:53], col="red")

# Graficamos los intervalos de confianza hasta la fila 53
lines(age, fitci$confint[1:53,"upr"], col="red")
lines(age, fitci$confint[1:53,"lwr"], col="red")

# Accedemos a los coeficiente del modelo de los no expuestos, a partir de la fila 54 y graficamos la recta
lines(age, coef(fitE)[54:106], col="blue")

# Graficamos los intervalos de confianza a partir de la fila 54
lines(age, fitci$confint[54:106,"upr"], col="blue")
lines(age, fitci$confint[54:106,"lwr"], col="blue")

Estas bandas de confianza o intervalos de confianza simultáneos sirven para identificar las diferencias más evidentes con los datos, pero son muy conservadoras.

  1. Un poco más difícil de interpretar y corresponde a incluir el intervalo de confianza directamente de la diferencia de esperanzas para una edad fija entre expuestos y no expuestos, es decir

\[\begin{align*} E(Y;group=High (low), age)-E(Y;group=No, age)&= [\beta_0 + \beta_1 age] - [(\beta_0 + \beta_2) + (\beta_1 + \beta_3) age]\\ &= -\beta_2-\beta_3 age \end{align*}\]

Notar que ahora nos interesan aquellos valorescuyos intervalos estén por abajo de 0 pues con eso conluiremos que para esas edades se observan en promedio menor capacidad vital en los expuestos.

Creamos una malla de valores asociados a la edad

# Repetimos lo hecho anteriormente
age <- seq(from = 39, to = 65, by = .5)
K <- cbind(0, 0, -1, -age)

fitE <- glht(fitred2, linfct = K)
fitci <- confint(fitE, level = 0.90)

# Diferencia entre los expuestos y los no expuestos (es una recta)
plot(age, coef(fitE), col="black", type="l", main="Diferencia de E(vitcap) entre expuestos y no", ylab="E(Y;High-Low) - E(Y;No)")


lines(age, fitci$confint[,"upr"], col="black")
lines(age, fitci$confint[,"lwr"], col="black")

# Recta horizonal sobre 0
abline(h=0, col="blue")

A partir de los 54 años se observa que la esperanza de la capacidad vital es menor en las personas expuestas a cadmio en comparación con las no expuestas. Esta diferencia aumenta con la edad