library(tibble)
library(dplyr)
library(ggplot2)

Referències

  • L. Cayuela & M. De la Cruz, Análisis de datos ecológicos en R, 2022.
  • A. García Pérez, Cuadernos de estadística aplicada: biología y ciencias ambientales, 2023.

Introducció

L’anàlisi de la variança (ANOVA) i de la covariança (ANCOVA) són casos específics dels models lineals on tenim variables categòriques.

En ANOVA totes les variables són categòriques, mentre que en ANCOVA hi ha variables tant categòriques com contínues. Ambdós models no difereixen massa dels altres models lineals, però presenten algunes peculiaritats que caldrà tenir en compte.


Anàlisi de la variança (ANOVA)

En ANOVA, la hipòtesi nul·la no serà que la pendent de la recta serà 0, car no hi ha recta. En aquest tipue de model el que farem serà comparar mitjanes \(\mu_j\) entre diferents grups (nivells del facto). En aquest cas donc la \(H_0\) serà la igualtat de les mitjanes dels diferents grups \[ \small H_0: \mu_a=\mu_b=...=\mu_k=\mu \] contra una \(H_1\) que no totes les mitjanes són iguals \[ \small H_1: \exists~ \mu_j\neq\mu~~~j=1,2,...,k \]

Exemple Chalcosoma atlas

Para demostra com s’ajusta ANOVA a les dades farem servir les dades d’un estudi alomètric sobre l’escarbat atles, Chalcosoma atlas L. En aquest cas concret s’estudia com canvia la grandària de l’aparat genital en mascles de diferents poblacions de l’espècies.

Observació ràpida de les dades

Las dades es troben al data frame chalco de la libreria ADER

library(ADER)

data("chalco")
str(chalco)
## 'data.frame':    156 obs. of  3 variables:
##  $ Localidad     : Factor w/ 4 levels "C. Malasia","E. Tailandia",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Long_cuerpo   : num  65.2 55 60.5 59.5 64.9 63.9 59.6 55.1 58 60.6 ...
##  $ Long_genitalia: num  12.2 12 12.6 12.3 12.6 12.3 12.1 12.2 12 12.3 ...

Veiem que les dades provenen d’una mostra de 156 individus i presenten 3 variables: contínues (Long_cuerpo i Long_genitalia) i una categòrica (Localidad) que divideix la mostra en 4 poblacions.

levels(chalco$Localidad)
## [1] "C. Malasia"   "E. Tailandia" "Nias"         "S. Sumatra"

En este caso posarem Población com a variable explicativa de Long_genitalia, que implica que la mitjana de dita variable és igual en les 4 poblacions. És a dir, tal com hem viast abans: \[ \small H_0: \mu_1=\mu_2=\mu_3=\mu_4=\mu \]

La \(H_1\), en canvi, serà \[ \small H_1: \exists~ \mu_j\neq\mu~~~j=1,2,3,4 \]

Anàlisi exploratòria de la relació entre ambdues variables

Això es pot fer molt bé amb un gràfic de caixes:

boxplot(Long_genitalia~Localidad, data = chalco)

El gràfic sembla al·ludir a certes diferències que en tot cas no són tan grans. Per tant per saber si existeix alguna diferenència significativa haurem de fer servir un test d’hipòtesi específic que compari les variacions entre grups i dins de cada grup generant els valors residuals:

\[ \small y_{ij}= \bar y + (\bar y_j-\bar y)+(y_{ij}-\bar y) \] ## Ajustament del model

Per al realitzar l’ajustament del model a les dades utilitzarem la funció lm(), tenit present que cal recordar de transformar la variable explicativa en categòrica, si no ho és. En aquest cas, no cal fer res, ya que Localidad es categòrica.

lm.chalco <- lm(Long_genitalia~Localidad, data = chalco)

Anàlisi de la variança

De manera anàloga el test s’analitza amb la funció anova():

anova(lm.chalco)
## Analysis of Variance Table
## 
## Response: Long_genitalia
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Localidad   3 10.506  3.5020  10.815 1.757e-06 ***
## Residuals 152 49.221  0.3238                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El p valor de la columna Pr(>F) de la taula es realment petit. Que confirma la exitència de almenys una diferència entre mitjanes.

Exploració dels coeficients estimats

