[cc'ing back to r-sig-mixed-models, with apologies]
On 04/15/2011 02:38 PM, Kevin Wright wrote:
Ben, Even Doug has "Variety" as a fixed effect and on the right side of the bar
print(Om1 <- lmer(yield ~ nitro + Variety + (1 |
Block/Variety), + Oats), corr = FALSE) Source: cran.r-project.org/web/packages/ <http://cran.r-project.org/web/packages/>*lme4*/vignettes/Implementation.pdf I'm confused by your comment. Kevin
Interesting. I'm a bit confused myself. I guess these are two
different models (I think):
* if we use (1|Block/Variety), we are estimating a random intercept
across blocks and a random intercept across varieties within blocks;
all random effect values are independent of each other (e.g. eps_B(i) vs
eps_B(j) for i \neq j, eps_V(i) vs eps_V(j) for i \neq j, eps_B(i) vs
eps_V(j) for *all* i,j). Parameters are sigma^2_B and sigma^2_V.
* if we use (Variety|Block), we are estimating variance of the
intercept across blocks, estimates of the variance of two contrasts
('Marvellous' vs 'Golden Rain' and 'Victory' vs 'Golden Rain') across
blocks, as well as the correlations among the contrasts -- parameters
are sigma^2(GR,M vs GR,V vs GR), rho_{all 3 pairwise contrasts}.
The first way is certainly the classic way to do it (it's also the way
it's presented in MASS). It's also more parsimonious, although in
principle the information you get out from the new way could be
interesting (allows e.g. for non-sphericity -- different correlations
among levels within blocks). I guess the second way reduces to the
first when there is a single value of sigma^2 and a single value of rho ...
cheers
Ben Bolker
On Fri, Apr 15, 2011 at 8:45 AM, Ben Bolker <bbolker at gmail.com
<mailto:bbolker at gmail.com>> wrote:
Kevin Wright <kw.stat <at> gmail.com <http://gmail.com>> writes:
>
> > I am trying to teach myself R and replicate some previous SAS
analysis.
> > Could someone please help me translate the following SAS code
into R.
> >
> > Proc mixed method=ml
> > Class Group Treatment Stream Time Year;
> > Model Logrpk=Treatment Time Treatment*Time;
> > Random Group Stream (Group Treatment) Year(Time);
> >
>
> Assuming you have a data frame "dat" with these factors: Group
Treatment
> Stream Time Year
> And continuous response: logrpk
>
> This code is a starting point: (I'm not sure exactly what the SAS
syntax
> means).
>
> require(lme4)
> m1 = lmer(logrpk ~ treatment*time + (1|Group) +
(1|Stream:Group:Treatment) +
> (1|Year:Time), data=dat)
>
Can I please suggest that (Treatment|Stream:Group) or something
like it is more appropriate than (1|Stream:Group:Treatment)? In
general, what goes on the LEFT of the bar is an intercept or fixed
effect (i.e. something that varies between groups); what goes on
the RIGHT of the bar is a grouping variable. Thus if a fixed effect
terms ends up on the right of the bar, something funny is going on.
Ben Bolker
______________________________________________
R-help at r-project.org <mailto: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.