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

Introducció

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.
M.Y. Cabrero, A. García, Análisis estadístico de datos espaciales con QGIS y R, 2020.
S.N. Wood, Generalized Additive Models: An Introduction with R, 2017.

Models binomials o models de regressió logística

Un cas particular de model linial generalitzat són els models binomials. Les dues princiapls variants d’aquest model difereixen per la forma amb què es mesura la variable de resposta. Aquesta pot

Agafar valors binaris, representant que un esdeveniment ha tingut lloc (èxit) o no (fracàs).
Agafar valors proporcionals, on cada observació representaria el nombre d’èxits en un nombre fix d’assaigs.

La principal característica d’aquest model és que el fenòmen de la variable de resposta es parametritza com una probabilitat, amb uns valors acotats entre 0 i 1. Aquest tipus de model en la literatura es coneixen com models de regressió logística.

En els GLM binomials la resposta consisteix en dos valors: nombre d’èxits (\(k\)) i nombre total d’assaigs (\(n\)). Ambdós tipus de model s’assumeix que les observacions procedeixen d’una distribució binomial amb distribució de probabilitat: \[ \small P(k)=(_k^n)p^k(1-p)^{(n-k)} \] on: \[ \small (_k^n)=\frac{n!}{k!(n-k)!} \] i \(p\) és la probabilitat d’un èxit.

Aquí doncs presentarem ambdós tipus del model binomial mitjançant casos d’estudi, proposant com ajustar-los, avaluar-los i interpretar-los.

GLM de resposta binària

Moltes varibles de resposta són de tipus binari, com ara - mortalitat: les observacions en els individus només poden agafar dos valors (viu, mort). - epidemiologia: individus infectats o no (indexs de morbiditat). - etc.

En tots aquests casos podem investigar quines variables explicatives relacionades amb l’assignació d’una classe o de l’altra amb els models GLM, on la variable de resposta hauria de ser codificada amb dos valors, 0 i 1. La codificació dependrà de la hipòtesi de partida i de la natura del procès a investigar.

Per a justar un model binomial en R cal utilitzar la funció glm(), especificant l’argument family=binomial, que per defecte té associada una funció de tipus logit, que es defineix com \[ \small logit(y)=\log\bigg(\frac{p}{1-p}\bigg) \] L’avantatge de fer servir el cocient de probabilitat com a variable de resposta en compte de \(p\) és que pot agafar valors per sobre de 1, permetent així l’ús del predictor linial, \(\eta\), per relacionar-lo amb el logit. \[ \small logit(y)=\log\bigg(\frac{p}{1-p}\bigg)= \eta=\beta_0+\beta_1x_1+...+\beta_qx_q \]

Això implica que, per a obtenir prediccions sobre la variable de resposta hem d’aplicar la inversa a la funció de vincle \[ \small p=\frac{e^\eta}{1+e^\eta}=\frac{e^{(\beta_0+\beta_1x_1+...+\beta_qx_q)}}{1+e^{(\beta_0+\beta_1x_1+...+\beta_qx_q)}} \] La forma d’aquesta funció és una corba sigmoïdal, acotada entre els valors 0 i 1.

Cas en estudi: l’estratègia d’una alga invasora

S’estudia l’estratègia reproductiva d’una espècie d’alga verda Codium fragile, en relació dues espècies del mateix gènere, C. tomentosum i C. vermilara amb què cohabita. Les dedades son freuit de mostreig aleatori a diversos punts de la Costa Cantàbrica. En total hi ha 359 individus de les tres espècies. Les dedades estan disponibles amb el paquet ADER.

library(ADER)
data("algas")

str(algas)

## 'data.frame':    359 obs. of  3 variables:
##  $ Sp    : Factor w/ 3 levels "Fra","Tom","Ver": 3 3 3 3 3 3 3 3 3 3 ...
##  $ Long  : num  1 2 4 4 4 4.5 4.5 4.5 5 5 ...
##  $ Estado: int  0 0 0 0 0 0 0 0 0 0 ...

Estado: és la variable resposta. Només pot agafar dos valors, 0 i 1, que representen l’absència o la presència d’estructures reproductores.
Long: variable explicativa que indica la llargada de l’individu.
Sp: variable explicativa que indica l’espècie de l’individu (Fra, C. fragile, Tom, C. Tometosum, Ver, C. vermilara).

La hipòtesi seria que la grandària mínima per reproduir-se seria menor en la espècie invasora (Fra, Codium fragile) que en les altres dues.

Anàlisi exploratòria

