Konfidenzintervalle für Vorhersagen für ein nichtlineares gemischtes Modell (nlme)

Ich möchte 95% -Konfidenzintervalle für die Vorhersagen eines nichtlinearen gemischten nlmeModells erhalten. Da dies innerhalb von nichts Standardmäßigem vorgesehen ist, habe nlmeich mich gefragt, ob es richtig ist, die Methode der "Bevölkerungsvorhersageintervalle" zu verwenden, die in Ben Bolkers Buchkapitel im Kontext von Modellen beschrieben wird , die auf der Idee von mit maximaler Wahrscheinlichkeit passen Resampling von festen Effektparametern basierend auf der Varianz-Kovarianz-Matrix des angepassten Modells, Simulation von darauf basierenden Vorhersagen und anschließende Ermittlung der 95% -Perzentile dieser Vorhersagen, um die 95% -Konfidenzintervalle zu erhalten?

Der Code dafür sieht folgendermaßen aus: (Ich verwende hier die 'Loblolly'-Daten aus der nlmeHilfedatei.)

library(effects)
library(nlme)
library(MASS)

fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
    data = Loblolly,
    fixed = Asym + R0 + lrc ~ 1,
    random = Asym ~ 1,
    start = c(Asym = 103, R0 = -8.5, lrc = -3.3))

xvals=seq(min(Loblolly$age),max(Loblolly$age),length.out=100)
nresamp=1000
pars.picked = mvrnorm(nresamp, mu = fixef(fm1), Sigma = vcov(fm1)) # pick new parameter values by sampling from multivariate normal distribution based on fit
yvals = matrix(0, nrow = nresamp, ncol = length(xvals))

for (i in 1:nresamp) 
{
    yvals[i,] = sapply(xvals,function (x) SSasymp(x,pars.picked[i,1], pars.picked[i,2], pars.picked[i,3]))
} 

quant = function(col) quantile(col, c(0.025,0.975)) # 95% percentiles
conflims = apply(yvals,2,quant) # 95% confidence intervals

Jetzt, da ich meine Vertrauensgrenzen habe, erstelle ich ein Diagramm:

meany = sapply(xvals,function (x) SSasymp(x,fixef(fm1)[[1]], fixef(fm1)[[2]], fixef(fm1)[[3]]))

par(cex.axis = 2.0, cex.lab=2.0)
plot(0, type='n', xlim=c(3,25), ylim=c(0,65), axes=F, xlab="age", ylab="height");
axis(1, at=c(3,1:5 * 5), labels=c(3,1:5 * 5)) 
axis(2, at=0:6 * 10, labels=0:6 * 10)   

for(i in 1:14)
{
    data = subset(Loblolly, Loblolly$Seed == unique(Loblolly$Seed)[i])   
    lines(data$age, data$height, col = "red", lty=3)
}

lines(xvals,meany, lwd=3)
lines(xvals,conflims[1,])
lines(xvals,conflims[2,])

Hier ist die grafische Darstellung mit den so erhaltenen 95% -Konfidenzintervallen:

Ist dieser Ansatz gültig, oder gibt es andere oder bessere Ansätze, um 95% -Konfidenzintervalle für die Vorhersagen eines nichtlinearen gemischten Modells zu berechnen? Ich bin mir nicht ganz sicher, wie ich mit der zufälligen Effektstruktur des Modells umgehen soll ... Sollte man vielleicht über zufällige Effektniveaus mitteln? Oder wäre es in Ordnung, Konfidenzintervalle für ein durchschnittliches Thema zu haben, die näher an dem zu liegen scheinen, was ich jetzt habe?

Was Sie hier getan haben, sieht vernünftig aus. Die kurze Antwort lautet, dass die Probleme bei der Vorhersage von Konfidenzintervallen aus gemischten Modellen und aus nichtlinearen Modellen größtenteils mehr oder weniger orthogonal sind auf seltsame Weise interagieren.

• Gemischte Modellprobleme : Versuchen Sie, auf Bevölkerungs- oder Gruppenebene Vorhersagen zu treffen? Wie berücksichtigen Sie die Variabilität der Zufallseffektparameter? Konditionieren Sie Beobachtungen auf Gruppenebene oder nicht?
• Nichtlineare Modellprobleme : Ist die Stichprobenverteilung der Parameter normal? Wie berücksichtige ich die Nichtlinearität bei der Fehlerweitergabe?

Ich gehe davon aus, dass Sie auf Bevölkerungsebene Vorhersagen treffen und Konfidenzintervalle als Bevölkerungsebene konstruieren - mit anderen Worten, Sie versuchen, die vorhergesagten Werte einer typischen Gruppe zu zeichnen , ohne die Unterschiede zwischen den Gruppen in Ihrem Vertrauen einzubeziehen Intervalle. Dies vereinfacht die Probleme mit gemischten Modellen. In den folgenden Darstellungen werden drei Ansätze verglichen (siehe unten für Codedump):

