El famós cas dels ous de verat: una aproximació pròpia

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

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

Referències

• 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.
• D.I, Walton, Eco-Stats: Data Analysis in Ecology: From t-tests to Multivariate Abundances, 2023.
• S.N. Wood, Generalized Additive Models: An Introduction with R, 2017.

 

El cas: Estimació de l’abundància d’una espècie de Scombridae mitjançant mostreig i recompte d’ous

Aquest és un cas molt conegut, tractat de manera molt completa per Wood i reprès i adaptat per García Perez. És interessant per la gran complexitat d’elemnts que hi entren en joc i per la possibilitat d’analitzar-la des de diferents aproximacions.

L’estimació del nombre de peixos dels bancs de pesca és de gran importància avui dia atès que alguns estan gairebé esgotats i diversos organismes fixen les quotes de captura depenent d’aquestes estimacions.

La dificultat fonamental en l’estimació de tal nombre rau en la impossibilitat de la obtenció de dades, és a dir, de peixos vius. Per aquesta raó el recompte d’ous de les espècies és el mètode habitual d’estimació del nombre d’adults (o més sovint massa d’adults) ja que la taxa de producció d’ous per quilogram d’adult és fàcil d’obtenir a partir dels peixos capturats a les xarxes d’arrossegament.

Més concretament, s’estableix una cadena d’estacions de mostreig a l’àrea de recerca i s’anota el nombre d’ous de l’espècie en estudi que apareixen en unes malles extretes del fons del mar, coneixent-ne el volum d’aigua mostrejada. S’anota així el nombre d’ous de l’espècie a l’estudi per metre cúbic d’aigua.

Les dades mack de la llibreria gamair corresponen a un estudi d’aquestes característiques sobre ous de verat (Scomber scombrus) dut a terme per Borchers et al. (1997). En aquestes dades apareixen, entre d’altres, les variables egg.count, nombre d’ous en estat I al mostreig realitzat (és a dir, a la malla extreta); b.depth, profunditat del llit marí a la qual es va fer el mostreig; c.dist, distància des del lloc de mostreig al contorn de llit marí de 200 metres, distància mesurada en graus; lon, longitud de l’estació de mostreig en graus est; lat, latitud de l’estació de mostreig a graus nord; salinity, salinitat de l’aigua a l’estació de mostreig; temp.surf, temperatura de la superfície de l’aigua a l’estació de mostreig; temp.20m, temperatura de l’aigua a 20 metres de profunditat a l’estació de mostreig.

library(gamair)

data("mack")
str(mack)

## 'data.frame':    634 obs. of  16 variables:
##  $ egg.count  : num  0 0 0 1 4 3 0 1 0 0 ...
##  $ egg.dens   : num  0 0 0 18.3 90.2 ...
##  $ b.depth    : num  4342 4334 4286 1438 166 ...
##  $ lat        : num  44.6 44.6 44.6 44 44 ...
##  $ lon        : num  -4.65 -4.48 -4.3 -2.87 -2.07 -2.13 -2.27 -2.35 -2.42 -2.48 ...
##  $ time       : num  8.23 9.68 10.9 19.33 8.78 ...
##  $ salinity   : num  35.7 35.7 35.7 35.7 NA ...
##  $ flow       : num  417 405 377 420 354 373 375 364 375 362 ...
##  $ s.depth    : num  104 98 101 98 101 100 100 102 100 100 ...
##  $ temp.surf  : num  15 15.4 15.9 16.6 16.7 16.6 17.1 17.1 17.3 16.9 ...
##  $ temp.20m   : num  15 15.4 15.9 16.6 16.7 16.6 17.1 17.1 17.3 16.9 ...
##  $ net.area   : num  0.242 0.242 0.242 0.242 0.242 0.242 0.242 0.242 0.242 0.242 ...
##  $ country    : Factor w/ 4 levels "EN","fr","IR",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ vessel     : Factor w/ 4 levels "CIRO","COSA",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ vessel.haul: num  22 23 24 93 178 179 181 182 183 184 ...
##  $ c.dist     : num  0.84 0.859 0.893 0.396 0.04 ...

El propòsit és determinar un model que expliqui prou bé a la variable dependent egg.count en funció de les altres covariables independents.

