Comment tester l'effet d'une variable de regroupement avec un modèle non linéaire?

15

J'ai une question concernant l'utilisation d'une variable de regroupement dans un modèle non linéaire. Étant donné que la fonction nls () ne permet pas les variables de facteur, j'ai du mal à déterminer si l'on peut tester l'effet d'un facteur sur l'ajustement du modèle. J'ai inclus un exemple ci-dessous où je souhaite adapter un modèle de croissance «von Bertalanffy» saisonnier à différents traitements de croissance (le plus souvent appliqué à la croissance des poissons). Je voudrais tester l'effet du lac où les poissons ont grandi ainsi que la nourriture donnée (juste un exemple artificiel). Je connais une solution de contournement à ce problème - en appliquant un test F comparant les modèles ajustés aux données regroupées et les ajustements séparés, comme indiqué par Chen et al. (1992) (ARSS - "Analyse de la somme résiduelle des carrés"). En d'autres termes, pour l'exemple ci-dessous,

entrez la description de l'image ici

J'imagine qu'il existe un moyen plus simple de le faire dans R en utilisant nlme (), mais je rencontre des problèmes. Tout d'abord, en utilisant une variable de regroupement, les degrés de liberté sont plus élevés que ceux que j'obtiens avec mon ajustement de modèles séparés. Deuxièmement, je ne peux pas imbriquer les variables de regroupement - je ne vois pas où est mon problème. Toute aide en utilisant nlme ou d'autres méthodes est grandement appréciée. Voici le code de mon exemple artificiel:

###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
)

Référence: Chen, Y., Jackson, DA et Harvey, HH, 1992. Une comparaison de von Bertalanffy et des fonctions polynomiales dans la modélisation des données de croissance des poissons. 49, 6: 1228-1235.

Marc dans la boîte
la source

Réponses:

6

X1,...,XpOuiF

Oui=F(X1,...,Xp)+ε

εN(0,σ2)FBmBL1L0

Le modèle non stratifié est clairement un sous-modèle du modèle stratifié, donc le test du rapport de vraisemblance est approprié pour voir si le modèle plus grand vaut la complexité supplémentaire - la statistique du test est

λ=2(L1-L0)

λχ2mp-p=p(m-1)pχ2

Macro
la source
Suggérez-vous d'adapter m modèles distincts, d'additionner la vraisemblance logarithmique de chaque L1 = SUM (LL_i, i de 1 à m), puis de procéder à la vraisemblance? De plus, L0 est-il un modèle avec le prédicteur catégoriel en question inclus (avec m-1 variables muettes par exemple)?
B_Miner
L0BB
Merci pour votre suggestion Macro. Cela semble aller dans le sens de ce que j'ai déjà fait - bien que vous suggériez la comparaison de la probabilité plutôt que le test F. Dans mon exemple, le test F compare également les résidus d'ajustement unique à la somme des résidus de plusieurs ajustements appliqués à chaque niveau de prédicteur catégorique. Je suppose que je me demandais si on pouvait le faire dans un modèle mixte en une seule étape plutôt que de monter plusieurs modèles. De plus, une telle stratégie permettrait-elle de tester les facteurs imbriqués?
Marc dans la case du
Je ne pense pas que vous serez en mesure de contourner le montage de plusieurs modèles afin de comparer les modèles. De plus, oui, le test du rapport de vraisemblance peut être utilisé pour tester les facteurs imbriqués.
Macro
2

J'ai trouvé qu'il est possible de coder des variables catégorielles avec nls (), simplement en multipliant les vecteurs vrai / faux dans votre équation. Exemple:

# 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
la source