summary(lm.chalco)
## 
## Call:
## lm(formula = Long_genitalia ~ Localidad, data = chalco)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.68136 -0.36201  0.06985  0.38591  1.01864 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           12.03529    0.09759 123.322  < 2e-16 ***
## LocalidadE. Tailandia  0.68971    0.13609   5.068 1.15e-06 ***
## LocalidadNias          0.24606    0.12253   2.008   0.0464 *  
## LocalidadS. Sumatra    0.02026    0.14669   0.138   0.8903    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5691 on 152 degrees of freedom
## Multiple R-squared:  0.1759, Adjusted R-squared:  0.1596 
## F-statistic: 10.81 on 3 and 152 DF,  p-value: 1.757e-06

Interpretació dels resultats

Quan la variable explicativa és categòrica el model hauria de quedar representat per una equació amb tants estimadors com nivells del factor i on cada estimador multiplica una variable que pot prendre valor 0 1. Aquest tipus de variable és una variable dummy o indicadora.

El model seria doncs \[ \small y_{ij}=\beta_1x_1+\beta_2x_2+\beta_3x_3+\beta_4x_4+\epsilon_i~~~ \epsilon \sim N(0,\sigma^2) \] Cada una d’aquestes noves variables indica la pertenència de cada idndividu a la variable, agafant valor 1 si és que sí i 0 si és que no.

Tanmateix, tal com es pot observar, al nostre summary(), veiem que un dels 4 coeficients estimats és l’ordenada a l’origen (Intercept). En els models de tipus ANOVA l’ordenada a l’origen resumeix un dels nivells del factor, que passa a ser el nivell de referència (en aquest cas serà C. Malasia), mentre que els altres coeficients indiquen quant es modifica aquest valor de referència en les diferents poblacions: \[ \small y_{ij}=\beta_0+\beta_1x_1+\beta_2x_2+\beta_3x_3+\epsilon_{ij}~~~ \epsilon_{ij} \sim N(0,\sigma^2) \]

contrasts() i model.matrix()

Dues funcions de R interessants per explorar les dades en aquests casos son contrasts() i model.matrix(). El primero explica la codificación de R para un factor concreto:

contrasts(chalco$Localidad)
##              E. Tailandia Nias S. Sumatra
## C. Malasia              0    0          0
## E. Tailandia            1    0          0
## Nias                    0    1          0
## S. Sumatra              0    0          1

El segon ens dona la taula completa de dummy variables que ha fet servir R per ajustar el model

head(model.matrix(lm.chalco))
##   (Intercept) LocalidadE. Tailandia LocalidadNias LocalidadS. Sumatra
## 1           1                     0             0                   1
## 2           1                     0             0                   1
## 3           1                     0             0                   1
## 4           1                     0             0                   1
## 5           1                     0             0                   1
## 6           1                     0             0                   1

Equació de la regressió

L’equació serà \[ y_i=12.035+0.690x_1+0.246x_2+0.020x_3 \] amb les variables en el mateix ordre que en el summary().

Com els coeficients \(\beta_1,\beta_2,\beta_3\) només donen la diferència amb l’ordenada a l’origen, haurem de fer el càlcu manual o utilitzar predict():

predict(lm.chalco, list(Localidad=levels(chalco$Localidad)))
##        1        2        3        4 
## 12.03529 12.72500 12.28136 12.05556

o amb tapply()

tapply(chalco$Long_genitalia,chalco$Localidad, mean)
##   C. Malasia E. Tailandia         Nias   S. Sumatra 
##     12.03529     12.72500     12.28136     12.05556

o amb dplyr

chalco |>
  dplyr::group_by(Localidad) |>
  dplyr::summarise(mean = mean(Long_genitalia))
## # A tibble: 4 × 2
##   Localidad     mean
##   <fct>        <dbl>
## 1 C. Malasia    12.0
## 2 E. Tailandia  12.7
## 3 Nias          12.3
## 4 S. Sumatra    12.1

Si, com en aquest cas, el test F rebutja la hipòtesi nula \(\small H_0:\mu_1=\mu_2=\mu_3=\mu_4=\mu\), el pas següent serà averiguar les pobleacions que es deferencien de la resta.

Ajustar un model sense nivell de referència

