I am trying to use lmer or lme to fit a mixed-effects model. My data looks
something like this:

y: response
factor: two levels of a treatment...coded as "Non-Restored"/"Restored"
site: nested within factor (specifically I have 5 sites which are controls
and 5 which are treatments). This is simply coded as 1...10
x2: continuous centered covariate, the same for all measurements within a
site
measurements: nested within sites (ranging from 8-12 per site)

I would like to estimate the effect of the treatment (a fixed effect), the
covariate (another fixed effect). I would also like to estimate the random
effect for sites, but would like this to possibly be different for each
level of the factor. Finally I would like the error term to also possibly
be different across each level of the factor.

What I have tried:

lme(fixed = y ~ factor + x2, random =~ (1|site),weights =
varIdent(form=~1|factor))
lme(fixed = y ~ factor + x2, random =~ (1 + factor|site),weights =
varIdent(form=~1|factor))
##Both of these seem wrong since they don't account for the nested aspect,
and the second estimates random effects by all ten sites for both the
intercept (i.e. "Non-Restored") and the factor (i.e. Restored) making them
highly correlated,

So, I then found here: http://glmm.wikidot.com/faq#modelspec, that I should
represent the model as:

lme(fixed= y~factor + x2,random=~(1|factor:site), weights =
varIdent(form=~1|factor))

However, I got this error:

Error in getGroups.data.frame(dataMix, groups) :
  invalid formula for groups

So then I tried using lmer (thinking that the : specification was the
issue that lme didn't like). So I tried:

lmer(y~factor + x2+ (1|factor:site), weights = varIdent(form=~1|factor))


but got this error:

Error in model.frame.default(data = data, weights = varIdent(form = ~1 |  :
  variable lengths differ (found for '(weights)')


which I assumed because lmer doesn't support varClasses (is that true?)

So then i just tried this (even though I was pretty sure I will get the
same error term across levels of my factor):

lmer(y~factor + x2+ (1|factor:site))

and I got this error:

Error in validObject(.Object) :
  invalid class “mer” object: Slot Zt must by dims['q']  by dims['n']*dims['s']
In addition: Warning messages:
1: In Restoration:Site.Number :
  numerical expression has 140 elements: only the first used
2: In Restoration:Site.Number :
  numerical expression has 140 elements: only the first used

So, now I'm here asking for help because I'm not sure how to do the
things I want. (and this seems like a fairly straightforward
issue...but I'm stymied as to how to code it)

By the way I fit these models on the two subsets of the data defined
by the factor (i.e. Restored vs. Not)

lme(fixed = y ~ 1 + x2, random=~1|site, subset=Restoration=="Non-Restored")

lme(fixed = y ~ 1 + x2, random=~1|site, subset=Restoration=="Restored")

and the estimates made sense, but clearly I should be able to include
this in one model especially for estimating the fixed effect of the
covariate as well as the effect of restoration.

Thanks for advice,

Kris

	[[alternative HTML version deleted]]

