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

Question

Me gustaría obtener el 95% de intervalos de confianza en las predicciones de un no-lineal mixto nlme modelo. Como nada estándar se proporciona para ello, dentro de nlme, me pregunto si es correcto usar el método de "población intervalos de predicción", como se describe en Ben Bolker del capítulo del libro en el contexto de los modelos de ajuste con la máxima probabilidad, basado en la idea de repetición de muestreo de los parámetros de efecto fijo basado en la varianza-covarianza de la matriz del modelo ajustado, la simulación de las predicciones basadas en el este, y luego tomar el 95% de los percentiles de estas predicciones para obtener el 95% los intervalos de confianza?

El código para hacer esto se ve de la siguiente manera : (He aquí el uso de la 'Incienso' datos de la nlme archivo 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 mi límites de confianza puedo crear 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 parcela con el 95% de intervalos de confianza obtenidos de esta manera:

Es este enfoque válido, o hay alguna otra o de mejores métodos para calcular el 95% de intervalos de confianza en las predicciones de una relación no lineal modelo mixto? No estoy del todo seguro de cómo lidiar con el efecto aleatorio estructura de modelo... en caso de un promedio de quizás más de efecto aleatorio de los niveles? O sería ACEPTAR para que los intervalos de confianza para una media de tema, que parece estar más cerca de lo que tengo ahora?

Answer 1

Lo que has hecho aquí parece razonable. La respuesta corta es que la mayor parte de los problemas de la predicción intervalos de confianza de los modelos mixtos y de modelos no lineales son más o menos ortogonal, es decir, lo que usted necesita preocuparse acerca de los dos conjuntos de problemas, pero no (que yo sepa) interactuar de alguna manera extraña.

• Modelo mixto cuestiones: ¿estás tratando de predecir en la población o el nivel de grupo? ¿Cómo se cuenta la variabilidad en el de efectos aleatorios parámetros? Estás acondicionado en el nivel del grupo de observaciones o no?
• No lineales del modelo cuestiones: es la distribución de muestreo de los parámetros Normales? ¿Cómo se cuenta la no linealidad cuando la propagación de error?

Todo, voy a suponer que usted pronostica que en el nivel de la población y la construcción de intervalos de confianza como el nivel de la población - en otras palabras, usted está tratando de graficar los valores de la predicción de un típico grupo, y no como el entre-grupo variación en los intervalos de confianza. Esto simplifica el modelo mixto de problemas. Los siguientes grácos comparar los tres enfoques (ver más abajo para código de volcado de memoria):

• población en intervalos de predicción: este es el enfoque que se trató anteriormente. Se asume que el modelo es correcto y que la distribución de muestras de la fijo-parámetros de efecto multivariante Normal; también se ignora la incertidumbre en el de efectos aleatorios parámetros
• bootstrapping: he implementado jerárquica de arranque; nos remuestrear tanto a nivel de los grupos y dentro de grupos. Dentro del grupo de muestreo las muestras de los residuos y se agregan a las predicciones. Este enfoque hace que el menor número de suposiciones.
• método delta: se supone que ambos Normalidad multivariante de muestreo de distribuciones y que la no linealidad es suficientemente débil para permitir un segundo orden de la aproximación.

También podríamos hacer paramétrico de arranque ...

Aquí están los CIs se trazan a lo largo con los datos ...

... pero apenas podemos ver las diferencias.

Zoom en restando fuera de los valores pronosticados (rojo=bootstrap, azul=PPI, cian=método delta)

En este caso, el bootstrap intervalos son en realidad más estrecha (por ejemplo, es de suponer que el muestreo de las distribuciones de los parámetros son en realidad un poco más delgada de cola de lo Normal), mientras que el PPI y delta-método de los intervalos 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)