"Steeno, Gregory S" <[EMAIL PROTECTED]> 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 [EMAIL PROTECTED]
Statistics Department 608/262-2598
University of Wisconsin - Madison http://www.stat.wisc.edu/~bates/
______________________________________________
[EMAIL PROTECTED] mailing list
http://www.stat.math.ethz.ch/mailman/listinfo/r-help
