superponer - Mejor ajuste para un modelo lineal

superponer graficas en r

Question

Estoy ajustando algunas líneas y siento que le estoy diciendo exactamente a R cómo encajarlas, pero siento que hay algo (algún factor o efecto) que desconozco que impide un buen ajuste.

Mi unidad experimental es "trama" como en la trama de campo, lo que siento es confuso.

Los datos se pueden encontrar: https://www.dropbox.com/s/a0tplyvs8lxu1d0/rootmeansv2.csv . con

df$plot.f<-as.factor(df$plot)
dfG<-groupedData(mass ~ year|plot.f, data=df)
dfG30<-dfG[dfG$depth == 30,]

Simplemente, tengo masa en el tiempo y la ajusto a cada unidad experimental con el modelo:

fit <- lme(mass ~ year , random = ~ 1 | plot, data = df)

y con la plot (augPred(fit)) obtengo estos ajustes para cada unidad experimental ("trama"):

¿Qué debo hacer para permitir que la pendiente varíe más entre las unidades experimentales? No me interesa esto desde una perspectiva estadística, pero desde una perspectiva predictiva, por lo que cualquier cosa en el modelo puede manipularse para que esas líneas se muevan.

Answer 1

Esta respuesta va a ser un poco enciclopédica, porque hay varios puntos importantes acerca de su pregunta. tl; dr agregar el year*plot como un efecto aleatorio es el primer paso, pero el ajuste es en realidad un poco problemático, aunque al principio no parece así: al centrar la variable del year se soluciona el problema, pero se transforma el registro los datos pueden ser una idea aún mejor.

ACTUALIZACIÓN : OP estaba haciendo un análisis en subset(df,Depth==30) . Accidentalmente hice el análisis en el conjunto completo de datos, lo que puede haber dificultado las cosas (la heterocedasticidad documentada a continuación probablemente no es tan mala para un subconjunto de los datos), pero la voy a dejar como está para fines pedagógicos. (y por pereza).

Obtenga los datos:

df <- read.csv("rootmeansv2.csv")
library(nlme)
gdf <- groupedData( mass ~ year | plot, data=df)

Agregar la interacción año por argumento al modelo como un efecto aleatorio:

fit0 <- lme(mass ~ year , random = ~ year | plot, data = gdf)

Sin embargo, si miramos los resultados vemos que esto no ayuda mucho en absoluto: el término del year (varianza entre parcelas en la pendiente) es muy pequeño.

##             Variance     StdDev       Corr  
## (Intercept) 9.522880e-12 3.085916e-06 (Intr)
## year        7.110616e-08 2.666574e-04 0.32  
## Residual    3.885373e-01 6.233276e-01       

Esto sucede a veces (un ajuste singular), pero, de lo contrario, el resumen de lme no indica de ninguna otra manera que podría haber algo mal. Sin embargo, si tratamos de obtener intervalos de confianza, vemos que puede haber problemas:

 intervals(fit0)
 ##   Error in intervals.lme(fit0) : 
 ##   cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance

La otra forma de verificar dos lme4 es ajustar el mismo modelo en lme4 .

library(lme4)
VarCorr(fit1 <- lmer(mass ~ year +(year | plot), data = gdf))

##  Groups   Name        Std.Dev.   Corr  
##  plot     (Intercept) 0.66159674       
##           year        0.00059471 -1.000
##  Residual             0.62324403       
## Warning messages:
## 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
##   Model failed to converge with max|grad| = 0.132739 (tol = 0.002, component 1)
## 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
##   Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

Aparentemente recibimos respuestas y advertencias aparentemente diferentes sobre la convergencia (esto es en la versión de desarrollo, 1.1-7, que no sufre las advertencias falsas positivas notificadas ampliamente para 1.1-6). Parece que el lme4 de lme4 es ligeramente mejor:

c(logLik(fit1)-logLik(fit0))
## 0.1775161

Uno de los posibles problemas de datos para ajustes complejos tales como los modelos mixtos es un centrado y escalado pobre de las variables de predicción. En este caso, como el year es un predictor continuo que comienza en 2008, la intersección y las pendientes están extremadamente correlacionadas (¡la intersección es el valor predicho en el año 0!). Podemos solucionar el problema en este caso centrando: también sería razonable restar el mínimo (es decir, el año de inicio es 0 en lugar de 2008)

c. <- function(x) scale(x,center=TRUE,scale=FALSE)
VarCorr(fit2 <- update(fit1,.~ c.(year) +(c.(year) | plot)))

##  Groups   Name        Std.Dev. Corr 
##  plot     (Intercept) 0.53798       
##           c.(year)    0.10515  1.000
##  Residual             0.59634       

Obtenemos respuestas más razonables (y ninguna advertencia), aunque todavía tenemos intersecciones y pendientes perfectamente correlacionadas (ahora positivamente), esto solo dice que el modelo está sobreajustando ligeramente los datos, pero no hay mucho que podamos hacer al respecto (podríamos forzar la correlación a cero ajustando el modelo ~c.(year)+(c.(year)||plot)) , pero esto tiene sus propios problemas).

Ahora pruébalo en lme :

VarCorr(fit3 <- update(fit0,
                       fixed.=~c.(year),
                       random=~c.(year)|plot,
         control=lmeControl(opt="optim")))
## plot = pdLogChol(c.(year)) 
##             Variance   StdDev    Corr  
## (Intercept) 0.28899909 0.5375864 (Intr)
## c.(year)    0.01122497 0.1059479 0.991 
## Residual    0.35559015 0.5963138       

Los resultados son similares (aunque la correlación es solo de 0.991 en lugar de 1.0): la diferencia en logaritmo-verosimilitudes es en realidad un poco mayor, pero aún pequeña. (El ajuste todavía es un tanto problemático: tuve que cambiar los optimizadores como se muestra en el argumento lmeControl para obtener aquí).

Las cosas ahora se ven mejor:

plot(augPred(fit3))

Sin embargo, la gráfica residual versus ajustada le muestra un problema por el que debe preocuparse:

plot(fit3)

La forma del embudo muestra que tienes un problema de heterocedasticidad. El gráfico de ubicación de la escala lo muestra aún más claro:

plot(fit3,sqrt(abs(resid(.)))~fitted(.),type=c("p","smooth"))

La solución más obvia es la de la transformación de registros de datos:

df$logmass <- log10(df$mass)  ## log10() for interpretability
gdfL <- groupedData( logmass ~ year | plot, data=df)
VarCorr(fitL1 <- lme(logmass ~ c.(year) , random = ~ c.(year) | plot,
             data = gdfL))
## plot = pdLogChol(c.(year)) 
##             Variance     StdDev     Corr  
## (Intercept) 0.1842905872 0.42929080 (Intr)
## c.(year)    0.0009702893 0.03114947 -0.607
## Residual    0.1052946948 0.32449144       

La variación entre años es pequeña nuevamente, pero esta vez es probablemente porque no es necesaria. No hay advertencias de lmer cuando lmer el modelo equivalente, y obtenemos resultados idénticos; el ajuste no es singular (no varianzas cero o correlaciones perfectas); y los intervals(fitL1) funcionan.

Las predicciones parecen razonables

plot(augPred(fitL1))

Lo mismo ocurre con el diagrama de diagnóstico:

plot(fitL1)

Aunque parece haber algo un poco gracioso sobre 2008 (los ejes se voltean porque lme prefiere trazar lme horizontalmente - factor(year) le dice a lme que use un diagrama de caja en lugar de un diagrama de dispersión).

plot(fitL1,factor(year)~resid(.))