Ejemplo Ancova: Parte I

Problema

Se desea analizar si la exposición a cadmio, un metal pesado, produce una alteración (disminución) en la capacidad respiratoria. Se tiene información de la capacidad vital (vitcap: es la cantidad máxima de aire que una persona puede expulsar de los pulmones tras una inhalación máxima) para un conjunto de 84 hombres que trabajan en la industria del cadmio. Cada individuo fue clasificado de acuerdo con su exposición a vapores de cadmio en uno de tres grupos (group): 1 - exposición por más de 10 años; 2 - exposición por menos de 10 años, y 3 - sin exposición.

También se espera que por diseño los hombres en el grupo de exposición por más de 10 años tengan una edad mayor al resto. Además se sabe que la capacidad respiratoria tiende a ser menor con la edad, por lo que se decide considerar a la edad (age) en el modelo. Con un análisis exploratorio se puede observar que no hay individuos del grupo con exposición por más de 10 años que tengan una edad menor a 39, por lo que limitaremos el estudio a sólo los individuos de edad 39 y más.

Recordemos que los problemas de tipo anova y ancova son utilizados para problemas del tipo de causalidad. Por ejemplo, para el problema que queremos abordar, será que si nosotros impactamos la exposición al cadmino, podríamos mejorar o modificar la capacidad vital. Para ello debemos estar seguros de que estamos comparando variables de la misma naturaleza y que no tenemos variables confusoras que puedan distorsionar posibles relaciones.

Además, en estos tipos de problemas nos centramos en las pruebas de hipótesis, pues para el caso de la \(R^{2}\) aún cuando no sea alta, en realidad nos interesa ver si hay alguna relación entre las variables involucradas, aunque dicha relación sea muy pequeña.

Así, comenzamos entonces por cargar los datos

CADdata <- read.table("cadmium.txt", header=TRUE, sep=" ", dec=".")
str(CADdata)
## 'data.frame':    84 obs. of  3 variables:
##  $ group : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ age   : int  39 40 41 41 45 49 52 47 61 65 ...
##  $ vitcap: num  4.62 5.29 5.52 3.71 4.02 5.09 2.7 4.31 2.7 3.03 ...

por la información que tenemos del problema notamos que la variable group debe ser convertida a factor

CADdata$group <- factor(CADdata$group, levels=c(1, 2, 3), labels=c("High", "Low", "No") )

str(CADdata)
## 'data.frame':    84 obs. of  3 variables:
##  $ group : Factor w/ 3 levels "High","Low","No": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age   : int  39 40 41 41 45 49 52 47 61 65 ...
##  $ vitcap: num  4.62 5.29 5.52 3.71 4.02 5.09 2.7 4.31 2.7 3.03 ...

Podemos observar una gráfica global que representan pares de variables

ggpairs(data=CADdata, title="Datos", aes(colour = group))

notamos en el scatterplot de age vs vitcap que la relación no parece ser lineal; vemos también que los grupos tienen diferentes medias. Tenenemos mayor frecuencia para el grupo de personas que no han tenido exposición al cadmio respecto al resto de grupos.

Asimismo, parace ser que la capacidad vital disminuye conforme se tiene más edad. Consideremos la gráfica que se encuentra en la esquina superior derecha, en la cual no está involucrada la variable edad. Notamos que la capacidad vital disminuye con las personas que han tenido una exposición alta al cadmio, pero vemos también un detalle: la capacidad vital de las personas que han tenido una poca exposición al cadmio es mayor a la capacidad vital de las personas que no han tenido exposición al cadmio, lo cual nos causa ruida en las posibles conclusiones, siendo así que es necesario analizar la información incluyendo la variable edad, pues también observamos que la edad tiene un efecto sobre la variable vitcap.

Para complementar veamos las siguientes gráficas:

# Numero de graficas| margenes de las graficas
par(mfrow=c(1,2)); par(mar=c(4,4,2,1.5))

# Boxplot vitcap vs group
boxplot(vitcap ~ group, data = CADdata, col = "white", outline=FALSE)
# Agregamos los puntos de los datos
stripchart(vitcap ~ group, data = CADdata,
           # agregamos ruido a los puntos
           method = "jitter",
           pch = 19,
           col = 2:4,
           vertical = TRUE,
           # Agregamos los puntos al boxplot
           add = TRUE)

# boxplot age vs group
boxplot(age ~ group, data = CADdata, col = "white", outline=FALSE)
# Agregamos los puntos de los datos
stripchart(age ~ group, data = CADdata,
           method = "jitter",
           pch = 19,
           col = 2:4,
           vertical = TRUE,
           add = TRUE)

vemos que hay una dispersión grande para el grupo que tiene una alta exposición en el gráfico de la izquierda.

Luego, consideremos el diagrama de dispersión que habíamos visto antes sobre age vs vitcap:

plot(CADdata$age, CADdata$vitcap, col=c("green","red","blue")[CADdata$group], xlab = "Age", 
     ylab = "Vital Capacity (L)", xlim=c(17,68), ylim=c(2.5,6.1), 
     pch=c(16,16,16), main = "Exposure to cadmium")

# Legendas para cada grupo
legend(60, 5.8, levels(CADdata$group), col=c("green","red","blue"), pch = c(16,16,16) )

Notemos que tenemos pocos datos de color verde (exposición alta). Si nos centramos en dichos datos verdes notamos cierto decrecimiento de la capacidad vital al aumentar la edad. Para los puntos rojos notamos que, al parecer, a partir de los 40 (+/-) la tendencia es a la baja. Podemos ver también del gráfico que parece ser que la pérdida de la capacidad respiratoria para las personas menores a los 40 (+/-) es muy lenta, pero llega una cierta edad en la cual la pérdida de la capacidad comienza a ser más rápida. Para el grupo en verde no tenemos datos para edades menores a los 40 (+/-), pero si notamos que se nota claramente la tendencia a la baja en la capacidad respiratoria para dicho grupo. De ahí que sea necesario realizar una estandarización en nuestros datos, para ello filtraremos primero los datos:

# Filtramos los datos por edad. Consideraremos a la
# persona mas joven del grupo de exposicion alta

summary(CADdata[CADdata$group=="High",])
##   group         age         vitcap   
##  High:12   Min.   :39   Min.   :2.7  
##  Low : 0   1st Qu.:41   1st Qu.:3.0  
##  No  : 0   Median :48   Median :3.9  
##            Mean   :50   Mean   :3.9  
##            3rd Qu.:58   3rd Qu.:4.7  
##            Max.   :65   Max.   :5.5
CADdata <- CADdata[CADdata$age >= 39, ]

lo cual nos permite tener grupos más homogéneos.

Volvemos a graficar pero 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)

Trabajo sobre el modelo

Importante analizar el metadato de niveles, pues el primero será el nivel de referencia en los modelos

levels(CADdata$group)
## [1] "High" "Low"  "No"

Aquí la referencia es High. Procedemos a realizar el ajuste del modelo con interacciones:

\[ 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) \]

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.000592

donde podemos ver 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*}\]

en las cuales tenemos tres rectas distintas, lo cual parece lógico de acuerdo a las gráficas que vimos antes sobre los datos filtrados.

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)\).