Selección y comparación de modelos en GAMM utilizando R (mgcv)

Question

Estoy ajustando un GAMM con estructura de correlación, utilizando una familia no gaussiana. Aquí hay un ejemplo de mi modelo global: M0 <- gamm(response ~ var1*var2 + var3 + s(var4) + s(var5) + s(var6,var7), random=list(placeID= ~1), correlation= corAR1(form= ~ year | placeID), data=data, family=quasipoisson)

Me gustaría hacer una selección de términos en el modelo global para eliminar cualquier variable no relevante, pero no estoy seguro de cómo debería hacerse. Basándome en el texto de ayuda de gam.selection - ¡y esperando que se aplique también a GAMM! - Hasta ahora he utilizado la selección escalonada hacia atrás. El texto dice "Es perfectamente posible realizar una selección hacia atrás utilizando los valores p de la forma habitual" . Así, he empezado por los suavizadores (eliminando los suavizadores de cualquier variable lineal y/o no significativa) y luego he pasado a las variables lineales. Me interesaría saber si cree que este método es útil en el GAMM o no. Por favor, diga también si tiene algo más que sugerir (por ejemplo, el GCV mencionado más abajo).

También me hubiera gustado comparar modelos que difieren en su estructura respecto a las coordenadas (suavizadas o no, interacción incluida o no). Sospecho que s(latitude,longitude) podría estar sobreajustado, y me gustaría comprobarlo de alguna manera. Dado que el modelo utiliza PQL, deduzco que el AIC no se recomienda (aunque se incluye en la parte lme).

El gam.selection El texto menciona diferentes puntuaciones: "En general, el método más coherente desde el punto de vista lógico para decidir qué términos incluir en el modelo es comparar las puntuaciones GCV/UBRE/ML de los modelos con y sin el término". ¿Este método sólo sirve para comparar modelos que difieren por una sola variable? Si es así, ¿tiene alguna sugerencia sobre cómo comparar modelos como:

• M1 <- gamm(response ~ var1*var2 + s(lat,long)+ s(var5), random=list(placeID= ~1), correlation= corAR1(form= ~ year | placeID), data=data, family=quasipoisson)
• M2 <- gamm(response ~ var1*var2 + lat + s(var5), random=list(placeID= ~1), correlation= corAR1(form= ~ year | placeID), data=data, family=quasipoisson)

Tenía la impresión de que mi modelo GAMM global M0 utiliza el GCV (por defecto), pero cuando llamo a gam.check(M0$gam) Sólo obtengo los cuatro gráficos y ninguna puntuación GCV u otra salida de texto. ¿Hay alguna manera de obtener la puntuación GCV de un objeto GAMM?

Answer 1

No estás adaptando un GAM sino un GAMM. Si mira ?gamm verá que hay dos opciones para el ajuste de un modelo AMM y sólo una cuando se ajusta un GAMM:

  method: Which of ‘"ML"’ or ‘"REML"’ to use in the Gaussian additive
          mixed model case when ‘lme’ is called directly. Ignored in
          the generalized case (or if the model has an offset), in
          which case ‘gammPQL’ is used.

Como se está equipando con un quasipoisson familia, la jerarquía de ajuste del modelo es:

gamm() ---> MASS::glmmPQL() ---> nlme::lme()

y la probabilidad que se optimiza es una cuasi-verosimilitud penalizada (PQL). Como los splines se tratan como parte de los efectos fijos y parte de los efectos aleatorios en esta formulación del GAMM, la selección de la suavidad se reduce al uso de PQL y no tengo conocimiento de los resultados que pueden aconsejar sobre lo bien que funciona PQL en ese sentido (a diferencia de la selección REML con GAMs regulares).

Si sólo necesita efectos aleatorios simples, podría ser posible trabajar en el lado del GAM:

• Sur ?magic hay un ejemplo de ajuste del GAM con datos correlacionados (pero en el caso gaussiano.
• Sur ?bam hay una opción para especificar un parámetro AR(1), pero de nuevo, en el caso gaussiano.

Tendría que utilizar bs = "re" en un s() para el efecto aleatorio.