R help - Feb 2004 - nlme_crossed AND nested random effects

Dear R-help group,

How can I define a lme with 3 factors(a,b,c), where c is nested in b, 
and a is crossed with b/c?
I think that:

lme(response ~ ..., data = Data,
      random = pdBlocked(list(pdIdent(~ a - 1), pdIdent(~ b - 1))))

is one part of the answer and:

lme(response~...,  data=Data, random=~1|b/c)

is the other part of the answer but how can I combine them??

Could anybody please help me?

Many greetings,
Thomas

PS: Is there good literatur available which could help working with 
mixed random effects (I just know Pinheiro and Bates 2000), perhaps 
even on the internet?
	[[alternative HTML version deleted]]

Fabbro Thomas <thomas.fabbro at unifr.ch> writes:
> How can I define a lme with 3 factors(a,b,c), where c is nested in b, 
> and a is crossed with b/c?
> I think that:
> 
> lme(response ~ ..., data = Data,
>       random = pdBlocked(list(pdIdent(~ a - 1), pdIdent(~ b - 1))))
> 
> is one part of the answer and:
> 
> lme(response~...,  data=Data, random=~1|b/c)
> 
> is the other part of the answer but how can I combine them??
It's messy and inefficient with the current lme syntax, which was
specifically designed for nested random effects and does not handle
crossed random effects elegantly.

(I am working on other methods for lme that will make this much
easier.  however, it will be some time before I am able to release
even experimental versions of that code.  The changes involve tearing
apart all the lme code and rebuilding the underlying data structures
from scratch.  This is the fourth time that I have done that and I
keep thinking of a line from an old blues song about "keep doing it
wrong till you do it right".)

You need to generate a factor that corresponds to b/c.  Depending upon
how c is coded you may need to create it as

Data$bc = factor(paste(Data$b, Data$c, sep = "/"))

You will also need a factor with only one level

Data$const = factor(rep(1, nrow(Data)))

and you just generate the random effects model matrix from the various
sections

lme(response ~ ..., data = Data,
  random = list(const = pdBlocked(pdIdent(~a - 1), pdIdent(~ b - 1), 
    pdIdent(~ bc - 1)))

As I said, this is ugly and inefficient but it isn't the type of model
for which the current lme was designed.