[R-sig-eco] Can I weight LS means estimates by elevation in lme4/lmerTest packages?

[Philippi, Tom]

Matthew--
Do you get what you want in terms of estimated effect of ecosystem type by
dropping Elevation from your model?  That difference would include the
component of the response to the difference in mean elevation.

If not, from a pure parameter estimation perspective, if your simulated
data are representative, you have elevation ranging from 1200:1700 (m) in
grassland ecosystems, but 1800:3000 in shrublands.  Therefore, you have 2
df for elevation: a linear effect (Elevation) plus a step effect
(Ecosystem) split at 1750m, so forced parallel lines that don't overlap.
Your estimated effect for ecosystem is the offset of those parallel lines.
However, any nonlinearities in grass cover as a function of elevation will
be forced into that step (ecosystem) parameter, which could be misleading
or masking.

Is a linear effect of elevation reasonable for your range of elevations, or
are there likely to be curved or unimodal relationships where too low gets
too dry and too high gets too cold for grasses?  Before I did any model
fitting, I would plot 6 gam or loess lines of grass cover by elevation for
the 6 ecosystem * quality combinations.

I hope that this gives you food for thought, or alternatively that I'm way
off base and someone provides you with exactly the narrow answer you
requested.

Tom 2


On Fri, Nov 28, 2014 at 5:35 PM, Matthew Van Scoyoc <scoyoc at gmail.com>
wrote:

> Good afternoon r-sig-ecology,
>
> I'm running linear mixed models using the lme4 and lmerTest packages to
> examine ecosystem structure in grasslands and shrublands. The grasslands
> are located at lower elevations than the shrublands, and I would like to
> weight the estimates from LS means to reflect the differences in
> elevations. My colleague says there is an easy way to do it in SAS, but I
> haven't found a way to do it in R.
>
> There is an example dataset and workflow. Here I would be examining grass
> cover differences between the two ecosystems, the quality of the
> ecosystems, and examining the interaction between ecosystem and quality. I
> am also interested in systematic changes throughout the study area (not
> represented in this example) so I don't want to run separate analyses on
> for each ecosystem. I just want to adjust the LS means estimates to reflect
> the differences in elevation and not a mean elevation of the sampled plots.
>
> >library("lme4")
> >library("lmerTest")
> >
> > df = data.frame(PlotID = rep(c(paste0("G", 1:30), paste0("S", 1:30)), 2),
> +                 SamplePeriod = as.factor(c(rep(2012, 30), rep(2014,
> 30))),
> +                 Ecosystem = rep(c(rep("Grassland", 30), rep("Shrubland",
> 30)), 2),
> +                 Quality = rep(rep(c(rep("Good", 10), rep("Moderate", 10),
> rep("Poor", 10)), 2), 2),
> +                 GrassCover = c(runif(10, min = 0.50, max = 0.85), # 2012
> Grassland Good
> +                                runif(10, min = 0.50, max = 0.60), # 2012
> Grassland Moderate
> +                                runif(10, min = 0.30, max = 0.40), # 2012
> Grassland Poor
> +                                runif(10, min = 0.25, max = 0.60), # 2012
> Shrubland Good
> +                                runif(10, min = 0.20, max = 0.45), # 2012
> Shrubland Moderate
> +                                runif(10, min = 0.05, max = 0.25), # 2012
> Shrubland Poor
> +                                runif(10, min = 0.50, max = 0.90), # 2014
> Grassland Good
> +                                runif(10, min = 0.50, max = 0.55), # 2014
> Grassland Moderate
> +                                runif(10, min = 0.30, max = 0.30), # 2014
> Grassland Poor
> +                                runif(10, min = 0.25, max = 0.60), # 2014
> Shrubland Good
> +                                runif(10, min = 0.20, max = 0.30), # 2014
> Shrubland Moderate
> +                                runif(10, min = 0.05, max = 0.15))) # 2014
> Shrubland Poor
> > Elevation = c(sample(1200:1700, size = 30, replace = T),
> sample(1800:3000, size = 30,
> >                    replace = T))
> > df$Elevation = c(Elevation, Elevation); rm(Elevation)
> >
> > lmm = lmer(GrassCover ~ Ecosystem*Quality + Elevation + (1|PlotID), data
> = df, REML = T)
> > anova(lmm,  ddf = "Satterthwaite", type = 3, method.grad = "Richardson")
> Analysis of Variance Table of type 3  with  Satterthwaite
> approximation for degrees of freedom
>                          Sum Sq      Mean Sq      NumDF      DenDF
>  F.value    Pr(>F)
> Ecosystem           1.68172     1.68172     1              113
>  51.239    8.79e-11 ***
> Quality               1.97119     0.98560     2               113
> 152.611   < 2.2e-16 ***
> Elevation             0.00311     0.00311     1              113
>  0.016      0.89978
> Ecosystem:Quality 0.04201    0.02101     2              113          3.604
>      0.03039 *
> ---
> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> > lsms = lsmeans(lmm, method.grad = "Richardson")
> >estimates = lsms[[1]]$Estimate
>
> At this point I would be plotting the estimates of the significant response
> variables and interactions to look at the differences. As I understand, the
> ecosystem effect and the interaction between ecosystem and quality are with
> an average elevation for all plots, and this could exaggerate LS means
> estimate.  Correcting for elevation would provide more accurate
> estimations. Right?
>
> Okay, thanks for you help.
> Cheers,
> MVS
> =====
> Matthew Van Scoyoc
>
> <
> https://mail.google.com/mail/?view=cm&fs=1&tf=1&to=mvanscoyoc@aggiemail.usu.edu
> >
> https://sites.google.com/site/scoyoc/
> =====
> Think SNOW!
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>

	[[alternative HTML version deleted]]