R: modelo aditivo generalizado sobre datos proporcionales

Question

Introducción

Estoy analizando datos poblacionales temporales del anfípodo Orchestia gammarellus . En varios momentos de cada año, se recogieron todos los animales en un pequeño punto, y se midieron varios rasgos del ciclo vital. Para estos rasgos, quiero saber si hay estacionalidad y/o tendencias. Tengo datos desde mediados de 2014.

En este ejemplo, me centro en la proporción de hembras reproductivamente activas en cada muestra. Esto alcanza su punto máximo cada año alrededor de julio, por lo que esperaría encontrar estacionalidad.

Los datos se han recogido de forma irregular, por lo que pensé que el análisis de series temporales no era apropiado. En cambio, he leído cosas buenas sobre los modelos aditivos generalizados (GAM), y estoy intentando utilizarlos en R. Hasta ahora no he tenido suerte (véase más abajo).

Métodos

Sistema de estudio

Los animales se recogen en las marismas de pastoreo de la isla de Schiermonnikoog, una isla barrera en el norte de los Países Bajos. Se trata de la parte occidental de la marisma salada, que es pastoreada por el ganado. El ganado se alimenta principalmente de Elytrigia atherica pero evita las más duras Juncus gerardii . Esto crea parches de "alta vegetación" J. gerardii hábitat, en el que el suelo está suelto, en medio de un mar de "vegetación baja" E. atherica hábitat, donde el suelo está compactado por el pisoteo.

En comparación, la parte oriental de la marisma salada no está pastoreada y no hay mosaicos de tipos de vegetación. Aquí el suelo es mucho más suelto, y las densidades de anfípodos mucho más altas.

Animales

El anfípodo Orchestia gammarellus es una especie clave en las marismas europeas. Al excavar pequeños túneles, modifica los niveles de nutrientes del suelo, estimula su desarrollo, favorece las condiciones de oxígeno del suelo y repercute en la sucesión de la vegetación. Sin embargo, en el suelo compactado y pastoreado, este comportamiento se limita a parches de "vegetación alta".

Al igual que los isópodos, los anfípodos tienen una bolsa de cría en la que se encuentran los juveniles, y son portadores de vida.

Diseño del estudio

En una parcela de unos 50 m2 en la marisma pastoreada, los animales se recogen idealmente una vez cada cuatro (en verano) o seis (en invierno). Se coloca un cilindro (de 15 cm de diámetro) sobre un trozo de estiércol de vaca en la vegetación baja, que es un punto caliente en las manchas de vegetación baja, y se recogen todos los animales en etanol al 70%.

Todos los animales se inspeccionan con un microscopio. Se registra el sexo (o los juveniles, para los que no es posible), el número de segmentos antenales (un indicador de la edad, ya que se añade uno en cada muda) y la longitud del primer segmento corporal (pereon 1, como indicador del tamaño corporal). En el caso de las hembras, registramos el estadio de la bolsa de cría (podemos identificar cuatro niveles, de los cuales sólo el cuarto es reproductivamente activo) y, en el caso de las hembras plenamente reproductoras, el tamaño de la cría (si está presente) y el estadio de la bolsa de cría (podemos identificar tres niveles).

Modelos estadísticos

Una de las variables que quiero investigar es la proporción de hembras plenamente reproductoras en las muestras, para la que quiero saber si hay estacionalidad y tendencia. Aquí están los datos:

brood_pouch <- data.frame(
  date = c("2015-07-08","2015-08-20","2015-10-07","2015-12-10","2016-02-26",
       "2016-04-29","2016-06-26","2016-07-12","2016-08-13","2016-09-29",
       "2016-11-29","2017-03-03","2017-04-25","2017-06-17","2017-07-11",
       "2017-08-13","2017-10-03","2017-11-29","2018-03-14","2018-04-29",
       "2018-06-17","2018-12-07"),
  n = c(101, 57, 75, 95, 118, 203, 197, 134, 77, 175, 314, 236, 64, 171, 103, 
    288, 49, 61, 133, 51, 75, 154),
  proportion = c(0.900990099009901, 0.210526315789474, 0, 0, 0, 
             0.108374384236453, 0.903553299492386, 0.985074626865672, 
             0.246753246753247, 0.0171428571428571, 0, 0, 0.03125, 
             0.906432748538012, 0.300970873786408, 0.197916666666667, 
             0.0612244897959184, 0, 0, 0.176470588235294, 0.986666666666667, 
             0.00649350649350649)
)

Como puede verse en la figura siguiente, las hembras reproductivamente activas son abundantes en verano pero están ausentes en invierno.

Si fuera un modelo lineal generalizado, crearía variables separadas de year et month y crear el siguiente modelo:

glm(proportion ~ year * month, data = brood_pouch, 
                             weight = n_total, family = binomial)

donde el peso es el número de individuos que componen cada proporción.

Con un modelo aditivo generalizado, esperaría crear un modelo similar, donde en este caso habría una clara estacionalidad pero no una clara tendencia.

He estado mirando el gam y el paquete prophet pero parecen operar con puntos de datos individuales por fecha.

Answer 1

Una forma de hacerlo con un GAM es establecer un modelo que descomponga el tiempo en diferentes componentes. En este caso, se quiere un componente estacional, potencialmente una tendencia aunque los datos no parezcan tenerla, y luego un resto, que supondremos por ahora que es ruido independiente, de manera que dado el modelo tenemos extracciones aleatorias de una distribución binomial con una expectativa dada por los valores ajustados del modelo.

Esta es una forma de conseguirlo utilizando mgcv .

Primero hay que convertir el date variable en algo útil; inicialmente un Date y luego lo procesamos para obtener información sobre el día del año (la semana del año también funcionaría), el año, que también lanzamos como un factor para uno de los modelos que veremos.

brood_pouch <- transform(brood_pouch,
                         date = as.Date(as.character(date)))
brood_pouch <- transform(brood_pouch,
                         doy   = as.numeric(format(date, '%j')),
                         year  = as.numeric(format(date, '%Y')),
                         fyear = as.factor(format(date, '%Y')))

El modelo razonable más sencillo podría ser

m1 <- gam(proportion ~ fyear + s(doy), data = brood_pouch, method = 'REML',
          weights = n, family = binomial)

que tiene una suavidad estacional más un efecto medio diferente para cada año. Extendiendo este modelo para incluir suavizaciones estacionales separadas para cada año, tenemos:

m2 <- gam(proportion ~ fyear + s(doy, by = fyear, k = 5), data = brood_pouch, 
          method = 'REML', weights = n, family = binomial)

Y el modelo completo con la tendencia a largo plazo modelada suavemente y la suavidad estacional permitida para variar con la tendencia estaría dada por

m3 <- gam(proportion ~ te(year, doy, k = c(3,6)), data = brood_pouch,
          method = 'REML', weights = n, family = binomial)

De estos tres, hay una clara preferencia por m2 :

> AIC(m1, m2 ,m3)
         df      AIC
m1 11.23595 169.6490
m2 14.99039 114.4300
m3 14.43998 130.7556

pero hay algunos problemas con ese modelo que un análisis adecuado podría tratar de abordar.

También se puede modelar la tendencia como un s(Date) pero eso no descompondrá los datos en efectos estacionales más el cambio entre años, y tendrá que aumentar k para obtener una spline lo suficientemente flexible como para modelar todos los datos en una sola toma.