Dos métodos para agregar efectos aleatorios a un GAM dan resultados muy diferentes. ¿Por qué es esto y cuál debe usarse?

Question

Una sección particular de la mgcv documentación le ofrece varios métodos de incorporación de efectos aleatorios en un modelo aditivo generalizado. Dos métodos son: 1) para agregar un suave plazo en las etiquetas de clase utilizando bs="re" en gam; 2) Uso de la función gamm, que incluye instalaciones similares a lme, combinado con las funciones existentes para gam. Sin embargo, en los datos simulados, los dos dan bastante diferente modelo se ajusta. ¿Por qué es eso y que uno debe usar?

x <- rnorm(1000) 
ID <- rep(1:200,each=5)
y <- x 
for(i in 1:200) y[which(ID==i)] <- y[which(ID==i)] + rnorm(1)
y <- y + rnorm(1000)
ID <- as.factor(ID)
m1 <- gam(y ~ x + s(ID,bs="re"))
m2 <- gamm(y ~ x, random=list(ID=~1) )
mean( (fitted(m1)-fitted(m2$gam))^2 ) 

Answer 1

Sospecho que es la diferencia en términos de lo que los valores ajustados se están obteniendo. Si nos fijamos en lo que yo llamaría el ajuste del modelo, el coeficiente de las estimaciones de la varianza de los términos, los modelos son idénticos. Compare summary(m2$lme) con summary(m1) y gam.vcomp(m1).

> summary(m1)

Family: gaussian 
Link function: identity 

Formula:
y ~ x + s(ID, bs = "re")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.05234    0.07932    0.66     0.51    
x            1.01375    0.03535   28.68   <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(ID) 167.1    199 5.243  <2e-16 ***
---
Signif. codes:  0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1

R-sq.(adj) =  0.674   Deviance explained = 72.9%
GCV = 1.2133  Scale est. = 1.0082    n = 1000
> summary(m2$lme)
Linear mixed-effects model fit by maximum likelihood
 Data: strip.offset(mf) 
       AIC     BIC    logLik
  3218.329 3237.96 -1605.165

Random effects:
 Formula: ~1 | ID
        (Intercept) Residual
StdDev:    1.025306 1.003452

Fixed effects: y.0 ~ X - 1 
                 Value  Std.Error  DF   t-value p-value
X(Intercept) 0.0523358 0.07922717 799  0.660578  0.5091
Xx           1.0137531 0.03535887 799 28.670404  0.0000
 Correlation: 
   X(Int)
Xx 0.014 

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.80375873 -0.67702485  0.04245145  0.64026891  2.59257295 

Number of Observations: 1000
Number of Groups: 200 

Vemos que las estimaciones para el $\hat{\beta}_x$ son 1.01375 y $\hat{\beta}_0$ ~0,05 en ambos modelos. Tenga en cuenta también la estimación de entre sujeto/grupo de varianza (desviación estándar), $\hat{\sigma}_{\mathrm{ID}}$ es dado como 1.025 en la salida de summary(m2$lme). La misma información puede ser calculada a partir de la gam modelo de uso de la gam.vcomp(), lo que da por m1

> gam.vcomp(m1)
   s(ID) 
1.027795

que está lo suficientemente cerca de un partido para nosotros, no te preocupes por eso.

Por lo tanto el fitted métodos se deben devolver distintos valores ajustados; si generamos ajustar los valores de m2$lme, entonces se obtienen los mismos valores como la producida por fitted(m1):

> mean((fitted(m1)-fitted(m2$lme))^2) 
[1] 2.966927e-07

que es para todos los intentos y propósitos 0.

fitted.lme está documentado que la devolución de las contribuciones de la población (un promedio de más de ID) y para el tema específico de los componentes. Esto es lo fitted.gam va a hacer para m1 porque representa el efecto aleatorio como un spline "fijo" efecto. En el caso de la gamm modelo, fitted.gam está volviendo a los valores ajustados para el "fijo" efecto parte del modelo, lo que explicaría la diferencia. (Estoy escribiendo "fijo" porque con ellas los términos "fijo" y "aleatorio" efectos de obtener un poco borrosa.)

En su libro, Simon Madera menciona este tema en el ejemplo que se ajuste a la Suela de los datos utilizando gamm(). Él habla sobre el uso de resid() y fitted() sobre el $lme componente de la gamm modelo excluyendo e incluyendo el de efectos aleatorios, respectivamente. Él dice que esto es apropiado para el modelo de diagnóstico en el ejemplo específico que aquí se está utilizando.

Lo que usted necesita dependerá del contexto de su uso específico/pregunta de investigación.

Si todo lo que usted necesita son simples efectos aleatorios como este y que están familiarizados con GAMs y mgcv entonces podría ser más simple de todos los de la ronda de usar sólo el efecto aleatorio de la spline con base gam() , en lugar de tener que lidiar con la extraña salida del híbrido que es un GAMM modelo ajustado a través de la gamm(). Como he mostrado anteriormente, los dos modelos son de eficacia equivalente, y la diferencia que usted informe es sólo por si los valores ajustados para incluir o excluir el tema (o ID) efectos específicos.