[R-sig-ME] random effect variance per treatment group in lmer
Alan Cobo-Lewis
alanc at umit.maine.edu
Thu Jul 12 00:39:44 CEST 2007
Dave,
How about using stratifying variance on level of drug using
( 0 + as.numeric(drug=="D") | Patient ) + ( 0 + as.numeric(drug=="P") | Patient )
Here's some code (whose sim also builds in a fixed effect of time that applies only to the Drug condition).
set.seed(500)
n.timepoints <- 8
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 <- factor(paste("Time-", rep(1:n.timepoints, n.subj.per.tx*2), sep=""))
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(time, drug, dv, Patient)
xyplot( dv~time|drug, group=Patient, type="l", data=dat.new )
# fit model treats time as a quantitative predictor
( fm.het <- lmer( dv ~ rep(1:n.timepoints, n.subj.per.tx*2)*drug + ( 0 + as.numeric(drug=="D") | Patient ) + ( 0 + as.numeric(drug=="P") | Patient ), data=dat.new ) )
HTH,
alan
>From: "Afshartous, David" <afshart at exchange.sba.miami.edu>
asked:
>
>All,
>I didn't receive a response to the query below sent to the general
>R-help mailing list so figured I'd try this mailing list. Apologies
>in advance if this is an overly simplistic question for this list; I
>just started
>w/ lmer after not using lme for awhile.
>Cheers,
>Dave
>
>
>__________________________________________________________
>
>All,
>
>How does one specify a model in lmer such that say the random effect for
>
>the intercept has a different variance per treatment group?
>Thus, in the model equation, we'd have say b_ij represent the random
>effect
>for patient j in treatment group i, with variance depending on i, i.e,
>var(b_ij) = tau_i.
>
>Didn't see this in the docs or Pinherio & Bates (section 5.2 is specific
>for
>modelling within group errors). Sample repeated measures code below is
>for
>a single random effect variance, where the random effect corresponds to
>patient.
>cheers,
>dave
>
>
>z <- rnorm(24, mean=0, sd=1)
>time <- factor(paste("Time-", rep(1:6, 4), sep=""))
>Patient <- rep(1:4, each = 6)
>drug <- factor(rep(c("D", "P"), each = 6, times = 2)) ## P = placebo, D
>= Drug
>dat.new <- data.frame(time, drug, z, Patient)
>fm = lmer(z ~ drug + time + (1 | Patient), data = dat.new )
--
Alan B. Cobo-Lewis, Ph.D. (207) 581-3840 tel
Department of Psychology (207) 581-6128 fax
University of Maine
Orono, ME 04469-5742 alanc at maine.edu
http://www.umaine.edu/visualperception
More information about the R-sig-mixed-models
mailing list