algas %>% 
  group_by(Sp, Estado) %>% 
  summarise(M = mean(Long), D = sd(Long), n = n(), 
            upper= M+1.96*D/sqrt(n), lower=M-1.96*D/sqrt(n)) %>%  # Càlcul de l'I.C.
  ggplot(aes(x = Sp, y=M, fill= factor(Estado))) +
  geom_crossbar(aes(ymin = lower, ymax = upper))+
  coord_flip() +
  labs(title = "Llargades mitjanes i intervals de confiança al 95% per a individus \nno reproductius (0) i reproductius(1) de tres espècies d'algues",
       y= "Llargada",
       x= "Espècie") +
  theme_bw()

  algas %>% 
  group_by(Sp, Estado) %>%
  ggplot(aes(x= Sp, y=Long, fill = factor(Estado))) +
  geom_boxplot()+
  coord_flip() +
  labs(title = "Gràfic de caixes per a individus no reproductius (0) i reproductius(1) \nde tres espècies d'algues",
       y= "Llargada",
       x= "Espècie")+
  theme_bw()

D’ambdues gràfiques es desprèn que hi ha clara separació en llargada entre individus en fase reproductiva i no. A més a més l’espècie C. fragile sembla tenir menor talla en ambdues fases.

Si existeix una interacció entre espècie i llargada, aquesta no queda molt evident a primera vista.

Ajustament i resolució d’hipòtesi

Ena aquest cas, farem servir, per al test d’hipòtesi farem servir la suma de quadrats de tipus I. No hi ha un efecte bloc a compensar, però la relació entre llargada i presència d’estructures reproductores i bastant òbvia. Per tant, una suma de quadrats de tipos I on introduïm primer la llargada i després l’espècie sembla ser més adequada.

glm.algas <- glm(Estado ~ Long+Sp+Long:Sp, data = algas, family = binomial)
anova(glm.algas)

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Estado
## 
## Terms added sequentially (first to last)
## 
## 
##         Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      358     476.38              
## Long     1  103.010       357     373.38 < 2.2e-16 ***
## Sp       2   31.771       355     341.60 1.262e-07 ***
## Long:Sp  2    0.459       353     341.15    0.7949    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Veiem com la interacció no és significativa. Caldrà per tant tornar a ajustar el model apartant la interacció

glm.algas2 <- glm(Estado ~ Long+Sp, data = algas, family = binomial)
anova(glm.algas2)

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Estado
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                   358     476.38              
## Long  1  103.010       357     373.38 < 2.2e-16 ***
## Sp    2   31.771       355     341.60 1.262e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretació del model

La funció `summary()`

La funció summary() ens donarà els coeficients del model

summary(glm.algas2)

## 
## Call:
## glm(formula = Estado ~ Long + Sp, family = binomial, data = algas)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.80703    0.42701  -6.574 4.91e-11 ***
## Long         0.26505    0.03174   8.351  < 2e-16 ***
## SpTom       -0.57476    0.34528  -1.665    0.096 .  
## SpVer       -1.84877    0.36933  -5.006 5.57e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 476.38  on 358  degrees of freedom
## Residual deviance: 341.60  on 355  degrees of freedom
## AIC: 349.6
## 
## Number of Fisher Scoring iterations: 5

Aquests coeficients faran referència al logaritme natural del cocient de probabilitats, és a dir al \(\small logit(y)=\log(\frac{p}{1-p})\).

Predicció de valors

Possibles valors de resposta segons canvis variables explicatives

nuevas_combinaciones <- expand.grid(
  Long = c(10, 11,20,21),
  Sp = c("Fra", "Tom","Ver")
)

# 3. Predicción
nuevas_combinaciones$probabilidad <- round(predict(glm.algas2,
                                 newdata=nuevas_combinaciones,
                                 type="response"),3)


nuevas_combinaciones

##    Long  Sp probabilidad
## 1    10 Fra        0.461
## 2    11 Fra        0.527
## 3    20 Fra        0.924
## 4    21 Fra        0.940
## 5    10 Tom        0.325
## 6    11 Tom        0.386
## 7    20 Tom        0.872
## 8    21 Tom        0.899
## 9    10 Ver        0.119
## 10   11 Ver        0.149
## 11   20 Ver        0.656
## 12   21 Ver        0.713

Nivell efectiu de la mediana

Un estadístic molt útil en el cas de regressions logístiques és el Nivell efectiu de la mediana, calculat com el negatiu del quocient entre intercepte i pendent, \(\small-\beta_0/\beta_1\). Aquest quocient ens indica el valor de la variable explicativa per al qual els dos possibles resultats són igual de probables.

