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

Introducció

Quan una variable de resposta és de tipus recompte, la distribució d’errors que se sol utilitzar és la distribució d’error de Poisson. Freqüentment s’esdevé que la variança creix més ràpid que la mitjana: fenomen que s’anomena sobredispersió. La família d’errors binomial negativa permet poder gestionar estructures d’errors sobredisperses per a variables tipus recompte.

GLM amb distribució d’errors Poisson

Aquest exemple presenta dades d’un estudi de l’efecte de la ramaderia i el règim de tenència de la terra, sobre la regeneració de l’araucària (Araucaria araucana). El treball mostreja les plàntules d’araucària a 26 parcel·les de 2000 m\(^2\), petits propietaris que feien un ús més intensiu de la terra, i grans empreses d’explotació forestal que deixaven a l’interior algun remanent de bosc nadiu.

La hipòtesi principal és que la ramaderia afecta negativament a la regeneració de l’araucària, però que aquest efecte podria diferir in funció de l’ús que es faci de la terra.

Aquest efecte diferencial de la ganaderia sobre la regeneració, depenent del tipus de tinència de la terra, quedaria representat per la interacció entre variable i factor. Les dades estan disponibles al dataset ara del paquet ADER:

library(ADER)
data(ara)
str(ara)

## 'data.frame':    26 obs. of  3 variables:
##  $ seedlings: int  7 10 0 5 54 0 30 14 1 2 ...
##  $ dungs    : int  5 0 330 960 115 2330 475 10 2000 130 ...
##  $ property : Factor w/ 2 levels "Campesino","Empresa Forestal": 2 2 2 2 2 2 2 2 2 1 ...

Les tres hipòtesis nul·les serien:

\(H_0:\beta=0\): la pendent de la relació entre la regeneració (és a dir la quantitat de plàntules, seedlings) i la intensitat de la ramaderia, mesurada pel nombre de deposicions, dungs.
\(H_0: \mu_1=\mu_2\): la regeneració mitjana de plàntules (seedlings) és similar en parcel·les gestionades per petits propietaris (Campesino) o per empreses forestals (Empresa forestal), property.
\(H_0: \beta_1=\beta_2\): la relació entre regeneració (seedlings) i intensivitat de la ramaderia (dungs) és igual per a petits propietaris (Campesino) i empresa forestal (Empresa forestal).

Anàlisi exploratòria

Com ja sabem, és important, abans d’ajustar qualsevol model a les dades, explorar-les gràficament.

library(patchwork)

a <- ara %>% 
  ggplot(aes(x=dungs,y=seedlings, colour = property)) +
  geom_point(size= 2) +
  facet_wrap(~ property) +
  labs(title = "a) Sense funció de vincle") +
  theme_bw() +
  theme(legend.position = "none")

b <- ara %>% 
  ggplot(aes(x=dungs,y=log(seedlings), colour = property)) +
  geom_point(size= 2) +
  facet_wrap(~ property) +
  labs(title = "b) Amb funció de vincle logarítmica") +
  theme_bw() +
  theme(legend.position = "none")

a / b +
  plot_annotation(title = "Relació entre intensivitat ramadera i regeneració de plàntules \ncondicionada pel règim de propietat",
                  theme = theme(plot.title = element_text(
                    colour = "dodgerblue3", hjust = 0.5, size = 14)))

Allò primer que podem observar en a) és que la relació entre variable de resposta i les explicatives no apareix linial; més aviat és una exponencial negativa. Això implicarà utilitzar una funció de vincle logarítmica, ja que ajudarà a linialitzar la relació, tal com es pot intuir observant b).

Ajustament del model i inferència en GLM

El model que es faria servir seria un ANCOVA, ja que hi conviuen una variable explicativa continua i una que es factorial. El model podria ser així:

glm.ara.interacc <- glm(seedlings ~ dungs + property + dungs:property,
                        data = ara, family = poisson)

Inferència: ratio de versemblança \(\chi^2\)**

La funció Anova() del paquet car utilitza el LR Chisq, que significa Likelihood-Ratio Chi-Square (ratio de versemblança \(\chi^2\)). Aquest estadístic representa la diferència en l’ajust del model entre el model complet i un model restringit que elimina un predictor específic, cosa que permet comprovar si la inclusió d’aquest predictor millora significativament el model.