• Bevölkerungsvorhersageintervalle : Dies ist der Ansatz, den Sie oben versucht haben. Es wird davon ausgegangen, dass das Modell korrekt ist und dass die Stichprobenverteilungen der Parameter mit festem Effekt multivariate Normalverteilungen sind. Außerdem werden Unsicherheiten bei den Zufallseffektparametern ignoriert
• Bootstrapping : Ich habe hierarchisches Bootstrapping implementiert. wir resampling sowohl auf gruppenebene als auch innerhalb von gruppen. Die gruppeninterne Stichprobe tastet die Residuen ab und fügt sie den Vorhersagen hinzu. Dieser Ansatz macht die wenigsten Annahmen.
• Delta-Methode : Hierbei wird sowohl eine multivariate Normalität der Stichprobenverteilungen als auch eine Nichtlinearität vorausgesetzt, die schwach genug ist, um eine Approximation zweiter Ordnung zu ermöglichen.

Wir könnten auch parametrisches Bootstrapping durchführen ...

Hier sind die CIs zusammen mit den Daten aufgetragen ...

... aber wir können kaum Unterschiede erkennen.

Vergrößern durch Abziehen der vorhergesagten Werte (Rot = Bootstrap, Blau = PPI, Cyan = Delta-Methode)

In diesem Fall sind die Bootstrap-Intervalle tatsächlich am engsten (z. B. sind die Stichprobenverteilungen der Parameter vermutlich tatsächlich etwas dünner als normal), während die Intervalle für die PPI- und Delta-Methode einander sehr ähnlich sind.

library(nlme)
library(MASS)

fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
            data = Loblolly,
            fixed = Asym + R0 + lrc ~ 1,
            random = Asym ~ 1,
            start = c(Asym = 103, R0 = -8.5, lrc = -3.3))

xvals <-  with(Loblolly,seq(min(age),max(age),length.out=100))
nresamp <- 1000
## pick new parameter values by sampling from multivariate normal distribution based on fit
pars.picked <- mvrnorm(nresamp, mu = fixef(fm1), Sigma = vcov(fm1))

## predicted values: useful below
pframe <- with(Loblolly,data.frame(age=xvals))
pframe$height <- predict(fm1,newdata=pframe,level=0)

## utility function
get_CI <- function(y,pref="") {
    r1 <- t(apply(y,1,quantile,c(0.025,0.975)))
    setNames(as.data.frame(r1),paste0(pref,c("lwr","upr")))
}

set.seed(101)
yvals <- apply(pars.picked,1,
               function(x) { SSasymp(xvals,x[1], x[2], x[3]) }
)
c1 <- get_CI(yvals)

## bootstrapping
sampfun <- function(fitted,data,idvar="Seed") {
    pp <- predict(fitted,levels=1)
    rr <- residuals(fitted)
    dd <- data.frame(data,pred=pp,res=rr)
    ## sample groups with replacement
    iv <- levels(data[[idvar]])
    bsamp1 <- sample(iv,size=length(iv),replace=TRUE)
    bsamp2 <- lapply(bsamp1,
        function(x) {
        ## within groups, sample *residuals* with replacement
        ddb <- dd[dd[[idvar]]==x,]
        ## bootstrapped response = pred + bootstrapped residual
        ddb$height <- ddb$pred +
            sample(ddb$res,size=nrow(ddb),replace=TRUE)
        return(ddb)
    })
    res <- do.call(rbind,bsamp2)  ## collect results
    if (is(data,"groupedData"))
        res <- groupedData(res,formula=formula(data))
    return(res)
}

pfun <- function(fm) {
    predict(fm,newdata=pframe,level=0)
}

set.seed(101)
yvals2 <- replicate(nresamp,
                    pfun(update(fm1,data=sampfun(fm1,Loblolly,"Seed"))))
c2 <- get_CI(yvals2,"boot_")

## delta method
ss0 <- with(as.list(fixef(fm1)),SSasymp(xvals,Asym,R0,lrc))
gg <- attr(ss0,"gradient")
V <- vcov(fm1)
delta_sd <- sqrt(diag(gg %*% V %*% t(gg)))
c3 <- with(pframe,data.frame(delta_lwr=height-1.96*delta_sd,
                             delta_upr=height+1.96*delta_sd))

pframe <- data.frame(pframe,c1,c2,c3)

library(ggplot2); theme_set(theme_bw())
ggplot(Loblolly,aes(age,height))+
    geom_line(alpha=0.2,aes(group=Seed))+
    geom_line(data=pframe,col="red")+
    geom_ribbon(data=pframe,aes(ymin=lwr,ymax=upr),colour=NA,alpha=0.3,
                fill="blue")+
    geom_ribbon(data=pframe,aes(ymin=boot_lwr,ymax=boot_upr),
                colour=NA,alpha=0.3,
                fill="red")+
    geom_ribbon(data=pframe,aes(ymin=delta_lwr,ymax=delta_upr),
                colour=NA,alpha=0.3,
                fill="cyan")


ggplot(Loblolly,aes(age))+
    geom_hline(yintercept=0,lty=2)+
    geom_ribbon(data=pframe,aes(ymin=lwr-height,ymax=upr-height),
                colour="blue",
                fill=NA)+
    geom_ribbon(data=pframe,aes(ymin=boot_lwr-height,ymax=boot_upr-height),
                colour="red",
                fill=NA)+
    geom_ribbon(data=pframe,aes(ymin=delta_lwr-height,ymax=delta_upr-height),
                colour="cyan",
                fill=NA)