Wie teste ich den Effekt einer Gruppierungsvariablen mit einem nichtlinearen Modell?

Ich habe eine Frage zur Verwendung einer Gruppierungsvariablen in einem nichtlinearen Modell. Da die Funktion nls () keine Faktorvariablen zulässt, habe ich mich bemüht, herauszufinden, ob man die Auswirkung eines Faktors auf die Modellanpassung testen kann. Im Folgenden habe ich ein Beispiel aufgeführt, in dem ich ein "saisonalisiertes von Bertalanffy" -Wachstumsmodell für verschiedene Wachstumsbehandlungen (die am häufigsten für das Fischwachstum angewendet werden) verwenden möchte. Ich möchte die Wirkung des Sees, in dem der Fisch gewachsen ist, sowie das angegebene Futter testen (nur ein künstliches Beispiel). Ich kenne eine Problemumgehung für dieses Problem - Anwenden eines F-Tests, der Modelle vergleicht, die an gepoolte Daten angepasst sind, im Vergleich zu separaten Anpassungen, wie von Chen et al. (1992) (ARSS - "Analysis of residual sum of squares"). Mit anderen Worten, für das Beispiel unten,

Bildbeschreibung hier eingeben

Ich stelle mir vor, dass es einen einfacheren Weg gibt, dies in R mit nlme () zu tun, aber ich habe Probleme. Erstens sind die Freiheitsgrade bei Verwendung einer Gruppierungsvariablen höher als bei der Anpassung separater Modelle. Zweitens kann ich keine Gruppierungsvariablen verschachteln - ich sehe nicht, wo mein Problem liegt. Jede Hilfe mit nlme oder anderen Methoden wird sehr geschätzt. Unten ist Code für mein künstliches Beispiel:

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

Referenz: Chen, Y., Jackson, DA und Harvey, HH, 1992. Ein Vergleich von von Bertalanffy- und Polynomfunktionen bei der Modellierung von Fischwachstumsdaten. 49, 6: 1228 & ndash; 1235.

r mixed-model nls Marc in der Kiste
quelle

Antworten:

Sie können nach den Werten des kategorialen Prädiktors schichten und die Übereinstimmungen vergleichen. Zum Beispiel : Angenommen , Sie kontinuierlich Prädiktoren haben und abhängige Variable . Ich glaube, nls () gibt die maximale Wahrscheinlichkeitsschätzung von so an, dass $X_{1}, ..., X_{p}$ $Y$ $f$

Y = f (X_{1}, . . ., X_{p}) + ε

$Y = f(X_1, ..., X_p) + \varepsilon$

Dabei werden und nichtlinear parametrisiert (siehe nls-Hilfedatei). Angenommen, Sie haben einen kategorialen Prädiktor mit Ebenen und schichten die Daten basierend auf den Werten von und passen das Modell in jede Schicht ein. Da es sich um disjunkte Teilmengen der Daten handelt, ist die Log-Wahrscheinlichkeit für das geschichtete Modell die Summe der Wahrscheinlichkeit innerhalb jeder Schicht, die mit der logLik-Funktion aus einem nls-Modell extrahiert werden kann (Sie können auch die Log-Wahrscheinlichkeit abrufen) aus dem nicht geschichteten Modell mit logLik). $\varepsilon \sim N(0,\sigma^2)$ $f$ $B$ $m$ $B$ $L_1$ $L_0$

Das nicht geschichtete Modell ist eindeutig ein Untermodell des geschichteten Modells, daher ist der Likelihood-Ratio-Test geeignet, um festzustellen, ob das größere Modell die zusätzliche Komplexität wert ist - die Teststatistik ist es

λ = 2 (L_{1} - L_{0})

$\lambda = 2(L_1-L_0)$

If the categorical predictor truly has no effect, $\lambda$ has a $\chi^2$ distribution with degrees of freedom equal to $mp - p = p(m-1)$ , where $p$ is the number of parameters underlying your non-linear regression function. You can use pchisq() to calculate $\chi^2$ p-values.

Macro
quelle

Are you suggesting to fit m separate models, sum the log likelihood from each L1= SUM(LL_i, i from 1 to m) and then proceed with the likelihood? Also, is L0 a model with the categorical predictor in question included (with m-1 dummy variables for example)?

B_Miner

Yes, I am suggesting that.

L_{0}

$L_0$ is the maximized likelihood when you've left

B

$B$ out of the model entirely (i.e. you're fitting the regression function for the entire data set, not stratifying by values of

B

$B$ ).

Macro

Thanks for your suggestion Macro. This seems to be in the direction of what I have already done - although you suggest comparison of likelihood rather than the F-test. In my example, the F-test also compares the single fit residuals to the summation of residuals from several fits applied to each categorical predictor level. I guess I was wondering if one can do this within a mixed model in one step rather than fitting several models. Also, would such a strategy allow for nested factor testing?

Marc in the box

I don't think you'll be able to get around fitting several models in order to compare models. Also, yes, the likelihood ratio test can be used to test for nested factors.

Macro

I found that it is possible to code categorical variables with nls(), simply by multiplying true/false vectors into your equation. Example:

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