Thanks Thierry,

I wasn't aware that lme understood the nesting structure implicit in the
data. However, I'm still trying to understand what having a random effect
on the intercept means for sites that are Restored (i.e., factor = 1). In
other words, what if I think the variance across sites (i.e. the variance
in the random effect) is different for each level of my fixed factor? So,
for example as I mentioned in my post, I ran:

lme(fixed = y ~ 1 + x2, random=~1|site,subset=factor=="Restored")
I get an estimate for the random effect sd as: 5.47 (I interpret this as
the variation across sites in mean response for Restored sites)

lme(fixed = y ~ 1 + x2, random=~1|site,subset=factor=="Non-Restored")
I get an estimate for the random effect sd as: 27.4 (I interpret this as
the variation across sites in mean response for Non-Restored sites)

So when I run the model I first ran (i.e. the one you said is correct), I
get a random effect sd on the intercept of 7.7. So to me this means the
random effect variance would be the same for both groups since the only
difference between Restored and Non-Restored is the fixed effect constant
estimated in the model. Yet, looking at the results from these two models
above, the random effect estimate seems to be quite different.

Thanks in advance for professional and courteous dialogue through this user
group.

Kris



On Tue, Jul 30, 2013 at 8:22 AM, ONKELINX, Thierry <Thierry.ONKELINX@inbo.be
> wrote:

> 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@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@r-project.org [mailto:
> r-sig-mixed-models-bounces@r-project.org] Namens Kristofor Voss
> Verzonden: dinsdag 30 juli 2013 5:11
> Aan: r-sig-mixed-models@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.
>



-- 
Kristofor A. Voss
Ph.D. Candidate
Bernhardt Lab
Nicholas School of the Environment
Duke University
P.O. Box 90328
Durham, NC 27708
619-249-3614
kristofor.voss@duke.edu
Office: FFSC 3304

	[[alternative HTML version deleted]]