Tanmateix, això seria només en el cas del grup de referència (C. fragile).

Per a qualsevol grup addicional (nivell \(i\) del Factor \(F\)) seria (assumint que el efecte dosi sigui el mateix per a tots els grups): \[ \small x_{ED50}= -\frac{\beta_0+\beta_{F_i}}{\beta_{1}} \]

Tal com veiem aquí sota:

coef.alg <- glm.algas2$coefficients
coef.alg

## (Intercept)        Long       SpTom       SpVer 
##  -2.8070319   0.2650514  -0.5747628  -1.8487745

C. Fragile (nivell de referència)

-coef.alg[1]/coef.alg[2]

## (Intercept) 
##    10.59052

C. Tomentosum

-(coef.alg[1]+coef.alg[4])/coef.alg[2]

## (Intercept) 
##    17.56568

C. Vermilara

-(coef.alg[1]+coef.alg[4])/coef.alg[2]

## (Intercept) 
##    17.56568

Deviança explicada pel model i capacitat predictiva

La deviança, \(\small D^2=1-D/D_{nul·la}\), seria

1 - glm.algas2$deviance / glm.algas2$null.deviance

## [1] 0.2829237

En models de resposta binària, \(D^2\) de 0.2-0.3 són la tònica. En general, per sobre de 0.3, lajustament es jutja bastant bo.

Taula de classificació

En tot cas, per a models amb resposta discreta, com aquest, es poden utilitzar estadístics que mesuren el poder predictiu del model: és a dir la capacitat de predir el resultat a del model ajustat. Això se sol resumir en una taula de classificació en què es creuen els valors observats i els seus corresponents valors observats. Com aquests són probabilitats es fa servir una probabilitat llindar (p.ex. \(p_0=0.50\)) per convertir-los en uns i zeros.

table(algas$Estado, fitted(glm.algas2)>0.5)

##    
##     FALSE TRUE
##   0   192   31
##   1    48   88

És a dir de 136 individus reproductors el model prediu correctament 88, i de 223 no reproductors el model prediu 192.

Dues mesures útils del poder predictiu son la sensibilitat, o la proporció d’èxits correctament classificats i la especificitat, o proporció de fracassos correctament classificats. En aquest exemple la sensibilitat seria 88/136=0.65, i la especificitat seria 192/223=0.86.

La corba ROC

El problema de les tables de classificació i dels estadístics que se’n deriven és que són molt sensibles a l’elecció de la probabilitat llindar (\(p_0\)). Una alternativa fora utilitzada és la corba ROC (característica operativa del receptor): un gràfic que mostra la sensibilitat i l’especificitat per a tots els possibles valors de probabilitat llindar \(p_0\).

library(pROC)

rocplot <- roc(Estado ~ fitted(glm.algas2), data=algas)
plot(rocplot, xlab= "Especificidad", ylab= "Sensibilidad")

La corba ROC es pot obtenir amb la funció roc() del paquet pROC. L’àrea sota la corba ROC (AUC) proporciona un sesum del valor predictiu del model. L’AUC estima la probabilitat qur i les prediccions del model coincideixin. Un AUC de 0.5 (la diagonal del gràfic) indicaria que les prediccions del model no són millors que predictors completament aleatoris. En general, s’assumeix que un model amb bona capacitat discriminant té un AUC per sobre de 0.80.

auc(rocplot)

## Area under the curve: 0.8473

En aquest cas, l’AUC arriba gairebé a 0.85, que indica una bona capacitat discriminant del model.

Gràfics de predicció

Per obtenir les prediccions de probabilitat de tenir estructures reproductores (\(p\)) serà necessari aplicar la funció de vincle sobre el predictor linial \(\eta\). Per això generarem un data frame amb un vector de valors de la llargada de l’alga per a cada una de les 3 espècies. Sobre aquest data frame generarem les prediccions, amb la funció predict(), amb l’argument type=response, important per obtenir \(p\) (si no ho fem així, predict() ens donarà els valors de \(\tiny \log\bigg(\frac{p}{1+p}\bigg)\), sense aplicar la funció inversa).

new_algas <-  data.frame(Sp = rep(c("Fra","Tom","Ver"), each=50),
                               Long = rep(seq(range(algas$Long)[1],
                                   range(algas$Long)[2], 
                                   length=50), length=150)
                               )

# 'type = "response"' és per a aplicar la funció inversa al logit.
new_algas$predict <- predict(glm.algas2,
                            newdata = new_algas,
                            type = "response")

