# Non-linear mixed model (nlme) with nested random effect, do not know how to include nested random effect when bootstrapping

I have fitted a nonlinear model with nested random effect, and I'm trying to estimate confidence intervals on predictions by nested bootstrapping, by following this explanation. I have edited my orginal question as I first had a typos in the model specification.

I'm unsure how to modify the bootstrapping script to include nested random effect, and I hope I could get some help. I get one error message, see below.

library(nlme)
library(MASS)
fm1 <- nlme(Beat_freq ~ SSasymp(WL, Asym, R0, lrc),
data = WBF,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1|Species/ID,
start = c(Asym = 27, R0 = 7000, lrc = 0.22))
xvals <-  with(WBF,seq(min(WL),max(WL),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(WBF,data.frame(WL=xvals))
pframe$WingBeatFreq <- 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="Species/ID") { 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$WingBeatFreq <- 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,WBF,"Species/ID"))))

c2 <- get_CI(yvals2,"boot_")

pframe <- data.frame(pframe,c1,c2)
write.csv(pframe, file = "pframe.csv")


I've adapted the sampfun() function (I think it works, but not carefully tested). If nlmer were better developed, or if nlme had a useful simulate method (sigh), it would be easier to do this by parametric bootstrapping.

WBF <- read.csv("WBF.csv")
library(nlme)
fm1 <- nlme(Beat_freq ~ SSasymp(WL, Asym, R0, lrc),
data = WBF,
fixed = Asym + R0 + lrc ~ 1,
random = Asym ~ 1|Species/ID,
start = c(Asym = 27, R0 = 7000, lrc = 0.22))


Nested bootstrapping:

sampfun <- function(fitted,data,idvar=c("Species","ID"),
resp="Beat_freq") {
pp <- predict(fitted,levels=length(idvar))
rr <- residuals(fitted)
dd <- data.frame(data,pred=pp,res=rr)
## sample top-level groups with replacement
iv1 <- levels(dd[[idvar[1]]])
bsamp1 <- sample(iv1,size=length(iv1),replace=TRUE)
## sample next level, *within* top-level groups
bsamp2 <- lapply(bsamp1,
function(id1) {
iv2 <- unique(as.character(dd[dd[[idvar[1]]]==id1,idvar[2]]))
sample(iv2,size=length(iv2),replace=TRUE)
})
## sample at lowest level:
## for loop rather than nested lapply() to reduce(?) confusion
bsamp3 <- list()
for (i in seq_along(bsamp2)) {
bsamp3[[i]] <- lapply(bsamp2[[i]],
function(x) {
ddb <- dd[dd[[idvar[1]]]==bsamp1[i] &
dd[[idvar[2]]]==as.character(x),]
ddb[[resp]] <- ddb$pred + sample(ddb$res,size=nrow(ddb),replace=TRUE)
return(ddb)
})
}
## flatten everything
res <- do.call(rbind,lapply(bsamp3,
function(x) do.call(rbind,x)))
if (is(data,"groupedData"))
res <- groupedData(res,formula=formula(data))
return(res)
}

sampfun(fm1,WBF,c("Species","ID"))

• aha. It would help to be a little bit more explicit in your question that you are running into problems at the bootstrapping stage (I didn't follow the link to see the whole thing)! Yes, the bootstrapping process will have to be adjusted ... it's probably going to get considerably more complex, as you have to implement bootstrapping nested at both levels. Sep 1 '16 at 23:11
• Thanks. I will change the question and be more more explicit.
– Kes
Sep 2 '16 at 6:52
• Thanks. The adapted sampfun() function works great. Although, I'm unsure how to adapt rest of the code to get lwr and upr bootstrap set.seed(101) yvals2 <- replicate(nresamp, pfun(update(fm1,data=sampfun(fm1,WBF,"Species/ID")))) c2 <- get_CI(yvals2,"boot_")
– Kes
Sep 3 '16 at 16:35