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

Introducció

En els models lineals, les variables explicatives es combinen linealment per modelar el seu efecte sobre la variable de resposta. De vegades en tot cas la linealitat és insuficient per a capturar l’estructura de les dades. Una alternativa que ja coneixem és la transformació de la variable de resposta, o l’ús de funcions de vincle, incorporades als models lineals generalitzats per tal de linealitzar relacions que lineals no són. Una altra alternativa pot ser representada per la transformació d’una o més variables explicatives (per exemple, termes quadràtics o cúbics dins d’aquestes).

Però en certes ocasions aquestes transformacions resulten insuficients i es necessitarà d’un marc de modelació més flexible. Certs mètodes com els mètodes additius generalitzats, la regressió no lineal, els arbres de classificació i regressió (CART) o les xarxes neuronals, aporten alternatives molt versàtils per a la modelació de relacions no lineals.

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.

Dades

En aquest capítol farem servir el data frame dry de la llibreria ADER que mostra la relació entre richness, riquesa d’arbres, i NDVI (Normalized Difference Vegetation Index), que podem definir com a productivitat. Dades obtingudes amb imatges de satèl·lit , preses entre 2001 i 2007, sobre parcel·les de 0.1 ha de bosc tropical a Mèxic. El valor de NDVI oscil·la entre 0 i 1, indicant els valor més alts més alta productivitat.

library(ADER)
data('dry')
str(dry)
## 'data.frame':    96 obs. of  2 variables:
##  $ richness: int  18 22 20 24 19 13 11 10 16 15 ...
##  $ ndvi    : num  0.802 0.792 0.792 0.756 0.695 ...

Models additius generalitzats (GAM)

Descripció teòrica general

Quan les dades no mostren una clara relació lineal, i els models polinòmics són insuficients, els models additius generalitzats poden ser una bona alternativa als LM i els GLM. Els GAM utilitzen corbes de suavització per modelar la relació resposta/explicativa, tot permetent relacions no lineals.

La suavització

La suavització es pot considerar un entremig entre la recta de regressió i una corba que passi per tots el punts del gràfic de dispersió.

library(patchwork)

a <- dry %>% 
  ggplot(aes(x=ndvi, richness)) +
  geom_point()+
  geom_smooth(method = "lm", se = F, color="blue", linewidth= 0.5) +
  labs(title = "a.") +
  theme_bw()

b <- dry %>% 
  ggplot(aes(x=ndvi, richness)) +
  geom_point()+
  geom_line(color="blue") +
  labs(title = "b.") +
  theme_bw()


a+b

Mentre la recta de regressió (a.) pot arribar a ser exageradament simplista per a representar la relació entre riquesa d’espècies i NDVI, la segona opció (b.) tampoc ens aporta informació. Entremig d’aquestes dues aproximacions, podríem traçar un ajustament no lineal que descrigués com varia \(y\) en front d’\(x\).

La suavització permet obtenir una funció no lineal a través d’estimacions de la variable de resposta, \(y\), per a rangs parcials de valors de \(x\). Els valors estimats d’\(y\) per a cada rang s’uneixen després per formar una recta no lineal.

Existeixen diferents mètodes de suavització. Els dos que se solen fer servir més són les Splines i la regressió local.

Splines

Una spline és una regressió polinòmica per trams, on es divideixen les dades d’acord als valors de la variable explicativa i a cada tram de dades s’ajusta una regressió polinòmica normalment de grau baix, posant com a condició que els extrems de cada tram successiu coincideixin. Normalment es fan servir polinomis de grau 3.

Regressió local o regressió local ponderada

La regressió local (LOESS) ajusta múltiples corbes de regressió lineal tenint en compte a cada punt només un subconjunt de dades que siguin pròximes a la dada focal. També aquí caldrà definir quants punts pròxims tindrem en compte per fer l’ajustament.

També es pot assignar més pes als punts més pròxims i menys pes als més allunyats mitjançant ponderació lineal: mètode que es coneix com regressió local ponderada (LOWESS).

El model additiu

Els models GAM poden considerar-se extensions de models GLM en què s’adscriguin funcions no lineals per a variables explicatives individuals, tot mantenint l’additivitat dels efectes. És a dir, el predictor lineal \(\eta\) \[ \small \eta_i =\beta_0+\beta_{1}x_{i1}+\beta_2x_{i2}+...+\beta_p x_{ip} + \epsilon_i \] es converteix ara en \[ \small \eta_i = \beta_0 + f_1(x_{i1}) + f_2(x_{i2})+...+ f_p(x_{ip}) + \epsilon_i \] on cada component lineal \(\small \beta_{ij}\) és substituït per una funció no lineal de suavització, \(\small f_j(x_{ij})\). Aquest model es defineix com additiu perquè s’ajusta una funció individual per a cada variable i, posteriorment, se’n sumen els efectes.

