library(readxl)
algas=read_excel("datos/algae.xlsx")
library(ggplot2)
ggplot(algas,aes(x=log(PO4),y=log(Chla),col=season)) + geom_point()+geom_smooth(method="lm")

algas$season=factor(algas$season)
modelFull=lm(log(Chla)~log(PO4)*season,data=algas)
modelParallel=lm(log(Chla)~log(PO4)+season,data=algas)
modelSimple=lm(log(Chla)~log(PO4),data=algas)
anova(modelFull,modelParallel)
## Analysis of Variance Table
## 
## Model 1: log(Chla) ~ log(PO4) * season
## Model 2: log(Chla) ~ log(PO4) + season
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    179 233.04                           
## 2    182 236.33 -3   -3.2947 0.8436 0.4717
anova(modelParallel,modelSimple)
## Analysis of Variance Table
## 
## Model 1: log(Chla) ~ log(PO4) + season
## Model 2: log(Chla) ~ log(PO4)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)  
## 1    182 236.33                             
## 2    185 248.24 -3   -11.906 3.0563 0.0297 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nos quedamos con el modelo de lineas paralelas:

summary(modelParallel)
## 
## Call:
## lm(formula = log(Chla) ~ log(PO4) + season, data = algas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.75374 -0.90381 -0.05106  0.79809  3.12540 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -1.1997     0.3672  -3.267   0.0013 ** 
## log(PO4)       0.7469     0.0715  10.446   <2e-16 ***
## seasonspring  -0.7137     0.2493  -2.862   0.0047 ** 
## seasonsummer  -0.5531     0.2574  -2.149   0.0330 *  
## seasonwinter  -0.3282     0.2420  -1.356   0.1768    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.14 on 182 degrees of freedom
##   (13 observations deleted due to missingness)
## Multiple R-squared:  0.3919, Adjusted R-squared:  0.3785 
## F-statistic: 29.32 on 4 and 182 DF,  p-value: < 2.2e-16

El modelo estima la distancia de cada recta a la recta ausente (autumn), y sólo nos indica que la ordenada en invierno coincide con la ordenada en otoño, y que a su vez difiere de primavera y verano, pero nos queda la duda ¿serán iguales las ordenadas de primavera y verano?.

Estamos nuevamente ante un problema de comparaciones múltiples.

library(multcomp)
posthoc<-glht(modelParallel,linfct=mcp(season="Tukey"))
summary(posthoc)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = log(Chla) ~ log(PO4) + season, data = algas)
## 
## Linear Hypotheses:
##                      Estimate Std. Error t value Pr(>|t|)  
## spring - autumn == 0  -0.7137     0.2493  -2.862   0.0239 *
## summer - autumn == 0  -0.5531     0.2574  -2.149   0.1412  
## winter - autumn == 0  -0.3282     0.2420  -1.356   0.5276  
## summer - spring == 0   0.1606     0.2373   0.677   0.9055  
## winter - spring == 0   0.3854     0.2199   1.753   0.2988  
## winter - summer == 0   0.2249     0.2296   0.979   0.7609  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(posthoc)
## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = log(Chla) ~ log(PO4) + season, data = algas)
## 
## Quantile = 2.5929
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                      Estimate lwr      upr     
## spring - autumn == 0 -0.71368 -1.36015 -0.06721
## summer - autumn == 0 -0.55310 -1.22059  0.11439
## winter - autumn == 0 -0.32824 -0.95586  0.29938
## summer - spring == 0  0.16058 -0.45471  0.77587
## winter - spring == 0  0.38544 -0.18476  0.95564
## winter - summer == 0  0.22486 -0.37060  0.82031

Datos simulados

Se muestra a continuación una situación algo más compleja:

modelFull=lm(y~x*grupo,datos)
modelParallel=lm(y~x+grupo,datos)
modelSimple=lm(y~x,datos)
anova(modelFull,modelParallel)
## Analysis of Variance Table
## 
## Model 1: y ~ x * grupo
## Model 2: y ~ x + grupo
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    392 211.40                                  
## 2    395 223.32 -3   -11.926 7.3715 8.128e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(modelFull,modelSimple)
## Analysis of Variance Table
## 
## Model 1: y ~ x * grupo
## Model 2: y ~ x
##   Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
## 1    392 211.4                                  
## 2    398 536.5 -6   -325.11 100.48 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El mejor modelo es el que requiere pendientes y ordenadas distintas:

summary(modelFull)
## 
## Call:
## lm(formula = y ~ x * grupo, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.71134 -0.38016  0.00098  0.40750  2.25627 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.01526    0.07421   0.206 0.837230    
## x            0.72918    0.08353   8.729  < 2e-16 ***
## grupo2       1.99687    0.10441  19.126  < 2e-16 ***
## grupo3      -0.01468    0.10544  -0.139 0.889370    
## grupo4       1.49667    0.10497  14.258  < 2e-16 ***
## x:grupo2     0.02885    0.11303   0.255 0.798677    
## x:grupo3    -0.29563    0.11512  -2.568 0.010597 *  
## x:grupo4    -0.37381    0.10851  -3.445 0.000633 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7344 on 392 degrees of freedom
## Multiple R-squared:  0.6972, Adjusted R-squared:  0.6918 
## F-statistic: 128.9 on 7 and 392 DF,  p-value: < 2.2e-16

El resumen anterior nos informa de que la ordenada de la recta ajustada a los datos del grupo 3 no se diferencia significativamente de la ordenada de la recta ajustada a los datos del grupo 1 (lo cual es obvio pues ambas rectas se cruzan en el origen de coordenadas).

