[R-sig-ME] [R] Conveting SAS Proc mixed to R code

[Ben Bolker]

  [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.
> 
>