Les funcions poden ser de qualsevol tipus (splines, LOESS, LOWESS, regressions polinòmiques o lineals). Dins del mateix model es poden plantejar funcions diferents per a cada una de las variables explicatives. Per exemple, es poden plantejar relacions lineals amb algunes i no lineals amb d’altres mitjançant l’ajustament de funcions, per exemple \[ \small \eta_i= \beta_0 + \beta_1 x_{1i} + f_2 (x_{2i}) + \epsilon_i \]

NOTA: Les funcions de suavització només es poden aplicar a variable contínues, mai a les categòriques.

a. Ajustament amb la llibreria gam

Amb mètode splines

El paquet gam disposa de funcions per ajustar i treballar models additius generalitzats. La funció gam() d’aquest model permet executar suavitzacions pel mètode splines i per LOESS.

Per a implementar una suavització amb splines podem fer servir la funció s() que ajusta splines suavitzadores sobre la variable explicativa dins la fórmula de gam().

library(gam)
gam1 <- gam::gam(richness ~ gam::s(ndvi), data = dry, family = poisson)

summary(gam1)
## 
## Call: gam::gam(formula = richness ~ gam::s(ndvi), family = poisson, 
##     data = dry)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8649 -1.8516  0.2694  1.4236  3.1329 
## 
## (Dispersion Parameter for poisson family taken to be 1)
## 
##     Null Deviance: 353.1647 on 95 degrees of freedom
## Residual Deviance: 350.4069 on 94 degrees of freedom
## AIC: 762.348 
## 
## Number of Local Scoring Iterations: 4 
## 
## Anova for Parametric Effects
##              Df  Sum Sq Mean Sq F value Pr(>F)
## gam::s(ndvi)  1   2.768  2.7679  0.8243 0.3662
## Residuals    94 315.638  3.3579

La funció summary ofereix un resum de l’ajustament, incloent-hi els tests de significació que obtindríem amb anova().

Interpretació del resultat

La funció summary.Gam() retorna dues taules ANOVA: una per a la part lineal i una altra per a la part aproximada de manera no paramètrica. A la primera s’inclou també la part lineal del termes no lineals, amb un grau de llibertat, mentre que la resta de graus de llibertat s’inclouen a la segona taula. So no s’han especificat, s() aplica per defecte 4 graus de llibertat. D’acord amb el test del quocient de versemblança (Anova for Nonparametric Effects) existe un efecto de l’NDVI sobre la riquesa d’arbres en boscos tropicals secs.

Però, com és aquest efecte? El model que ajustat és \[ \small \eta_i=\beta_0 + s(NDVI_i) + \epsilon_i \] sent \(\small s(NDVI_i)\) una funció no lineal (amb un nombre variable de nusos) de la variable explicativa NDVI.

Tot i que podem conèixer els coeficients del model

coefficients(gam1)
##  (Intercept) gam::s(ndvi) 
##    2.1616869    0.6209708

l’efecte de l’NDVI no es pot interpretar en funció del signe i de la magnitud del coeficient, com faríem en el cas de LMs o GLMs: és necessari en canvi representar gràficament les prediccions del model.

La funció plot.Gam per defecte genera un gràfic que representa la spline en funció de la variable explicativa (gràfic de l’esquerra), que equival a representar les prediccions del model, però a diferent escala (veure gràfic de la dreta que, a més a més, s’ha dibuixat amb ggplot2). La spline està en l’escala del predictor lineal. Per tant, si volem conèixer el valor esperat de riquesa per a un valor d’NDVI concret, haurem d’agafar el valor corresponent de la spline, sumar-lo a l’ordenada a l’origen i aplicar la funció de vincle inversa (en aquest cas la funció exponencial) a aquesta suma.

Aquestes operacions es poden obviar aplicant el paràmetre type = "response" quan executem la funció predict.

library(grid)

pred <- as.data.frame(predict(gam1, type = "response", se.fit = T)[1:2]) %>% 
  mutate(upper = fit + 1.96*se.fit,
         lower= fit - 1.96*se.fit)

# Gràfic de la dreta
 P <- dry %>% 
  ggplot(aes(x=ndvi,y=richness))+ 
  geom_point()+
  geom_line(data=pred,
        aes(x=dry$ndvi,y=fit)) +
  geom_line(data=pred,
        aes(x=dry$ndvi,y=upper), linetype=3) +
  geom_line(data=pred,
        aes(x=dry$ndvi,y=lower), linetype=3) +
  theme_bw()

 vp.BottomRight <- viewport(height=unit(.85, "npc"), width=unit(0.47, "npc"), 
                           just=c("left","top"), 
                           y=0.96, x=0.53)
 