En primer lloc podem geolocalitzar les densitats, egg.dens (una altra de les variables subministrades amb les dades encara que no considerada al model a ajustar) aprofitant les variables lon i lat . Això, ens donarà una idea de la grandària de l’àrea de mostreig: important per entendre les diferències que existeixen entre variables.

Per fer això geo localitzarem els punts de mostreig en un mapa d’ESRI creant una capa de visualització de densitat dels ous als diferents punts de mostreig.

Mapa interactiu de gelocalització dels punts de mostreig amb grandària del punts proporcional a la densitat de la presència d’ous de verat al punt de mostreig

library(sf)
library(tmap)

mack_sf <- st_as_sf(mack[,c(5,4,2,1,3,16,7,10,11,12)], 
                    coords = c("lon", "lat"), crs = 4326)

tmap_mode("view")


tm_shape(mack_sf) + 
  tm_dots(size = 'egg.dens', 
          fill = "yellow2", fill_alpha = .5,
          popup.vars=c("Densitat ous en fase I (ous/m2): "= "egg.dens",
                       "Nombre d'ous en fase I: "= "egg.count",
                       "Profunditat del fons marí: "= "b.depth",
                       "Distància de la isòbata de 200m en graus: "= "c.dist",
                       "Salinitat de l'aigua: "= "salinity",
                       "Temperatura de l'aigua en superfície: "= "temp.surf",
                       "Temperatura de l'aigua a 20m de profund.: "= "temp.20m",
                       "Àrea de mostreig en m2: "= "net.area"))+
    tmap_options(basemaps = "Esri.WorldImagery")

Models GLM

Començarem l’anàlisi amb una aproximació gradual, mirant primer com funcionen els models lineals generalitzats amb aquestes dades.

Ajustament amb una GLM Poisson

Com que egg.count, el recompte d’ous, és la variable de resposta, com a primer pas sembla raonable ajustar un model GLM de Poisson.

fit.glm <- glm(formula = egg.count~ lon+lat+b.depth+c.dist+salinity+temp.surf+temp.20m,
    data = mack,
    family = poisson)
summary(fit.glm)

## 
## Call:
## glm(formula = egg.count ~ lon + lat + b.depth + c.dist + salinity + 
##     temp.surf + temp.20m, family = poisson, data = mack)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.319e+01  6.021e+00  -5.513 3.53e-08 ***
## lon         -2.919e-01  1.257e-02 -23.224  < 2e-16 ***
## lat         -8.325e-02  1.985e-02  -4.194 2.74e-05 ***
## b.depth      1.365e-04  2.011e-05   6.787 1.14e-11 ***
## c.dist      -5.233e-01  3.890e-02 -13.454  < 2e-16 ***
## salinity     1.198e+00  1.606e-01   7.460 8.65e-14 ***
## temp.surf   -9.571e-02  3.614e-02  -2.648  0.00809 ** 
## temp.20m    -2.536e-01  3.001e-02  -8.451  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 9656.4  on 329  degrees of freedom
## Residual deviance: 3737.4  on 322  degrees of freedom
##   (304 observations deleted due to missingness)
## AIC: 4679.4
## 
## Number of Fisher Scoring iterations: 6

Aparentment, l’ajustament és correcte, amb unes les variables explicatives seleccionades que tenen efecte significatiu sobre la resposta.

Avaluarem ara la distribució dels valors residuals,

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

• Es veu, als dos gràfics de l’esquerra de valors predits contra valors residuals \(\small (Pearson ~Residuals\) i \(\small \sqrt{Pearson ~Residuals})\), que no hi ha homocedasticitat, amb el típic patró d’embut: al començament hi ha uns quants punts concentrats al voltant de 0 però al final hi ha una gran variació.
• El QQ-plot dels valors residuals de dalt a la dreta indica una manca de normalitat en els valors de la dreta.
• El gràfic Residual vs Leverage de baix a la dreta, mostra que hi ha dades anòmales per al model considerat.

Possible problemàtica de sobredispersió

Cal afegir que el quocient Residual deviance / Degrees of freedom és 11.58, molt alt, indicant una problemàtica de sobredispersió.

Provant amb un model GLM de distribició binomial quasi-Poisson

