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.
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")

 A: 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"))