par(mar= c(5,4,1.2,18))
# Gràfic de l'esquerra
plot(gam1,se=T)

print(P, vp=vp.BottomRight)

En tot cas, el model sembla suggerir que, en ambient tropical sec, el valor de la riquesa d’arbres sigui màxim per a valor extrems de l’escala de de NDVI.

Variabilitat explicada pel model

Aplicant la fórmula \(\small D^2=1-D/D_{nul}\) que ja coneixem i, basant-nos en els valors del summary, obtenim

1-289.1961/353.1647
## [1] 0.1811297

és a dir un 18% aproximat de la variabilitat observada. Cayuela i De la Cruz consideren aquest ordre de magnitud com acceptable, atesa la baixa resolució de les dades espacials. Tot i així, és possible que que incloent més variables al model es pugui arribar a explicar més variabilitat de la variable de resposta.

Els graus de llibertat en el procés de suavització

En la suavització per splines, el nombre de graus de llibertat efectius (edf) per defecte en la funció és 4. Es pot canviar el nombre de graus de llibertat especificant-lo dins de la funció s() en la fórmula del model. Cal destacar que el nombre de graus de llibertat no ha de ser necessàriament un nombre sencer.

Exemples:

gam_1df <- gam::gam(richness ~ gam::s(ndvi, df=1), data = dry, family = poisson)
gam_2df <- gam::gam(richness ~ gam::s(ndvi, df=2), data = dry, family = poisson)
gam_4df <- gam::gam(richness ~ gam::s(ndvi, df=4), data = dry, family = poisson)
gam_8df <- gam::gam(richness ~ gam::s(ndvi, df=8), data = dry, family = poisson)
par(mfrow= c(2,2))
plot(gam_1df, main= "Splines amb 1 edf")
plot(gam_2df, main= "Splines amb 2 edf")
plot(gam_4df, main= "Splines amb 4 edf")
plot(gam_8df, main= "Splines amb 8 edf")

Amb mètode LOESS

Ajustament

Si volem passar d’una suavització amb Splines a una amb regressió local (LOESS), només hauríem de canviar la funció s() per la funció lo() dins la fórmula.

gam2 <- gam::gam(richness ~ gam::lo(ndvi), data = dry, family = poisson)
summary(gam2)
## 
## Call: gam::gam(formula = richness ~ gam::lo(ndvi), family = poisson, 
##     data = dry)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.8649 -1.8516  0.2694  1.4236  3.1329 
## 
## (Dispersion Parameter for poisson family taken to be 1)
## 
##     Null Deviance: 353.1647 on 95 degrees of freedom
## Residual Deviance: 350.4069 on 94 degrees of freedom
## AIC: 762.348 
## 
## Number of Local Scoring Iterations: 4 
## 
## Anova for Parametric Effects
##               Df  Sum Sq Mean Sq F value Pr(>F)
## gam::lo(ndvi)  1   2.768  2.7679  0.8243 0.3662
## Residuals     94 315.638  3.3579
coefficients(gam2)
##   (Intercept) gam::lo(ndvi) 
##     2.1616869     0.6209708

Veiem que les dades del summary els coeficients són molt similars als que hem obtingut amb la suavització per splines, tot i que la variabilitat \(D^2\) explicada pel model és lleugerament superior.

1 - 277.0453/353.1647
## [1] 0.2155351

Gràfics

El gràfic també sembla força similar tot i que es nota una tendència de la spline a obrir més l’interval de confiança en els valors extrems.

par(mfrow= c(1,2))
plot(gam1, se=T, main= "Suavització SPLINE", ylim= c(-0.45,0.7))
plot(gam2, se=T, main= "Suavització LOESS", ylim= c(-0.45,0.7))

Regulació de l’ajustament utilitzant span

En la suavització per regressió local, el paràmetre de suavització no és el nombre de nusos (knots), sinó el nombre de punts propers que s’utilitzen per a la regressió local que es denomina span. Per defecte, span agafa el valor de 0.5, és a dir que es pren en consideració el 50% de punts repartits a la dreta i a l’esquerra del punt focal. Tot i així, aquest valor es pot modificar, com es pot veure a sota

gam_1s <- gam::gam(richness ~ gam::lo(ndvi, span = 0.1), data = dry, family = poisson)
gam_2s <- gam::gam(richness ~ gam::lo(ndvi, span = 0.2), data = dry, family = poisson)
gam_5s <- gam::gam(richness ~ gam::lo(ndvi, span = 0.5), data = dry, family = poisson)
gam_9s <- gam::gam(richness ~ gam::lo(ndvi, span = 0.9), data = dry, family = poisson)
par(mfrow= c(2,2))
plot(gam_1s, main= "span = 0.1")
plot(gam_2s, main= "span = 0.2")
plot(gam_5s, main= "span = 0.5")
plot(gam_9s, main= "span = 0.9")

