[R] Nested Design Coding Question
Douglas Bates
bates at stat.wisc.edu
Wed Feb 19 19:43:29 CET 2003
"Steeno, Gregory S" <gregory_s_steeno at groton.pfizer.com> writes:
> I'm a SAS user who is slowly but surely migrating over to R. I'm trying to
> find the proper code to analyze a nested design. I have four
> classification variables, L (fixed), A (random within L), D (random within
> L), and I (random within L). The model I'm interested in is
>
> L A(L) D(L) I(L) A:D:I(L),
>
> where the interaction is interpreted as the lack-of-fit term. I've tried
> variants of the lme function similar to these,
>
> lme(response~L, data, random=~Lab/(A+L+I+A:D:I),
> lme(response~1, data, random=~Lab/(A+L+I+A:D:I),
> lme(response~L, data, random=~1/(A+L+I+A:D:I).
>
> All give results different from SAS, and all give warning messages regarding
> either false- or non-convergence.
>
> For reference, the abbreviated SAS code is,
>
> model response = L;
> random A(L) D(L) I(L) A:D:I(L);
>
> Can anyone shed some light? I'd be very appreciative.
You can use lme to fit models with crossed random effects like this
but not easily. The algorithms for lme are tuned for nested random
effects. If you have A random within L, you must first establish a
unique set of levels for the factor
data$aL = getGroups(data, ~ 1 | L/A, level = 2)
data$dL = getGroups(data, ~ 1 | L/D, level = 2)
data$iL = getGroups(data, ~ 1 | L/I, level = 2)
Then you must establish an artificial grouping factor with only one level
data$grp = as.factor(rep(1, nrow(data)))
Once you have the artifical grouping factor you create a blocked
variance-covariance matrix whose blocks are multiples of the identity
applied to indicator variables. In R the indicators for a factor are
generated from a formula like ~ aL - 1. (By the way, D and I are poor
choices for variable names - see Venables and Ripley (Springer, 2002,
p. 13).)
Putting it all together provides the highly unintuitive function call
fm = lme(response ~ L, data, random = list(
grp = pdBlocked(list(pdIdent(~ aL - 1),
pdIdent(~ dL - 1),
pdIdent(~ iL - 1),
pdIdent(~ aL:dL:iL - 1))))
--
Douglas Bates bates at stat.wisc.edu
Statistics Department 608/262-2598
University of Wisconsin - Madison http://www.stat.wisc.edu/~bates/
More information about the R-help
mailing list