R

非線性混合模型 (nlme) 預測的置信區間

  • August 22, 2016

我想獲得非線性混合nlme模型預測的 95% 置信區間。由於沒有提供任何標準來做到這一點nlme，我想知道使用“人口預測區間”的方法是否正確，正如Ben Bolker 的書一章中在模型擬合最大似然的上下文中概述的那樣，基於以下思想根據擬合模型的方差-協方差矩陣對固定效應參數進行重採樣，基於此模擬預測，然後取這些預測的 95% 的百分位數以獲得 95% 的置信區間？

執行此操作的代碼如下所示：（我在這裡使用nlme幫助文件中的“Lobolly”數據）

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% 置信區間？我不完全確定如何處理模型的隨機效應結構……是否應該平均超過隨機效應水平？或者是否可以為一個普通主題設置置信區間，這似乎更接近我現在的置信區間？

你在這裡所做的看起來很合理。簡短的回答是，在大多數情況下，從混合模型和非線性模型預測置信區間的問題或多或少是正交的，也就是說，你需要擔心這兩組問題，但他們不需要（我知道of) 以任何奇怪的方式進行交互。

  • 混合模型問題：您是要在總體還是群體層面進行預測？您如何解釋隨機效應參數的可變性？您是否以小組級別的觀察為條件？
  • 非線性模型問題：參數的抽樣分佈是否正常？傳播誤差時如何考慮非線性？

在整個過程中，我將假設您在總體水平上進行預測並將置信區間構建為總體水平 - 換句話說，您正在嘗試繪製典型組的預測值，而不包括您的置信度中的組間變化間隔。這簡化了混合模型問題。下圖比較了三種方法（代碼轉儲見下文）：

  • 人口預測區間：這是您在上面嘗試的方法。它假設模型是正確的，並且固定效應參數的抽樣分佈是多元正態分佈；它還忽略了隨機效應參數的不確定性
  • bootstrapping：我實現了分層引導；我們在組級別和組內重新採樣。組內抽樣對殘差進行採樣並將它們添加回預測。這種方法做出的假設最少。
  • delta方法：這假設採樣分佈的多元正態性和非線性足夠弱以允許二階近似。

我們也可以做參數引導…

這是與數據一起繪製的 CI…

在此處輸入圖像描述

…但我們幾乎看不到差異。

通過減去預測值來放大（red=bootstrap，blue=PPI，cyan=delta 方法）

在此處輸入圖像描述

在這種情況下，bootstrap 間隔實際上是最窄的（例如，可能參數的採樣分佈實際上比 Normal 稍微細尾），而 PPI 和 delta 方法間隔彼此非常相似。

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)

引用自：https://stats.stackexchange.com/questions/231074

comments powered by Disqus

相關問答

Confidence-Interval

當您擁有全部人口時，是否適合在數據上放置“誤差線”？

  • January 21, 2022
查看問答
Hypothesis-Testing

關於文章“拋棄 p 值。改用 Bootstrap 置信區間”的三個問題

  • November 9, 2021
查看問答
Normal-Distribution

估計參數函數周圍的置信區間

  • November 8, 2021
查看問答
Confidence-Interval

計算兩個獨立比例之差的標準誤

  • October 10, 2021
查看問答
R

如何從二元邏輯回歸模型中獲得兩個概率之間差異的置信區間？

  • August 10, 2021
查看問答
Self-Study

Z值可以被認為是標準偏差的數量嗎？

  • July 29, 2021
查看問答