GAMM con datos espaciales fragmentados - respuesta de presencia/ausencia

Question

Bastante atascado así que espero que alguien pueda ayudarme en la dirección correcta.

Algunos antecedentes de los datos : Durante los últimos 15 años se han recogido datos sobre una especie específica durante el período comprendido entre el inicio del verano y el final del otoño en tres zonas geográficas diferentes. Esto ha acumulado un conjunto de datos bastante amplio. Un resumen del número de individuos investigados por semana en dos de las zonas:

table(data$year, data$week, data$area)

Área 1:

         26   27   28   29   30   31   32   33   34   35   36   37   38   39   40   41   42   43   44
  2002    0    0    0    0    0    0    0    0  931 1147 1142 1275 1238 1164 1295 1388 1270  707  138
  2003    0    0    0    0    0    0  595  378  499  782 1283 1319  977  820  874 1167  793  762  271
  2004    0    0    0    0    0    0    0    0    0    0  597  646  660  609  552  113   62    0    0
  2007    0    0   58   57   95  261  170  202  249  266  155  252  173  178   94   76   47   92   55
  2008    0    0    0    0    0    0  145  178  169  200  209  171  154  104  144   35   45    0    0
  2009    0    0    0    0    0    0   43   53   44   69   34   56   44   30   62   59    0    0    0
  2010    0    0   72   55  171  136  132  122  103   78  124   93   23   86    0    0    0    0    0
  2012   67    0   78   41   93   87   51  216  278  328  177  235  274  300  306  282  318  125   99
  2013    0   27   89  109   73   76  160  184  188  208  243  188  204  176  170  208  199  156  164
  2014    0    0    0    0    0    0    0    0   52   78  108  109  132  142  125  110  129   84   65
  2015    0    0    0    0   59   78  125  130  134  144  130  127  128  192   44   89   53   33   41

Área 2:

      26   27   28   29   30   31   32   33   34   35   36   37   38   39   40   41   42   43   44
  2002    0    0    0    0    0  216  417  480  549  668  758  811  625  827  603  104    0    0    0
  2003   88  101  681  320  362  110  512   80  504    0  331   54  289  237    0    0    0    0    0
  2004    0    0    0    0    0    0    0    0    0    0  320  352  435  420  392    0    0    0    0
  2007    0  137  167  189  277  227  222  241  262  264  157  205   45   81   44   75   70    0    0
  2008    0    0   45   50   81   97   88   82  113   75   79    0    0    0    0    0    0    0    0
  2009    0    0  131  144  221  168  189  190  184  222  153  204  150   56    0    0    0    0    0
  2010    0    0   71  139  123  197  155  149  133  129  158  190   49   61   53    0    0    0    0
  2012    0    0    0    0    0   59  109  202  140  122  157  107  197  192   45    0    0    0    0
  2013    0    0    0    0   82   27  151   91   77   85  122   99  140  108    0    0    0    0    0
  2014    0    0    0  104   90   95   91  113   53  155  154  125  176   68   65   60   56   56   82
  2015    0    0    0    0  171  187  173   57   74   77  170  174  173  173   67   79  200    0    0

Como puede ver, mis datos tienen bastantes valores perdidos tanto en la escala temporal como en la espacial. Y también se ha investigado un número variable de individuos en diferentes años y en diferentes áreas (debido a la variación anual de personal, financiación, clima, etc.). En la tabla he eliminado algunos años en los que sólo se realizaron muestreos en una o dos de las zonas.

Uno de los parámetros investigados para cada (y cada) individuo es una variable de presencia/ausencia ( parásito ). Espero poder montar un GAM con semana como término de suavización y zona como factor con el binario parásito como respuesta. Es decir, quiero ver si el momento en el que hay un pico de parásitos y si esto varía entre zonas. (La experiencia nos dice que sí). Una versión simple de mi modelo (sin todas las variables explicativas) es:

model<-gamm4(parasite~s(week, by=area)+area+s(size), family=binomial, data=data)

(He incluido el parámetro tamaño sólo para destacar que también tengo variables explicativas que se relacionan con el individuo)

Así que mi pregunta principal

¿Cómo puedo incluir el año en mi modelo? En realidad, no me interesa la variación anual, sino principalmente la tendencia general de la respuesta en una zona y la diferencia entre zonas. Por eso he "agrupado" los años. Pero sé que debería incluir el año como algún tipo de variable aleatoria (factor), pero no estoy seguro de cómo..

También estoy preparado para escuchar que mis datos están demasiado fragmentados en la escala temporal y/o espacial para utilizar un GAM, pero también estoy muy agradecido por cualquier aportación sobre el enfoque estadístico correcto para mi pregunta.

Answer 1

Podrías añadir s(Year) + s(Year, by = area, m = 1) que generará una tendencia común ( s(Year) ) y area -suaves desviaciones de esta tendencia común. El m = 1 en el area -El término específico utiliza una penalización sobre la primera derivada del spline y, por lo tanto, penaliza las desviaciones de una línea plana y horizontal, que cuando se acopla con los modelos de tendencia común cómo el Year efecto en cada área se aparta del global/encuesta Year efecto.

Si quieres tratar Year como suave pero también verdaderamente aleatorio, podrías hacer uno de

1. s(Year) + s(Year, area, bs = "fs", m = 1) o
2. s(Year, area, bs = "fs")

La primera forma es sólo el "verdadero efecto aleatorio" del by versión que mostré anteriormente (efecto común/población más area -salidas específicas). La segunda forma ajusta un intercepto aleatorio y suave para cada area directamente.

Si incluyes Year como efecto aleatorio estándar (intercepción aleatoria para Year ), esto sería muy parecido a ajustar una spline en Year dado que de todos modos estás representando los splines como efectos aleatorios en este modelo.

(Supongo que quieres gamm4() porque tienes algunos efectos aleatorios complejos que necesitan lme4? Si no es así, mira los splines de efectos aleatorios en gam() a través de la bs = "re" base y con method = "REML" .)