fit.glm.q <- glm(egg.count~ lon+lat+b.depth+c.dist+salinity+temp.surf+temp.20m,
    data = mack, family = quasipoisson())

summary(fit.glm.q)

## 
## Call:
## glm(formula = egg.count ~ lon + lat + b.depth + c.dist + salinity + 
##     temp.surf + temp.20m, family = quasipoisson(), data = mack)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.319e+01  2.445e+01  -1.358  0.17547    
## lon         -2.919e-01  5.103e-02  -5.720 2.43e-08 ***
## lat         -8.325e-02  8.059e-02  -1.033  0.30234    
## b.depth      1.365e-04  8.164e-05   1.672  0.09553 .  
## c.dist      -5.233e-01  1.579e-01  -3.314  0.00103 ** 
## salinity     1.198e+00  6.519e-01   1.837  0.06707 .  
## temp.surf   -9.571e-02  1.467e-01  -0.652  0.51469    
## temp.20m    -2.536e-01  1.218e-01  -2.082  0.03816 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 16.48351)
## 
##     Null deviance: 9656.4  on 329  degrees of freedom
## Residual deviance: 3737.4  on 322  degrees of freedom
##   (304 observations deleted due to missingness)
## AIC: NA
## 
## Number of Fisher Scoring iterations: 6

Aquest model no sembra presentar resultats acceptables, amb tan sols 3 variables significatives i sense resoldre les problemàtiques de sobre dispersió (Residual deviance / Degrees of freedom = 3717/321)

Provant amb un model GLM de distribició binomial negativa

fit.bn <- MASS::glm.nb(egg.count~ lon+lat+b.depth+c.dist+salinity+temp.surf+temp.20m,
    data = mack)
summary(fit.bn)

## 
## Call:
## MASS::glm.nb(formula = egg.count ~ lon + lat + b.depth + c.dist + 
##     salinity + temp.surf + temp.20m, data = mack, init.theta = 0.8161460249, 
##     link = log)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.289e+01  1.775e+01  -0.727   0.4675    
## lon         -3.178e-01  6.487e-02  -4.898 9.66e-07 ***
## lat         -7.266e-02  1.113e-01  -0.653   0.5140    
## b.depth      1.335e-04  9.277e-05   1.439   0.1503    
## c.dist      -9.002e-01  1.916e-01  -4.697 2.64e-06 ***
## salinity     6.593e-01  4.495e-01   1.467   0.1425    
## temp.surf   -2.257e-01  1.455e-01  -1.551   0.1209    
## temp.20m    -2.435e-01  1.023e-01  -2.380   0.0173 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.8161) family taken to be 1)
## 
##     Null deviance: 841.85  on 329  degrees of freedom
## Residual deviance: 352.87  on 322  degrees of freedom
##   (304 observations deleted due to missingness)
## AIC: 1934.7
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.8161 
##           Std. Err.:  0.0859 
## 
##  2 x log-likelihood:  -1916.6930

par(mfrow= c(2,2))
plot(fit.bn)

En aquest cas, sí sembla que s’han corregit les problemàtiques de sobredispersió, amb un quocient Residual deviance/Degrees of freedom= 1.084, però seguim tenint tan sols 3 variables significatives.

Models GAM

Model GAM family=poisson complet

Hem de provar doncs amb un model més general que els GLM: un model GAM. Com que les variables longitud i latitud poden anar juntes a l’ajustament, atès que són covariables de localització, el primer ajustament l’executem a baix, on hem afegit l’argument scale=-1 per indicar a l’algoritme que el paràmetre d’escala de la Poisson és desconegut, i l’argument gamma=1.4 per fixar problema del sobreajustament.

fit.gam1 <- mgcv::gam(egg.count~s(lon,lat,bs="ts",k=100)+ s(b.depth,bs="ts")+
                 s(c.dist,bs="ts")+s(salinity,bs="ts")+ s(temp.surf,bs="ts")+
                 s(temp.20m,bs="ts"),
               data=mack,
               family=poisson,
               scale=-1,
               gamma=1.4)

Interpretació dels graus de llibertat efectius en un GAM

Imprimint el model GAM que acabem d’ajustar veiem que mgcv ens dona uns valors de graus de llibertat:

