5
$\begingroup$

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.

You can download the data set here.

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")))) ## Error in eval(expr, envir, enclos) : object 'ID' not found c2 <- get_CI(yvals2,"boot_") pframe <- data.frame(pframe,c1,c2) write.csv(pframe, file = "pframe.csv")
Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

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"))
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ Commented Sep 1, 2016 at 23:11
  • $\begingroup$ Thanks. I will change the question and be more more explicit. $\endgroup$ Commented Sep 2, 2016 at 6:52
  • $\begingroup$ 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_") $\endgroup$ Commented Sep 3, 2016 at 16:35

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.