Taula ANOVA

La taula ANOVA que genera la funció conté:

L’estadístic: El valor LR Chisq (sovint anomenat desviació o valor \(\chi^{2}\)) és l’estadístic de la prova. Com més alt sigui el valor, més desviació explica el predictor i més forta serà l’evidència que la variable és important.
Df (graus de llibertat): Representa el nombre de paràmetres associats amb aquest predictor específic que s’estan provant.
Pr(>Chisq): Aquest és el valor p. Respon si el valor LR Chisq és prou gran per ser estadísticament significatiu. Si aquest valor és inferior al nivell de significació escollit (per exemple, 0.05), significa que el predictor contribueix significativament al model.

Una altra qüestió important és sobre el tipus de sumes de quadrats a utilitzar en aquesta prova. Atesa certa previsible colinealitat i el disseny desequilibrat (només 9 observacions de Empresa forestal sobre 26) serà més apropiat utilitzar les sumes de quadrats de tipus II.

car::Anova(glm.ara.interacc, type=2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: seedlings
##                LR Chisq Df Pr(>Chisq)    
## dungs            68.799  1  < 2.2e-16 ***
## property         80.598  1  < 2.2e-16 ***
## dungs:property   37.974  1   7.17e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Tenir en compte en tot cas que aquesta funció executa el type=2 per defecte.

Per al nostre model, tant les variables explicatives com llur interacció tenen un efecte significatiu sobre la variable de resposta. Això significa que l’efecte de la menor o major intensivitat ramadera té un efecte sobre la regeneració de l’araucària que canvia entre petites propietats i empreses forestals.

Interpretació del model

La funció del model

L’output de summary() en glm() és bastant semblant al de lm().

summary(glm.ara.interacc)

## 
## Call:
## glm(formula = seedlings ~ dungs + property + dungs:property, 
##     family = poisson, data = ara)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     2.497423   0.135408  18.444  < 2e-16 ***
## dungs                          -0.010211   0.001789  -5.706 1.15e-08 ***
## propertyEmpresa Forestal        0.576261   0.170467   3.380 0.000724 ***
## dungs:propertyEmpresa Forestal  0.009041   0.001803   5.015 5.31e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 347.56  on 25  degrees of freedom
## Residual deviance: 198.02  on 22  degrees of freedom
## AIC: 282.3
## 
## Number of Fisher Scoring iterations: 5

La part determinista del model s’expressaria amb l’equació \[ \small\log (\text{regeneració})=2.497+0.576~\text{empr.for}+ (-0.010+0.009\text{ empr.for})\text{ depos}. \]

L’ordenada en l’origen (2.497) resumeix la mitjana del logaritme de la regeneració en terres de petits propietaris en absència de ramaderia.
L’efecte de la ramaderia és negatiu, pero la pendent disminueix (-0.10+0.09), tornant-se menys negativa en terres d’empreses forestals.
Aquests coeficients fan referència al predictor linial \(\eta\): per predir la resposta haurem d’emprar la funció de vincle inversa.
A diferència dels models linials, els tests de Wald (el coeficient entre cada coficient l’error estàndard) es rotulen z test, ja que en aquest cas es comparen el valor es compar amb els de la distribució normal.

Deviança aplicada pel model

Ja sabem (Nota 18) que la deviança és \[ \small D^2 =1-\frac{D}{D_{null}} \]

on \(D\) és la deviança residual del model i \(D_{null}\) és la deviança d’un model on tan sol s’ajusta l’ordenada d’origen. Al summary les trobem com Residual deviance i Null deviance, respectivament.

Aquest valors es poden obtenir directament amb la funció deviance()

1 - deviance(glm.ara.interacc) / deviance(update(glm.ara.interacc, .~1))

## [1] 0.4302522

Aquesta explicaria el 43% de la variabilitat total observada en el nombre de plàntules.

Gràfics de predicció

Quan la variable explicativa és continua, se sol mostrar la recta de regressió que representa els valors esperats pel model sobre els gràfics de dispersió de les ambdues variables, podentnt-hi afegir bandes d’intervals de confiança sobre la/les línea/es de predicció. Naturalment tot això es complica quan tenim vàries variables explicatives.

L’estratègia seria veure com canvia la variable de resposta amb cada una de les variables explicatives mantenint la resta constants.

Sense intervals de confiança, amb `ggplot2`

new_data <- list()

for (i in 1:length(levels(ara$property))) {
  new_data[[i]] <-  data.frame(seedlings = seq(range(ara$seedlings)[1],
                                   range(ara$seedlings)[2], length=50),
                               dungs= seq(range(ara$dungs)[1],
                                   range(ara$dungs)[2], length=50),
                               property = rep(levels(ara$property)[i], 50)
  )
  }

new_data <- bind_rows(new_data)

# Important calcular la inversa de la predicció ja que és log(y)=...
# L'alternativa és posar el paràmetre 'type = "response"'.
new_data$predict <- exp(predict(glm.ara.interacc,
                            newdata = new_data))

ara %>% 
  ggplot(aes(x=dungs, y= seedlings, colour = property, linetype = property)) +
  geom_point() +
  geom_line(data = new_data, aes(x=dungs, y=predict, colour = property)) +
  labs(title = "Predicció de l'efecte de la ramaderia en la regeneració \nd'araucària segons el règim de tenència de la terra",
       x= "Nombre de deposicions",
       y= "Nombre de plàntules") +
  theme_bw() +
  theme(plot.title = element_text(size = 12, hjust = 0.5))

Amb intervals de confiança, amb `visreg`

scale: Per defecte, el model es representa en l’escala del predictor lineal. Si scale=‘response’ per a un glm, s’aplicarà la funció d’enllaç inversa (exp, que hem emprat abans per representar el gràfic manualment) de manera que el model es representi en l’escala de la resposta original.

library(visreg)
visreg(fit=glm.ara.interacc, 
       xvar = "dungs", 
       by= "property", 
       overlay = T,
       scale = "response")

Amb intervals de confiança amb `ggplot2`

nd <- list()

for (i in 1:length(levels(ara$property))) {
  nd[[i]] <-  data.frame(seedlings = seq(range(ara$seedlings)[1],
                                   range(ara$seedlings)[2], length=100),
                               dungs= seq(range(ara$dungs)[1],
                                   range(ara$dungs)[2], length=100),
                               property = rep(levels(ara$property)[i], 100)
  )
  }

nd <- bind_rows(nd)

# No cal calcular la inversa de la predicció, ja ho fa predict(tipe="response")

nd.camp.predict <- as.data.frame(predict(glm.ara.interacc,
        newdata = nd[nd$property=="Campesino",],
        type = "response",
        se= T)) %>% 
  mutate(property= "Campesino") %>% 
  select(1,2,4)
nd.empr.predict <- as.data.frame(predict(glm.ara.interacc,
        newdata = nd[nd$property=="Empresa Forestal",],
        type = "response",
        se= T)) %>% 
  mutate(property= "Empresa Forestal") %>% 
  select(1,2,4)

nd.predict <- bind_rows(nd.camp.predict,nd.empr.predict)


ara %>% 
  ggplot(aes(x=dungs, y= seedlings, colour = property)) +
  geom_point() +
  geom_line(data = nd.predict, aes(x=nd$dungs, y=fit), linewidth=0.9) +
  geom_ribbon(data = nd.predict, 
              mapping = aes(x= nd$dungs,
                            y=fit,
                            ymin= (fit - 1.96*se.fit), 
                            ymax = (fit + 1.96*se.fit)), 
              alpha =0.2, fill = "dodgerblue", outline.type = 'full')+
  labs(title = "Predicció de l'efecte de la ramaderia en la regeneració \nd'araucària segons el règim de tenència de la terra",
       x= "Nombre de deposicions",
       y= "Nombre de plàntules") +
  theme_bw() +
  theme(plot.title = element_text(size = 12, hjust = 0.5),
        legend.position = "bottom")

Revisió de supòsits

par(mfrow=c(2,2))
plot(glm.ara.interacc)

En aquest cas, tractant-se d’una distribució de Poisson, més que anar a controlar normalitat i homoscedasticitat hem de veure com creix la dispersió respecte a la mitjana. Aixó es pot veure amb Residuals vs Fitted, no es pot observar que la dispersió aumenta més que la mitjana: la qual cosa ens suggereix sobredispersió.

Sobredispersió i subdispersió

Quan es defineix una distribució d’errors Poisson, assumeix que la vanriança és igual que la mitjana. Però es pot esdevenir que la variança sigui major o menor que la mitjana: aquests fenòmens s’anomenen sobredispersió i subdispersió.

La sobredispersió pot aparèixer com a conseqüència de no aver considerat en el model o o més variables explicatives.
La subdispersió pot aparèixer, entre altres circumstàncies, quan existeixen estructures d’autocorrelació no definides per les dades, és a dir grups d’observacions correlacionades entre si.

La sobredispersió

Sembla ser un problema força comú trballant amb models Poisson de dades ecològiques.

Si no podem definir alguna variable explicativa nova que que redueixin l’increment de la variança residual,la solució és especificar una estructura de variança que inclogui la sobredispersió.

Existexen 2 alternatives:

1. Utilitzar una distribució quasi-Poisson

Tot i que en realitat no existixi una família d’errors amb aquest nom, sinó que es tracta d’una manera d’ajustar la variança, al cap i a la la fi, una distribució d’errors. És una distribució d’errors de Poisson on es calcula a posteriori la dispersió \(\phi\) que corregeix la variança residual de la forma: \[ \small Var(y|x)=\phi\mu \]

Corregeix els errors estàndard dels coeficients estimats, però no els coeficients.

El primer problema que presenta el màtode és que en R no presenta valors d’\(AIC\). Això perquè la quasi-Poisson no és ben bé una família de distribució.

En R, aquest model s’obté especificant l’argument family=quasipoisson dins de la funció glm().

glm.ara.q <- glm(seedlings ~ dungs + property + dungs:property,
                        data = ara, family = quasipoisson())

summary(glm.ara.q)$dispersion

## [1] 11.29489

El paràmetre de dispersió és doncs \(\phi=11.295\) i indica que per a cada canvi unitari de la mitjana, la variança creix més d’11 unitats.

Si mirem el sesum del model executant summary(),

summary(glm.ara.q)

## 
## Call:
## glm(formula = seedlings ~ dungs + property + dungs:property, 
##     family = quasipoisson(), data = ara)
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     2.497423   0.455076   5.488 1.63e-05 ***
## dungs                          -0.010211   0.006014  -1.698    0.104    
## propertyEmpresa Forestal        0.576261   0.572902   1.006    0.325    
## dungs:propertyEmpresa Forestal  0.009041   0.006059   1.492    0.150    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.29489)
## 
##     Null deviance: 347.56  on 25  degrees of freedom
## Residual deviance: 198.02  on 22  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5

Es pot observar que el quocient entre “Residual deviance” i graus de llibertat no canvia. Tampoc canvien els gràfics dels residuals, com tampoc els coeficients estimats. Però sí, canvien els errors estàndard associats als coeficients que augmenten considerablement.

També aumenten considerablement els tests de \(\chi^2\) de l’ANOVA test (perquè es divideix la deviança explicada per de cada terme pel paràmetre \(\phi\), de dispersió, en aques cas, 11.29489). La qual cosa pot tenir implicacions en la inferència, com en aquest cas, on la interacció dungs:property deixaria de ser significativa.

car::Anova(glm.ara.q, type=2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: seedlings
##                LR Chisq Df Pr(>Chisq)   
## dungs            6.0911  1   0.013586 * 
## property         7.1358  1   0.007556 **
## dungs:property   3.3620  1   0.066716 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. Utilitzar una distribució d’errors binomial negativa (cap. següent)

Aquest familia de distribució d’errors compta també amb un paràmetre de dispersió o escalat \(\phi\) que corregeix la variança residual. La binomial negativa és una família de distribució veritable, \(\phi\) forma part del model i no es calcula ad hoc. Per la mateixa raó, permet calcular l’AIC, per si l’hem d’emprar en el context de la selecció de models.

Els valors residuals en la funció `glm()`

La funció glm() utilitza deiversos tipus de càlcul de valors residuals:

Residual brut

És la diferència entre valor observat i valor i valor esperat \[ \small residual~brut=y_i-\hat y_i \]

residuals(glm.ara.interacc, type = "response")

Residual temporal

Els residuals temporals estandarditzen el residual brut dividint-lo pel valor esperat, \(\hat y\) \[ \small residual~temporal= \frac{y_i-\hat y_i}{\hat y_i} \]

residuals(glm.ara.interacc, type = "working")

També per defecte, és

glm.ara.interacc$residuals

Residual de Pearson

Els residuals de Pearson són una versió reescalada dels residuals temporals i estandarditzen el residu brut dividint-lo per l’arrel quadrada de \(\hat y_i\): \[ \small residual~de~Pearson= \frac{y_i-\hat y_i}{\sqrt{\hat y_i}} \]

residuals(glm.ara.interacc, type = "pearson")

Residuals de deviança

Els residuals de deviança escalen els residuals bruts en funció de la contribució de cada observació a la deviança total són els residuals que s’utilitzen més sovint en l’avaluació dels supòsits de GLM. En models de Poisson, aquests residus es poden calcular de la següent manera \[ \small residual~de~deviança= \sqrt{2·(y_i·\log(\frac{y_i}{\hat y_i})-(y_i-\hat y_i))} \]

Per defecte, la funció residuals processa aquesta opció per a ajustos GLM

residuals(glm.ara.interacc)

La subdispersió

És un problema menys comú en biologia.

Les dues formes més comunes per gestionar aquest problema són:

Intentar determinar la causa de la subdispersió (normalment algun factor de correlació entre dades) i seleccionar l’estructura del model més adequada.
Utilitzar models de Poisson generalitzats (GP). Igual que els binomials negatius disposen d’un paràmetre de dispersió o d’escala que redueix la variança. La funció genpoisson() de la liberia VGAM permet ajustar models GP.

Com reconèixer la sobredispersió i la subdispersió

La sobredispersió i la subdispersió es poden detectar calculant el quocient entre deviança residual (Residual deviance) i el nombre de graus de llibertat. Aquest quocient hauria de ser a prop de 1 si la família de distribució d’errors representa bé la variança residual.

També existeixen algunes funcions específiques de liberies diverses que poden ajudar en aquest labor d’avaluació.

GLM amb distribució d’errors binomial negativa

Una de les alternatives més recomanades per a gestionar els problemes de sobredispersió és la distribució d’errors binomial negativa.

La funció de R que executa aquesta distribució es glm.nb() de MASS.Hi ha diverses fórmules que relacionen la varianza amb la mitjana. La funció glm.nb() utilitza: \[ \small Var(y|x)=\mu+\mu^2/\phi \]

Ajustament i `summary`

Reproduim doncs l’ajustament amb la variable de resposta seedlings, el factor property i la covariable dungs.

glm.ara.nb <- MASS::glm.nb(seedlings ~ dungs + property + dungs:property, ara)

Observem després el compendi que ens proporciona summary. Veiem que tant els coeficients estimats com els errors estàndard ja han canviat respecte a l’ajustament quasiPoisson.

summary(glm.ara.nb)

## 
## Call:
## MASS::glm.nb(formula = seedlings ~ dungs + property + dungs:property, 
##     data = ara, init.theta = 1.13715532, link = log)
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                     2.152593   0.362295   5.942 2.82e-09 ***
## dungs                          -0.005663   0.002468  -2.294   0.0218 *  
## propertyEmpresa Forestal        1.113371   0.557993   1.995   0.0460 *  
## dungs:propertyEmpresa Forestal  0.003945   0.002528   1.561   0.1186    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.1372) family taken to be 1)
## 
##     Null deviance: 49.506  on 25  degrees of freedom
## Residual deviance: 28.040  on 22  degrees of freedom
## AIC: 155.22
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.137 
##           Std. Err.:  0.384 
## 
##  2 x log-likelihood:  -145.218

El paràmetre de dispersió \(\phi\) (que summary indica com theta) és \(\phi=1.137\), mentre el quocient entre Residual deviance i graus de llibertat ha disminuït molt: ara és 1.27. Tot això indica que hem corregit gran part de la sobredispersió.

Anàlisi ANOVA

car::Anova(glm.ara.nb, type = 2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: seedlings
##                LR Chisq Df Pr(>Chisq)    
## dungs           12.9047  1  0.0003278 ***
## property        11.6580  1  0.0006393 ***
## dungs:property   3.5106  1  0.0609766 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Veiem que en aquest model la interacció entre covariable i factor ha perdut significació. Caldrà per tant reduir el model, tornant a ajustar-lo sense la interacció covariable:factor.

glm.ara.nb2 <- MASS::glm.nb(seedlings ~ dungs + property, ara)
summary(glm.ara.nb2)

## 
## Call:
## MASS::glm.nb(formula = seedlings ~ dungs + property, data = ara, 
##     init.theta = 0.9538647699, link = log)
## 
## Coefficients:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               1.8440677  0.2836331   6.502 7.95e-11 ***
## dungs                    -0.0020671  0.0006192  -3.338 0.000842 ***
## propertyEmpresa Forestal  1.5951773  0.5003644   3.188 0.001432 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.9539) family taken to be 1)
## 
##     Null deviance: 43.340  on 25  degrees of freedom
## Residual deviance: 27.748  on 23  degrees of freedom
## AIC: 156.43
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.954 
##           Std. Err.:  0.306 
## 
##  2 x log-likelihood:  -148.426

On veiem que les dades de dispersió milloren encara més i en l’anàlisi de la variança ara tots els coeficients presenten valor \(p\) significatiu.

car::Anova(glm.ara.nb2, type = 2)

## Analysis of Deviance Table (Type II tests)
## 
## Response: seedlings
##          LR Chisq Df Pr(>Chisq)    
## dungs      11.271  1  0.0007875 ***
## property   10.063  1  0.0015130 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Gràfics i tests de dispersió

Els gràfics ara mostren una millores substancial de la dispersió que disminueix clarament: fixar-se com disminueix el rang dels gràfics. També han desaparegut els valors palanca que havia al gràfic Residual vs. Leverage.

par(mfrow= c(2,2))
plot(glm.ara.nb2)

Podem també dur a terme el test sobredispersió amb la funció gof() de aods3.

aods3::gof(glm.ara.nb2)

## D  = 27.7484, df = 23, P(>D) = 0.2254439
## X2 = 27.9001, df = 23, P(>X2) = 0.2195832

El test és positiu, amb uns valors \(p\) molt alts, que confirmen que aquesta distribució captura més correctament l’estructura de dispersió de les nostres dades al voltant de les mitjanes.

GLM per a recomptes en ecologia

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

Stefano Sanfilippo

Gener 2026

Referències

Introducció

GLM amb distribució d’errors Poisson

Anàlisi exploratòria

Ajustament del model i inferència en GLM

Inferència: ratio de versemblança \(\chi^2\)**

Taula ANOVA

Interpretació del model

La funció del model

Deviança aplicada pel model

Gràfics de predicció

Sense intervals de confiança, amb `ggplot2`

Amb intervals de confiança, amb `visreg`

Amb intervals de confiança amb `ggplot2`

Revisió de supòsits

Sobredispersió i subdispersió

La sobredispersió

1. Utilitzar una distribució quasi-Poisson

2. Utilitzar una distribució d’errors binomial negativa (cap. següent)

Els valors residuals en la funció `glm()`

Residual brut

Residual temporal

Residual de Pearson

Residuals de deviança

La subdispersió

Com reconèixer la sobredispersió i la subdispersió

GLM amb distribució d’errors binomial negativa

Ajustament i `summary`

Anàlisi ANOVA

Gràfics i tests de dispersió

GLM per a recomptes en ecologia

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

Stefano Sanfilippo

Gener 2026

Referències

Introducció

GLM amb distribució d’errors Poisson

Anàlisi exploratòria

Ajustament del model i inferència en GLM

Inferència: ratio de versemblança \(\chi^2\)**

Taula ANOVA

Interpretació del model

La funció del model

Deviança aplicada pel model

Gràfics de predicció

Sense intervals de confiança, amb ggplot2

Amb intervals de confiança, amb visreg

Amb intervals de confiança amb ggplot2

Revisió de supòsits

Sobredispersió i subdispersió

La sobredispersió

1. Utilitzar una distribució quasi-Poisson

2. Utilitzar una distribució d’errors binomial negativa (cap. següent)

Els valors residuals en la funció glm()

Residual brut

Residual temporal

Residual de Pearson

Residuals de deviança

La subdispersió

Com reconèixer la sobredispersió i la subdispersió

GLM amb distribució d’errors binomial negativa

Ajustament i summary

Anàlisi ANOVA

Gràfics i tests de dispersió

Sense intervals de confiança, amb `ggplot2`

Amb intervals de confiança, amb `visreg`

Amb intervals de confiança amb `ggplot2`

Els valors residuals en la funció `glm()`

Ajustament i `summary`