fit.gam1

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## egg.count ~ s(lon, lat, bs = "ts", k = 100) + s(b.depth, bs = "ts") + 
##     s(c.dist, bs = "ts") + s(salinity, bs = "ts") + s(temp.surf, 
##     bs = "ts") + s(temp.20m, bs = "ts")
## 
## Estimated degrees of freedom:
## 62.46  1.86  0.00  1.06  0.00  6.19  total = 72.57 
## 
## GCV score: 7.374066

Aquest valor són importants per poder valorar la participació de les diferents variables al model

En un Model Additiu Generalitzat (GAM), els graus de llibertat efectius (\(edf\), és així com es diuen exactament) d’una variable determinen la flexibilitat de la funció de suavitzat (spline). No hi ha un únic valor “correcte”, sinó que els graus de llibertat se seleccionen automàticament o es limiten mitjançant un paràmetre de penalització segons la complexitat del patró real en les dades.

Els graus de llibertat efectius (\(edf\))

En teoria, un cop ajustat el model, el valor de \(edf\) final d’una variable indica la forma de la relació:

Valor de \(edf\)	Tipus de relació	Significat Visual
\(edf = 1\)	Lineal	Una línia recta (equivalent a una regressió lineal clàssica)
\(1 < edf \leq 2\)	Curvatura simple	Una corba suau amb una sola inflexió (ex. quadràtica o parabòlica)
\(edf > 2)\)	No lineal complexa	Corbes amb múltiples pics, valls o patrons ondulats
\(edf\) proper a \(k-1\)	Restringit / Saturat	El model volia més flexibilitat però es va quedar sense espai (k és molt baix).

A la llum dels graus de llibertat estimats (\(edf\)) que proporciona la sortida, el model ha pres decisions molt clares i dràstiques sobre cada variable. S’empra bs = "ts" (thin plate regression splines with shrinkage) perquè aquest mètode pot reduir l’edf a exactament 0, eliminant completament les variables que no aporten al model.

Interpretació del model fit.gam1 a la llum dels \(edf\)

1. (lon, lat, bs = “ts”, k = 100) → edf = 62.46

   • Diagnòstic: Aquesta interacció espacial bidimensional és altament complexa. Com que l’\(edf\) està en 62.46, es troba relativament lluny del límit teòric màxim (k-1 = 99). Significa que el model va tenir prou “espai” per capturar l’estructura espacial.
   • Què fer: Per seguretat, es pot executar mgcv::gam.check(). Si el valor p per a aquest terme és superior a 0,05, el k = 100 actual és perfecte. Si és menor a 0.05, és possible que necessiteu pujar k a 120 o 150. En aquest cas el valor p és 1.00. És a dir, k és perfect

2. s(b.depth, bs = “ts”) → edf = 1.86

   • Diagnòstic: La profunditat de fons (bottom depth) té una relació no lineal molt simple (probablement una lleugera corba en forma de U o paràbola). El valor per defecte de k (que acostuma a ser 10) és més que suficient.

3. s(c.dist, bs = “ts”) → edf = 0.00

   • Diagnòstic: La distància a la costa (coast distance) ha estat eliminada del model. Gràcies al terme “ts”, la penalització va reduir els graus de llibertat a zero absolut. No té cap efecte lineal ni no lineal sobre el recompte d’ous.
   • Què fer: Pots deixar-la al model (ja que no consumeix graus de llibertat) o eliminar el terme del tot per simplificar la fórmula.

4. s(salinity, bs = “ts”) → edf = 1.06

   • Diagnòstic: La salinitat té un efecte purament lineal (edf ≈ 1). No requereix una corba matemàtica per ser representada.

5. s(temp.surf, bs = “ts”) → edf = 0.00

   • Diagnòstic: Igual que la distància a la costa, la temperatura superficial no té cap impacte al model i ha estat penalitzada fins a zero.

6. s(temp.20m, bs = “ts”) → edf = 6.19

   • Diagnòstic: La temperatura a 20 metres té una relació fortament no lineal i complexa (una corba amb múltiples fluctuacions, becs o valls). El valor per defecte de k és suficient per capturar aquesta manera sense problemes.

Model GAM family=poisson amb reducció de variables a ran d’anàlisi de l’\(\small edf\)