b. Ajustament amb la llibreria mgcv

Un altre mètode per ajustar models additius generalitzats està disponible amb la llibreria mgcv. La funció principal del paquet és gam() (com la del paquet gam, atenció!). Aquesta funció selecciona automàticament el grau de suavització utilitzant un algoritme específic (Penalized regression splines) que optimitza la relació entre complexitat (nombre de nusos i de graus efectius de llibertat) i el grau d’ajustament a les dades (log-likelihood).

Ajustament i summary

mgcv.gam1 <- mgcv::gam(formula = richness ~ s(ndvi), data = dry, family=poisson)
summary(mgcv.gam1)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## richness ~ s(ndvi)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.53575    0.02917   86.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##          edf Ref.df Chi.sq p-value    
## s(ndvi) 8.18   8.81   83.4  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.199   Deviance explained = 24.9%
## UBRE = 1.9527  Scale est. = 1         n = 96

Els resultats s’interpreten de la mateixa manera que els del paquet gam, tot i que les taules de sortida són bastant més clares i usables. Entre altres coses interessants, hi podem trobar el coeficient estimat de l’ordenada a l’origen, la percentual de deviança explicada (\(D^2\)) i l’\(R^2\) ajustat, tot i que, igual que en els GLM, en els GAM és més adequat utilitzar la deviana explicada que l’\(R^2\), com a mesura de la bondat de l’ajustament a les dades.

Representació gràfica

Funciona igual que el paquet gam, fent servir plot()

library(grid)

pred.mgcv <- as.data.frame(predict(mgcv.gam1, type = "response", se.fit = T)[1:2]) %>% 
  mutate(upper = fit + 1.96*se.fit,
         lower= fit - 1.96*se.fit)

# Gràfic de la dreta
 P <- dry %>% 
  ggplot(aes(x=ndvi,y=richness))+ 
  geom_point()+
  geom_line(data=pred.mgcv,
        aes(x=dry$ndvi,y=fit)) +
  geom_line(data=pred.mgcv,
        aes(x=dry$ndvi,y=upper), linetype=3) +
  geom_line(data=pred.mgcv,
        aes(x=dry$ndvi,y=lower), linetype=3) +
  theme_bw()

 vp.BottomRight <- viewport(height=unit(.85, "npc"), width=unit(0.47, "npc"), 
                           just=c("left","top"), 
                           y=0.96, x=0.53)
 
par(mar= c(5,4,1.2,18))
# Gràfic de l'esquerra
plot(mgcv.gam1)

print(P, vp=vp.BottomRight)

L’ajustament que se selecciona automàticament és similar al que obtenim amb el paquet gam seleccionant 8 graus de llibertat efectius: de fet al summary, a l’apartat Approximate significance of smooth terms, els graus de llibertat efectius (edf) són 8.18. Si volguéssim fixar els graus de llibertat, hauríem de de fer-ho especificant spline (k=..., fx=TRUE).

Avaluació del model

Atès que utilitzem funcions de suavització justament per incorporar relacions no lineals, la falta de linealitat no hauria de ser mai un problema d’aquests models. D’altra banda estem assumint una funció de distribució específica per a capturar dels errors. Ens hem de preguntar doncs si la funció és adequada. Hauríem de mirar llavors als gràfics dels residus.

Amb la funció gam.check()

La funció gam.check() dibuixa alguns gràfics de diagnòstic que podrien ser útils. Es pot elegir entre alguns tipus de residuals.

par(mfrow= c(2,2))
mgcv::gam.check(mgcv.gam1, type = "pearson", pch=19)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 3 iterations.
## Gradient range [6.350042e-07,6.350042e-07]
## (score 1.952694 & scale 1).
## Hessian positive definite, eigenvalue range [0.008173242,0.008173242].
## Model rank =  10 / 10 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##           k'  edf k-index p-value
## s(ndvi) 9.00 8.18    0.92    0.21

Aquí, tot i no ser visibles pautes clares, es pot entreveure cera irregularitat en la distribució dels residus.

Amb DHARMa

També amb GAM podem utilitzar per avaluar l’adequació de la funció de distribució d’errors seleccionada. Això ho podem fer indiferentment amb objectes creats amb gam i amb mgcv. Atenció que per poder utilitzar aquest llibreria és necessària la instal·lació també del paquet mgcViz.

