¿Cómo probar el efecto de una variable de agrupación con un modelo no lineal?

15

Tengo una pregunta sobre el uso de una variable de agrupación en un modelo no lineal. Dado que la función nls () no permite variables de factor, he estado luchando por descubrir si uno puede probar el efecto de un factor en el ajuste del modelo. He incluido un ejemplo a continuación donde deseo ajustar un modelo de crecimiento "estacionalizado de von Bertalanffy" a diferentes tratamientos de crecimiento (más comúnmente aplicado al crecimiento de peces). Me gustaría probar el efecto del lago donde crecieron los peces y la comida que se les dio (solo un ejemplo artificial). Estoy familiarizado con una solución alternativa a este problema: la aplicación de una prueba F que compara modelos ajustados a datos agrupados frente a ajustes separados como se describe por Chen et al. (1992) (ARSS - "Análisis de la suma residual de cuadrados"). En otras palabras, para el ejemplo a continuación,

ingrese la descripción de la imagen aquí

Me imagino que hay una manera más simple de hacer esto en R usando nlme (), pero estoy teniendo problemas. En primer lugar, al usar una variable de agrupación, los grados de libertad son más altos que los que obtengo con mi ajuste de modelos separados. En segundo lugar, no puedo anidar las variables de agrupación: no veo dónde está mi problema. Cualquier ayuda usando nlme u otros métodos es muy apreciada. A continuación hay un código para mi ejemplo artificial:

###seasonalized von Bertalanffy growth model
soVBGF <- function(S.inf, k, age, age.0, age.s, c){
    S.inf * (1-exp(-k*((age-age.0)+(c*sin(2*pi*(age-age.s))/2*pi)-(c*sin(2*pi*(age.0-age.s))/2*pi))))
}

###Make artificial data
food <- c("corn", "corn", "wheat", "wheat")
lake <- c("king", "queen", "king", "queen")

#cornking, cornqueen, wheatking, wheatqueen
S.inf <- c(140, 140, 130, 130)
k <- c(0.5, 0.6, 0.8, 0.9)
age.0 <- c(-0.1, -0.05, -0.12, -0.052)
age.s <- c(0.5, 0.5, 0.5, 0.5)
cs <- c(0.05, 0.1, 0.05, 0.1)

PARS <- data.frame(food=food, lake=lake, S.inf=S.inf, k=k, age.0=age.0, age.s=age.s, c=cs)

#make data
set.seed(3)
db <- c()
PCH <- NaN*seq(4)
COL <- NaN*seq(4)
for(i in seq(4)){
    age <- runif(min=0.2, max=5, 100)
    age <- age[order(age)]
    size <- soVBGF(PARS$S.inf[i], PARS$k[i], age, PARS$age.0[i], PARS$age.s[i], PARS$c[i]) + rnorm(length(age), sd=3)
	PCH[i] <- c(1,2)[which(levels(PARS$food) == PARS$food[i])]
	COL[i] <- c(2,3)[which(levels(PARS$lake) == PARS$lake[i])]
	db <- rbind(db, data.frame(age=age, size=size, food=PARS$food[i], lake=PARS$lake[i], pch=PCH[i], col=COL[i]))
}

#visualize data
plot(db$size ~ db$age, col=db$col, pch=db$pch)
legend("bottomright", legend=paste(PARS$food, PARS$lake), col=COL, pch=PCH)


###fit growth model
library(nlme)

starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)

#fit to pooled data ("small model")
fit0 <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  start=starting.values
)
summary(fit0)

#fit to each lake separatly ("large model")
fit.king <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  start=starting.values,
  subset=db$lake=="king"
)
summary(fit.king)

fit.queen <- nls(size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  start=starting.values,
  subset=db$lake=="queen"
)
summary(fit.queen)


#analysis of residual sum of squares (F-test)
resid.small <- resid(fit0)
resid.big <- c(resid(fit.king),resid(fit.queen))
df.small <- summary(fit0)$df
df.big <- summary(fit.king)$df+summary(fit.queen)$df

F.value <- ((sum(resid.small^2)-sum(resid.big^2))/(df.big[1]-df.small[1])) / (sum(resid.big^2)/(df.big[2]))
P.value <- pf(F.value , (df.big[1]-df.small[1]), df.big[2], lower.tail = FALSE)
F.value; P.value


###plot models
plot(db$size ~ db$age, col=db$col, pch=db$pch)
legend("bottomright", legend=paste(PARS$food, PARS$lake), col=COL, pch=PCH)
legend("topleft", legend=c("soVGBF pooled", "soVGBF king", "soVGBF queen"), col=c(1,2,3), lwd=2)

