Przedziały ufności dotyczące prognoz dla nieliniowego modelu mieszanego (NLM)

Chciałbym uzyskać 95% przedziały ufności dla prognoz nieliniowego nlmemodelu mieszanego . Ponieważ nie ma w tym celu standardu nlme, zastanawiałem się, czy słuszne jest zastosowanie metody „przedziałów prognozowania populacji”, jak opisano w rozdziale książki Bena Bolkera w kontekście modeli pasujących z najwyższym prawdopodobieństwem , w oparciu o ideę ponownie próbkować parametry ustalonego efektu w oparciu o macierz wariancji-kowariancji dopasowanego modelu, symulując oparte na tym prognozy, a następnie biorąc 95% percentyle tych prognoz, aby uzyskać 95% przedziały ufności?

Kod do wykonania tego wygląda następująco: (Używam tutaj danych „Loblolly” z nlmepliku pomocy)

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

Teraz, gdy mam już granice pewności, tworzę wykres:

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,])

Oto wykres z 95% przedziałami ufności uzyskanymi w ten sposób:

Czy to podejście jest poprawne, czy istnieją jakieś inne lub lepsze podejścia do obliczania 95% przedziałów ufności dla prognoz nieliniowego modelu mieszanego? Nie jestem całkowicie pewien, jak poradzić sobie ze strukturą efektu losowego modelu ... Czy powinna być średnia ponad poziomy efektów losowych? Czy też byłoby OK mieć przedziały ufności dla przeciętnego przedmiotu, które wydawałyby się bliższe temu, co mam teraz?

To, co tu zrobiłeś, wygląda rozsądnie. Krótka odpowiedź jest taka, że ​​w przeważającej części kwestie przewidywania przedziałów ufności z modeli mieszanych i modeli nieliniowych są mniej więcej ortogonalne , to znaczy, że musisz się martwić o oba zestawy problemów, ale tak nie jest (wiem o tym ) wchodzą w interakcje w jakikolwiek dziwny sposób.

• Problemy modelu mieszanego : czy próbujesz przewidzieć na poziomie populacji czy grupy? Jak uwzględnić zmienność parametrów efektów losowych? Czy uzależniasz się od obserwacji na poziomie grupy, czy nie?
• Problemy z modelem nieliniowym : czy rozkład próbkowania parametrów jest normalny? Jak uwzględnić nieliniowość podczas propagowania błędu?

W całym założeniu zakładam, że przewidujesz na poziomie populacji i konstruujesz przedziały ufności jako poziom populacji - innymi słowy, próbujesz wykreślić przewidywane wartości typowej grupy, a nie uwzględniać zmienności między grupami w zakresie Twojej pewności interwały. Upraszcza to problemy z modelem mieszanym. Poniższe wykresy porównują trzy podejścia (zrzut kodu poniżej):

• przedziały przewidywania populacji : to podejście wypróbowałeś powyżej. Zakłada się, że model jest poprawny i że rozkłady próbkowania parametrów o stałym efekcie są wielowymiarowe Normalne; ignoruje również niepewność parametrów efektów losowych
• ładowanie : zaimplementowałem hierarchiczne ładowanie; dokonujemy ponownego próbkowania zarówno na poziomie grup, jak i wewnątrz grup. Próbkowanie wewnątrz grupy pobiera próbki pozostałości i dodaje je z powrotem do prognoz. Takie podejście przyjmuje najmniej założeń.
• metoda delta : zakłada to zarówno wielowymiarową normalność rozkładów próbkowania, jak i to, że nieliniowość jest wystarczająco słaba, aby umożliwić przybliżenie drugiego rzędu.

Możemy również zrobić ładowanie parametryczne ...

Oto wykresy CI wraz z danymi ...

... ale trudno nam dostrzec różnice.

Powiększanie odejmując przewidywane wartości (czerwony = bootstrap, niebieski = PPI, cyan = delta)

W tym przypadku interwały ładowania są w rzeczywistości najwęższe (np. Przypuszczalnie rozkłady próbkowania parametrów są w rzeczywistości nieco cieńsze niż normalne), podczas gdy interwały PPI i metody delta są bardzo do siebie podobne.

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)