I am running a mixed effects model using the following code: m1<- lmer(DV~TMT*TMT2+(1|Block/TMT1),verbose=T) An example dataset is here: http://pastebin.com/bHug5kTt The two explanatory variables are both factors with two levels (treatment and control). Treatment 2 is split within Treatment 1 which is in turn within block. The model output I get is as follows (sorry for the formatting!) Linear mixed model fit by REML Formula: Richness ~ NT * WT + (1 | Block/WT) Data: rich2013 AIC BIC logLik deviance REMLdev 93.04 101.3 -39.52 81.87 79.04 Random effects: Groups Name Variance Std.Dev. WT:Block (Intercept) 3.3667e+00 1.8348e+00 Block (Intercept) 8.2334e-19 9.0738e-10 Residual 6.1667e-01 7.8528e-01 Number of obs: 24, groups: WT:Block, 12; Block, 6 Fixed effects: Estimate Std. Error t value (Intercept) 12.0000 0.8148 14.728 NTNX -0.5000 0.4534 -1.103 WTS -3.6667 1.1523 -3.182 NTNX:WTS 3.3333 0.6412 5.199 Correlation of Fixed Effects: (Intr) NTNX WTS NTNX -0.278 WTS -0.707 0.197 NTNX:WTS 0.197 -0.707 -0.278 I have a few questions - i'm sorry if this is not the right place to ask them (they are quite simple!): In interpreting this model, am I correct in thinking that the output I get tells me the effect of one level of each factor compared with the other level of that factor, and also the effect of the interaction in terms of the effect of one combination of treatments compared to another combination of treatments? My main concern lies in how to report this model in a paper or thesis. It seems to be common practice to the "main effect" of each factor in a model (and perhaps the main effect of an interaction?), and then to discuss the difference between factor levels later. Is this correct? If the above is correct, in an lm, I would simple use summary.aov(model) to give me the summary of the main effects of each explanatory factor and the interaction and report these, along with their test statistic and p value (subjective, I know). However I do not know how to do this in lme4, or indeed if this is even a correct approach. Specifically then, I am wondering: a) Is the extraction of these "main effects" possible in lme4 and if so, how is it done? b) How best to present the results of the two treatments... c) How best to present the results of the interaction between the two treatments at all levels of this interaction (so for all combinations of factor levels that make up the 2x2 factorial experiment). I have struggled to find much information relating to how to present such models online... Many thanks in advance! Sarah
Presenting results of a mixed model containing factors
9 messages · Ben Bolker, Henrik Singmann, Sarah Dryhurst
2 days later
Sarah Dryhurst <s.dryhurst at ...> writes:
I am running a mixed effects model using the following code: m1<- lmer(DV~TMT*TMT2+(1|Block/TMT1),verbose=T) An example dataset is here: http://pastebin.com/bHug5kTt The two explanatory variables are both factors with two levels (treatment and control). Treatment 2 is split within Treatment 1 which is in turn within block. The model output I get is as follows (sorry for the formatting!) Linear mixed model fit by REML Formula: Richness ~ NT * WT + (1 | Block/WT) Data: rich2013 AIC BIC logLik deviance REMLdev 93.04 101.3 -39.52 81.87 79.04 Random effects: Groups Name Variance Std.Dev. WT:Block (Intercept) 3.3667e+00 1.8348e+00 Block (Intercept) 8.2334e-19 9.0738e-10 Residual 6.1667e-01 7.8528e-01 Number of obs: 24, groups: WT:Block, 12; Block, 6
Note here that your among-block variance is effectively zero ...
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.0000 0.8148 14.728
NTNX -0.5000 0.4534 -1.103
WTS -3.6667 1.1523 -3.182
NTNX:WTS 3.3333 0.6412 5.199
Correlation of Fixed Effects:
(Intr) NTNX WTS
NTNX -0.278
WTS -0.707 0.197
NTNX:WTS 0.197 -0.707 -0.278
I have a few questions - i'm sorry if this is not the right place to
ask them (they are quite simple!):
In interpreting this model, am I correct in thinking that the output I
get tells me the effect of one level of each factor compared with the
other level of that factor, and also the effect of the interaction in
terms of the effect of one combination of treatments compared to
another combination of treatments?
As is generally true in R modeling (with the default treatment contrasts), the main effects represent the difference between the treatment and the control at the baseline level of the other effect. The interaction represents the difference between the double-treatment effect and the combined (additive) effect of the two treatments.
My main concern lies in how to report this model in a paper or thesis. It seems to be common practice to the "main effect" of each factor in a model (and perhaps the main effect of an interaction?), and then to discuss the difference between factor levels later. Is this correct? If the above is correct, in an lm, I would simple use summary.aov(model) to give me the summary of the main effects of each explanatory factor and the interaction and report these, along with their test statistic and p value (subjective, I know). However I do not know how to do this in lme4, or indeed if this is even a correct approach.
Specifically then, I am wondering: a) Is the extraction of these "main effects" possible in lme4 and if so, how is it done?
The "main effects" are the effects you see, but they have to be interpreted carefully.
b) How best to present the results of the two treatments... c) How best to present the results of the interaction between the two treatments at all levels of this interaction (so for all combinations of factor levels that make up the 2x2 factorial experiment). I have struggled to find much information relating to how to present such models online...
Most of these questions aren't specific to lme4, it's just that the standard suite of answers get more complicated in _any_ situation (glm(), lmer(), etc.) where there isn't a simple, unique additive decomposition of effects. There is a fair amount of controversy even about how to handle main effects in the presence of interactions (the "SAS type III SSQ/marginality" argument): e.g. google for "venables exegeses linear models". There is a strong effect of the interaction in your case. If you want to get 'average' effects of the main effects, which may or may not make sense, you can use sum-to-zero contrasts. Ben Bolker
Hi Ben, Thank you for your reply! I don't really want to calculate the main effects as it doesn't make much biological sense (to me!). I just wasn't sure whether this was "required" in terms of statistical reporting. That interaction effect is what I am interested in largely, as it's the combined effect of the different treatments that is my focus. With regards to the lack of variance at the Block level, would you recommend dropping this level here? It doesn't seem to make too much sense to keep it there... Thank you once again Sarah
On Sun, Sep 8, 2013 at 10:47 PM, Ben Bolker <bbolker at gmail.com> wrote:
Sarah Dryhurst <s.dryhurst at ...> writes:
I am running a mixed effects model using the following code: m1<- lmer(DV~TMT*TMT2+(1|Block/TMT1),verbose=T) An example dataset is here: http://pastebin.com/bHug5kTt The two explanatory variables are both factors with two levels (treatment and control). Treatment 2 is split within Treatment 1 which is in turn within block. The model output I get is as follows (sorry for the formatting!) Linear mixed model fit by REML Formula: Richness ~ NT * WT + (1 | Block/WT) Data: rich2013 AIC BIC logLik deviance REMLdev 93.04 101.3 -39.52 81.87 79.04 Random effects: Groups Name Variance Std.Dev. WT:Block (Intercept) 3.3667e+00 1.8348e+00 Block (Intercept) 8.2334e-19 9.0738e-10 Residual 6.1667e-01 7.8528e-01 Number of obs: 24, groups: WT:Block, 12; Block, 6
Note here that your among-block variance is effectively zero ...
Fixed effects:
Estimate Std. Error t value
(Intercept) 12.0000 0.8148 14.728
NTNX -0.5000 0.4534 -1.103
WTS -3.6667 1.1523 -3.182
NTNX:WTS 3.3333 0.6412 5.199
Correlation of Fixed Effects:
(Intr) NTNX WTS
NTNX -0.278
WTS -0.707 0.197
NTNX:WTS 0.197 -0.707 -0.278
I have a few questions - i'm sorry if this is not the right place to
ask them (they are quite simple!):
In interpreting this model, am I correct in thinking that the output I
get tells me the effect of one level of each factor compared with the
other level of that factor, and also the effect of the interaction in
terms of the effect of one combination of treatments compared to
another combination of treatments?
As is generally true in R modeling (with the default treatment contrasts), the main effects represent the difference between the treatment and the control at the baseline level of the other effect. The interaction represents the difference between the double-treatment effect and the combined (additive) effect of the two treatments.
My main concern lies in how to report this model in a paper or thesis. It seems to be common practice to the "main effect" of each factor in a model (and perhaps the main effect of an interaction?), and then to discuss the difference between factor levels later. Is this correct? If the above is correct, in an lm, I would simple use summary.aov(model) to give me the summary of the main effects of each explanatory factor and the interaction and report these, along with their test statistic and p value (subjective, I know). However I do not know how to do this in lme4, or indeed if this is even a correct approach.
Specifically then, I am wondering: a) Is the extraction of these "main effects" possible in lme4 and if so, how is it done?
The "main effects" are the effects you see, but they have to be interpreted carefully.
b) How best to present the results of the two treatments... c) How best to present the results of the interaction between the two treatments at all levels of this interaction (so for all combinations of factor levels that make up the 2x2 factorial experiment). I have struggled to find much information relating to how to present such models online...
Most of these questions aren't specific to lme4, it's just that the standard suite of answers get more complicated in _any_ situation (glm(), lmer(), etc.) where there isn't a simple, unique additive decomposition of effects. There is a fair amount of controversy even about how to handle main effects in the presence of interactions (the "SAS type III SSQ/marginality" argument): e.g. google for "venables exegeses linear models". There is a strong effect of the interaction in your case. If you want to get 'average' effects of the main effects, which may or may not make sense, you can use sum-to-zero contrasts. Ben Bolker
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
NERC PhD Student Community Ecology and Global Change Department of Biology Imperial College, London email: sarah.dryhurst08 at imperial.ac.uk
Sarah Dryhurst <s.dryhurst at ...> writes:
Hi Ben, Thank you for your reply! I don't really want to calculate the main effects as it doesn't make much biological sense (to me!). I just wasn't sure whether this was "required" in terms of statistical reporting. That interaction effect is what I am interested in largely, as it's the combined effect of the different treatments that is my focus.
I think it would still be worth reporting the main effects, as their size puts the size of the interaction in perspective (e.g. I would generally like to be able to judge the size of the interaction *relative to* the main effects, not just its magnitude/t statistic/ p value ...
With regards to the lack of variance at the Block level, would you recommend dropping this level here? It doesn't seem to make too much sense to keep it there...
Doesn't matter too much, since the results will be almost identical. May be worth checking without, to double-check that the zero-variance result hasn't thrown off the optimization, but I would personally probably err on the side of reporting it (or say that it was in the original model but estimated as being effectively zero).
8 days later
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20130918/4c8de521/attachment.pl>
Hi Sarah,
One way to globally set sum to zero contrasts is by setting the contrasts via options:
options(contrasts=c('contr.sum', 'contr.poly'))
This is automatically done when loading afex which also provides you with the overall effect you are interested in (using the Kenward-Rogers approximation):
require(afex)
dat <- read.table("http://pastebin.com/raw.php?i=bHug5kTt", header = TRUE)
mixed(DV~TMT1*TMT2+(1|Block/TMT1), dat)
## Effect stat ndf ddf F.scaling p.value
## 1 (Intercept) 377.3464 1 5 1 0.0000
## 2 TMT1 3.2653 1 5 1 0.1306
## 3 TMT2 13.2433 1 10 1 0.0045
## 4 TMT1:TMT2 27.0271 1 10 1 0.0004
Alternatively you can get the same results from car::Anova:
Anova(m1)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: DV
## Chisq Df Pr(>Chisq)
## TMT1 3.2653 1 0.0707600 .
## TMT2 13.2433 1 0.0002736 ***
## TMT1:TMT2 27.0271 1 2.006e-07 ***
## ---
## Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1
As your design is fully balanced this is equivalent to type three tests (the default for afex):
Anova(m1, type = 3)
Hope that helps.
Cheers,
Henrik
Am 18.09.2013 18:20, schrieb Sarah Dryhurst:
Hi Ben, I'm sorry to dredge up this post again - I thought my reply had sent but I checked and it was sat in my drafts folder. My question was on the sum to zero contrasts you suggested. Is there a recommended way to calculate these in R, or any references I can read about them (I am a complete novice). I was looking into using anova or drop1 on my models to get these "average" effects of the main effects and interaction, however these two only seem to take into account fixed effects in their calculations (do they use Wald statistics?) which seems pointless to me. What is your view on these two? My design is a balanced one... I know the glmm wiki recommends using a mcmc approach: "Tests of effects (i.e. testing that several parameters are simultaneously zero)
From worst to best:
- Wald chi-square tests (e.g. car::Anova)
- Likelihood ratio test (via anova or drop1)
- *For balanced, nested LMMs* where df can be computed: conditional
F-tests
- *For LMMs*: conditional F-tests with df correction (e.g. Kenward-Roger
in pbkrtest package)
- MCMC or parametric, or nonparametric, bootstrap comparisons
(nonparametric bootstrapping must be implemented carefully to account for
grouping factors)"
I have done this via pvals.fnc(m1) in languageR but it still only gives the
a similar output to my original (with no way of knowing the average
effects). My concern is that the average effects seem to be what
supervisors/reviewers want reported i.e. the effect of Treatment1, rather
than the effect of Treatment1 level A compared to Treatment 1 level B
etc....
Any thoughts would be much appreciated! I'm finding it hard to find a
consensus anywhere. It is difficult to track down examples of reporting
these things - most focus seems to be on interpretation. Thank you once
again for your advice.
Sarah
On Mon, Sep 9, 2013 at 10:54 PM, Ben Bolker <bbolker at gmail.com> wrote:
Sarah Dryhurst <s.dryhurst at ...> writes:
Hi Ben, Thank you for your reply! I don't really want to calculate the main effects as it doesn't make much biological sense (to me!). I just wasn't sure whether this was "required" in terms of statistical reporting. That interaction effect is what I am interested in largely, as it's the combined effect of the different treatments that is my focus.
I think it would still be worth reporting the main effects, as their size puts the size of the interaction in perspective (e.g. I would generally like to be able to judge the size of the interaction *relative to* the main effects, not just its magnitude/t statistic/ p value ...
With regards to the lack of variance at the Block level, would you recommend dropping this level here? It doesn't seem to make too much sense to keep it there...
Doesn't matter too much, since the results will be almost identical. May be worth checking without, to double-check that the zero-variance result hasn't thrown off the optimization, but I would personally probably err on the side of reporting it (or say that it was in the original model but estimated as being effectively zero).
_______________________________________________ R-sig-mixed-models at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
Dipl. Psych. Henrik Singmann PhD Student Albert-Ludwigs-Universit?t Freiburg, Germany http://www.psychologie.uni-freiburg.de/Members/singmann
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20130918/e52695ac/attachment.pl>
Am 18.09.2013 19:29, schrieb Sarah Dryhurst:
Can I ask, do the sum to zero contrasts/type 3 anova tests take into account the random effects in these estimations?
Yes.
[snip]
Dipl. Psych. Henrik Singmann PhD Student Albert-Ludwigs-Universit?t Freiburg, Germany http://www.psychologie.uni-freiburg.de/Members/singmann
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20130919/058c6d22/attachment.pl>