1
$\begingroup$

I am fitting a linear mixed effect models with two factors (mPair with 6 levels, and spd_des with 3 levels) and their interaction. I obtain inconsistent results depending on the contrasts that I choose, and I would like to understand why and how to deal with it.

If I use treatment contrasts, I obtain the following (I only copy the relevant info of the results)

> options(contrasts = c("contr.treatment","contr.poly"))
> linM1 <- lme(cc_marg ~ mPair*spd_des , random = ~mPair|ratID, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )
> summary(linM1)

Fixed effects: cc_marg ~ mPair * spd_des 
                         Value  Std.Error DF   t-value p-value
(Intercept)          1.4628761 0.09618167 94 15.209511  0.0000
mPairRFVI           -0.8180718 0.10454920 94 -7.824754  0.0000
mPairVLRF           -0.7990828 0.13193991 94 -6.056415  0.0000
mPairVLVI           -0.6077804 0.13734253 94 -4.425289  0.0000
mPairVMRF           -0.7444267 0.13294167 94 -5.599649  0.0000
mPairVMVI           -0.4799995 0.12194383 94 -3.936234  0.0002
spd_des15           -0.0830016 0.07990370 94 -1.038771  0.3016
spd_des20           -0.0856339 0.08321984 94 -1.029008  0.3061
mPairRFVI:spd_des15 -0.1576193 0.13500809 94 -1.167481  0.2460
mPairVLRF:spd_des15  0.0866510 0.11385875 94  0.761039  0.4485
mPairVLVI:spd_des15  0.0083311 0.13500809 94  0.061708  0.9509
mPairVMRF:spd_des15  0.0184844 0.11385875 94  0.162345  0.8714
mPairVMVI:spd_des15 -0.0672286 0.13500809 94 -0.497960  0.6197
mPairRFVI:spd_des20 -0.1705514 0.14201095 94 -1.200973  0.2328
mPairVLRF:spd_des20  0.0899629 0.11949193 94  0.752879  0.4534
mPairVLVI:spd_des20 -0.0626845 0.14359174 94 -0.436547  0.6634
mPairVMRF:spd_des20 -0.0106400 0.11969131 94 -0.088895  0.9294
mPairVMVI:spd_des20 -0.0608750 0.14286017 94 -0.426116  0.6710

I interpret this as follows: Since the interaction terms are all non-significant, then the factor spd_des (also non-significant) does not influence the data at any level of the factor mPair.

On the other hand, using sum contrasts I obtain the following results.

> options(contrasts = c("contr.sum","contr.poly"))
> linM2 <- lme(cc_marg ~ mPair*spd_des , random = ~mPair|ratID, data=dat_trf, na.action=na.omit, method = "ML", control=lCtr )
> summary(linM2)

Fixed effects: cc_marg ~ mPair * spd_des 
                     Value  Std.Error DF   t-value p-value
(Intercept)      0.8137433 0.04791890 94 16.981678  0.0000
mPair1           0.5929117 0.06609665 94  8.970373  0.0000
mPair2          -0.3341386 0.04969616 94 -6.723629  0.0000
mPair3          -0.1472874 0.07260892 94 -2.028503  0.0453
mPair4          -0.0328631 0.08993236 94 -0.365421  0.7156
mPair5          -0.1488959 0.06991733 94 -2.129600  0.0358
spd_des1         0.0743293 0.02315254 94  3.210416  0.0018
spd_des2        -0.0272358 0.02325774 94 -1.171043  0.2445
mPair1:spd_des1 -0.0181081 0.04414399 94 -0.410206  0.6826
mPair2:spd_des1  0.0912334 0.05726538 94  1.593168  0.1145
mPair3:spd_des1 -0.0769866 0.04518813 94 -1.703691  0.0917
mPair4:spd_des1  0.0000066 0.05743544 94  0.000114  0.9999
mPair5:spd_des1 -0.0207337 0.04518548 94 -0.458859  0.6474
mPair1:spd_des2  0.0004559 0.04558473 94  0.010002  0.9920
mPair2:spd_des2 -0.0478225 0.05730295 94 -0.834555  0.4061
mPair3:spd_des2  0.0282279 0.04525282 94  0.623781  0.5343
mPair4:spd_des2  0.0269011 0.05747689 94  0.468034  0.6408
mPair5:spd_des2  0.0163141 0.04525367 94  0.360503  0.7193

> anova.lme(linM2,type="marginal")
              numDF denDF   F-value p-value
(Intercept)       1    94 288.37740  <.0001
mPair             5    94  35.30799  <.0001
spd_des           2    94   5.17279  0.0074
mPair:spd_des    10    94   0.52159  0.8710

The results are now telling me that the first level of the factor spd_des is significant; i.e. the mean of the data at that level of spd_des is significantly different from the grand mean (Intercept), and since the interactions are all non-significant, this is true at all levels of mPair.

So, with treatment contrasts spd_des does not influence the data at any level of mPair, and with sum contrast spd_des influence the data at all level of mPair. How do I deal with this? What result should I trust? Thanks in advance for your help.

You can find this question also at this link. Any help is highly appreciated.