algas %>% 
  ggplot(aes(x=Long, y= Estado)) +
  geom_point() +
  geom_line(data = new_algas, aes(x=Long, y=predict, colour = Sp)) +
  labs(title = "Relació entre probabilitat d'estructures reproductives \ni llargada segons espècie d'alga",
       x= "Llargada (cm)",
       y= "Prob.(Estat reproductiu)",
       colour= "Especie") +
  theme_bw() +
  scale_color_manual(labels = 
                       c("C. fragile", "C. tomentosum", "C. vermilara"), 
                     values = c("dodgerblue3","darkred","darkgreen")) +
  theme(plot.title = element_text(size = 12, hjust = 0.5),
        legend.position = c(.95, .1),
        legend.justification = c('right', 'bottom'))

S’observa com la probabilitat de trobar estructures reproductives augmenta amb amb la llargada de l’individu. Tanmateix amb igualtat de llargada l’espècie invasora Codium fragile té més probabilitat de presentar estructures reproductores respecte a les altres.

Tests post hoc

En algun cas, pot ser intressant realitzar un test post hoc per identificar quins grups són diferents dels altres. Per això es pot emprar la funció glt() de la llibreria multcomp.

Com que la interacció entre covariable i factor no es significativa, compararem simplement els nivells del factor Sp.

summary(multcomp::glht(glm.algas2, linfct = multcomp::mcp(Sp="Tukey")))

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: glm(formula = Estado ~ Long + Sp, family = binomial, data = algas)
## 
## Linear Hypotheses:
##                Estimate Std. Error z value Pr(>|z|)    
## Tom - Fra == 0  -0.5748     0.3453  -1.665 0.217806    
## Ver - Fra == 0  -1.8488     0.3693  -5.006  < 1e-04 ***
## Ver - Tom == 0  -1.2740     0.3126  -4.076 0.000134 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

El summary mostra que, per a una llargada determinada, no hi ha diferència entre C. Tomentosum i C. Fragile, però sí, entre aquests dos i C. vermilara.

En el cas d’haver tingut interacció, una altrernativa interessant podria haver estat comparar les pendents de les rectes de dues en dues, que es pot fer amb lstrends() de emmeans.

lst <- emmeans::lstrends(glm.algas2, ~ Sp, var= "Long")
pairs(lst)

Revisió dels supòsits del model

La funció genèrica plot() mostra els principals gràfics de residus quan l’argument és un objecte glm.

par(mfrow= c(2,2))
plot(glm.algas2)

Aquests gràfics solen tenir formes peculiars però ens poden donar algunes informació útil.

Gràfic Residuals vs Fitted: residuals estandarditzats (Pearson) vs. valors esperats

plot(glm.algas2,1)

Tenim dues línes de punts perquè prediem la probabilitat d’una variable que agafa valors 0 i 1. Si el valor observat és 0, llavors el valors predit serà més gran i el residu negatiu, mentre que si el valor observat és 1, el valor predit serà menor i el residu positiu.

En general, com més expliqui el model, més propers seran els dos extrems de les línies, positiva i negativa, respectivament.

Gràfic Normal Q-Q

plot(glm.algas2,2)

Es valora com sempre, l’adherència dels punts a la línia.

La llibreria DHARMa

En GLM, els gràfic dels residus no sempre no sempre ens permeten d’avaluar correctament la idonitat del model. És veritat que els residus re-escalats (de Pearson o de deviança) permeten millorar-ne la interpretació.En tot cas, aquests gràfics no són sempre fàcils d’interpretar directament.

En resposta a això, es va publicar el paquet DHARMa. Aquest paquet disposa d’una sèrie de funcions que generen i avaluen residus estandarditzats amb calors acotats a 0-1, que poden ser interpretats de manera tan intuïtiva com els lineals. Aquests residus s’obtenen mitjançant simulacions com en el bootstrapping.

Els residus estandarditzats es poden calcular amb la funció simulateResiduals(). Si el nombre d’observacions és baix podem establir un nombre de simulacions alt, per exemple n=1000. Si el model és complex i el nombre d’observacions alt, es pot establir n=250.

simul <- DHARMa::simulateResiduals(glm.algas2, n=1000)

par(mfrow= c(2,2),
    mar= c(4,4,3,2))
DHARMa::plotQQunif(simul,
     testUniformity = FALSE, 
     testOutliers = FALSE,
     testDispersion = FALSE,
     xaxs = "i", yaxs = "i")
DHARMa::plotResiduals(simul,xaxs = "i", yaxs = "i")

En principi, ambdós gràfics mostren que el model és l’adequat.