fit.gam2 <- mgcv::gam(egg.count~s(lon,lat,bs="ts",k=100)+ 
                        s(b.depth,bs="ts")+
                        salinity+
                        s(temp.20m,bs="ts"),
               data=mack,
               family=poisson,
               scale=-1,
               gamma=1.4)

summary(fit.gam2)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## egg.count ~ s(lon, lat, bs = "ts", k = 100) + s(b.depth, bs = "ts") + 
##     salinity + s(temp.20m, bs = "ts")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  66.7402    23.4615   2.845  0.00480 **
## salinity     -1.8468     0.6628  -2.786  0.00573 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                edf Ref.df     F  p-value    
## s(lon,lat)  62.462     99 3.745  < 2e-16 ***
## s(b.depth)   1.854      9 0.412 0.065007 .  
## s(temp.20m)  6.210      9 2.467 0.000201 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.807   Deviance explained = 87.9%
## GCV = 7.3681  Scale est. = 5.023     n = 330

El summary ara ens dona dues taules: una de coeficients paramètric de la qual deduïm que existeix una relació negativa entre recompes d’ous i salinitat: el recompte d’ous disminueix a mesura creix la salinitat.

Cal destacar que tant l’ordenada a l’origen com la variable salinity són significatives.

La taula de termes suavitzats ens indica que tant el paràmetre de localització (lat,lon) com el de temperatura a 20 m de profunditat són significatius per a explicar el recompte d’ous, mentre que veiem que b.depth no arriba a tenir un valor p significatiu. A més a més, tenint en compte quela profunditat del fons (b.depth) sol estar fortament lligada a la posició geogràfica, és molt probable que el terme espacial s(lon, lat) (que té un \(edf\) enorme de ~62) ja estigui capturant de forma implícita l’estructura de la profunditat de la costa, fent que la variable individual b.depth es torni redundant.

La idea és comparar posteriorment els dos models directament mitjanant el Criteri d’Informació d’Akaike (AIC)

Model GAM family=poisson sense la variable b.depth

Ajustem per tant el model sense b.depth

fit.gam3 <- mgcv::gam(egg.count~s(lon,lat,bs="ts",k=100)+ 
                        # s(b.depth,bs="ts")+
                        salinity+
                        s(temp.20m,bs="ts"),
               data=mack,
               family=poisson,
               scale=-1,
               gamma=1.4)

summary(fit.gam3)

## 
## Family: poisson 
## Link function: log 
## 
## Formula:
## egg.count ~ s(lon, lat, bs = "ts", k = 100) + salinity + s(temp.20m, 
##     bs = "ts")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  67.6694    23.4989   2.880  0.00431 **
## salinity     -1.8732     0.6639  -2.822  0.00515 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                edf Ref.df     F  p-value    
## s(lon,lat)  63.495     99 4.535  < 2e-16 ***
## s(temp.20m)  6.572      9 2.665 0.000129 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.805   Deviance explained = 87.9%
## GCV = 7.3428  Scale est. = 5.0464    n = 330

Veiem, com era previsible que tots els termes implicats en el model fit.gam3 són significatius. A més a més, comparant l’AIC d’ambdós models,

AIC(fit.gam2, fit.gam3) %>% 
  mutate(diff= AIC-lead(AIC))

##                df      AIC      diff
## fit.gam2 73.52511 2236.466 -1.617587
## fit.gam3 73.06699 2238.083        NA

veiem que la diferència d’AIC, perden la variable profunditat és mínima (1.617). Amb la qual cosa podem condsidera perfectament raonable apartarla del model

La variable profunditat de fons es retira doncs del model final pel fet que no assoleix el nivell de significació estadística estàndard (\(\alpha = 0.05, p = 0.065\)) i la seva exclusió resulta en un model més parsimoniós (amb menor o equivalent AIC).

Avaluació del model

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

## 
## Method: GCV   Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [1.010947e-11,1.237273e-09]
## (score 7.342831 & scale 5.046431).
## Hessian positive definite, eigenvalue range [0.0348293,0.410599].
## Model rank =  110 / 110 
## 
## 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(lon,lat)  99.00 63.50    1.10    1.00
## s(temp.20m)  9.00  6.57    0.92    0.23