En principi, es pot ajustar un model que no tingui un nivell de referència: és una operació molt senzilla que consiteix en afegur un \(-1\) a la formula de lm():

summary(lm(Long_genitalia ~ Localidad - 1, data = chalco))
## 
## Call:
## lm(formula = Long_genitalia ~ Localidad - 1, data = chalco)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.68136 -0.36201  0.06985  0.38591  1.01864 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## LocalidadC. Malasia   12.03529    0.09759   123.3   <2e-16 ***
## LocalidadE. Tailandia 12.72500    0.09484   134.2   <2e-16 ***
## LocalidadNias         12.28136    0.07408   165.8   <2e-16 ***
## LocalidadS. Sumatra   12.05556    0.10951   110.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5691 on 152 degrees of freedom
## Multiple R-squared:  0.9979, Adjusted R-squared:  0.9979 
## F-statistic: 1.82e+04 on 4 and 152 DF,  p-value: < 2.2e-16

Tanmateix cal destacar que està bastant desaconsellat suprimir l’ordenada a l’origen ajustant un model i, en el cas que es faci, es recomana atenció i prudència en la interpretació dels coeficients i dels resultats dels tests.

Test post hoc i comparacions múltiples

Com hem vist, existeixen diferències entre els nivells del factor en la variable de resposta. Per verificar quines d’aquestes diferències són estadísticament significatives tenim a disposició uns tests específics. Dos que són molt utilitzats són el del test de Bonferroni i el de Tukey.

Amb R podem dur a terme el test de Bonferroni amb la funció pairwise.t.test() que genera una matru de valors p que ens indiquen si les mitjanes de cada parell de poblacions són iguals (\(\text{p-value} > 0.05\)) o diverses entre si (\(\text{p-value} < 0.05\)).

pairwise.t.test(chalco$Long_genitalia, chalco$Localidad,
                p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  chalco$Long_genitalia and chalco$Localidad 
## 
##              C. Malasia E. Tailandia Nias  
## E. Tailandia 6.9e-06    -            -     
## Nias         0.2783     0.0019       -     
## S. Sumatra   1.0000     4.9e-05      0.5383
## 
## P value adjustment method: bonferroni

Com que aquesta funció, malgrat disposar de diferent alternativas, no pot executar el mètode Tukey, en aquest cas haurem d’utilitzar una funció específica: TukeyHSD(), l’argument de la qual és un objecte aov() que recupera lm()

TukeyHSD(aov(lm.chalco))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = lm.chalco)
## 
## $Localidad
##                                diff         lwr        upr     p adj
## E. Tailandia-C. Malasia  0.68970588  0.33619841  1.0432133 0.0000068
## Nias-C. Malasia          0.24606181 -0.07222327  0.5643469 0.1894836
## S. Sumatra-C. Malasia    0.02026144 -0.36079022  0.4013131 0.9990557
## Nias-E. Tailandia       -0.44364407 -0.75626988 -0.1310183 0.0017799
## S. Sumatra-E. Tailandia -0.66944444 -1.04578188 -0.2931070 0.0000476
## S. Sumatra-Nias         -0.22580038 -0.56926469  0.1176639 0.3232117

Ambdós mètodes donen el mateix resultat, identificant el nivell E. Tailandia com l’únic que és significativament diferent dels altres.

El segon mètode permet fer se representat gràficament tamé

par(mar=c(3,6,3,0))

plot(TukeyHSD(aov(lm.chalco)),
     cex.axis=.5, las=1)

Avaluació dels supòsits del model

Podem avaluar els supòsits que ja hem vist. Cal destacar que el suòsit de linialitat no té gaire sentit en el cas de variables categòriques.

par(mfrow= c(2,2))
plot(lm.chalco)

A les gràfiques no s’observa cap problema aparent. Ara, les prediccions no es fan al llarg d’un continu, sinó sobre nivells predefinits. Per tant, no s’observarà un núvol de punts. Tanmateix no representa un problema, sempre i quan els punts es distribueixin de manera constant.

Vamos a testar la normalitat dels valors residuals con el Test de Shapiro-Wilk i l’homoscedasticitat amb el test de Levene:

shapiro.test(residuals(lm.chalco))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lm.chalco)
## W = 0.97446, p-value = 0.00538
car::leveneTest(Long_genitalia ~ Localidad, data=chalco)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   3  2.0797 0.1052
##       152