DHARMa::testUniformity(simul, plot = F)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.02817, p-value = 0.9382
## alternative hypothesis: two-sided

GLM de resposta de tipus proporció

En molts casos (mortalitar, infecció, resposta clínica, germinació, etc.), cada observació correspon amb una proporció. En aquests casos és necessari conèixer el nombre de vegades que s’esdevé un fenomen i el nombre de vegades que no.

Els processos que impliquen una resposta proporcional es poden es poden modelar en R amb un GLM d’una família de tipus binomial. La diferència respecte al cas de dades amb resposta binària, és que, a més d’especificar la família, haurem d’especificar també el nombre d’èxits i fracassos que estarem modelant en una única resposta. Per això, haurem d’unir dos vectors en un sol objecte \(y\), mitjançant el manament cbind() que representarà la variable de resposta.

Naturalment això presenta quatre problemes:

Els errors no estan normalment distribuïts.
La variança no és constant.
Les prediccions del model no queden acotades als valors entre 0 i 1.
Perdem informació sobre la grandària de la mostra, \(n\), sobre la qual s’estima la proporció.

Cas d’estudi: la germinació de l’ailant

La base de dades ailanto presenta dades de l’estudi de germinació de plàntules d’ailant. Són 48 safates amb 40 alvèols cada safata que es van posar en condicions diferents de llum i humitat: 100%, 70%, 40%, 10% d’exposició a la llum i 90%, 70%, 50% d’humitat. Per a cada safata es va anotar el nombre de llavors germinades.

data("ailanto")
str(ailanto)

## 'data.frame':    48 obs. of  4 variables:
##  $ Semillas   : int  40 40 40 40 40 40 40 40 40 40 ...
##  $ Luz        : num  0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
##  $ Agua       : num  0.5 0.5 0.5 0.5 0.7 0.7 0.7 0.7 0.9 0.9 ...
##  $ Germinacion: int  29 31 27 27 34 35 27 40 25 30 ...

El primer pas serà transformar las variables Luz i Agua a categòricas. Cal dir que aquestes variables haurien pogut ser tractades com a contínues si hi hagués hagut una relació mínimament lineal o linealitzable amb la percentual de germinació.

ailanto <- ailanto %>% 
  mutate(Luz = factor(scales::percent(Luz), 
                      levels = c("10%","40%","70%", "100%")),
         Agua = factor(scales::percent(Agua), levels = c("50%","70%","90%"))
         ) %>% 
  mutate(yprop = Germinacion / Semillas)

ailanto %>% 
  ggplot(aes(x=Luz,y=yprop, fill = Agua)) +
  geom_boxplot() +
  facet_wrap(~ Agua) +
  labs(title = 
         "Ratio de germinació de l'ailant segons nivells de llum i d'aigua",
       x= "Nivells de llum",
       y= "Prop. Germinació",
       subtitle = "Nivells d'aigua") +
  theme_bw()+
  theme(plot.subtitle = element_text(hjust = 0.5),
        legend.position = "none")

Ajustament del model

Per ajustat el model a les dades amb la funció glm(), es poden seguir 2 vies: una que utilitza la proporció d’èxit que acabem de calcular. Necessita que s’afegeixi el paràmetre weights a la funció glm(). El segon mètode crea un vector doble d’èxits i fracassos que col·loca com a variable de resposta a la funció.

Mètode amb `weights`

glm.ailant <- glm(yprop ~ Luz * Agua, 
                  weights = Semillas, 
                  data = ailanto,
                  family = binomial)

Mètode unint èxit i fracàs amb `rbind()`

y <- rbind(ailanto$Germinacion,ailanto$Semillas- ailanto$Germinacion)

glm.ailant2 <- glm(y ~ Luz * Agua, data = ailanto, family = binomial)

Contrast d’hipòtesi

car::Anova(glm.ailant, type=3)

## Analysis of Deviance Table (Type III tests)
## 
## Response: yprop
##          LR Chisq Df Pr(>Chisq)    
## Luz        57.176  3  2.357e-12 ***
## Agua       10.522  2   0.005189 ** 
## Luz:Agua   51.020  6  2.935e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Es tria utilitzar una suma de quadrats de tipus III per exclusió:

Descartem el tipus I ja que coneixem poc els efectes de ambdues variables i cap d’aquestes és un bloc que generi diferències en la mostra.
Podria haver-hi l’alternativa del tipus II, però la mostra està compensada, com es pot veure fàcilment agrupant les safates segons les condicions d’aigua i llum:

ailanto %>% 
  summarise(n = n(), .by = c(Agua, Luz))

