1
$\begingroup$

I have a (lme) mixed model of this type:

model=lme(trait1~tl_s+year+month_s+I(month_s^2)+I(month_s^3),random=~1|ID,correlation=corAR1(form=~tim),data=depm3,method="REML") 

I'm interested in estimating the repeatability of the trait I'm modelling (trait1 in this case). I get a point estimate with this code:

R=as.numeric(VarCorr(mod))[1]/(as.numeric(VarCorr(mod))[1]+as.numeric(VarCorr(mod))[2])

But what would be the code to extract confidence interval for R estimate? I read that parametric bootstrapping would be an option here, but I don't really know how to implement it. I have check some recent posts on the topic like this one...

How to calculate estimated proportions and their confidence intervals from a mixed model?

...but I haven't really found the answer to my question.

Any advise would be appreciated

Thanks

$\endgroup$
1
  • $\begingroup$ There is no package to deal with lme repeatability when it includes autocorrelation and variance functions. I tried to make the function myself using the boot function for bootstrapping in order to get the confidence intervals. However, the question here is: does the autocorrelation and variance function remain valid for the new dataset after bootstrapping?? $\endgroup$
    – HFS
    Sep 22, 2021 at 10:49

1 Answer 1

0
$\begingroup$

There is no package to deal with lme repeatability when it includes autocorrelation and/or variance functions. I tried to code the function myself using the boot() function for bootstrapping in order to get the confidence intervals.

library(boot)

rsq <- function(formula,cor.b,weights.b,random.b, data,indices) {

    d <- data[indices,] # allows boot to select sample

    fit <- lme(formula, data=d,correlation=cor.b,weights=weights.b,random=random.b)

    VarComps=VarCorr(fit)

    var_G <-as.numeric(as.character(VarComps[1,1]))

    var_R <- as.numeric(as.character(VarComps[2,1]))

    var_P <- sum(as.numeric(VarComps [,1]))

    R=var_G / var_P

    return(R)

} 
results <- boot(data=depm3, statistic=rsq, R=2, formula=trait1~tl_s+year+month_s+I(month_s^2)+I(month_s^3), random.b=~1|ID, cor.b=corAR1(form=~tim))

boot.ci(results, type="bca") #to get confidence intervals

Usually, R should be around 1000 but it will take too long to be done.

However, the question here is: does the autocorrelation and variance function remain valid for the new dataset after bootstrapping?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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