Es detecta certa falta de normalitat. Però no es veuen problemes aparents d’heteroscedasticitat.


Anàlisi de la covariança (ANCOVA)

L’anàlisi ANCOVA és equivalent a un ANOVA en què s’inclogui a més a més una variable explicativa contínua (denominada normalment covariable), per corregir l’hipotètic efecte d’alguna mena de gradient a l’hora de comparar mitjanes entre nivells d’un factor.

També es pot considerar equivalent a una regressió en la qual s’inclou un factor per a corregir les possibles diferències entre grups observacionals o experimentals.

Exemple firefly

Al data set firefly de ADER es presenten els resultats d’un experiment en el qual s’intenta comprovar si el nombre d’aparellaments de les femelles de Photinus ignitus té impacte en la grandària de les postes d’ous.

data("firefly")
str(firefly)
## 'data.frame':    40 obs. of  3 variables:
##  $ weight: num  9.69 11.52 19.53 31.38 30.42 ...
##  $ eggs  : num  50 56 77 68 58 53 51 56 43 36 ...
##  $ mating: Factor w/ 2 levels "single","triple": 1 1 1 1 1 1 1 1 1 1 ...

El nompre d’aparellaments podria tenir efecte tant positiu1, com negatiu2 en la posta d’ous.

A més a més (i aquí intervé la covariable), atès que la capacitat reproductiva dels insectes pot dependre també del seu pes corporal (per exemple, per les reserves acumulades en la fase de larva) pot ser raonable raonable assumir que les femelles més grans ponguin postes més grans.

En aquest esperiment doncs, dos grups de femelles es van sotmetre o bé a un únic aparellament o bé a tres secuencials amb diferents mascles.

Hi ha dues hipòtesis nul·les a comprovar:

  • \(H_0:\) No hi ha relació entre la gràndària de la posta (nombre d’ous postos) i pes corporal de l’insecte.
  • \(H_0:\) Desconnectant l’efecte del pes corporal, no hi ha relació entre la grandària de la posta i el nombre d’aparellaments.

La pregunta de la investigació està relacionada amb la segona hipòtesis nul·la. Tanmateix, habrá que testar també la primera com a part del procés inferencial.

Utilitzarem la suma de quadrats de tipus I o secuencial.

Homogeneïtat de pendents

Als supòsits que ja coneixem per als models linials (indipendència, linialitat, normalitat, heteroscedasticitat), n’hem d’afegir un cinquè: que la relació entre variable de resposta i covariable (és a dir la pendent de la regressió) és la mateixa en tots el nivells del factor.

Una manera d’explorar aquest relació és a través d’un gràfic amb panells. Nosaltres en aquest cas proposem utilitzar ggplot2, encara que existeix un model específic també en lattice:

firefly %>% 
  ggplot(aes(x=weight,y=eggs))+
  geom_point() +
  facet_wrap(~ mating)+
  geom_smooth(method = "lm", se = F) +
  theme_bw()

Aquí basicament s’observa com es comporten la relació linial Pes-Posta, separant les dades segons el tipus d’aparellament.

Observem que es generen 2 rectes amb pendents lleugerament diferents. Però no sabem si significativament diferents.

Test d’hipòtesi d’homogeneïtat de pendents

Cal incloure al model un terme d’interacció entre el factor i la covariable.

firefly0 <- lm(eggs ~ weight + mating + weight:mating, data = firefly)
anova(firefly0)
## Analysis of Variance Table
## 
## Response: eggs
##               Df Sum Sq Mean Sq F value    Pr(>F)    
## weight         1 7181.9  7181.9 36.1871 6.625e-07 ***
## mating         1 1793.2  1793.2  9.0353  0.004803 ** 
## weight:mating  1  245.2   245.2  1.2353  0.273741    
## Residuals     36 7144.8   198.5                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ara doncs obtenim una taula de l’anàlisi de la variança amb les sumes dels quadrats, els estadístics F i els nivells de significància de cada terme del model.

Veiem que tant la covariable (weight) com el factor (mating) tenen efecte significatiu en la variable de resposta, però no el terme d’interacció wight:mating. En aquest model es compleix la presumpció d’homogeneïtat de les pendents. Això significa que entre el pes de la femella i els ous postos és constant, indipendentment del nombre d’aparellaments.

