[R-sig-ME] Nested random effect within fixed effect factor?

[ONKELINX, Thierry]

Dear Kristofor,

Since 'site' is nested in 'factor' it is nonsense to have a random effect of 1 + factor|site or 1|factor:site. The first one cannot be estimate since you have only one factor per site. The second one is (in the case of your dataset) equivalent to 1|site.

lmer() does not handle correlation structure nor variance structures.

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

is the correct specification of your model and takes the implicit nesting into account.

Best regards,

Thierry


ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx op inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey


-----Oorspronkelijk bericht-----
Van: r-sig-mixed-models-bounces op r-project.org [mailto:r-sig-mixed-models-bounces op r-project.org] Namens Kristofor Voss
Verzonden: dinsdag 30 juli 2013 5:11
Aan: r-sig-mixed-models op r-project.org
Onderwerp: [R-sig-ME] Nested random effect within fixed effect factor?

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]]

* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.