Доверительные интервалы прогнозов для нелинейной смешанной модели (nlme)

Я хотел бы получить 95% доверительные интервалы на предсказаниях нелинейной смешанной nlmeмодели. Поскольку для этого не предусмотрено ничего стандартного nlme, мне хотелось бы знать, правильно ли использовать метод «интервалов прогнозирования населения», как описано в главе книги Бена Болкера в контексте моделей, подходящих с максимальной вероятностью , основанных на идее пересобрать параметры фиксированного эффекта на основе матрицы дисперсии-ковариации подобранной модели, смоделировать прогнозы на основе этого, а затем взять 95% -ный процентиль этих прогнозов, чтобы получить 95% доверительные интервалы?

Код для этого выглядит следующим образом: (здесь я использую данные «Loblolly» из nlmeфайла справки)

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

Теперь, когда у меня есть пределы доверия, я создаю график:

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

Вот график с 95% доверительными интервалами, полученными таким образом:

Является ли этот подход допустимым, или существуют ли другие или более эффективные подходы для расчета 95% доверительных интервалов для прогнозов нелинейной смешанной модели? Я не совсем уверен, что делать со структурой случайных эффектов модели ... Следует ли усреднять, возможно, уровни случайных эффектов? Или было бы нормально иметь доверительные интервалы для среднего субъекта, которые, казалось бы, были ближе к тому, что у меня сейчас?

То, что вы сделали здесь, выглядит разумным. Короткий ответ заключается в том, что по большей части проблемы прогнозирования доверительных интервалов из смешанных моделей и из нелинейных моделей являются более или менее ортогональными , то есть вам нужно беспокоиться об обоих наборах проблем, но они этого не делают (я знаю о) взаимодействовать любым странным образом.

• Проблемы смешанной модели : вы пытаетесь прогнозировать на уровне населения или группы? Как вы учитываете изменчивость параметров случайных эффектов? Вы обусловливаете наблюдения на уровне группы или нет?
• Проблемы нелинейной модели : является ли выборочное распределение параметров нормальным? Как мне учесть нелинейность при распространении ошибки?

Везде, я буду предполагать, что вы прогнозируете на уровне населения и строите доверительные интервалы как уровень населения - другими словами, вы пытаетесь построить прогнозируемые значения типичной группы, не включая различия между группами в своей уверенности интервалы. Это упрощает проблемы смешанной модели. На следующих графиках сравниваются три подхода (см. Ниже дамп кода):

• интервалы прогнозирования численности населения : этот подход вы опробовали выше. Предполагается, что модель верна и что выборочные распределения параметров с фиксированным эффектом являются многомерными Normal; он также игнорирует неопределенность в параметрах случайных эффектов
• начальная загрузка : я реализовал иерархическую загрузку; мы повторяем выборку как на уровне групп, так и внутри групп. Внутригрупповая выборка выборки остатков и добавляет их обратно в прогнозы. Этот подход делает наименьшее количество предположений.
• дельта-метод : он предполагает как многомерную нормальность распределений выборки, так и нелинейность, достаточно слабую, чтобы допустить приближение второго порядка.

Мы могли бы также сделать параметрическую загрузку ...

Вот графики, нанесенные вместе с данными ...

... но мы вряд ли можем увидеть различия.

Увеличение масштаба путем вычитания прогнозируемых значений (красный = начальная загрузка, синий = ИЦП, голубой = дельта-метод)

В этом случае интервалы начальной загрузки на самом деле самые узкие (например, предположительно, распределения выборки параметров на самом деле немного тоньше, чем нормальные), в то время как интервалы PPI и дельта-метода очень похожи друг на друга.

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)