-----Original Message-----
From: Afshartous, David
Sent: Friday, July 13, 2007 4:32 PM
To: Alan Cobo-Lewis
Cc: r-sig-mixed-models at r-project.org; Afshartous, David;
Andrew Robinson
Subject: RE: random effect variance per treatment group in lmer
Alan,
Thanks again for the help. I defined time.num, Dind, and
Pind outside the lmer call as per Doug Bates' suggestion and
now it works and I get the same results as your output w/ the
suggested lmer call:
( fm.het <- lmer( dv ~ time.num*drug + ( 0 + Dind | Patient )
+ ( 0 + Pind | Patient ), data=dat.new ) )
More importantly, this achieves the goal I was seeking
originally. Looks like the key to the issue is not defining
a random effect via "(1 | Patient)"
as once this is done it doesn't seem possible to stratify;
but rather splitting the "1" into the two pieces as per the
drug dichotomy, and getting separate random effects variances
for each while mainting drug as a fixed effect.
Thanks again to you and everyone else for the help!
Cheers,
Dave
-----Original Message-----
From: Alan Cobo-Lewis [mailto:alanc at umit.maine.edu]
Sent: Friday, July 13, 2007 12:27 PM
To: Afshartous, David
Cc: r-sig-mixed-models at r-project.org; "
"Afshartous at basalt.its.maine.edu; Andrew Robinson
Subject: Re: random effect variance per treatment group in lmer
Hmm, could it be a word-wrap issue? I just verified that the code
works on R 2.5 with lme4 and lattice packages installed.
See http://www.umaine.edu/visualperception/lme4het
The model I posted assumes that the subj-to-subj variability in
baseline score has higher variability in the "D" condition
"P" condition. You are correct that there's no random
I don't think you *want* there to be a random effect for
drug doesn't have levels sampled from some population (it's
fixed effect).
Instead, the code I posted has the random effect of Patient with a
higher variability in one level of drug than in the other level of
drug.
If you're looking some the drug to have a different time course for
some patients than for other patients I think that'd be a random
effect of Time. (To avoid overparameterization I think you'd almost
surely want to treat time as a quantitative predictor if
to model it as having a random effect.)
One of the reasons that I constructed a full example with a new sim
was to get a big-enough data set together so that the lmer
obviously matched the parameters of the sim (to verify the
of the sim and the suggested model formula). Another reason is that
the code of your sim didn't actually have different subj-to-subj
variability in the "D"
condition than it did in the "P" condition, so a model with such an
effect was overparameterized for your original simulated
this be why the lmer call I suggested failed for you?
BTW, nice to see someone from the UM School of Business (I didn't
notice until just now what your email address said). I worked there
about 20 years ago on the Applied AI Reporter when I was an
undergraduate psychology major. But, looking at the SBA web site, I
recognize barely anyone, and those I do recognize probably wouldn't
recognize me.
alan
"Afshartous, David" <afshart at exchange.sba.miami.edu> on
13, 2007 at 10:44 AM -0500 wrote:
Alan,
Thanks for the suggestion. I noticed the following error
( fm.het <- lmer( dv ~ rep(1:n.timepoints,
n.subj.per.tx*2)*drug + (
+ as.numeric(drug=="D") | Patient ) + ( 0 + as.numeric(drug=="P") |
Patient ), data=dat.new ) )
Error: length(term) == 3 is not TRUE
I tried a few changes to the model but the error still
keep checking. I assume the rationale for the structure of
call, where you use as.numeric as opposed to just drug
insure that you do not introduce a random effect for drug into the
model?
Regards,
Dave
-----Original Message-----
From: Alan Cobo-Lewis [mailto:alanc at umit.maine.edu]
Sent: Wednesday, July 11, 2007 6:40 PM
To: r-sig-mixed-models at r-project.org
Cc: " "Afshartous at basalt.its.maine.edu; Afshartous,
Robinson
Subject: Re: random effect variance per treatment group in lmer
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
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 <-
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 ) #
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>
All,
I didn't receive a response to the query below sent to
R-help mailing list so figured I'd try this mailing list.
advance if this is an overly simplistic question for this
started w/ lmer after not using lme for awhile.
Cheers,
Dave