Análisis de la covarianza.

Para este tema utilizaremos los datos del archivo coast.xlsx.

# Lectura de datos
library(readxl)
coast=read_excel("datos/coast.xlsx")
algae=read_excel("datos/algae.xlsx")

Asociación NO2-NH4

Para explicar la asociación entre estas dos variables, podemos representarlas gráficamente:

ggplot(coast,aes(x=NO2,y=NH4))+geom_point()+geom_smooth(method="lm")

Vemos aquí que el valor de NH4 crece con el valor de NO2 de manera lineal. La ecuación de esta recta es de la forma:

\[NH_4 = \beta_0 + \beta_1 NO_2\] que en R podemos estimar mediante:

lm(NH4~NO2,data=coast)
## 
## Call:
## lm(formula = NH4 ~ NO2, data = coast)
## 
## Coefficients:
## (Intercept)          NO2  
##      10.851        3.395

Podemos visualizar el posible efecto la estación del año en que se miden estas variables representando los puntos de color distinto según se haya realizado la medición en noviembre o en mayo:

ggplot(coast,aes(x=NO2,y=NH4, col=month))+geom_point()+geom_smooth(method="lm")

Este gráfico sugiere que la relación entre ambas variables queda mejor representada por dos rectas:

\[Noviembre:\;\;\; NH_4 = \beta_0 + \beta_1 NO_2\] \[Mayo:\;\;\; NH_4 = \beta'_0 + \beta'_1 NO_2\] que podemos estimar mediante:

lm(NH4~NO2,data=subset(coast,month=="may"))
## 
## Call:
## lm(formula = NH4 ~ NO2, data = subset(coast, month == "may"))
## 
## Coefficients:
## (Intercept)          NO2  
##      10.541        3.126
lm(NH4~NO2,data=subset(coast,month=="nov"))
## 
## Call:
## lm(formula = NH4 ~ NO2, data = subset(coast, month == "nov"))
## 
## Coefficients:
## (Intercept)          NO2  
##      12.357        3.344

Como podemos apreciar las pendientes son muy parecidas, y las ordenadas se diferencian en aproximadamente 2 unidades; ello significaría que para un mismo valor de \(NO_2\), el valor de NH4 en noviembre es dos unidades mayor que en junio.

Podemos preguntarnos si para representar la asociación entre estas dos variables es mejor el primer modelo simple, o el modelo que comprende dos rectas, una para cada mes. Esta pregunta se traduce en que debemos decidir si las diferencias entre los coeficientes estimados en cada una de las dos rectas son significativas (atribuíbles al azar) o no. Si las diferencias en pendiente u ordenada fuesen atribuíbles al azar, podríamos concluir que no hemos detectado diferencias entre noviembre y mayo, y basta con una única recta para describir la relación entre \(NO_2\) y \(NH_4\).

Para llevar a cabo este contraste, debemos integrar ambos modelos en una única ecuación; ello puede hacerse fácilmente si definimos la variable \(I_{may}\) como:

\[I_{may}=\begin{cases} 0 & \textrm{Si la medida se toma en noviembre}\\ 1 & \textrm{Si la medida se toma en mayo} \end{cases}\]

Si consideramos entonces el modelo:

\[NH_{4}=\beta_{0}+\gamma_{0}I_{may}+\beta_{1}NO_{2}+\gamma_{1}I_{may}\cdot NO_{2}\] tenemos que, para los datos medidos en noviembre, el valor de \(I_{may}\) es 0, por lo que el modelo se reduce a:

\[NH_{4}=\beta_{0}+\beta_{1}NO_{2}\] es decir, la primera de las dos rectas anteriores; y cuando los datos corresponden a mayo, entonces \(I_{may}=1\), y el modelo puede escribirse como:

\[NH_{4}=(\beta_{0}+\gamma_{0})+(\beta_{1}+\gamma_{1}) NO_{2}\] Considerando \(\beta'_0=\beta_0+\gamma_0\) y \(\beta'_1=\beta_1+\gamma_1\), esta ecuación corresponde a la segunda de las rectas anteriores.

Por tanto, contrastar si existen diferencias significativas entre las dos rectas es equivalente a contrastar si los coeficientes \(\gamma_0\) y \(\gamma_1\) son cero o no:

\[\begin{cases} H_{0}: & \gamma_{0}=\gamma_{1}=0\\ H_{1}: & \gamma_{0}\neq0\,\,\,\textrm{ó}\,\,\,\gamma_{1}\neq0 \end{cases}\]

De acuerdo a lo señalado en el capítulo Modelos lineales el contraste anterior puede expresarse como:

\[H_{0}:\left(\begin{array}{cccc} 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 \end{array}\right)\left(\begin{array}{c} \beta_{0}\\ \gamma_{0}\\ \beta_{1}\\ \gamma_{1} \end{array}\right)=\left(\begin{array}{c} 0\\ 0 \end{array}\right)\]

Llamando

\[A=\left(\begin{array}{cccc} 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 \end{array}\right)\]

\[\beta=\left(\begin{array}{c} \beta_{0}\\ \gamma_{0}\\ \beta_{1}\\ \gamma_{1} \end{array}\right)\]

El contraste se expresa como \(H_0: A\beta=0\). Cuando \(H_0\) es cierta el test estadístico:

\[F=\frac{1}{q\hat{\sigma}^{2}}\hat{\beta}^{t}\mathbf{A}^{t}\left(\mathbf{A}\left(\mathbf{X}^{t}\mathbf{X}\right)^{-1}\mathbf{A}^{t}\right)^{-1}\mathbf{A\hat{\beta}}\]

(siendo \(q\) el número de filas de la matriz \(\mathbf{A}\)) sigue una distribución \(F\) de Fisher con \(q\) y \(n-p-1\) grados de libertad (siendo \(p\) el número de variables explicativas). El p-valor del contraste es

\[p-valor=P\left(F_{q,n-p}>F_{obs}\right)=1-P\left(F_{q,n-p} \le F_{obs}\right)\]

siendo \(F_{obs}\) el valor observado del test estadístico \(F\). Un p-valor grande (>0.05) indica que puede aceptarse \(H_0\). Un p-valor pequeño indica que la muestra contiene información suficiente para rechazar dicha hipótesis.

Para el caso particular que nos ocupa, pdemos construir las matrices \(X\) e \(Y\):

coast$Imay=ifelse(coast$month=="may",1,0)
coast$NO2Imay=coast$Imay*coast$NO2
X=cbind(1,coast$Imay,coast$NO2,coast$NO2Imay)
Y=coast$NH4

Estimamos el modelo y su varianza residual:

beta=solve(t(X)%*%X)%*%t(X)%*%Y
beta
##            [,1]
## [1,] 12.3571765
## [2,] -1.8165053
## [3,]  3.3442094
## [4,] -0.2179961
p=ncol(X)-1
n=nrow(X)
yXB=Y-X%*%beta
sigma2=(t(yXB)%*%yXB)/(n-p-1)
sigma2
##          [,1]
## [1,] 30.22049

Por último, podemos construir el test estadístico F y calcular su p-valor de manera matricial:

A=matrix(rbind(c(0,1,0,0),c(0,0,0,1)),nrow=2)
Ab=A%*%beta
numF=t(Ab)%*%solve(A%*%solve(t(X)%*%X)%*%t(A))%*%Ab
q=nrow(A)
Fexp=numF/(q*sigma2)
Fexp
##          [,1]
## [1,] 35.73226
pvalor=1-pf(Fexp,q,n-p-1)
pvalor
##              [,1]
## [1,] 8.881784e-16

En R estos cálculos pueden hacerse de manera más sencilla estimando el modelo completo y el modelo sin los parámetros que se quiere contrastar:

modelFull=lm(NH4~NO2+Imay+NO2Imay,data=coast)
modelSimple=lm(NH4~NO2,data=coast)
anova(modelFull,modelSimple)
## Analysis of Variance Table
## 
## Model 1: NH4 ~ NO2 + Imay + NO2Imay
## Model 2: NH4 ~ NO2
##   Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
## 1   1236 37353                                  
## 2   1238 39512 -2   -2159.7 35.732 8.198e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El p-valor tan pequeño indica que se puede rechazar la hipótesis nula; por tanto \(\gamma_0\) ó \(\gamma_1\) son distintos de cero, esto es el modelo requiere la presencia de dos rectas.