Cuan un terme d’interacció no és significatiu cal reduir la taula, tornant a executar-la sense aquest:

firefly1 <- lm(eggs ~ weight + mating, data = firefly)
anova(firefly1)
## Analysis of Variance Table
## 
## Response: eggs
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## weight     1 7181.9  7181.9 35.9584 6.347e-07 ***
## mating     1 1793.2  1793.2  8.9782  0.004856 ** 
## Residuals 37 7390.0   199.7                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(firefly1)
## 
## Call:
## lm(formula = eggs ~ weight + mating, data = firefly)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.146  -8.685   0.319   7.429  36.613 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    9.1650     5.9631   1.537  0.13282    
## weight         1.5987     0.3141   5.089 1.07e-05 ***
## matingtriple -13.8046     4.6071  -2.996  0.00486 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.13 on 37 degrees of freedom
## Multiple R-squared:  0.5484, Adjusted R-squared:  0.524 
## F-statistic: 22.47 on 2 and 37 DF,  p-value: 4.096e-07

Equació general

\[ \small\text{eggs}=9.165+1.599~(\text{weight})-13.805~(\text{matingtriple}) \]

Revisió dels supòsits

par(mfrow=c(2,2))

plot(firefly1)

Els supòsit de linialitat i normalitat semblen respectar-se. Algun dubte tenim pel que fa a l’homogeneïtat de la variança (veure gràfic 3).

El tests de Shapiro i Levene confirmen el que acabem di dir sobre normalitat i heteroscedasticitat:

shapiro.test(residuals(firefly1))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(firefly1)
## W = 0.96585, p-value = 0.2639
car::leveneTest(residuals(firefly1), firefly$mating)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value   Pr(>F)   
## group  1  9.2361 0.004279 **
##       38                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Què passa quan la interacció es positiva?

L’homogenïtat entre pendents no hauria de ser considerada només un supòsit del model ANCOVA. Les interaccions entre factors i covariables poden representar efectes de considerable interès biològic. Les diferències entre pendents indiquen que els tractament afecta la relació entre covariable (\(x\)) i resposta \(y\). Explicar això pot ser tan important com explicar els efectes principals.

Exemple arrecifes

Un exemple d’interacció significadiva, el podem veure en aquest data set on s’analitza la competència entre algues i invertebrats bentònics sèssils en tres esculls californians. Les dades es van recollir en setze transectes de 40 m x 2 m durant 8 anys. Es tracta d’un data frame amb 69 observacions sobre les 3 variables següents:

  • sitio: Ubicació del transecte. Un factor amb tres nivells: *Carpinteria, Mohawk, Naples**.
  • algas: Percentatge de cobertura d’algues.
  • invertebrados: Percentatge de cobertura d’invertebrats sèssils.
data("arrecifes")
str(arrecifes)
## 'data.frame':    69 obs. of  3 variables:
##  $ sitio        : Factor w/ 3 levels "Carpinteria",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ algas        : num  66.7 50.6 43.4 65.7 14.3 47.2 46.2 6.9 23.1 13.9 ...
##  $ invertebrados: num  23.1 44.2 36.8 25.7 50 18.1 43.8 50 35.9 39.2 ...

Test d’hipòtesi

S’introdueix primer sitio, utilizant les sumes de quadrats de tipus I o seqüencial, i després la covariable.

arrecife0 <- lm(invertebrados ~ sitio + algas + sitio:algas, 
   data = arrecifes)
anova(arrecife0)
## Analysis of Variance Table
## 
## Response: invertebrados
##             Df Sum Sq Mean Sq  F value    Pr(>F)    
## sitio        2 8433.4  4216.7  46.9017 3.354e-13 ***
## algas        1 9109.7  9109.7 101.3262 9.681e-15 ***
## sitio:algas  2 1240.4   620.2   6.8985  0.001954 ** 
## Residuals   63 5664.0    89.9                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El test d’hipòtesi rebutja clarament l’homogeneïtat de pendents. Els resultats mostren que l’efecte de cobertura de les algues en la cobertura d’invertebrats no és idèntic en les tres eculleres.