##    Agua  Luz n
## 1   50%  10% 4
## 2   70%  10% 4
## 3   90%  10% 4
## 4   50%  40% 4
## 5   70%  40% 4
## 6   90%  40% 4
## 7   50%  70% 4
## 8   70%  70% 4
## 9   90%  70% 4
## 10  50% 100% 4
## 11  70% 100% 4
## 12  90% 100% 4

D’altra banda, triant una opció o una altra d’aquestes 3, les diferències en aquest serien petites.

Interpretació del model

La funció `summary`

summary(glm.ailant)

## 
## Call:
## glm(formula = yprop ~ Luz * Agua, family = binomial, data = ailanto, 
##     weights = Semillas)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      0.90756    0.17467   5.196 2.04e-07 ***
## Luz40%          -0.42331    0.23876  -1.773 0.076234 .  
## Luz70%           1.80049    0.37037   4.861 1.17e-06 ***
## Luz100%          0.65059    0.27156   2.396 0.016586 *  
## Agua70%          0.82704    0.28201   2.933 0.003361 ** 
## Agua90%          0.08483    0.24978   0.340 0.734129    
## Luz40%:Agua70%   0.37511    0.39170   0.958 0.338248    
## Luz70%:Agua70%  -0.45057    0.57930  -0.778 0.436701    
## Luz100%:Agua70% -0.95936    0.40345  -2.378 0.017412 *  
## Luz40%:Agua90%   1.27675    0.37950   3.364 0.000767 ***
## Luz70%:Agua90%  -1.23669    0.46250  -2.674 0.007497 ** 
## Luz100%:Agua90%  0.77568    0.43519   1.782 0.074687 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 218.38  on 47  degrees of freedom
## Residual deviance: 101.00  on 36  degrees of freedom
## AIC: 274.56
## 
## Number of Fisher Scoring iterations: 5

Com que hi ha una interacció significativa, es fa molt difícil entendre els efectes de cada factor sobre la variable de resposta mirant exclusivament als coeficients estimats.

Representació del model

Seria convenient treure gràfics amb les respostes mitjanes per a cada combinació dels nivells dels factors.

emmeans::emmip(glm.ailant, Agua ~ Luz, 
               type = "response", CIs=T)+
  theme_bw()

Deviança explicada pel model

Com ja hem vist, la deviança \(D^2\) explicada del model seria:

D2 <- 1 - glm.ailant$deviance / glm.ailant$null.deviance
D2

## [1] 0.537519

Però, al contrari d’allò que hem vist per als models de resposta binària (on una \(D^2\) de 0.2-0.3 és la tònica), en els models GLM de resposta tipus proporció, la deviança explicada pot apropar-se a 1 si l’ajustament és bo. En aquest cas, explicant el model un 53.7% de deviança total, es pot considerar un ajustament raonablement bo.

Prediccions

Taula de valors mitjans

nuevas_combinaciones <- expand.grid(
  Luz = unique(ailanto$Luz),
  Agua = unique(ailanto$Agua)
)

# 3. Predicción
nuevas_combinaciones$p.media <- round(predict(glm.ailant,
                                 newdata=nuevas_combinaciones,
                                 type="response"),3)

nuevas_combinaciones %>% 
  flextable::flextable()

Luz	Agua	p.media
10%	50%	0.713
40%	50%	0.619
70%	50%	0.937
100%	50%	0.826
10%	70%	0.850
40%	70%	0.844
70%	70%	0.956
100%	70%	0.806
10%	90%	0.730
40%	90%	0.864
70%	90%	0.826
100%	90%	0.918

Taula de valors mitjans i intervals de confiança