Podemos decidir si las dos rectas tienen o no la misma pendiente. Ello es equivalente a contrastar \(H_0: \gamma_1=0\). Del mismo modo que antes, podemos construir ahora la matriz \(A=\left(\begin{array}{cccc} 0 & 0 & 0 & 1\end{array}\right)\), por lo que el contraste anterior puede escribirse como \(H_0: A\beta=0\)

A=matrix(c(0,0,0,1),nrow=1)
Ab=A%*%beta
numF=t(Ab)%*%solve(A%*%solve(t(X)%*%X)%*%t(A))%*%Ab
q=nrow(A)
Fexp=numF/(q*sigma2)
Fexp
##          [,1]
## [1,] 0.886524
pvalor=1-pf(Fexp,q,n-p-1)
pvalor
##          [,1]
## [1,] 0.346605

que en R puede resolverse como:

modelParallel=lm(NH4~NO2+Imay,data=coast)
anova(modelFull,modelParallel)
## Analysis of Variance Table
## 
## Model 1: NH4 ~ NO2 + Imay + NO2Imay
## Model 2: NH4 ~ NO2 + Imay
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1   1236 37353                           
## 2   1237 37379 -1   -26.791 0.8865 0.3466

El p-valor indica que puede aceptarse que \(\gamma_1\) es cero y por tanto que el mejor modelo es el formado por dos rectas paralelas con la misma pendiente:

summary(modelParallel)
## 
## Call:
## lm(formula = NH4 ~ NO2 + Imay, data = coast)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.2405  -3.7959   0.0095   3.6387  18.5533 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  12.7111     0.5116  24.845   <2e-16 ***
## NO2           3.2573     0.1134  28.734   <2e-16 ***
## Imay         -2.6511     0.3156  -8.401   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.497 on 1237 degrees of freedom
## Multiple R-squared:  0.4438, Adjusted R-squared:  0.4429 
## F-statistic: 493.5 on 2 and 1237 DF,  p-value: < 2.2e-16

La distancia estimada entre las rectas es el valor del coeficiente \(\gamma_0=-2.6511\).

Podemos obtener intervalos de confianza para estos parámetros mediante:

confint(modelParallel)
##                 2.5 %    97.5 %
## (Intercept) 11.707347 13.714797
## NO2          3.034888  3.479685
## Imay        -3.270178 -2.032024

 

En realidad, para llevar a cabo estas estimaciones con R no es necesario construir la variable \(I_{may}\). Los modelos anteriores puede construirse igualmente mediante:

modelFull=lm(NH4~NO2*month,data=coast)
modelParallel=lm(NH4~NO2+month,data=coast)
modelSimple=lm(NH4~NO2,data=coast)
anova(modelFull,modelSimple)
## Analysis of Variance Table
## 
## Model 1: NH4 ~ NO2 * month
## Model 2: NH4 ~ NO2
##   Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
## 1   1236 37353                                  
## 2   1238 39512 -2   -2159.7 35.732 8.198e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(modelFull,modelParallel)
## Analysis of Variance Table
## 
## Model 1: NH4 ~ NO2 * month
## Model 2: NH4 ~ NO2 + month
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1   1236 37353                           
## 2   1237 37379 -1   -26.791 0.8865 0.3466
anova(modelParallel,modelSimple)
## Analysis of Variance Table
## 
## Model 1: NH4 ~ NO2 + month
## Model 2: NH4 ~ NO2
##   Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
## 1   1237 37379                                  
## 2   1238 39512 -1   -2132.9 70.585 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 

 

Ejercicios.

Ajustar modelos de análisis de la covarianza a los siguientes casos:

ggplot(coast,aes(x=Temperature,y=NH4, col=month))+geom_point()+geom_smooth(method="lm")

ggplot(coast,aes(x=Hg,y=As, col=month))+geom_point()+geom_smooth(method="lm")

ggplot(coast,aes(x=DIN,y=COD))+geom_point()+geom_smooth(method=lm)

ggplot(coast,aes(x=DIN,y=COD,col=month))+geom_point()+
  geom_smooth(method=lm)

ggplot(coast,aes(x=Cu,y=Pb,col=month))+geom_point()+
  geom_smooth(method=lm)

algae1=subset(algae,size!="small")
ggplot(algae1,aes(x=mxPH,y=mnO2,col=size)) + geom_point() + geom_smooth(method=lm)