使用非线性混合效应模型 (nlme) 的预测变量进行预测的置信区间

Question

我正在尝试使用 Ben Bolker here 描述的引导方法为非线性混合效应模型生成 95% 的置信区间。我可能误解了所使用的一些功能。我的目标是模拟预测变量每个级别（在本示例中为年份）的 95% 置信区间。

下面是使用

FlexParamCurve

包的数据集

penguin.data

的可重现代码。需要明确的是，提供的代码是 Bolker 博士在回答上述链接问题时提供的代码的修改版本。

library(FlexParamCurve) #also loads package 'nlme' which is needed.
library(ggplot2)

set.seed(1234)

##creating model
fm2 <- nlme(weight ~ SSlogis(ckage, Asym, R0, lrc),
            data = penguin.data,
            fixed= list(Asym ~ year,
                        R0 ~ year,
                        lrc ~ year),
            random = Asym ~ 1,
            start = c(Asym = 1000, 0,
                      R0 = 21, 0,
                      lrc = 1, 0),
            control = list(maxIter = 100),
            na.action = na.pass)

#created simulated x ('ckage') values to use in prediction below
xvals.peng <- with(penguin.data,seq(min(ckage),max(ckage),length.out=100))

nresamp <- 100

## 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"))) #function to get CI 
}

##creating the data frame to use for predictions
pengframe <- with(penguin.data, data.frame(ckage = xvals.peng))

##Tried to use for weight predictions and it did not work
pengframe$weight <- predict(fm2,newdata=pengframe,level=0)

Error in eval(predvars, data, env) : object 'year' not found

这不起作用，因为模型中的固定效应“年份”在

pengframe

中缺失，因此我无法使用

predict()

函数。这完全有道理，所以我尝试使用

rbind()

函数解决方法：

###this is where I created the column (year) and its values separately before rbinding the data frames. There are only two levels in year.
pengframe1 <- with(penguin.data,data.frame(ckage=xvals.peng))
pengframe1$year <- as.factor('2000')
pengframe2 <- with(penguin.data,data.frame(ckage=xvals.peng))
pengframe2$year <- as.factor('2002')

pengframe <- rbind(pengframe1, pengframe2)

#predicting weight for each year now works
pengframe$weight <- predict(fm2,newdata=pengframe,level=0)

head(pengframe)


sampfun <- function(fitted,data,idvar="bandid") {
  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=pengframe,level=0)
}

yvals2 <- replicate(nresamp, 
                    pfun(update(fm2, 
                                data = sampfun(fm2, 
                                               penguin.data, 
                                               "bandid"))))

peng2 <- get_CI(yvals2,"boot_")
head(peng2)
pengframe <- data.frame(pengframe,peng2)
head(pengframe)

ggplot(pengframe, aes(ckage, weight, color = year)) + 
  geom_smooth() + #this is for simplicity purposes only. I use geom_func() in my real dataset
  geom_ribbon(pengframe, mapping = aes(x = ckage, ymin = boot_lwr, ymax =boot_upr, group=year, fill = year), alpha = 0.3)

此方法通过使用与其他方法相同的估计“ckage”，为我提供了每年 95% 的置信区间，这是我的预期结果。

我想确认这种方式在统计上是否合理？我怀疑这种方法还可以，但我才刚刚开始掌握非线性混合模型。我还想问是否有更直接的方法[即将其添加到

with()

函数中，在首次创建

xvals.peng

时最初模拟“ckage”以简化流程]。我将使用性别而不是年份，并且我将嵌套随机因素（1 |组/id），这可能完全是一个不同的问题。

Answer 1

我觉得基本上没问题。但是您可以使用

expand.grid()

以更便宜的价格获得新数据。此外，您还可以保持工作空间更干净：实际上您只需要一个可以扩展的新 data.frame

newdata.peng

。

> ## create newdata
> newdata.peng <- with(penguin.data, 
+                      expand.grid(
+                        ckage=seq(min(ckage), max(ckage), length.out=100), 
+                        year=unique(year)
+                      ))
> 
> 
> ## add predicted weight to newdata
> newdata.peng$weight <- predict(fm2, newdata=newdata.peng, level=0)
> 
> 
> ## bootstrap
> nresamp <- 100
> set.seed(1234)
> yvals2 <- 
+   replicate(nresamp, {
+     pfun(update(fm2, data=sampfun(fm2, penguin.data, "bandid")))
+   })

也许您需要更多数量级的复制。而且它会永远运行，如果您使用的是 Linux，请尝试

parallel::parSapply

或更好的

parallel::mclapply

。

> ## add CI weight to newdata
> newdata.peng <- cbind(newdata.peng, get_CI(yvals2, "boot_"))

不知道为什么你想在绘图时平滑某些东西或类似的东西。我认为你已经有了明确的数据，不再需要进行后处理。这是一种没有 ggplot 的方法。

> ## plot
> ci_cols <- c('boot_lwr', 'weight', 'boot_upr')
> years <- sort(as.integer(as.character(unique(newdata.peng$year))))
> 
> par(mar=c(4, 4, 2, 2) + .1)
> plot.new(); plot.window(xlim=range(newdata.peng$ckage), 
+                         ylim=range(newdata.peng[ci_cols]))
> for (i in 1:2) axis(i, axTicks(i))
> for (i in 1:2) mtext(c('chick age (days)', 'chick mass (g)')[i], i, 3)
> for (i in (seq_along(years))) {
+   matlines(
+     subset(newdata.peng, year == years[i], select=ci_cols), 
+     col=i + 1L, lty=2:1)
+ }
> box()
> legend('bottomright', legend=c(years, '95% CI'), col=c(2:3, 8), lty=c(1, 1, 2))

---

数据：

> data('penguin.data', package='FlexParamCurve')
> library(nlme)
> fm2 <- nlme(weight ~ SSlogis(ckage, Asym, R0, lrc), data=penguin.data, 
+             fixed=list(Asym ~ year, R0 ~ year, lrc ~ year), 
+             random=Asym ~ 1, start=c(Asym=1000, 0, R0=21, 0, lrc=1, 0), 
+             control=list(maxIter=100), na.action=na.pass)

功能取自@BenBolker的帖子。