# R linear mixed model results

I want to investigate whether Group1 affects Resp variable in my model. Group1 is an ordinal variable that can assume 4 values (1 2 3 4).

Group2 is also an ordinal variable taht can assume 5 values (0 1 2 3 4)

I create the null model (without Group1).

model.null =  lmer(Resp~Group2+Gender+Age+BMI+(1|Subject)+(1|Day_type),data=table_data,REML=FALSE)


I create the full model:

model.full = lmer(Resp~Group1+Group2+Gender+Age+BMI+(1|Subject)+(1|Day_type),data=table_data,REML=FALSE)


I run anova to see whether the 2 models are significantly different:

anova(model.null,model.full)
Data: table_data
Models:
model.null: Resp ~ Group2 + Gender + Age + BMI + (1 | Subject) + (1 | Day_type)
model.full: Resp ~ Group1 + Group2 + Gender + Age + BMI + (1 | Subject) + (1 |
model.full:     Day_type)
Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
model.null 11 18320 18391 -9148.8    18298
model.full 14 18314 18405 -9143.0    18286 11.588      3   0.008935 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


according to the output of anova I can conclude that are significantly different: i.e. Group1 affects the response.

Now I run

summary(model.full)


and get

Fixed effects:
Estimate Std. Error t value
(Intercept)  2.428580   0.462824   5.247
Group1   2   0.119003   0.163572   0.728
Group1   3   0.210836   0.171478   1.230
Group1   4   0.562406   0.196697   2.859
Group2   1   0.069754   0.139780   0.499
Group2   2   0.139545   0.148745   0.938
Group2   3   0.094811   0.162958   0.582
Group2   4   0.600394   0.214628   2.797
GenderMale   0.451459   0.095112   4.747
Age          0.005298   0.005544   0.956
BMI         -0.003777   0.008470  -0.446


Are all the 4 values of Group1 significantly different? What can I say looking at this output?

The anova test doesn't really say that the categories of Group1 are significantly associated with the outcome, but that that model provides a better fit for the data. The lmer() command doesn't output p-values, but there are a couple of workarounds: http://mindingthebrain.blogspot.se/2014/02/three-ways-to-get-parameter-specific-p.html

I use the lmerTest package so that I get the p-values (option #2 on the link above). You can then just install that package and rerun the models, and you'll get the p-values. You can also run this code on the p-values (option #1 in the link above):

t.values <- c(0.728, 1.230, 2.859)
2 * (1 - pnorm(abs(t.values)))
[1] 0.466613586 0.218697105 0.004249788


From the looks of it, it seems that Group1 category 4 is significant against the reference category 1, but groups 2 and 3 are not significant.

Since Group1 is an ordinal variable, and the estimated coefficients are increasing (though not exactly linear) you might try a model with Group1 as a continuous variable which would assume that the effect is linear. If you compare this model to the one you already have, and there is no significant difference, you have some evidence that the assumption of a linear affect may be reasonable so you can present that model instead.

• Thanks for your answer. I compared the 2 models (grouped,continuous) and I got a p=0.4555. I would like anyway to use the grouped values, so I can estimate the marginal mean response for each group and test whether they statistically differ.. Would it be a good approach? – gabboshow Sep 14 '15 at 10:02
• What I mean is: I use mixed models to assess an impact of the variable of interest (Group1) I stratify the responses according to Group1 categories and I compute the marginal mean for each group (adjusted for the other fixed effects) – gabboshow Sep 14 '15 at 10:04
• Yes, that's perfectly fine! – JonB Sep 14 '15 at 10:05