Tot i les millores que hem anat introduint al model, aquest segueix presentant una distribució de residuals insatisfactòria, sobretot en el gràfic de valors residuals per quantils teòrics (Q-Q plot). Si els punts del gràfic Q-Q formen una corba en forma de “S” o es disparen fortament cap amunt a l’extrem dret, hi ha un problema de sobredispersió (la variabilitat de les dades és molt més gran del que s’esperava) o un excés de zeros (moltes estacions de mostreig sense ous).

Per tant, hem de descartar també aquest model.

Origen de la problemàtica

L’ús de family = poisson amb una variable resposta de tipus recompte d’ous (egg.count) en ecologia marina és la causa directa de la sobredispersió extrema. Els ous de peixos (com el sorell o el verat del dataset mack) es distribueixen en agregacions massives a l’oceà (cardums de fresa), la qual cosa genera una quantitat enorme de zeros (llances sense ous) combinada amb llances amb desenes de milers d’ous. La distribució de Poisson assumeix que la variància és igual a la mitjana (Var(Y) = μ), cosa que mai ocorre en aquestes dades de fresa, on la variància és exponencialment més gran.

• Família Tweedie (family = tw()): És actualment l’estàndard d’or en ecologia pesquera**. Maneja de forma excel·lent una combinació de molts zeros continus i valors extremadament alts.
• Binomial Negativa (family = nb()): Ideal si les teves dades són estrictament nombres enters (0, 1, 2…) i hi ha molta variància. És molt superior a la distribució poisson.
• Famílies per a excés de zeros (family = ziap): Si el teu gràfic Q-Q falla principalment perquè tens una enorme quantitat de zeros reals (zones on mai hi ha fresa).

Solució alternativa la família de distribució Tweedie

En la funció gam(), per a aquesta situació que es plantejar és un canvi de família de distribució d’errors. Les alternatives serien (disposades en ordre decreixent d’eficàcia): - Família Tweedie (family = tw()): És actualment l’estàndard d’or en ecologia pesquera. Gestiona de forma excel·lent una combinació de molts zeros continus i valors extremadament alts. - Binomial Negativa (family = nb()): Ideal si les teves dades són estrictament nombres enters (0, 1, 2…) i hi ha molta variància. És molt superior a la distribució poisson. - Famílies per a excés de zeros (family = ziap)**: Si el teu gràfic Q-Q falla principalment perquè tens una enorme quantitat de zeros reals (zones on mai hi ha fresa).

En ecologia pesquera, les dades de Captura per Unitat d’Esforç (CPUE) o biomassa combinen dues naturaleses:

• Discreta: La probabilitat de llançar la xarxa i que surti completament buida (abundants zeros estructurats).
• Contínua: Si hi ha captura, la biomassa obtinguda és una variable contínua positiva (quilos o tones de peix) altament variable.