#plot "small" model (pooled data)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100))
pred <- predict(fit0, tmp)
lines(tmp$age, pred, col=1, lwd=2)

#plot "large" model (seperate fits)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100), lake="king")
pred <- predict(fit.king, tmp)
lines(tmp$age, pred, col=2, lwd=2)
tmp <- data.frame(age=seq(min(db$age), max(db$age),,100), lake="queen")
pred <- predict(fit.queen, tmp)
lines(tmp$age, pred, col=3, lwd=2)



###Can this be done in one step using a grouping variable?
#with "lake" as grouping variable
starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)
fit1 <- nlme(model = size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  fixed = S.inf + k + c + age.0 + age.s ~ 1,
  group = ~ lake,
  start=starting.values
)
summary(fit1)

#similar residuals to the seperatly fitted models
sum(resid(fit.king)^2+resid(fit.queen)^2)
sum(resid(fit1)^2)

#but different degrees of freedom? (10 vs. 21?)
summary(fit.king)$df+summary(fit.queen)$df
AIC(fit1, fit0)


###I would also like to nest my grouping factors. This doesn't work...
#with "lake" and "food" as grouping variables
starting.values <- c(S.inf=140, k=0.5, c=0.1, age.0=0, age.s=0)
fit2 <- nlme(model = size ~ soVBGF(S.inf, k, age, age.0, age.s, c), 
  data=db,
  fixed = S.inf + k + c + age.0 + age.s ~ 1,
  group = ~ lake/food,
  start=starting.values
)

Referencia: Chen, Y., Jackson, DA y Harvey, HH, 1992. Una comparación de von Bertalanffy y las funciones polinómicas en el modelado de datos de crecimiento de peces. 49, 6: 1228-1235.

Marc en la caja
fuente

Respuestas:

6

X1,...,XpagYF

Y=F(X1,...,Xpag)+ε

donde εnorte(0 0,σ2)FsimetrosiL1L0 0

El modelo no estratificado es claramente un submodelo del modelo estratificado, por lo que la prueba de razón de probabilidad es apropiada para ver si el modelo más grande vale la complejidad adicional: el estadístico de prueba es

λ=2(L1-L0 0)

Si el predictor categórico realmente no tiene efecto, λχ2metropag-pag=pag(metro-1)pagχ2

Macro
fuente
¿Está sugiriendo ajustar m modelos separados, sumar la probabilidad logarítmica de cada L1 = SUMA (LL_i, i de 1 a m) y luego proceder con la probabilidad? Además, ¿es L0 un modelo con el predictor categórico en cuestión incluido (con variables ficticias m-1, por ejemplo)?
B_Miner
Sí, estoy sugiriendo eso. L0 0sisi
Gracias por su sugerencia Macro. Esto parece estar en la dirección de lo que ya he hecho, aunque sugiere una comparación de probabilidad en lugar de la prueba F. En mi ejemplo, la prueba F también compara los residuos de ajuste único con la suma de los residuos de varios ajustes aplicados a cada nivel de predictores categóricos. Supongo que me preguntaba si se puede hacer esto dentro de un modelo mixto en un solo paso en lugar de ajustar varios modelos. Además, ¿tal estrategia permitiría la prueba de factores anidados?
Marc en la caja el
No creo que pueda moverse ajustando varios modelos para comparar modelos. Además, sí, la prueba de razón de verosimilitud se puede usar para evaluar los factores anidados.
Macro
2

Descubrí que es posible codificar variables categóricas con nls (), simplemente multiplicando los vectores verdadero / falso en su ecuación. Ejemplo:

# null model (no difference between groups; all have the same coefficients)
nls.null <- nls(formula = percent_on_cells ~ vmax*(Time/(Time+km)),
            data = mehg,
            start = list(vmax = 0.6, km = 10))

# alternative model (each group has different coefficients)
nls.alt <- nls(formula = percent_on_cells ~ 
              as.numeric(DOC==0)*(vmax1)*(Time/(Time+(km1))) 
            + as.numeric(DOC==1)*(vmax2)*(Time/(Time+(km2)))
            + as.numeric(DOC==10)*(vmax3)*(Time/(Time+(km3)))
            + as.numeric(DOC==100)*(vmax4)*(Time/(Time+(km4))),
            data = mehg, 
            start = list(vmax1=0.63, km1=3.6, 
                         vmax2=0.64, km2=3.6, 
                         vmax3=0.50, km3=3.2,
                         vmax4= 0.40, km4=9.7))
housetyrell
fuente