The *best* approach is to have good arguments for a particular model. Model comparison based only on a set of observed data is neccesarily suboptimal.

Is the linear relationship between the response and time reasonable?

Is the chosen error model (e.g. normal) reasonable?

I don't understand why you want to compare the fixed-effects models when you know that you have longitudinal data? These models both ignore that the errors are partially correlated!

For the random effect (Time|Subject) be aware that here intercept and slope are correlated. If you want the model allow to fit a random intercept that is uncorrelated to the random slope, then use (Time-1|Subject) + (1|Subject).

Fit the model you deem resonable and then interpret the fitted coefficients. If you think this model is too complicated, then you should not built it in the first plase. Instead make a simpler model and check it's performance (w.r.t. to aim of your model, what in this case typically is prediction or forecast). I don't see the logic behind building different models and comparing them on a given set of data.

If you really have a couple of different models, then you should better explicitely state the (relative) prior credibility of these models. You can then use the data to adjust these credibilities. That would be a "Bayesian model comparison", if you like. In this way you at least have some clear statements about your models, conditioned on the data, whereas your model comparisons only give you information about the data, conditioned on the models (which credability is nowhere stated).

The *best* approach is to have good arguments for a particular model. Model comparison based only on a set of observed data is neccesarily suboptimal.

Is the linear relationship between the response and time reasonable?

Is the chosen error model (e.g. normal) reasonable?

I don't understand why you want to compare the fixed-effects models when you know that you have longitudinal data? These models both ignore that the errors are partially correlated!

For the random effect (Time|Subject) be aware that here intercept and slope are correlated. If you want the model allow to fit a random intercept that is uncorrelated to the random slope, then use (Time-1|Subject) + (1|Subject).

Fit the model you deem resonable and then interpret the fitted coefficients. If you think this model is too complicated, then you should not built it in the first plase. Instead make a simpler model and check it's performance (w.r.t. to aim of your model, what in this case typically is prediction or forecast). I don't see the logic behind building different models and comparing them on a given set of data.

If you really have a couple of different models, then you should better explicitely state the (relative) prior credibility of these models. You can then use the data to adjust these credibilities. That would be a "Bayesian model comparison", if you like. In this way you at least have some clear statements about your models, conditioned on the data, whereas your model comparisons only give you information about the data, conditioned on the models (which credability is nowhere stated).

Hello Jochen, first of all thank you for your reply.

I have good reasons to believe that the grouping variable is an important factor to use as an input for the model. My goal is just to prove it formally.

The correlation between observations is not being ignored, given that it is being modeled through the random effects part. I just started building the model from the fixed-effects part but would inevitably include the random effects part.

My main doubt was whether start building the model from the fixed-effects or the random effects part.

I have found in a book (Mixed Effects Models and Extensions in Ecology with R) a strategy to model construction called top-down. It states (page 121):

1. Start with a beyond optimal model (maximum possible fixed-effects).

2. Find the optimal structure of the random component (through AIC or BIC) comparisons.

3. Find the optimal fixed structure (with t-tests or wald z for the fixed-effects).

The author also states a step-up strategy but it is not further explained in the book.

As I said, there is no really good strategy to build a model based only on data. You cannot substitute a scientific idea, an exogenous aim (what should the model be good for?) and subject knowledge by statistics (t-values, p-values, information criteria etc) calculated from the data at hand. That's why there are so many proposals - and it is impossible to decide which one is better (or "good" in a more absolute meaning).

As long as you clearly explain what you do (and ideally why you do it this way and not the other), everything should be fine, IMHO.

Hi Felipe,

I generally agree with Jochen’s answer and would like to outline the approach I’ve learned to investigate research questions like yours:

1. Start by building your full model including all fixed and random effects:

• lmer(response ~ time:groups + time + groups + (1 + time|Subject), REML = FALSE, data = yourdata)

If you want to exclude the correlation between random intercept and slope use the following (or what Jochen suggested):

• lmer(response ~ time:groups + time + groups + (0 + time|Subject) + (1|Subject), REML = FALSE, data = yourdata)

2. Then build your intercept-only or null-model. This model should include all random effects and all fixed effects that you consider as control variables rather than test variables (which has to be a theoretical decision based on your study question and what you already know about your study organism)

3. Compare your full model with your null-model using the anova-function:

• anova(full, null, test = “Chisq”)

4. Only in case the full model is significantly better than the null model, proceed and establish the p-values for the fixed effects using the drop1 function:

• drop1(full, test = “Chisq”)

I don’t think intermediate models (e.g., without random effects) are valid if you have repeated measurements.

I find that one of the best and up-to-date sources on mixed models in the internet this this FAQ:

https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html

Hope that helps,

Urs

I forgot to add:

drop1 will not give you any results for the variables included in the interaction. Thus, if the interaction is not significant, you have to exclude it from the model and use the drop1 function again (in case the interaction is significant, the p-values for the main effects are not interpretable because the two variables depend on each other).

I agree to all what @Jan said but with this point: "if the interaction is not significant, you have to exclude it from the model".

Why should one have to exlude "non-significant" terms from the model? The term is/should be in the model because because it makes sense. If the data is not sufficient to interpret the sign of the MLE for the respective coefficient, it does not mean that this term has to be excluded. This would be, btw, the interpretation of a non-significant result - what would be logically flawed! Additionally, removing the interaction can have a considerable impact on other coefficients in the model (especially on those which are also part of the interaction). So one might get a considerably wrong picture there (again, assuming the interaction would make sense).

This is very similar for all coefficients (simple or interactions). It is common practice to simplify models by removing those terms whch coefficients are non-significant. This is a very bad practice, often leading to bad models, with little or no theoretical justification and very poor performance in "creating insights" and in prediction. And again: removing a predictor because it was "not-significant" mans to interpret a non-significant result, saying that it means that the coefficient is zero (up to some negligible difference). But that's not a conclusion that is supported by a non-significant result. In fact, the coefficient may be really much different from zero, and this can be very relevant - it was just not enough data available to see this in front of the noise in the data.

I did not mean to drop the interaction and then re-calculate the p-value of the full-model vs. the null-model. This should be done with the inclusion of the interaction.

However, it's not possible to establish p-values for the main effects using the drop1 function as long as the interaction is in the model. Therefore, it's necessary to remove the interaction and use the drop1 function again if you are interested in the p-values for the main effects.

In my understanding, starting with the full-null model comparison, and only proceed with the establishment of p-values for single predictors if the full-null model comparison leads to a p-value below the chosen alpha value should prevent data dredging. Or how would you calculate the p-values for main effects in case the interaction is not significant?

The meaning of the coefficients in the model change when you remove the interaction. Even if the contrasts are build to have a comparable meaning of the "main effects" coefficients in the models with and without interaction, both models will get different estimates of the residual variance. Take-home message: be very clear about what you really want, and then find a way how to get it. Don't adjust your question to the tools you know (and then misinterpret the result)!

And to repeat myself: non-significance does NOT mean absence. Judging a term in a model as negligible because it's non-significant is the interpretation of non-significance. That's a no-go.

Hello Jochen, I'm not sure I understood when you say that non-significant terms must not be excluded from the model.

That means keeping the non-significant terms in the model in the case where there are sufficient evidence that these terms are relevant?

I see the logic where the non-significance does not mean absence of evidence, just not enough evidence (or data). However, how to proceed in a case where the model is being used to discover if these effects influence the overall response?

As to the way I see it, if your theory tells you to include X1, X2 and the interaction between X1 and X2, and the interaction is not significant, you should not drop it. First, because dropping it implies to impose the restriction that the effect of X1 # X2 is zero. This is usually false. Usually you estimate a non-zero effect. You may not be able to claim with, say, 95% confidence, that it is different from zero. But this is very different from saying that you can claim that it is zero! I think that this was one of the things @Jochen was hinting at.

Hi Felipe, I said that excluding a term just because it's non-significant is a no-go. If you had evidence or agruments that a term is (or can be) relevant, it should be part of the model - indendent if the data at hand is sufficient to give a "significant" p-value (for the null hypothesis that the coefficient is zero). So, I'd say "yes" to your first question.

To your second question: if your aim is to discover the effect of a term in a model, you need to plan the experiment to give enough data to get a sufficiently precise estimate of the coefficient (often indirectly controlled by setting a minimum relevant size and a power). After all, one would likely need to establish the performance of the model in comparison with an alternative (possibly simpler) model. In what range does it work good enough? Where does it fail? And eventually, it's the theory behind the models that should suggest (or reject) what terms are in a model. If there is no theory, the model is only descriptive, and you can only look how good the predcitions are.

I forgot to say that I believe there is a second reason for not doing so, connected to the fact that if you run and re-run the model as many times as you want for the same data, you are actually "tricking" with the thresholds related to given pvalues (if you do not adjust them to take into account that you have carried out repeated runs of the model).

Hello Jochen,

Thank you. This answer explains perfectly both of my questions.

Hej Urs,you said

" drop1 will not give you any results for the variables included in the interaction. Thus, if the interaction is not significant, you have to exclude it from the model and use the drop1 function again (in case the interaction is significant, the p-values for the main effects are not interpretable because the two variables depend on each other)"

Is this a particular problem of the drop1 algorythm? Because more generally, the p-values are of course interpretable, just not as main but conditional effects (which are often very fundamental to really understand the interaction). Also, is there still a good/strong reason to rely on drop1 anymore after the development of lmerTest-package which gives you p-value estimates directly?