Per obtenir l’error estàndard (EE) a l’escala de percentatge (probabilitat) d’una GLM (binomial a R, primer heu d’extreure les prediccions i els errors estàndard a l’escala de logit, construir els intervals de confiança i després transformar aquests valors mitjançant la funció logit invers.

Per què no podeu transformar l’EE directament simplement. La funció predict.glm() calcula els errors estàndard linealment a l’escala de probabilitats logarítmiques. Si sol·liciteu prediccions directament a l’escala de percentatge mitjançant type = "response", els errors estàndard faltaran o seran matemàticament defectuosos perquè la conversió no és lineal.

p2 <- predict(glm.ailant, newdata = nuevas_combinaciones, type = "link", se.fit = TRUE)

# Extreure valors mitjans i EE en escala logit
p2link <- p2$fit
p2se.link <- p2$se.fit 

# Calcular el límits de l'interval al 95% en escala logit
lo.link <- p2link - (1.96 * p2se.link)
up.link <- p2link + (1.96 * p2se.link)

# Funció per tornar a l'escala de la variable de resposta
inv_logit <- function(x) { exp(x) / (1 + exp(x)) }

# Convertir en l'escala original
pred2 <- data.frame(
  predicted = inv_logit(p2link),
  lower     = inv_logit(lo.link),
  upper     = inv_logit(up.link)
  )
df <- cbind(nuevas_combinaciones[,-3], pred2)

df %>% 
  flextable::flextable()

Luz	Agua	predicted	lower	upper
10%	50%	0.7125000	0.6376539	0.7772857
40%	50%	0.6187500	0.5412099	0.6906751
70%	50%	0.9375000	0.8877456	0.9660453
100%	50%	0.8260870	0.7596224	0.8771462
10%	70%	0.8500000	0.7859443	0.8973897
40%	70%	0.8437500	0.7789647	0.8921752
70%	70%	0.9562500	0.9110855	0.9790018
100%	70%	0.8062500	0.7376473	0.8603109
10%	90%	0.7295597	0.6553075	0.7928726
40%	90%	0.8636364	0.7998873	0.9093779
70%	90%	0.8258065	0.7578834	0.8777482
100%	90%	0.9182390	0.8642873	0.9519353

df %>% 
  ggplot(aes(x=Luz, y=predicted, fill = Luz))+
  geom_bar(stat = "identity") +
  facet_wrap(~Agua) +
  geom_errorbar(ymin= df$lower, ymax= df$upper, width=.3) +
  labs(x= "Nivells de llum",
       y= "Predicció (prop.)",
       subtitle = "Nivells d'aigua") +
  theme_bw() +
  theme(legend.position = "none",
        plot.subtitle = element_text(hjust = 0.5))

Test post hoc

En aquest tipus de models factorials, sol ser interessant calcular quines combinacions de variables són iguals o diferents entre si. Els tests post hoc serveixen per això. Molt útil en aquest sentit és el paquet emmeans mitjançant dues funcions emmeans() i contrast() que ens permeten analitzar l’efecte en la variable de resposta, comparant els nivells de dos en dos.

Fixar-se com quadren aquests valors \(p\) observant les barres d’error de la gràfica anterior.

seed.emm <- emmeans::emmeans(glm.ailant, ~ Luz*Agua)

emmeans::contrast(seed.emm, method= "pairwise", simple= "each", combine=F, adjust= "tukey")

## $`simple contrasts for Luz`
## Agua = 50%:
##  contrast   estimate    SE  df z.ratio p.value
##  10% - 40%    0.4233 0.239 Inf   1.773  0.2864
##  10% - 70%   -1.8005 0.370 Inf  -4.861 <0.0001
##  10% - 100%  -0.6506 0.272 Inf  -2.396  0.0779
##  40% - 70%   -2.2238 0.365 Inf  -6.094 <0.0001
##  40% - 100%  -1.0739 0.264 Inf  -4.067  0.0003
##  70% - 100%   1.1499 0.387 Inf   2.970  0.0158
## 
## Agua = 70%:
##  contrast   estimate    SE  df z.ratio p.value
##  10% - 40%    0.0482 0.311 Inf   0.155  0.9987
##  10% - 70%   -1.3499 0.445 Inf  -3.031  0.0130
##  10% - 100%   0.3088 0.298 Inf   1.035  0.7290
##  40% - 70%   -1.3981 0.444 Inf  -3.152  0.0088
##  40% - 100%   0.2606 0.296 Inf   0.881  0.8146
##  70% - 100%   1.6587 0.435 Inf   3.811  0.0008
## 
## Agua = 90%:
##  contrast   estimate    SE  df z.ratio p.value
##  10% - 40%   -0.8534 0.295 Inf  -2.893  0.0199
##  10% - 70%   -0.5638 0.277 Inf  -2.035  0.1750
##  10% - 100%  -1.4263 0.340 Inf  -4.194  0.0002
##  40% - 70%    0.2896 0.316 Inf   0.916  0.7963
##  40% - 100%  -0.5728 0.373 Inf  -1.537  0.4153
##  70% - 100%  -0.8625 0.359 Inf  -2.405  0.0761
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## 
## $`simple contrasts for Agua`
## Luz = 10%:
##  contrast  estimate    SE  df z.ratio p.value
##  50% - 70%  -0.8270 0.282 Inf  -2.933  0.0094
##  50% - 90%  -0.0848 0.250 Inf  -0.340  0.9384
##  70% - 90%   0.7422 0.284 Inf   2.610  0.0246
## 
## Luz = 40%:
##  contrast  estimate    SE  df z.ratio p.value
##  50% - 70%  -1.2022 0.272 Inf  -4.422 <0.0001
##  50% - 90%  -1.3616 0.286 Inf  -4.766 <0.0001
##  70% - 90%  -0.1594 0.320 Inf  -0.498  0.8723
## 
## Luz = 70%:
##  contrast  estimate    SE  df z.ratio p.value
##  50% - 70%  -0.3765 0.506 Inf  -0.744  0.7373
##  50% - 90%   1.1519 0.389 Inf   2.959  0.0087
##  70% - 90%   1.5283 0.441 Inf   3.468  0.0015
## 
## Luz = 100%:
##  contrast  estimate    SE  df z.ratio p.value
##  50% - 70%   0.1323 0.289 Inf   0.459  0.8906
##  50% - 90%  -0.8605 0.356 Inf  -2.415  0.0417
##  70% - 90%  -0.9928 0.352 Inf  -2.822  0.0132
## 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

Revisió dels supòsits del model

Igual que en exemples anteriors, es poden denerar gràfics de valors residuals amb la funció plot. En aquest cas, fo farem amb autoplot que s’executa amb ggplot2 i ggfortify. El resultat és gairebé idèntic al que s’obté amb plot.glm().

library(ggfortify)

glm.ailant %>% 
  autoplot() +
  theme_bw()

Els gràfics no suggereixen cap patró especial els valors residuals obtinguts. Tanmateix si observem el summary, veiem que Residual deviance és 101.00 amb 36 graus de llibertat, que significa una ratio de 2.80.Cosa que podria representar sobredispersió.

Realitzarem per tant el test de sobredispersió amb la funció gof() de aods3

aods3::gof(glm.ailant)

## D  = 100.9971, df = 36, P(>D) = 4.422451e-08
## X2 = 87.0328, df = 36, P(>X2) = 4.101207e-06

El test indica que, efectivament, hi ha sobredispersió en els residuals de l’ajustament. El problema podria ser degut a problemàtiques d’autocorrelació en les mostres o algun tipus colinealitat per algun factor que no hem sabut controlar en el procediment experimental.

L’alternativa és en aquest cas és corregir els errors estàndard dels coeficients estimats ajustant una distribució quasi-binomial. Recordem que això no canviarà els patrons observats en els gràfics de residus, ni la partició de la deviança, peò sí els errors estàndard dels coeficients i els test de significació de \(\chi^2\) sobre cada un dels termes del model.

glm.ailant.2 <- glm(yprop ~ Luz*Agua, weights = Semillas, data = ailanto, family = quasibinomial)

car::Anova(glm.ailant.2, type= 3)

## Analysis of Deviance Table (Type III tests)
## 
## Response: yprop
##          LR Chisq Df Pr(>Chisq)    
## Luz       23.6500  3  2.956e-05 ***
## Agua       4.3524  2   0.113472    
## Luz:Agua  21.1036  6   0.001758 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GLM amb distribucions d’errors binomials en ecologia

Notes d’estadística aplicada a les ciències ambientals: Nota Núm. 5

Stefano Sanfilippo

Abril 2026

Introducció

Referències

Models binomials o models de regressió logística

GLM de resposta binària

Cas en estudi: l’estratègia d’una alga invasora

Anàlisi exploratòria

Ajustament i resolució d’hipòtesi

Interpretació del model

La funció summary()

Predicció de valors

Possibles valors de resposta segons canvis variables explicatives

Nivell efectiu de la mediana

C. Fragile (nivell de referència)

C. Tomentosum

C. Vermilara

Deviança explicada pel model i capacitat predictiva

Taula de classificació

La corba ROC

Gràfics de predicció

Tests post hoc

Revisió dels supòsits del model

Gràfic Residuals vs Fitted: residuals estandarditzats (Pearson) vs. valors esperats

Gràfic Normal Q-Q

La llibreria DHARMa

GLM de resposta de tipus proporció

Cas d’estudi: la germinació de l’ailant

Ajustament del model

Mètode amb weights

Mètode unint èxit i fracàs amb rbind()

Contrast d’hipòtesi

Interpretació del model

La funció summary

Representació del model

Deviança explicada pel model

Prediccions

Taula de valors mitjans

Taula de valors mitjans i intervals de confiança

Test post hoc

Revisió dels supòsits del model

La funció `summary()`

Mètode amb `weights`

Mètode unint èxit i fracàs amb `rbind()`

La funció `summary`