La família Tweedie (específicament quan el paràmetre de potència \(1 < p < 2\) es defineix formalment com un procés compost de Poisson-Gamma:

• Suposa que els esdeveniments d’agregació de peixos (cardúmens) entren a la xarxa seguint una distribució de Poisson (procés discret).
• Suposa que el pes o la biomassa de cadascun d’aquests cardúmens segueix una distribució Gamma (procés continu).

Replantejament del codi

• Estimació del paràmetre de potència (p): En deixar family = tw(), mgcv estimarà automàticament per màxima versemblança un paràmetre de potència entre 1 i 2. Això modelarà les teves dades com un procés on la variància escala amb la mitjana de forma flexible.
• Ús de method = “REML”: El codi original s’utilitza scale = -1 per forçar l’estimació del paràmetre d’escala. En passar a Tweedie, és altament recomanable treure scale = -1 i afegir method = “REML” (Màxima versemblança Residual Penalitzada). REML és molt més robust per estimar simultàniament els suavitzats de l’espai (lon, lat), la temperatura i el paràmetre de potència de Tweedie, evitant el sobreajustament.
• També prescindirem del paràmetre gamma.

Retorn a les variables inicials ajustant un model GAM amb family=tw()

library(mgcv)
fit.gam.tw2 <- mgcv::gam(egg.count ~s(lon,lat,bs="ts",k=100)+
                           s(b.depth,bs="ts")+
                           s(c.dist,bs="ts")+
                           s(salinity,bs="ts")+
                           s(temp.surf,bs="ts")+
                           s(temp.20m,bs="ts"),
                        data = mack,
                        family = tw(),  # <--- Clau per absorbir la sobredisp.
                        method = "REML" # <--- Molt recomenable per a Tweedie
                        )

fit.gam.tw2

## 
## Family: Tweedie(p=1.352) 
## Link function: log 
## 
## Formula:
## egg.count ~ s(lon, lat, bs = "ts", k = 100) + s(b.depth, bs = "ts") + 
##     s(c.dist, bs = "ts") + s(salinity, bs = "ts") + s(temp.surf, 
##     bs = "ts") + s(temp.20m, bs = "ts")
## 
## Estimated degrees of freedom:
## 32.2253  2.7001  0.9388  0.0009  0.0006  4.0685  total = 40.93 
## 
## REML score: 941.3792

Observem que els termes s(salinity) i s(temp.surf) amb \(edf\) nul, són redundants per al model (\(edf\approx0\)), mentre el terme s(c.dist)¹ no és nul, però té un valor d’edf proper a 1, que ens fa dubtar si convertir-lo en terme lineal.

Per acabar de decidir representarem, el gràfic d’efectes parcials per veure com queda representat el terme c.dist

par(mfrow= c(2,2))
plot(fit.gam.tw2, select = 1)
plot(fit.gam.tw2, select = 2)
plot(fit.gam.tw2, select = 3)
plot(fit.gam.tw2, select = 6)

El terme c.dist queda representat per una línia perfecta de tendència negativa. Això significa que:

La doble penalització ja va fer la selecció: En activar select = TRUE o fer servir splines d’encongiment bs= "ts", com en el nostre cas, mgcv descompon el spline en una part lineal (l’espai nul) i una part no lineal. L’algorisme ha penalitzat la part no lineal fins a reduir-la a zero i, a més, ha encongit lleument el pendent lineal (per això l’edf va caure de 1.0 a 0.9). Si es canviés manualment a un terme lineal clàssic (+x), s’eliminaria aquesta petita penalització sobre el pendent, forçant un edf exacte de 1.0 i alterant el criteri d’optimització del model.

Cal destacar a més que l’enfocament modern de modelatge amb splines penalitzats defensa que és millor deixar que l’estimador REML o GCV decideixi la complexitat de manera contínua. Passar la variable a lineal de forma manual és fer “selecció de models a ull”, la qual cosa invalida en part els p-valors generals del model a causa del biaix de selecció.

Per tant, podem concloure que variable c.dist s’ha modelat inicialment mitjançant un spline suavitzat penalitzat (bs="ts"). El procediment de selecció per encongiment (shrinkage) de mgcv ha reduït els seus graus de llibertat efectius a \(edf = 0.94\)**, demostrant un efecte purament lineal. Aquest efecte va resultar ser altament significatiu \((p < 0.001)\), mostrant una relació constant i negativa respecte a la variable de resposta.

Per tant, no modificarem la fórmula de l’ajustament, passant c.dist com a terme lineal del model.

Model definitiu family=tw()

El model definitiu serà per tant el reduït

fit.gam.tw3 <- mgcv::gam(egg.count ~s(lon,lat,bs="ts",k=100)+
                           s(b.depth,bs="ts")+
                           s(c.dist,bs="ts")+
                           s(temp.20m,bs="ts", k=20), # k: Mejora intr. ex-post
                        data = mack,
                        family = tw(),  # <--- Clau per absorbir la sobredisp.
                        method = "REML" # <--- Molt recomenable per a Tweedie
                        )

summary(fit.gam.tw3)

## 
## Family: Tweedie(p=1.333) 
## Link function: log 
## 
## Formula:
## egg.count ~ s(lon, lat, bs = "ts", k = 100) + s(b.depth, bs = "ts") + 
##     s(c.dist, bs = "ts") + s(temp.20m, bs = "ts", k = 20)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.05517    0.06845   15.41   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                 edf Ref.df     F  p-value    
## s(lon,lat)  51.0391     99 2.563  < 2e-16 ***
## s(b.depth)   3.1808      9 2.355 4.63e-06 ***
## s(c.dist)    0.9941      9 1.138 1.74e-06 ***
## s(temp.20m)  6.2944     19 1.438 1.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.707   Deviance explained = 80.7%
## -REML = 1584.9  Scale est. = 3.1819    n = 634

Amb tots els temes amb molt alta significació i un nivell de deviància explicada del 80.7%.

Avaluació del model definitiu

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

## 
## Method: REML   Optimizer: outer newton
## full convergence after 7 iterations.
## Gradient range [-8.259072e-06,5.662976e-07]
## (score 1584.878 & scale 3.181924).
## Hessian positive definite, eigenvalue range [0.440546,443.7984].
## Model rank =  137 / 137 
## 
## 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(lon,lat)  99.000 51.039    0.85   0.035 *  
## s(b.depth)   9.000  3.181    0.93   0.850    
## s(c.dist)    9.000  0.994    0.89   0.365    
## s(temp.20m) 19.000  6.294    0.79  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Veiem com millora la dispersió: el gràfic és molt més acceptable. Amb un nivell de dispersió més baixa i homoscedasticitat acceptable.

El biaix a la dreta de l’histograma, s’explica per vàries raons:

El biaix a la dreta i els pics de deshoveig

La cua llarga cap a la dreta representa els residus dels llances on vas trobar “superagregacions” o pegats massius d’ous. Si el model va predir una densitat alta (per exemple, 500 ous) basats en la temperatura i geografia, però en el llanci real vas capturar 3000 ous, el residu serà un nombre positiu molt gran. Com que aquests esdeveniments extrems són pocs però intensos, dibuixen aquesta cua allargada cap a la dreta.

Per què el Tweedie tolera això a diferència d’una Gaussiana?

Els residus estandarditzats (deviància o Pearson) d’un model Tweedie no han de ser perfectament simètrics ni normals a l’histograma si les dades originals tenen un component discret (els zeros) i un component continu (els recomptes alts). El que és veritablement crucial és que al QQ-plot, els punts segueixin la línia diagonal general amb desviacions lleus als extrems.

A destacar de la taula de gam.check

El \(edf\) s’a plantat a 6.29 (molt lluny del límit de 19), però el test dóna tres asteriscs amb un k-index baix (0.79). L’explicació real: Quan l’\(edf\) no augmenta però el p-valor del test de bases cau a zero, el problema ja no és la manca de flexibilitat \((k)\). El que està passant és que hi ha autocorrelació residual o un problema de dependència espacial/temporal a la temperatura. En haver-hi llançaments de pesca molt propers entre si que tenen temperatures idèntiques i recomptes d’ous similars (o molts zeros junts), el test de residus independents de gam.check() es trenca i dóna un fals positiu (asteriscs).

Conclusió

Per tant podem concloure que el model definitiu de l’anàlisi de la posta d’ous de verat a les costes atèntiques d’Europa és

mgcv::gam(
  egg.count ~
    s(lon,lat,bs="ts",k=100)+
    s(b.depth,bs="ts")+
    s(c.dist,bs="ts")+
    s(temp.20m,bs="ts"),
  data = mack, family = tw(), method = "REML"
  )

Si haguéssim de filar extremadament fi, per netejar els asteriscs causats per l’autocorrelació, la solució teòrica seria migrar a un model GAMM (Models Additius Generalitzats Mixts) usant altres llibreries (veure Models mixts (LMM, GLMM, GAMM) en ecologia per afegir una estructura de correlació espacial.

Tot i això, és veritat que per a la majoria de publicacions i informes científics d’avaluació d’estocs, l’estructura actual amb mgcv::gam i method="REML és perfectament acceptable i robusta.

---

1. Detall ecològic important: c.dist és el el valor de la distància des de la ubicació de la mostra fins a la isòbata (és a dir, la línea de nivell de profunditat) de 200 m mesurada en graus com si els graus de latitud fossin iguals als graus de longitud, cosa que en realitat no és així. Context ecològic: Els verats utilitzen molt la plataforma continental europea per a la desova. La línia de contorn de 200 metres és un límit ecològic crític on la plataforma poc profunda s’endinsa a l’oceà profund. Mesurar la distància fins a aquest límit ajuda a modelar on és més probable que els verats ponguin ous.↩︎