En tot cas això no invalida el model, ja que explicar això pot ser interessant o més que explicar els efectes principals. Anirem per tant a examinar els resultats amb la funció summary()

summary(arrecife0)
## 
## Call:
## lm(formula = invertebrados ~ sitio + algas + sitio:algas, data = arrecifes)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -27.7363  -3.0847   0.6788   7.2087  14.1975 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        55.0834     4.8596  11.335  < 2e-16 ***
## sitioMohawk        27.5006     7.0323   3.911 0.000228 ***
## sitioNaples        33.1548     5.9730   5.551 6.07e-07 ***
## algas              -0.4509     0.1155  -3.905 0.000232 ***
## sitioMohawk:algas  -0.3504     0.1671  -2.097 0.040057 *  
## sitioNaples:algas  -0.7065     0.1931  -3.658 0.000520 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.482 on 63 degrees of freedom
## Multiple R-squared:  0.7683, Adjusted R-squared:  0.7499 
## F-statistic: 41.79 on 5 and 63 DF,  p-value: < 2.2e-16

Coeficients i equació del model

En aquest, cas el factor sitio té tres nivells. Per tant el nombre de coeficients estimats és de dos (Mowak i Naples), més un (Carpinteria) que queda resumit per l’ordenada en origen. El mateix passa amb el terme en interacció: l’efecte per la colonització d’algues (algas) resumeix la pendent de la recta de regressió en el nivell de referència (Carpinteria), mentre que els altres dos coeficientsmodificaran el valor de la pendent a cada una de les esculleres.

El model es podria per tant expressar amb la fórmula següent: \[ \small\text{invertebrados}=55.08+27.50~\text{Mohawk}+33.14~\text{Naples} -(0.45+0.35~\text{Mohawk}+0.71~\text{Naples})~\text{algas} \]

Representació del model

Com més complex sigui el model, més difícil serà interpretar els coeficients que s’obtinguin amb summary(). Una forma senzilla és generar prediccions i dibuixar totes les rectes en un únic gràfic per examinar visualment les diferències entre elles.

new_data <- list()

for (i in 1:length(levels(arrecifes$sitio))) {
  new_data[[i]] <-  data.frame(algas = seq(range(arrecifes$algas)[1],
                                   range(arrecifes$algas)[2], length=10),
                       sitio = rep(levels(arrecifes$sitio)[i], 10)
                       )
  
}

new_data <- bind_rows(new_data)
new_data$predict <- predict(arrecife0, newdata = new_data)

arrecifes %>% 
  ggplot(aes(x=algas, y= invertebrados,colour = sitio, linetype = sitio)) +
  geom_point() +
  geom_line(data = new_data, aes(x=algas,y=predict,colour = sitio)) +
  labs(title = "Predicció de l'efecte de la cobertura d'algues \nen la cobertura d'invertebrats segons localització",
       x= "Cobertura d'algues (%)",
       y= "Cobertura d'invertebrats (%)") +
  theme_bw() +
  theme(plot.title = element_text(size = 12, hjust = 0.5))

Alternativament, i de manera més immediata però potser no tant eficaç, podem fer servir la funció allEffects() que ja hem vist del paquet effects.

plot(effects::allEffects(arrecife0))

Test post hoc de comparació múltiple de pendents

A més a més podrem dur a terme un test post hoc de comparació múltiple, ja que rebutjant la hipòtesi nul·la de d’homogenïtat de pendents, sabem que, almenys una és distinta, però no sabem quina o quines.

Per a això podem fer servir lstrends() del paquet emmeans.

library(emmeans)
lst <- emmeans::lstrends(arrecife0, ~ sitio, var= "algas")
lst
##  sitio       algas.trend    SE df lower.CL upper.CL
##  Carpinteria      -0.451 0.115 63   -0.682   -0.220
##  Mohawk           -0.801 0.121 63   -1.043   -0.560
##  Naples           -1.157 0.155 63   -1.467   -0.848
## 
## Confidence level used: 0.95
pairs(lst)
##  contrast             estimate    SE df t.ratio p.value
##  Carpinteria - Mohawk    0.350 0.167 63   2.097  0.0987
##  Carpinteria - Naples    0.707 0.193 63   3.658  0.0015
##  Mohawk - Naples         0.356 0.196 63   1.813  0.1735
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Això ens indica que les pendents significativamente diferents són les de Carpinteria i Naples amb un valor p de 0.0015. La resta de pendents poden ser diferents pero no són significativament diferents.