Asimismo el p-valor asociado a la pendiente de la recta ajustada al grupo 2 nos informa de que dicha pendiente no difiere significativamente de la pendiente ajustada al grupo 1; en otras palabras, las rectas 1 y 2 son paralelas; la tabla anterior nos indica además que la pendiente de la recta 1 difiere significativamente de las pendientes de las rectas 3 y 4. Ahora bien, el gráfico sugiere que las rectas 3 y 4 pueden ser paralelas. ¿Como contrastar si \(\beta_3=\beta_4\)

Los contrastes de Tukey ahora no son adecuados, debido a la presencia de interacciones:

posthoc<-glht(modelFull,linfct=mcp(grupo="Tukey"))
## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found
## -- default contrast might be inappropriate
summary(posthoc)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = y ~ x * grupo, data = datos)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)    
## 2 - 1 == 0  1.99687    0.10441  19.126   <1e-04 ***
## 3 - 1 == 0 -0.01468    0.10544  -0.139    0.999    
## 4 - 1 == 0  1.49667    0.10497  14.258   <1e-04 ***
## 3 - 2 == 0 -2.01155    0.10490 -19.175   <1e-04 ***
## 4 - 2 == 0 -0.50020    0.10443  -4.790   <1e-04 ***
## 4 - 3 == 0  1.51134    0.10546  14.331   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Comparaciones múltiples en el modelo lineal

Ver desarrollo teórico en (http://www.ievbras.ru/ecostat/Kiril/R/Biblio/R_eng/Bretz%20Multiple%20Comparisons.pdf) páginas 41-45.

Podemos expresar los contrastes de interés en forma matricial; así el contraste \(\beta_3=\beta_4\) sería de la forma:

K=matrix(c(0,0,0,0,0,0,1,-1),nrow=1)
colnames(K)=names(coef(modelFull))
K
##      (Intercept) x grupo2 grupo3 grupo4 x:grupo2 x:grupo3 x:grupo4
## [1,]           0 0      0      0      0        0        1       -1
summary(glht(modelFull, linfct = K))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ x * grupo, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)
## 1 == 0  0.07818    0.10523   0.743    0.458
## (Adjusted p values reported -- single-step method)

Podemos contrastar simultáneamente la igualdad de pendientes \(\beta_3=\beta_4\) y \(\beta1=\beta2\). Al mismo tiempo podemos calcular las distancias de la recta 3 a la 4 y de la 1 a la 2 (dado que estas rectas son paralelas, estas distancias son iguales a las diferencias entre sus ordenadas); lo único que tenemos que hacer es elegir bien las posiciones de los coeficientes a comparar en la matriz K:

K=rbind(c(0,0,0,0,0,0,1,-1),  # Diferencias en pendientes entre rectas 3 y 4
        c(0,0,0,0,0,1,0,0),  # Diferencia en pendiente entre rectas 1 y 2
        c(0,0,1,0,0,0,0,0),  # Distancia entre las rectas 1 y 2
        c(0,0,0,1,-1,0,0,0))  # Distancia entre las rectas 3 y 4
colnames(K)=names(coef(modelFull))
K
##      (Intercept) x grupo2 grupo3 grupo4 x:grupo2 x:grupo3 x:grupo4
## [1,]           0 0      0      0      0        0        1       -1
## [2,]           0 0      0      0      0        1        0        0
## [3,]           0 0      1      0      0        0        0        0
## [4,]           0 0      0      1     -1        0        0        0
diferencias=glht(modelFull, linfct = K)
summary(diferencias)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ x * grupo, data = datos)
## 
## Linear Hypotheses:
##        Estimate Std. Error t value Pr(>|t|)    
## 1 == 0  0.07818    0.10523   0.743    0.913    
## 2 == 0  0.02885    0.11303   0.255    0.998    
## 3 == 0  1.99687    0.10441  19.126   <1e-05 ***
## 4 == 0 -1.51134    0.10546 -14.331   <1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(diferencias)
## 
##   Simultaneous Confidence Intervals
## 
## Fit: lm(formula = y ~ x * grupo, data = datos)
## 
## Quantile = 2.5017
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##        Estimate lwr      upr     
## 1 == 0  0.07818 -0.18506  0.34141
## 2 == 0  0.02885 -0.25391  0.31160
## 3 == 0  1.99687  1.73568  2.25806
## 4 == 0 -1.51134 -1.77517 -1.24751

En lugar de definir los contrastes de interés especificando la matriz K, podemos especificar también los nombres de los parámetros a contrastar (mediante el nombre que reciben los coeficientes correspondientes en el modelo):

diferencias=glht(modelFull, linfct = c("x:grupo3-x:grupo4 = 0", "x:grupo2 = 0",
                                   " grupo2 = 0","grupo3 -grupo4 = 0"))
summary(diferencias)
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Fit: lm(formula = y ~ x * grupo, data = datos)
## 
## Linear Hypotheses:
##                          Estimate Std. Error t value Pr(>|t|)    
## x:grupo3 - x:grupo4 == 0  0.07818    0.10523   0.743    0.913    
## x:grupo2 == 0             0.02885    0.11303   0.255    0.998    
## grupo2 == 0               1.99687    0.10441  19.126   <1e-05 ***
## grupo3 - grupo4 == 0     -1.51134    0.10546 -14.331   <1e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
confint(diferencias)
## 
##   Simultaneous Confidence Intervals
## 
## Fit: lm(formula = y ~ x * grupo, data = datos)
## 
## Quantile = 2.5017
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##                          Estimate lwr      upr     
## x:grupo3 - x:grupo4 == 0  0.07818 -0.18507  0.34142
## x:grupo2 == 0             0.02885 -0.25391  0.31161
## grupo2 == 0               1.99687  1.73567  2.25807
## grupo3 - grupo4 == 0     -1.51134 -1.77518 -1.24751