[R] Extracting level-1 variance from lmer()
Doran, Harold
HDoran at air.org
Tue Feb 5 20:35:11 CET 2008
Dave:
This is an attribute. So, you can get it as follows:
library(lme4)
example(lmer)
attr(VarCorr(fm1), "sc")
> -----Original Message-----
> From: r-help-bounces at r-project.org
> [mailto:r-help-bounces at r-project.org] On Behalf Of David Afshartous
> Sent: Tuesday, February 05, 2008 2:37 PM
> To: r-help at r-project.org
> Subject: [R] Extracting level-1 variance from lmer()
>
>
> All,
>
> How does one extract the level-1 variance from a model fit via lmer()?
>
> In the code below the level-2 variance component may be
> obtained via subscripting, but what about the level-1
> variance, viz., the 3.215072 term?
> (actually this term squared) Didn't see anything in the
> archives on this.
>
> Cheers,
> David
>
>
>
> > fm <- lmer( dv ~ time.num*drug + (1 | Patient.new), data=dat.new )
>
> > VarCorr(fm)
> $Patient.new
> 1 x 1 Matrix of class "dpoMatrix"
> (Intercept)
> (Intercept) 8.519916
>
> attr(,"sc")
> scale
> 3.215072
>
> > VarCorr(fm)[[1]][1]
> [1] 8.519916
> > VarCorr(fm)[[2]][1]
> Error in VarCorr(fm)[[2]] : subscript out of bounds
>
>
>
>
>
>
>
> ##########################################################
> set.seed(500)
> n.timepoints <- 4
> n.subj.per.tx <- 20
> sd.d <- 5;
> sd.p <- 2;
> sd.res <- 1.3
> drug <- factor(rep(c("D", "P"), each = n.timepoints, times =
> n.subj.per.tx))
> drug.baseline <- rep( c(0,5), each=n.timepoints,
> times=n.subj.per.tx ) #Patient <- rep(1:(n.subj.per.tx*2),
> each = n.timepoints) Patient.baseline <- rep( rnorm(
> n.subj.per.tx*2, sd=c(sd.d, sd.p) ), each=n.timepoints )
> time.baseline <-
> rep(1:n.timepoints,n.subj.per.tx*2)*as.numeric(drug=="D")
> dv <- rnorm( n.subj.per.tx*n.timepoints*2,
> mean=time.baseline+Patient.baseline+drug.baseline, sd=sd.res )
>
> dat.new <- data.frame(drug, dv)
> dat.new$Dind <- as.numeric(dat.new$drug == "D") dat.new$Pind
> <- as.numeric(dat.new$drug == "P") dat.new$time.num =
> rep(1:n.timepoints, n.subj.per.tx*2) dat.new$Patient.new =
> rep(1:20, each=8)
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
More information about the R-help
mailing list