

	I hope I'm not beating a dead horse here, but in your original post
you  said 'treatment effects are likely to vary with altitude', which
implies that you'd like to make some kind of formal inference about
altitude. If you use the categorical variable 'Plot' rather than the
numerical 'altitude' as a fixed effect, it's not surprising that you
get smaller AIC, since the fit to each plot is unconstrained by an
external variable -- but then you don't get to make inferences on
altitude. (as an aside, if you do use altitude, it might be better to
transform or center it in some way, as per the example in
Pinheiro/Bates Chap 1.4)

	As I recall from my long-ago stats class, an interaction effect of
p=0.06 is still considered significant and gives you license to test
treatment effects at individual levels of plot/altitude, if that's
what you want to do. 

	I'm not following/understanding the issue about pseudoreplication

	--Seth

----- Original Message -----
From: leverkus 
To:
Cc:
Sent:Wed, 30 Jan 2013 07:20:31 +0100
Subject:Re: [R-sig-ME] Generalized randomized block design

 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 model1 anova(model,model1)
 Model df AIC BIC logLik Test L.Ratio [2] 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 model1a anova(modela,model1a)
 Model df AIC BIC logLik Test L.Ratio [4] 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@swbigelow.net [5] 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 
 >>
 >> To:"Robert Kushler" ,
 >> 
 >>
 >> 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 [6] [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@r-project.org [7] [2] mailing list
 >> >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models [8]
[3]
 >>
 >> _______________________________________________
 >> R-sig-mixed-models@r-project.org [9] [4] mailing list
 >> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models [10] [5]
 >
 >
 > Links:
 > ------
 > [1] http://Std.Error [11]
 > [2] mailto:R-sig-mixed-models@r-project.org [12]
 > [3] https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models [13]
 > [4] mailto:R-sig-mixed-models@r-project.org [14]
 > [5] https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models [15]



Links:
------
[1] http://sitemail.hostway.com/http:
[2] http://L.Ratio
[3] http://sitemail.hostway.com/http:
[4] http://L.Ratio
[5] mailto:seth@swbigelow.net
[6] http://Std.Error
[7] mailto:R-sig-mixed-models@r-project.org
[8] https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[9] mailto:R-sig-mixed-models@r-project.org
[10] https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[11] http://Std.Error
[12] mailto:R-sig-mixed-models@r-project.org
[13] https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[14] mailto:R-sig-mixed-models@r-project.org
[15] https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


	[[alternative HTML version deleted]]