Suma de quadrats: quin model seleccionar (de tipus I, II o III)

Si el disseny factorial és totalment ortogonal, l’ús de diferents sumes de quadrats donarà els mateixos resultats, que equival a allò que hem vist d’aquest concepte parlant de regressió múltiple (veure Nota núm. 16). El concepte d’ortogonalitat en un model factorial és equivalent al de falta de colinealitat en un model de regressió.

Sumes de quadrats de tipus I

Quan en canvi el disseny no està tan equilibrat, diferents tipus de seumes de quadrats poden donar resultats diferents. El tipos de suma de quadrats a utilitzar depèn de la hipòtesi que estiguem comparant: si el nostre objectiu és testar interaccions, el tipus de SS és indiferent, perquè si la interacció és significativa no té massa sentit examinar els efectes individuals. La cosa més senzilla en aquest cas seria les sumes de quadrats de tipus I amb la funció stats::anova(): primer es testa l’efecte principal d’\(A\), després el de \(B\) i, finalment, es testa la ineracció \(A:B\):

Models de la hipòtesi nula i alternativa en l’anàlis seqüencial de la variança, amb sumes de quadrats de tipus I, d’un model \(y\sim A*B\) amb la funció stats::anova()

Término \(H_0\) \(H_1\)
\(A\) \(Y\sim 1\) \(Y\sim A\)
\(B\) \(Y\sim A\) \(Y\sim A+B\)
\(A:B\) \(Y\sim A+B\) \(Y\sim A+B+A:B\)

Un altre dels casos en què cal utilitzar el tipus I és quan un dels factors és un bloc. Això vol dir que, tot i no ser d’interès per a l’inverstigador, el factor imposa una estructura d’autocorrelació a les dades. Un exemple de bloc, particularment important en biologia, és l’emplaçament que imposa uns condicionants a la constitució d’una mostra.

Amb la qual cosa tindrem un factor d’interès real i un segon factor, l’emplaçament, que haurem d’afegir per incorporar la possible correlació entre mostres recollides al mateix emplaçament.

En aquest cas doncs seria pertinent utilitzar les sumes de quadrats de tipus I, introduint primer el bloc (emplaçament) i després la resta de variables, un cop calculat l’efecte de l’emplaçament.

Sumes de quadrats de tipus III

Quan no hi hagi un factor de tipus bloc o una variable que coneguem ja amb anterioritat, en la majoria dels casos, seria més convenient fer servir les sumes de quadrats de Tipus III. Són les més utilitzades (per defecte en SPSS), també perquè no depenen de la grandària de la mostra.

Important: les sumes de quadrats de tipus III són recomanables només si el disseny és totalment ortogonal o no existeix desequilibri entre les variables.

Models de la hipòtesi nula i alternativa en l’anàlis seqüencial de la variança, amb sumes de quadrats de tipus III, d’un model \(y\sim A*B\) amb la funció car::Anova()

Término \(H_0\) \(H_1\)
\(A\) \(Y\sim B+A:B\) \(Y\sim A+B+A:B\)
\(B\) \(Y\sim A+A:B\) \(Y\sim A+B+A:B\)
\(A:B\) \(Y\sim A+B\) \(Y\sim A+B+A:B\)

Sumes de quadrats de tipus II

Quan el desequilibri més gran les sumes de quadrats recomanades serien les de Tipus II.

Models de la hipòtesi nula i alternativa en l’anàlis seqüencial de la variança, amb sumes de quadrats de tipus II, d’un model \(y\sim A*B\) amb la funció car::Anova()

Término \(H_0\) \(H_1\)
\(A\) \(Y\sim B\) \(Y\sim A+B\)
\(B\) \(Y\sim A\) \(Y\sim A+B\)
\(A:B\) \(Y\sim A+B\) \(Y\sim A+B+A:B\)

  1. Tenint en compte que aquests insectes no s’alimenten en la fase adulta, el nombre d’ous podria estar condicionat pels nutients aportat pels espermatòfors.↩︎

  2. L’aparellament pot representa una despesa energètica per a la femella.↩︎