Hi there. OSX/R2.0
We are trying to implement a bootstrap of the coeffecients of a mixed effect model. In particular, we are interested in the intercept and slope of the random effects.
Following from the basics for a linear model, we construct our lme models and a boot function:
library(nlme)
library(boot)
data<-read.csv("~/data.csv")bootcoef<-function(data,index){
dat<-data[index,]
mod<-lme(Frames~Man2+Manip+Strings+Date.+Cut.,random=~Man2|ID.,data=dat)
fixef(mod)
}
boot.out<-boot(data,bootcoef,99)
boot.ci(boot.out) # produces information via boot.ci() that suggests this is not necessarily successful
ORDINARY NONPARAMETRIC BOOTSTRAP Call: boot(data = data, statistic = bootcoef, R = 99)
Bootstrap Statistics :
original bias std. error
t1* 16.125904015 5.299827478 9.98818463
t2* 0.010901682 -0.004134585 0.01621935
t3* -0.038168126 -0.078833467 0.35778286
t4* 1.101486342 -0.021886290 0.45720400
t5* 0.005982241 -0.009140839 0.01563175
t6* 2.729537567 0.287663533 1.56150779
> boot.ci(boot.out)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 99 bootstrap replicates
CALL :
boot.ci(boot.out = boot.out)
Intervals :
Level Normal Basic Studentized
95% ( -8.75, 30.40 ) ( -8.40, 36.01 ) (-40.93, 35.15 )
Level Percentile BCa
95% (-3.75, 40.65 ) (-5.55, 28.65 )
Calculations and Intervals on Original Scale
Some basic intervals may be unstable
Some studentized intervals may be unstable
Some percentile intervals may be unstable
Warning : BCa Intervals used Extreme Quantiles
Some BCa intervals may be unstable
Warning messages:
1: NaNs produced in: sqrt(tv[, 2])
2: Extreme Order Statistics used as Endpoints in: norm.inter(t, adj.alpha)
As stated above, we are interested in the ranef(mod) components..... including this instead of fixef(mod) results in the error:
> bootcoef<-function(data,index){
+ dat<-data[index,]
+ mod<- lme(Frames~Man2+Manip+Strings+Date.+Cut.+ID.*Man2+ID.*Manip,random=~Man2 |ID.,data=dat)
+ ranef(mod)
+ }
>
> boot.out<-boot(data,bootcoef,99)
Error: incorrect number of subscripts on matrix
Suggesting that the setup of the ranef(mod) list is different (clearly)....
Any suggestions on any of this? I have a sneaking suspicion this is not a straightforward issue.
Cheers andrew
------------------------------------------------------------------------ ---------
Dr. Andrew Beckerman
Department of Animal and Plant Sciences, University of Sheffield,
Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK
ph +44 (0)114 222 0026; fx +44 (0)114 222 0002
http://www.shef.ac.uk/beckslab
------------------------------------------------------------------------ ----------
______________________________________________ [EMAIL PROTECTED] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
