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")
$\endgroup$
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"))
$\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$
    – Ben Bolker
    Sep 1 '16 at 23:11
  • $\begingroup$ Thanks. I will change the question and be more more explicit. $\endgroup$
    – Kes
    Sep 2 '16 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$
    – Kes
    Sep 3 '16 at 16:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.