[R-sig-ME] Fwd: help with a cross-classified random effects model code in R.

[Douglas Bates]

---------- Forwarded message ----------
From: Douglas Bates <bates at stat.wisc.edu>
Date: Mon, May 12, 2008 at 12:38 PM
Subject: Re: help with a cross-classified random effects model code in R.
To: "Violet(Shu) Xu" <shuxu at ucdavis.edu>


On Fri, May 9, 2008 at 5:29 PM, Violet(Shu) Xu <shuxu at ucdavis.edu> wrote:
> Dear Dr. Bates,

>  I have a question about fitting a cross-classified random effects model in
> R.

There are now two packages in R for fitting linear mixed effects
models using maximum likelihood or REML.  The more recent of these
packages, called lme4, is the better choice for fitting models with
crossed or partially crossed factors for the random effects.  I would
recommend using that package rather than trying to cobble together a
solution using lme.


> The model can be expressed in SAS code as follows (from
> http://www.ats.ucla.edu/stat/sas/examples/mlm_ma_hox/chapter7.htm):

Thank you for the URL to the SAS code and results.  I enclose analyses
that I believe reproduce the SAS results using the development version
of the lme4 package, which you can install with

install.packages("lme4", repos = "http://R-forge.R-project.org")

This script uses the pupcross.dta file from the Stata form of the data
files, available as

http://www.ats.ucla.edu/stat/stata/examples/mlm_ma_hox/ma_hox_stata.ZIP


> proc mixed data =pupcross covtest noclprint method=ml;
>   class pupil pschool sschool;
>   model achiev =   / solution ddfm =satterth;
>   random intercept / subject=sschool;
>   random intercept / subject=pschool;
> run;
>
> I tried to fit it in R, however I did not get the same estimates,
> especially the estimate of the random effects .  Would you help me to
> figure out what is the problem in my code?
>
> My R(version 2.5.1) code as follows
> pupcross = read.csv("pupcross.csv", header=T, na.strings = '.')
> head(pupcross)
> pupcross$cons = 1
>
> pupcrossg = groupedData(achiev ~ 1 |cons,  data = pupcross)
>
> pupcrossg.m1 = lme(achiev~1, random =
> pdBlocked(list(pdIdent(~pschool-1),pdIdent(~sschool-1))),
>                  data = pupcrossg)
>
> Thank you for your time and consideration. Have a nice weekend .
>
> --
> Shu (Violet) Xu
> --*-*--*-*--*-*--*-*--*-*--*-*--
> Graduate student
> Quantitative Psychology program
> Department of Psychology
> University of California, Davis