Intervalos de confianza en las predicciones para un modelo mixto no lineal (nlme)

Me gustaría obtener intervalos de confianza del 95% en las predicciones de un nlmemodelo mixto no lineal . Como no se proporciona nada estándar para hacer esto dentro nlme, me preguntaba si es correcto usar el método de "intervalos de predicción de población", como se describe en el capítulo del libro de Ben Bolker en el contexto de modelos ajustados con la máxima probabilidad , basados ​​en la idea de remuestreo de parámetros de efectos fijos basados ​​en la matriz de varianza-covarianza del modelo ajustado, simulando predicciones basadas en esto y luego tomando los percentiles del 95% de estas predicciones para obtener los intervalos de confianza del 95%.

El código para hacer esto tiene el siguiente aspecto: (Aquí uso los datos 'Loblolly' del nlmearchivo de ayuda)

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

Ahora que tengo mis límites de confianza, creo un gráfico:

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

Aquí está la gráfica con los intervalos de confianza del 95% obtenidos de esta manera:

¿Es válido este enfoque, o existen otros enfoques mejores para calcular los intervalos de confianza del 95% en las predicciones de un modelo mixto no lineal? No estoy del todo seguro de cómo lidiar con la estructura de efecto aleatorio del modelo ... ¿Debería uno promediar tal vez por encima de los niveles de efecto aleatorio? ¿O estaría bien tener intervalos de confianza para un sujeto promedio, que parece estar más cerca de lo que tengo ahora?

Lo que has hecho aquí parece razonable. La respuesta corta es que, en su mayor parte, los problemas de predicción de intervalos de confianza a partir de modelos mixtos y de modelos no lineales son más o menos ortogonales , es decir, debe preocuparse por ambos conjuntos de problemas, pero no es así (eso sé) de) interactuar de alguna manera extraña.

• Problemas de modelos mixtos : ¿está tratando de predecir a nivel de población o de grupo? ¿Cómo se explica la variabilidad en los parámetros de efectos aleatorios? ¿Estás condicionando las observaciones a nivel de grupo o no?
• Problemas del modelo no lineal : ¿es normal la distribución de muestreo de los parámetros? ¿Cómo tomo en cuenta la no linealidad cuando propago un error?

En todo momento, supondré que está prediciendo a nivel de población y construyendo intervalos de confianza como el nivel de población; en otras palabras, está tratando de trazar los valores pronosticados de un grupo típico , y no incluye la variación entre grupos en su confianza intervalos. Esto simplifica los problemas del modelo mixto. Las siguientes gráficas comparan tres enfoques (ver más abajo el volcado de código):

• intervalos de predicción de población : este es el enfoque que probó anteriormente. Se supone que el modelo es correcto y que las distribuciones de muestreo de los parámetros de efectos fijos son multivariadas Normal; también ignora la incertidumbre en los parámetros de efectos aleatorios
• bootstrapping : implementé bootstrapping jerárquico; remuestreamos tanto a nivel de grupos como dentro de grupos. El muestreo dentro del grupo muestrea los residuos y los agrega de nuevo a las predicciones. Este enfoque hace la menor cantidad de suposiciones.
• Método delta : esto supone tanto la normalidad multivariada de las distribuciones de muestreo como que la no linealidad es lo suficientemente débil como para permitir una aproximación de segundo orden.

También podríamos hacer bootstrapping paramétrico ...

Aquí están los IC trazados junto con los datos ...

... pero apenas podemos ver las diferencias.

Acercarse restando los valores predichos (rojo = bootstrap, azul = PPI, cian = método delta)

En este caso, los intervalos de arranque son realmente más estrechos (p. Ej., Presumiblemente, las distribuciones de muestreo de los parámetros son en realidad un poco más delgadas que las normales), mientras que los intervalos PPI y del método delta son muy similares entre sí.

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)