[R-sig-ME] Generalized randomized block design
leverkus
leverkus at ugr.es
Wed Jan 30 07:20:31 CET 2013
Thanks for clarifying that, Seth. I just tried using altitude instead
of plot as suggested, and only subplot as a random term. However,
substituting altitude by the categorical Plot seems to work way better
(now I´m using the cover of annuals as response):
>
> model<-lme(asin(sqrt(annualcov))~Treat*altitude,random=~1|Subplot,method="ML",data=na.omit(divers))
> model1<-update(model,~.-Treat:altitude)
> anova(model,model1)
Model df AIC BIC logLik Test L.Ratio p-value
model 1 8 -247.3299 -220.4396 131.665
model1 2 6 -245.9759 -225.8082 128.988 1 vs 2 5.353978 0.0688
# for the interaction using altitude
>
> modela<-lme(asin(sqrt(annualcov))~Treat*Plot,random=~1|Subplot,method="ML",data=na.omit(divers))
> model1a<-update(model,~.-Treat:Plot)
> anova(modela,model1a)
Model df AIC BIC logLik Test L.Ratio p-value
model 1 11 -255.4568 -218.4826 138.7284
model1 2 7 -244.6295 -221.1004 129.3147 1 vs 2 18.82735 8e-04
# for the interaction using plot
From plotting the data I can see there is definitely an interaction
here, so I would trust the version with Plot more than with altitude
(Residuals in both cases look centered around 0). But I have this issue
of not knowing whether it is wrong to pseudoreplicate in this way. In
any case, if my plots were not located at different altitudes there
would still have to be a way of testing this basic design, right? Could
my modela be considered correct?
Thank you,
alex
El 2013-01-29 13:52, seth at swbigelow.net escribió:
> Have you tried specifying altitude (ie., as numeric), rather than a
> plot identifier, as a main effect? This may be what Rob was
> suggesting. Since altitute is a numeric gradient, and you seem to be
> suspecting that something is going on that is related to that
> gradient, this would seem worth a try -- implemented as an
> interaction
> model, and as a non-interaction model, and then comparing with
> likelihood ratio test. I would try just specifying 'subplot' as a
> random effect, then examining the residuals as in the 'rails'l
> example
> in Pinheiro & Bates
>
> -Seth
>
>
>
>> ----- Original Message -----
>>
>> From: leverkus <leverkus at ugr.es>
>>
>> To:"Robert Kushler" <kushler at oakland.edu>,
>> <r-sig-mixed-models at r-project.org>
>>
>> Cc:
>>
>> Sent:Mon, 28 Jan 2013 21:08:32 +0100
>>
>> Subject:Re: [R-sig-ME] Generalized randomized block design
>>
>> Thanks for your reply, Rob,
>>
>> I guess you are right about not modeling plot as a random effect.
>> In
>> any case, if I formulate it this way (as I understand you suggest):
>>
>> lme(diversity~Treatment*Plot,random=~1|Plot/Subplot)
>>
>> I don´t have enough df to calculate a Plot (altitude) main effect
>> but
>> only treatment and the treatment*Plot interaction. The summary of
>> the
>> fixed effects looks like this:
>>
>> Fixed effects: diversity ~ Treatment * Plot
>> Value Std.Error [1] DF t-value p-value
>> (Intercept) 0.8332827 0.03153322 186 26.425551 0.0000
>> TreatPCL 0.0250449 0.04557570 18 0.549524 0.5894
>> TreatSL -0.1618297 0.04459471 18 -3.628898 0.0019
>> Plot2 0.1346471 0.04459471 0 3.019351 NaN # where
>> these results with 0 df look like they shouldn´t be in the model.
>> Plot3 0.0561054 0.04459471 0 1.258118 NaN
>> TreatPCL:Plot2 -0.0617449 0.06376388 18 -0968337 0.3457
>> TreatSL:Plot2 -0.0339678 0.06306644 18 -0.538603 0.5968
>> TreatPCL:Plot3 0.0217470 0.06376388 18 0.341054 07370
>> TreatSL:Plot3 0.1790523 0.06306644 18 2.839106 0.0109
>>
>> My questions here are: 1) is it ok to include a Plot main effect in
>> the
>> model (as above) even though I don´t have df for it? 2) Would it
>> be
>> "allowed" instead to use diversity~Treatment+Treatment:Plot as
>> fixed
>> effects, without a Plot main effect? Or otherwise, 3) How wrong
>> would it
>> be in the random term to place plot at the level of subplots, so
>> that
>> random=~1|Plot:Subplot? I understand in this latter way I would be
>> pseudoreplicating plot.
>>
>> I guess the main issue is that it annoys me to have a term in the
>> model
>> which tells me nothing, and not knowing which values to report for
>> altitude (the fixed effects with 0 df or the random term resulting
>> from
>> the specification of the experimental structure).
>>
>> Thanks again,
>>
>> alex
>>
>> El 2013-01-28 15:56, Robert Kushler escribió:
>> > Since Plot is confounded with "Altitude" I suggest you treat
>> Altitude
>> > as a fixed effect
>> > and give up on trying to estimate a Plot variance component (2 df
>> is
>> > not enough info
>> > for that).
>> >
>> > Regards, Rob Kushler
>> >
>> >
>> > On 1/28/2013 8:57 AM, leverkus wrote:
>> >> Dear R users,
>> >>
>> >> I am struggling with the formulation in lme of a generalized
>> >> randomized block design with subsampling, and I would very
>> >> much appreciate some help. The experiment consists of 3 plots
>> (of
>> >> ca. 20 ha each) located at different altitudes on a
>> >> mountain slope. In each plot there are 9 subplots, which are 3
>> >> replicates of 3 post-fire wood management treatments. In
>> >> each subplot we sampled 8 transects for plants (except in one
>> >> subplot, where only 5 transects were sampled), and my
>> >> response variable is species diversity. In order to take account
>> for
>> >> the experimental design and get the correct number
>> >> of denominator degrees of freedom, I am using (1|Plot/Subplot)
>> in
>> >> the random effects. Subplot is a categorical variable
>> >> which joins treatment names (treatments are "SL", "NI", "PCL")
>> and
>> >> replicates (1,2,3): SL1, SL2, SL3, NI1.. This gives
>> >> me the correct replication: 3 plots and 27 subplots. As for now,
>> my
>> >> model looks like this:
>> >>
>> >> lme(diversity~Treatment,random=~1|Plot/Subplot)
>> >>
>> >> However, treatment effects are likely to vary with altitude, so
>> I
>> >> wish to test for the treatment x plot interaction.
>> >> This is where I am stuck. By including plot as a fixed effect
>> >> (diversity~Treatment*Plot) I have no df to calculate the
>> >> plot effect and this looks weird to me. Besides, I want to have
>> plot
>> >> as a random effect. Could anyone give me some
>> >> suggestions? (I don´t mind using lmer instead.)
>> >>
>> >> Thanks in advance,
>> >>
>> >> alex
>> >>
>> >> _______________________________________________
>> >> R-sig-mixed-models at r-project.org [2] mailing list
>> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models [3]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org [4] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models [5]
>
>
> Links:
> ------
> [1] http://Std.Error
> [2] mailto:R-sig-mixed-models at r-project.org
> [3] https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> [4] mailto:R-sig-mixed-models at r-project.org
> [5] https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list