2
$\begingroup$

I'm having difficulty modeling something… I think it's fairly simple, but I'm just not confident enough in my ability.

I setup an experiment with two treatment groups, six replicates in each treatment group, and measured "diversity" at three time points. I want to see whether "treatment" has an effect and if so, whether it has an effect across all timepoints. I have labelled an id for each of the replicates so that I can tell the model that I've repeatedly measured that replicate.

I want to do a repeated measures model and this is what I think I have to do:

lmer(diversity~treatment*time+(1|id)+(1|time)

Could someone tell me whether this is correct, or tell me how to find the correct model?

$\endgroup$
3
$\begingroup$

I get the impression you're just trying to basically do a repeated measures ANOVA but using multi-level modelling. You'll generally need to look up information on multiple regression and understand that well in order to step into the world of multi-level modelling. You should look through the site at the many many questions on using lmer.

Your model is incorrect. You want:

m <- lmer( diversity ~ treatment * time + (1|id) )

Given your description, time is not a grouping variable. The parentheses is where you specify your random effects structure(s) but only your grouping or nesting goes after the vertical bar. In your case the data are grouped by id. The random effect here is id and the intercept (1) varies across across it but none of your effects like time or treatment do. You may want to know if there is a random effect of time and / or treatment, allowing those to vary by id as well. For example, letting time and the intercept vary by id would be:

m <- lmer( diversity ~ treatment * time + (1+time|id) )

but there are other options. To see a further explanation of various random effect structures search the site and you shall find.

You should plot your data and check to see if the effect of time is some simple function (like linear or log) and include it as that. And you could test the function the standard way comparing categorical time predictor model to function time predictor model to see if categorical had more explanatory power, anova(modelTimeLogline, modelTimeCat).

Also, with your model specification you're not checking if your effect is present at all time points in the coefficients. It will only tell you if it's present at time point 0 because you have an interaction term.

m <- lmer( diversity ~ treatment * time + (1|id) )
anova(m)

That will give you an assessment of both the main effects and interaction as F values in a sequential ANOVA. If you need the main effect coefficient for the model then you first need to run just this model:

lmer(diversity ~ treatment + time + (1|id) + (1|time))

Note that it does not matter at all that the interaction is significant when it is in the model in order for the main effect coefficients to be uninterpretable. There's more you should know about multiple regression in order to do more advanced analysis like this.

$\endgroup$
1
$\begingroup$

I believe this is correct yes. Unless you expect there to be a linear response with time, you should indeed specify Time as a categorical random effect, as you have done: (1|Time).

You could also combine your random effects thus: (1|Time/ID)

Remember that when simplifying your model, even if the main effect of Time is not significant, it should be retained in the model since it is a random effect too.

lmer allows you to specify error distributions, so if you plot the residuals of your models and they don't look normal:

m<-lmer(diversity~treatment*time+(1|id)+(1|time))
plot(resid(m) ~ fitted(m),main="residual plot")

then you could try specifying another error distribution, which might improve the residual distribution. If your response data is a count for example, you could try a poisson distribution. If your response data were binomial, you would use a binomial error distribution.

e.g.

m<-lmer(diversity~treatment*time+(1|id)+(1|time),poisson)
$\endgroup$
  • $\begingroup$ No, that approach will inflate type I error. $\endgroup$ – Frank Harrell Jul 23 '13 at 11:09
  • $\begingroup$ Anova is ok, it's just that you can't use estimated trends to compose a hypothesis test for differences, i.e., to lower the number of degrees of freedom. This is double-dipping. $\endgroup$ – Frank Harrell Jul 23 '13 at 11:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.