$\endgroup$
  • $\begingroup$ please don't cross-post (this was also posted to r-sig-mixed-models@r-project.org, or, if you must, at least give links connecting the two posts ... $\endgroup$ – Ben Bolker Feb 21 at 23:27
  • $\begingroup$ I am sorry! In the past, I was told to re-post to the other channel after I asked a question. That's why I posted to both this time. I can cancel if that bothers, and certainly not post to both in the future. $\endgroup$ – Cristiano Feb 21 at 23:52
  • $\begingroup$ OK this time. Maybe edit your question above with this link to the r-sig-mixed post ? $\endgroup$ – Ben Bolker Feb 21 at 23:57
  • 1
    $\begingroup$ I have added the link as suggested $\endgroup$ – Cristiano Feb 22 at 0:04
  • $\begingroup$ Although I understand that you can't show all of your data, could you at least say whether each of the 18 combinations of factor levels has the same number of observations? A quick look suggests that they don't, but I could be wrong. Also, how many ratID values are there? I haven't used the type="marginal" anova setting in R; is that supposed to be like the "Type III" tests in SAS? $\endgroup$ – EdM Feb 25 at 20:11
1
$\begingroup$

There is an interesting technical difference between the contrast codings for which I have some ideas without a definitive answer. But there are also some statistical issues that might make that difference less relevant to your application.

Technically, the standard errors of the fixed-effect coefficient estimates found with the sum contrasts are of substantially smaller magnitude than are those based on the treatment contrasts. This suggests that the random effects modeled with sum contrasts did a better job of controlling for inter-rat differences than did those developed with treatment contrasts.

Without seeing the data it is hard to say exactly why this might be.

First, there seems to be a large set of missing data, with 180 data points intended (18 combinations with 1 test each on 10 rats) but with many fewer data points collected. Perhaps the pattern of missing values allowed better estimation with the sum contrasts.

Second, you have allowed for random effects related to mPair but not related to spd_des; you clearly wouldn't be able to evaluate random effects for both factors with so few data points. It's not clear to me how that would play out in terms of the interaction terms or the missing-data patterns.

Third, you are getting close to trying to estimate as many values from the data as you have observations. There are 18 fixed effects including the intercepts and the interaction terms. For random effects with 6 levels of mPair and an intercept (as I understand this coding) and 10 rats, you are also trying to determine 70 random effects here for a total of 88 values you are trying to estimate with a little over 100 data points.

I think that you would have to look at the design matrices of both the treatment-contrast and sum-contrast formulations produced from your specific data set, and how this plays out in the modeling of the random effects, to identify the reasons for this difference.

In practice, however, statistical considerations might make it less important to grapple with this problem. You have no indication of significant interaction terms. Once you demonstrate that point you could proceed simply to analyze the model without interactions and probably have more power to define precisley the remaining fixed-effect coefficients. I suspect that the differences between the results with treatment and sum contrasts would diminish if the interaction terms were removed.

Sum contrasts are a good standard and they allow the Type-III ANOVA tests that you wish to perform despite their violation of the principle of marginality. So it might be safest to continue with the sum contrasts.

How much sense Type-III tests actually make, however, can be questioned. Type-II tests might be more defensible here and would be applicable to both the treatment- and the sum-contrast formulations. As John Fox has said about crossed two-way ANOVA like this (with factors named A and B):

The most powerful test for the A main effect is the type-II test: for A after B ignoring the AB interaction; (and continuing) for B after A ignoring AB; and for AB after A and B. These tests are independent of the contrasts chosen to represent A and B, which is generally the case when one restricts consideration to models that conform to the principle of marginality.

W.N. Venables has made similar points about Type-III tests.

One last thought: it looks like spd_des might have a natural ordering (with some reference value and further values of 15 and 20). If so, you should be treating it as an ordered factor that can capture similarly ordered responses. It would then be handled with polynomial contrasts under your selection of options.

$\endgroup$
  • $\begingroup$ Thanks! I am starting to think, as you suggested, that my problem may be related to the fact that the number of coefficients to be estimated is close enough to the number of observations. Would it make sense to use a random intercept only model (as opposed to random intercept and slope, as currently is) to reduce the number of parameters? I have the feeling that removing non-significant interactions is controversial in the statistics community; see here. What do you think? $\endgroup$ – Cristiano Feb 27 at 21:40
  • $\begingroup$ @Cristiano the inclusion of random slopes in addition to intercepts ideally should be based on your knowledge of the subject matter. If you don't have enough data fairly to assess the random slopes, then you can't. An advantage of Type-II ANOVA is that you test the main effects of both factors first in a principled way (regardless of contrast coding) and then see if their interaction still matters. If the interaction is insignificant, then you already have an estimate of the main-effect significance. If interactions are significant, then the main effects aren't easily interpretable anyway. $\endgroup$ – EdM Feb 27 at 22:03
  • $\begingroup$ @Cristiano with respect to removing non-significant interactions, see for example Frank Harrell's course notes, Section 4.12.1 on developing models: "Check additivity assumptions by testing pre-specified interaction terms. Use a global test and either keep all or delete all interactions." So if the global test for interactions is insignificant then it can be OK to remove them all to simplify the model. Note that your model may be already overfit with 18 predictors (with interactions) and only a bit >100 observations; see Section 4.4 of those course notes. $\endgroup$ – EdM Feb 27 at 22:25
  • $\begingroup$ I have never really grasped the logic of how to choose the random effects. It sounds obvious to me that there is subject specificity, and therefore fitting a random slope model renders better predictions. And so, it does not sound fair to me to use fitting measures to do model selection. How do I then choose the random effects in exploratory studies like mine (i.e. I do not have previous knowledge on whether the slope should change consistently across subjects)? $\endgroup$ – Cristiano Feb 27 at 22:28
  • 1
    $\begingroup$ @Cristiano in exploratory studies use statistical tests to protect yourself from following up on false positives. Random effects can be tricky. See this page for useful ideas. The scale of your study might limit how many predictors you can consider. Forward selection of predictors based on individual relations to outcome is poor practice. Starting with a full model based on subject-matter knowledge and then removing less useful predictors can, however, be OK. See the Harrell notes linked above. $\endgroup$ – EdM Feb 28 at 1:16

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.