Skip to content

Extracting level-1 variance from lmer()

4 messages · Doran, Harold, David Afshartous, Dieter Menne

Original
Doran, Harold
# 
Dave:

This is an attribute. So, you can get it as follows:

library(lme4)
example(lmer) 
attr(VarCorr(fm1), "sc")
David Afshartous
# 
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
$Patient.new
1 x 1 Matrix of class "dpoMatrix"
            (Intercept)
(Intercept)    8.519916

attr(,"sc")
   scale 
3.215072
[1] 8.519916
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)
David Afshartous
# 
All,

I've fit some models via lme() and now I'm trying to fit similar models with
lmer() for some simulations I'm running.

The model below (fm1) has an intercept variance that depends on treatment
group.  How would one accomplish a similar stratification for the level-1
variance, i.e., the within-group or patient variance?

With lme() I was able to do this via:
update(fm1,
weights = varIdent(form = ~ 1 | treatment.ind) )

However, this does not work for lmer().  Is there any way around this?

Model and simulated data below.

Cheers,
David



 
fm1 = lmer (y ~  treatment.ind + ( 0 + placebo.ind | person1) + (0 +
treatment.ind | person1), data = fake)

###################

Strat.mean.var.simple <- function (J, K){
 # function to generate simple longit
     time <- rep(seq(0,1, ,length=K), J) # K measurements
     person <- rep(1:(J/2), each = K)
     treatment <- rep(0:1, each = J/2)
     treatment.ind <- rep(0:1, each = (J/2)*K)
     person1 <- rep(1:J, ,each = K)
     placebo.ind.1 <- treatment.ind < 1
     placebo.ind = ifelse( placebo.ind.1, 1, 0)
 #
     mu.a.true.P = 4.8
     mu.a.true.T = 8.8
     sigma.a.true.P = 2.2
     sigma.a.true.T = 4.2
     sigma.y.true = 1.2
 #
     a.true.P = rnorm (J/2, mu.a.true.P, sigma.a.true.P)
     a.true.T = rnorm (J/2, mu.a.true.T, sigma.a.true.T)
 #
     y.P <- rnorm( (J/2)*K, a.true.P[person], sigma.y.true)
     y.T <- rnorm( (J/2)*K, a.true.T[person], sigma.y.true)
     y <- c(y.P, y.T)
     return ( data.frame (y, time, person1, treatment.ind, placebo.ind))
 }
 
 
 J = 10
 K = 4
 set.seed(500)
 fake = Strat.mean.var.simple (J,K)
Dieter Menne
# 
David Afshartous <dafshartous <at> med.miami.edu> writes:
This is not yet implemented. See Douglas Bates'

http://article.gmane.org/gmane.comp.lang.r.lme4.devel/518

Dieter