[R-sig-ME] lmer adjusted means

[Allan Carson]

Hello 

 

I am analysing plant biomass data from revegetation testplots which I constructed in a split-strip-split plot design using the function lmer. For the experiment, the three factors are amelioration (2 levels), fertilizer (2 levels) and treatment (7 levels). Each testplot (block) has a single replicate of each treatment (total of 8 testplots). The blocks were constructed of just a single source of topsoil. Each block was divided into 2 subplots in which half was ameliorated the other half was not (amelioration - yes or no). Each subplot was divided again into two subsubplots (fertilizer - yes or no). Then within each subsubplot, the area was divided into 7 seeding treatments (treatments). My response variable is abovegound biomass m/^2 (grams). I have transformed my data to meet model assumptions, simplified my model using stepwise elimination and used a likelihood ratio test to determine whether the lowest level of nesting should be removed.


Here is my minimal model: 

fm6<-lmer(sqrt(abovegroundbiomass.m.2)~amelioration+fertilizer+treatment +(1|block),na.action = na.omit)

 

 

Here is my anova:


anova(fm6)

 

Analysis of Variance Table
             Df Sum Sq Mean Sq F value
amelioration  1  28.88   28.88  2.9129
fertilizer    1 521.43  521.43 52.5936
treatment     6 279.41   46.57  4.6970


 

attributes(VarCorr(fm6))$sc
(fm6fstat<- anova(fm6)[, 3]/ 3.148707^2)
nrow(tvd1new) - sum(anova(fm6)[, 1]) - 1 # d.f=90
round(pf(fm6fstat, anova(fm6)[, 1], 90, lower.tail = FALSE), 4)

[1] 0.0913 0.0000 0.0003

 

Here is my summary:

 

tval <- attributes(summary(fm6))$coef[, 3]
pval <- 2 * pnorm(abs(tval), lower = FALSE)

cbind(attributes(summary(fm6))$coef, `p value` = round(pval,
+ + 4))
                           Estimate Std. Error   t value p value
(Intercept)               1.6196612  0.9165348 1.7671573  0.0772
ameliorationyes           1.0155789  0.5950497 1.7067129  0.0879
fertilizeryes             4.3153841  0.5950497 7.2521411  0.0000
treatmentelymgla ecovar   5.0125584  1.1132360 4.5026917  0.0000
treatmentelymgla local    2.7699072  1.1132360 2.4881582  0.0128
treatmentfestsax ecovar   0.6953266  1.1132360 0.6245995  0.5322
treatmentlupiarc local    1.1391840  1.1132360 1.0233087  0.3062
treatmentnativemix ecovar 2.1091683  1.1132360 1.8946282  0.0581
treatmentnativemix local  3.2103166  1.1132360 2.8837701  0.0039

 

I have used a tukeys pairwise comparisons to look for significant differences between the levels of treatment:

 

summary(fm6c, test=adjusted("none"))

         

Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lmer(formula = sqrt(abovegroundbiomass.m.2) ~ amelioration + 
    fertilizer + treatment + (1 | block), na.action = na.omit)

Linear Hypotheses:
                                        Estimate Std. Error z value Pr(>|z|)    
elymgla ecovar - bromcil ecovar == 0      5.0126     1.1132   4.503 6.71e-06 ***
elymgla local - bromcil ecovar == 0       2.7699     1.1132   2.488 0.012841 *  
festsax ecovar - bromcil ecovar == 0      0.6953     1.1132   0.625 0.532234    
lupiarc local - bromcil ecovar == 0       1.1392     1.1132   1.023 0.306162    
nativemix ecovar - bromcil ecovar == 0    2.1092     1.1132   1.895 0.058142 .  
nativemix local - bromcil ecovar == 0     3.2103     1.1132   2.884 0.003929 ** 
elymgla local - elymgla ecovar == 0      -2.2427     1.1132  -2.015 0.043954 *  
festsax ecovar - elymgla ecovar == 0     -4.3172     1.1132  -3.878 0.000105 ***
lupiarc local - elymgla ecovar == 0      -3.8734     1.1132  -3.479 0.000503 ***
nativemix ecovar - elymgla ecovar == 0   -2.9034     1.1132  -2.608 0.009106 ** 
nativemix local - elymgla ecovar == 0    -1.8022     1.1132  -1.619 0.105464    
festsax ecovar - elymgla local == 0      -2.0746     1.1132  -1.864 0.062384 .  
lupiarc local - elymgla local == 0       -1.6307     1.1132  -1.465 0.142962    
nativemix ecovar - elymgla local == 0    -0.6607     1.1132  -0.594 0.552827    
nativemix local - elymgla local == 0      0.4404     1.1132   0.396 0.692391    
lupiarc local - festsax ecovar == 0       0.4439     1.1132   0.399 0.690107    
nativemix ecovar - festsax ecovar == 0    1.4138     1.1132   1.270 0.204074    
nativemix local - festsax ecovar == 0     2.5150     1.1132   2.259 0.023873 *  
nativemix ecovar - lupiarc local == 0     0.9700     1.1132   0.871 0.383580    
nativemix local - lupiarc local == 0      2.0711     1.1132   1.860 0.062820 .  
nativemix local - nativemix ecovar == 0   1.1011     1.1132   0.989 0.322594    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
(Adjusted p values reported -- none method)
 

 

I would like to generate the adjusted means for reporting but I cannot seem to figure out how to do this. I know for aov you can use model.tables() but I do not know what the command is within lmer.If anyone could help me with this, it would be greatly appreciated. Thanks in advance.

 

Allan