[R] crossed random effects- clarified version

[Sarah Mclean]

Hi,

here is a clarified version of my problem.

I have a total of 4*74 observations on 74 different
mums in 5 different populations of mums, subject to 6
treatments (2 moisture levels*3 substrate types).

I want to know if mum interacts with moisture level or
substrate type.

Population, moisture and substrate are fixed effects
and mum is a random effect within population. Plot is
a random whole-plot error for a split-plot design. The
data set is called fm.

This is the formula I used and the error message I
got:

fm$pmu <- getGroups(fm, ~1|pop/mum, level=2)
fm$grp = as.factor(rep(1,nrow(fm)))
fm$pl <- getGroups(fm, ~1|plot)
fm$mo <- getGroups(fm, ~1|moist)
fm$su <- getGroups(fm, ~1|sub)
> fm1 <- lme(sqrt(mass) ~ iheight + moist*sub*pop,
data=fm, random=list(grp=pdBlocked(list(pdIdent(~pl -
1), pdIdent(~pmu - 1),  pdIdent(~pmu:su - 1),
pdIdent(~pmu:mo - 1)))))
Error in chol((value + t(value))/2) : non-positive
definite matrix in chol

The model works if the interaction terms:
pdIdent(~pmu:su - 1), pdIdent(~pmu:mo - 1), are
removed, so they are causing the problem. 

The model also works if I test mums from one
population at a time so that mum no longer needs to be
nested, i.e. if I replace pmu with mu:
fm$mu <- getGroups(fm, ~1|mum)
> fm1 <- lme(sqrt(mass) ~ iheight + moist*sub,
data=fm, random=list(grp=pdBlocked(list(pdIdent(~pl -
1), pdIdent(~pmu - 1),  pdIdent(~mu:su - 1),
pdIdent(~mu:mo - 1)))))

It would be a lot faster if I can test all of the
populations at once instead of individually.

Any help would be much appreciated.

Thanks,
Sarah

http://mobile.yahoo.com.au - Yahoo! Mobile
- Check & compose your email via SMS on your Telstra or Vodafone mobile.