library(DHARMa)
library(mgcViz)
simul <- DHARMa::simulateResiduals(mgcv.gam1, 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")

DHARMa::testUniformity(simul, plot = F)
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.18099, p-value = 0.003182
## alternative hypothesis: two-sided
DHARMa::testOutliers(simul, plot = F)
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simul
## outliers at both margin(s) = 2, observations = 96, p-value = 0.02
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01041667
## sample estimates:
## outlier frequency (expected: 0.00229166666666667 ) 
##                                         0.02083333
DHARMa::testDispersion(simul, plot = F)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 2.47, p-value < 2.2e-16
## alternative hypothesis: two.sided

Els diferents test confirmen els dubtes manifestats a dalt. El gràfic Quantil-Quantil que permet verificar si els residuals del model segueixen la distribució de probabilitat teòrica elegida, manifesta certa desviació. Es manifesten problemes també d’heteroscedasticitat (gràfics de valors predits i residuals i test de kolgomorov-Smirnov) i de valors atípics.

Les problemàtiques podrien ser degudes a problemes de sobredispersió, tal com ens podria indicar el darrer test.

c. Ajustament amb mgcv per a GAM binomial negatiu

La funció gam() de mgcv disposa de la família de distribució d’errors binomial negativa, però, atenció, no gam() de gam (disposa com a alternativa de la quasi-Poisson).

Ajustament i summary del model

mgcv.gam.bn <- mgcv::gam(formula = richness ~ s(ndvi), data = dry, family = nb)
summary(mgcv.gam.bn)
## 
## Family: Negative Binomial(5.334) 
## Link function: log 
## 
## Formula:
## richness ~ s(ndvi)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.54580    0.05275   48.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##           edf Ref.df Chi.sq p-value   
## s(ndvi) 3.814  4.756  15.51 0.00992 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.165   Deviance explained = 16.1%
## -REML = 318.07  Scale est. = 1         n = 96

Representació gràfica del model

library(grid)

pred.mgcv.bn <- as.data.frame(predict(mgcv.gam.bn, type = "response", se.fit = T)[1:2]) %>% 
  mutate(upper = fit + 1.96*se.fit,
         lower= fit - 1.96*se.fit)

# Gràfic de la dreta
 P <- dry %>% 
  ggplot(aes(x=ndvi,y=richness))+ 
  geom_point()+
  geom_line(data=pred.mgcv.bn,
        aes(x=dry$ndvi,y=fit)) +
  geom_line(data=pred.mgcv.bn,
        aes(x=dry$ndvi,y=upper), linetype=3) +
  geom_line(data=pred.mgcv.bn,
        aes(x=dry$ndvi,y=lower), linetype=3) +
  theme_bw()

 vp.BottomRight <- viewport(height=unit(.85, "npc"), width=unit(0.47, "npc"), 
                           just=c("left","top"), 
                           y=0.96, x=0.53)
 
par(mar= c(5,4,1.2,18))
# Gràfic de l'esquerra
plot(mgcv.gam.bn)

print(P, vp=vp.BottomRight)

L’aspecte que queda més evident comparant aquest gràfic amb els anteriors és que s’amplia de manera molt significativa l’interval de confiança.

Avaluació del model

La llibreria DHARMa, des de fa uns mesos, disposa de la possibilitat d’avaluar també els models GAM de família de distribució binomial negativa.

simul.bn <- DHARMa::simulateResiduals(mgcv.gam.bn, n=1000)
par(mfrow= c(2,2),
    mar= c(4,4,3,2))
DHARMa::plotQQunif(simul.bn,
     testUniformity = FALSE, 
     testOutliers = FALSE,
     testDispersion = FALSE,
     xaxs = "i", yaxs = "i")
DHARMa::plotResiduals(simul.bn,xaxs = "i", yaxs = "i")

DHARMa::testUniformity(simul.bn, plot = F)
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.14118, p-value = 0.03929
## alternative hypothesis: two-sided
DHARMa::testOutliers(simul.bn, plot = F)
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simul.bn
## outliers at both margin(s) = 0, observations = 96, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.01041667
## sample estimates:
## outlier frequency (expected: 0.00166666666666667 ) 
##                                                  0
DHARMa::testDispersion(simul.bn, plot = F)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.76316, p-value = 0.19
## alternative hypothesis: two.sided

El gràfic Quantil-Quantil manifesta, sense estar encara perfecte, indica més adherència dels valors residuals al model. Tant els gràfic (esquerra), com el tast Kolgomorov manifesten encara una certa heteroscedasticitat, tot i així cal destacar la millora del p-valor. D’altra banda ja no tenim problemàtiques de valors atípics.

L’aspecte a destacar és que, tal com veiem en el test de dispersió, ja no apareix la sobredispersió que s’havia manifestat amb el model GAM de la família de distribució d’errors Poisson.

En conclusió, sembla que en conjunt el model GAM amb distribució d’errors binomial negativa s’ajusti millor a l’estructura de les dades.

Exemple GAM: la contaminació de l’aire a Chicago

La relació entre contaminació de l’aire i salut és un tema força controvertit, amb una munió de treballs epidemiològics intentant de posar en focus en aquesta relació.

La present és una variant d’un mètode d’anàlisi en l’epidemiologia de la contaminació, que s’ha tornat fora popular avui en dia.

El conjunt de dades chicago del paquet gamair proporciona dades de sèries temporals diàries que analitzen la relació entre els factors ambientals i la mortalitat diària. Abasta des de l’1 de gener de 1987 fins al 31 de desembre de 2000.

  • variable de resposta death, i com a possibles covariables independents que l’expliquin:
  • els nivells d’ozó a parts per bilió, variable o3median;
  • els nivells de diòxid de sofre, variable so2median;
  • la temperatura mitjana diària en graus Fahrenheit, variable tmpd;
  • els nivells de partícules per metre cúbic d’aire, variable pm10median;
  • i la covariable time que mesura el temps en què es va dur a terme l’estudi, ja que es creu que el nombre de morts també té a veure amb el moment estudiat de l’any, independentment de les variables ambientals.

Es pretén trobar un model GAM que relacioni la pol·lució de l’aire amb la salut. Per això es va dur a terme un estudi per Peng i Welty (2004) on es va considerar com a variable dependent el nombre diari de morts a la ciutat de Chicago durant el període que va durar l’estudi.

library(gamair)

data(chicago)
str(chicago)
## 'data.frame':    5114 obs. of  7 variables:
##  $ death     : int  130 150 101 135 126 130 129 109 125 153 ...
##  $ pm10median: num  -7.434 NA -0.827 5.566 NA ...
##  $ pm25median: num  NA NA NA NA NA NA NA NA NA NA ...
##  $ o3median  : num  -19.6 -19 -20.2 -19.7 -19.2 ...
##  $ so2median : num  1.928 -0.986 -1.891 6.139 2.278 ...
##  $ time      : num  -2556 -2556 -2554 -2554 -2552 ...
##  $ tmpd      : num  31.5 33 33 29 32 40 34.5 29 26.5 32.5 ...

Per les característiques de la variable de resposta, sembla raonable considerar un model de regressió Poisson però, per ser més general, considerarem un model GAM de Poisson de manera que la mitjana diària de morts \(µ_i=E[death_i]\) la modelitzem mitjançant el model GAM semiparamètric de la forma \[ \small\log \mu_i= \beta_1 \texttt{o3median}_i+\beta_2\texttt{so2median}_i+\beta_3\texttt{tempd}_i+\beta_4\texttt{pm10median}_i+f(\texttt{time}_i) \] on \(\small\mu_i\) segueix una distribució Poisson i \(\small f\) és una funció de suavització.

Per efectuar l’ajust executarem la funció mgcv::gam() on hem utilitzat com a base per a la suavització una regressió spline cúbica (expressada amb cr a l’argument bs) perquè és el més senzill d’executar en ser moltes dades (\(n=5114\)). El valor de k suggerit seria \(k = 20· 5114^{2/9} = 133.4\); agafarem 100 que a més és el màxim valor utilitzat, de fet, per la funció.

contamin <- mgcv::gam(death ~ s(time,bs="cr",k=100)+o3median+so2median+tmpd+pm10median,
           data=chicago,family=poisson) 

contamin
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## death ~ s(time, bs = "cr", k = 100) + o3median + so2median + 
##     tmpd + pm10median
## 
## Estimated degrees of freedom:
## 89.8  total = 94.78 
## 
## UBRE score: 0.3073298
par(mfrow= c(1,2))
mgcv::qq.gam(contamin, cex=0.8, col= "blue", pch= 19)
plot(contamin$fitted.values, contamin$y, col= "blue",
     xlab = "Fitted values", ylab = "Resonse")

simul.contamin <- DHARMa::simulateResiduals(contamin, n=1000)
DHARMa::testUniformity(simul.contamin, plot = F)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.020865, p-value = 0.02955
## alternative hypothesis: two-sided
DHARMa::testOutliers(simul.contamin, plot = F)
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simul.contamin
## outliers at both margin(s) = 31, observations = 4841, p-value =
## 3.714e-08
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.004354985 0.009077236
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.006403636
DHARMa::testDispersion(simul.contamin, plot = F)
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.3636, p-value < 2.2e-16
## alternative hypothesis: two.sided

Com acabem de veure, aquest model presenta diferents problemes: una certa heteroscedasticitat i uns outliers molt contundents.

Reperesentant la spline estimada, amb i sense valors residuals, podem evidenciar els outliers que en aquest cas són prou importants.

par(mfrow= c(2,1))
plot(contamin)
plot(contamin, residuals = T, cex=4)

Si representem deaths en ordre d’ID, podem veure que els 4 valors més alts, a la sèrie, s’hi esdevenen durant 4 dies seguits.

library(patchwork)

ch1 <- as.data.frame(chicago) |>
  ggplot() +
  geom_point(aes(x= seq(1,nrow(chicago),by=1), y= death),
             size= 0.5) +
  ggrepel::geom_text_repel(aes(x= seq(1,nrow(chicago),by=1),
                           y= death,
                           label=seq(1,nrow(chicago),by=1)),
                           size=2, max.overlaps = 20) +
  labs(x="Index",
       y="death") +
  theme_bw()

ch2 <- as.data.frame(chicago) |>
  ggplot() +
  geom_point(aes(x= seq(1,nrow(chicago),by=1), y= o3median),
             size= 0.5) +
  ggrepel::geom_text_repel(aes(x= seq(1,nrow(chicago),by=1),
                           y= o3median,
                           label=seq(1,nrow(chicago),by=1)),
                           size=2) +
  labs(x="Index",
       y="Nivell d'ozó") +
  theme_bw()

ch1 / ch2

Si en creem una sèrie temporal i prenem en consideració les variables corresponents a les morts, els nivells de PM10, els nivells d’ozó i la temperatura, veiem que el gran pic de morts coincideix amb el pic més important del nivells d’ozó i amb temperatures altes. N.Wood hi posa molta èmfasi, però a la sèrie hi ha més pics importants d’ozó (encara que no tant), però ni de bon tros coincideixen amb pics similars de mortalitat.

library(forecast)

ts.chicago <- ts(chicago, start = 1)
a <- ts.chicago[,1] %>% 
  autoplot(colour= "blue") +
  labs(title = "Nombre de morts diaries")+
  theme_bw()
b <- ts.chicago[,2] %>% 
  autoplot(colour= "green") +
  labs(title = "Nivell de PM10")+
  theme_bw()
c <- ts.chicago[,4] %>% 
  autoplot(colour= "red") +
  labs(title = "Nivell d'Ozó")+
  theme_bw()
d <- ts.chicago[,7] %>% 
  autoplot(colour= "violet") +
  labs(title = "Temperatura")+
  theme_bw()

a/c/d

D’altra banda, és veritat que els diferents tests de Philips-Ouliaris apunten a l’existència de correlació entre la variable de resposta i les covariables. Aquí facilitem el resultat del test entre la variable death i o3median.

tseries::po.test(ts.chicago[,c(1,2)])
## 
##  Phillips-Ouliaris Cointegration Test
## 
## data:  ts.chicago[, c(1, 2)]
## Phillips-Ouliaris demeaned = -11816, Truncation lag parameter = 48,
## p-value = 0.01

García Pérez i Wood afirmen (probablement ho han estudiat més a fons…) que les dades anòmales suggereixen que les altes temperatures y uns grans nivells d’ozó poden ser la causa de mortalitat (és a dir, ser bones covariables predictores) però uns dies més tard, ja que els es donen 3 o 4 dies abans de aquests màxims en les morts.

Per aquesta raó, proposen avançar 3 dies els valors de les covariables o3median, so2median, tmpd, pm10median, sumant els seus valors, mantenint els de time i els de la variable dependent death, encara que, en perdre per aquest avançament 3 dades, cal ometre els primers tres d’aquestes dues. Així, considerarem les noves variables:

death <- chicago$death[4:5114]
time <- chicago$time[4:5114]

i, per avançar sumant les dades de les altres quatre utilitzarem la funció lag.suma definida com:

lag.suma<-function(x)
{ n <- length(x)
b<- rep(0,n-3)
for(i in 0:(3-0)) b <- b+x[(i+1):(n-3+i)] 
b
}

tot definint ara les noves covariables lagged, amb retard

lag.o3<-lag.suma(chicago$o3median)
lag.tmp<-lag.suma(chicago$tmpd)
lag.pm10<-lag.suma(chicago$pm10median)
lag.so2<-lag.suma(chicago$so2median)

En ajustar el nou model GAM chic1, utilitzem funcions de suavització de totes les covariables, ajuntant les dues que, segons l’anàlisi de les dades, semblen tenir un efecte combinat encara que hem omès en aquest cas la base cr perquè la spline cúbica només admet una covariable.

chic1 <- mgcv::gam(death ~ s(time,bs="cr",k=100)+
                     s(lag.o3,lag.tmp,k=100)+ # efecte combinat ozó-temperatura
             s(lag.pm10,bs="cr",k=100)+s(lag.so2,bs="cr",k=100),
           family=poisson) 
par(mfrow=c(2,2))
par(mfrow= c(1,2))
mgcv::qq.gam(chic1, cex=0.8, col= "blue", pch= 19)
plot(chic1$fitted.values, chic1$y, col= "blue",
     xlab = "Fitted values", ylab = "Response")

Observem un clar millorament en la dispersió dels valors residuals.

Del compendi de paràmetres que ens facilita el summary() cal destacar que l’efecte combinat d’ozó i temperatura és destacable i que el model explica gairebé un 45% de la variabilitat, que ja podem començar a considerar de manera positiva.

summary(chic1)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## death ~ s(time, bs = "cr", k = 100) + s(lag.o3, lag.tmp, k = 100) + 
##     s(lag.pm10, bs = "cr", k = 100) + s(lag.so2, bs = "cr", k = 100)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 4.742249   0.001422    3335   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df Chi.sq p-value    
## s(time)           92.092 97.614 954.50  <2e-16 ***
## s(lag.o3,lag.tmp) 81.500 93.889 697.83  <2e-16 ***
## s(lag.pm10)        6.194  7.695  23.52  0.0024 ** 
## s(lag.so2)        16.992 20.573  24.38  0.2610    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.433   Deviance explained = 44.8%
## UBRE = 0.18691  Scale est. = 1         n = 4325

D’altra banda, S.N. Wood suggereix fer servir el logaritmes naturals per a \(PM_{10}\) i \(SO_2\).

chic2 <- mgcv::gam(death ~ s(time,bs="cr",k=100)+
                     s(lag.o3,lag.tmp,k=100)+ # efecte combinat ozó-temperatura
             s(log(lag.pm10),bs="cr",k=100)+s(log(lag.so2),bs="cr",k=100),
           family=poisson) 
summary(chic2)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## death ~ s(time, bs = "cr", k = 100) + s(lag.o3, lag.tmp, k = 100) + 
##     s(log(lag.pm10), bs = "cr", k = 100) + s(log(lag.so2), bs = "cr", 
##     k = 100)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 4.749723   0.003319    1431   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df  Chi.sq p-value    
## s(time)           46.794 56.110 137.194  <2e-16 ***
## s(lag.o3,lag.tmp) 62.949 79.577 695.220  <2e-16 ***
## s(log(lag.pm10))   2.042  2.587   4.262  0.2102    
## s(log(lag.so2))    2.499  3.166   8.744  0.0422 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.599   Deviance explained = 60.4%
## UBRE = 0.28297  Scale est. = 1         n = 788

Provada la solució, veiem que només l’índex de \(SO_2\) segueix presentant p valors significatius.

Amb la qual cosa apartarem \(PM10 del model\), deixant només les variables lagged la combinació de o3 i tmp, a més de so2.

chic3 <- mgcv::gam(death ~ s(time,bs="cr",k=100)+
                     s(lag.o3,lag.tmp,k=100)+ # efecte combinat ozó-temperatura
             s(log(lag.so2),bs="cr",k=100),
           family=poisson) 
summary(chic3)
## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## death ~ s(time, bs = "cr", k = 100) + s(lag.o3, lag.tmp, k = 100) + 
##     s(log(lag.so2), bs = "cr", k = 100)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 4.771581   0.002376    2008   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df Chi.sq p-value    
## s(time)           84.471 92.944 429.76 < 2e-16 ***
## s(lag.o3,lag.tmp) 68.590 84.872 485.40 < 2e-16 ***
## s(log(lag.so2))    1.001  1.002   8.96 0.00277 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.527   Deviance explained = 54.3%
## UBRE = 0.25557  Scale est. = 1         n = 1505

On veiem que la combinació ozó-temperatura el logaritme del índex d’òxid de sofre (ambdós retardats) tenen un efecte clarament significatiu sobre la mortalitat. A més a més el model explica més d’un 54% de la variabilitat.

En un primer anàlisi veiem que el model sembla presentar, pel que fa a homoscedasticitat i adherència al model dels residuals, millor característiques.

par(mfrow= c(2,2))
mgcv::gam.check(chic3)

## 
## Method: UBRE   Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-7.878913e-07,6.023003e-08]
## (score 0.2555663 & scale 1).
## Hessian positive definite, eigenvalue range [7.868398e-07,0.01616814].
## Model rank =  298 / 298 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##                     k'  edf k-index p-value
## s(time)           99.0 84.5    1.06    0.97
## s(lag.o3,lag.tmp) 99.0 68.6    1.06    1.00
## s(log(lag.so2))   99.0  1.0    1.02    0.65
simul.chic3 <- DHARMa::simulateResiduals(chic3, n=1000)
plot(simul.chic3)

Com acabem de veure, aquest model ja no presenta heteroscedasticitat i es veu adherència dels valors residuals respecte als predits.

Resten certs problemes de outliers i de sobre dispersió.