R help - Feb 2003 - Nested Design Coding Question

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.

Thanks.

Greg Steeno









LEGAL NOTICE\ Unless expressly stated otherwise, this message is ... [[dropped]]

"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/