2
$\begingroup$

Suppose I have a data frame like this (four groups, for each n = 12):

> head(my_data,n=8)

Subject Group Condition    Value
1       1     A    Cond_1 12.40407
2       1     A    Cond_2 14.89856
3       2     B    Cond_1 13.87088
4       2     B    Cond_2 14.31305
5       3     C    Cond_1 13.40773
6       3     C    Cond_2 13.48016
7       4     D    Cond_1 13.76183
8       4     D    Cond_2 12.60769
...
...

I want to examine if there are differences in Value between the four groups A, B, C, D, considering also the effect of Condition. So I have to do a mixed model. Before, I set specific contrasts:

> AvsB <- c(1,-1,0,0)
> AvsC <- c(1,0,-1,0)
> A vsD <- c(1,0,0,-1)
> AvsAvg <- c(1,0,0,0)
> contrasts(my_data$Group) <- cbind(AvsB,AvsC,AvsC,AvsAvg)

And fit the full model (suppose I have already compared the various models with the intercept model):

> model <- lme(Value~Group+Condition+Group:Condition, random=~1|Subject/Condition,
               data=my_data, method="ML")

However, sometimes the following error appears:

Error in MEEM(object, conLin, control$niterEM) :
Singularity in backsolve at level 0, block 1

Other times, the error doesn't appear but the output doesn't show all the contrasts. For example, here AvsAvg is disappeared:

Fixed effects: Value ~ Group + Condition + Group:Condition 
                              Value Std.Error DF   t-value p-value
(Intercept)               21.793124 1.2435711 44 17.524631  0.0000
GroupAvsB                 -0.255563 2.1539283 44 -0.118650  0.9061
GroupAvsC                  0.176721 2.1539283 44  0.082046  0.9350
GroupAvsD                  0.492208 2.1539283 44  0.228516  0.8203
ConditionCond_2           -0.016164 0.2029039 44 -0.079661  0.9369
GroupAvsB:ConditionCond_2 -0.270565 0.3514399 44 -0.769874  0.4455
GroupAvsC:ConditionCond_2  0.217103 0.3514399 44  0.617754  0.5399
GroupAvsD:ConditionCond_2 -0.277043 0.3514399 44 -0.788309  0.4347

Generally, it appears to me as if a maximum of three contrasts are shown.

Am I mistaking something with the contrasts or with the model?

Thanks in advance.

$\endgroup$
3
$\begingroup$

If you include the intercept/constant term in the model (which is included implicitly by default), then what you are essentially doing here is trying to fit 4 group means with 5 parameters: 4 contrasts plus the intercept. Obviously this isn't going to work, because there are more "unknowns" (parameters) than there are "equations" (group means).

To get the model to behave properly, you need to have 4 parameters for your 4 group means. This means that you need to either remove one of your contrasts (although evidently lme is doing this for you without asking!), or you need to suppress the intercept term.

The first solution would be the more common way of handling this. Conventionally we model $k$ group means with $k-1$ contrasts or dummies, plus the intercept.

If you want to go the second route, you can suppress the intercept term by inserting either +0 or -1 in the right-hand side of the formula that you pass to the model fitting function. For example, lme(Value~0+Group+Condition+Group:Condition, ...).

In any case, you need to be careful working with non-orthogonal contrasts. It is often not at all obvious what is actually being tested by any of the contrasts when they are embedded in a non-orthogonal set. I note in particular that the set of codes that you posted do not test the questions that you seem to think they are testing, either with or without the intercept in place of AvsAvg. To see this, we can invert the matrix of contrasts to see what are the actual linear combinations of group means being computed by each contrast (these are the rows in the inverted matrix). Here is what this looks like for the set of codes you posted if we suppress the intercept:

> AvsB <- c(1,-1,0,0)
> AvsC <- c(1,0,-1,0)
> AvsD <- c(1,0,0,-1)
> AvsAvg <- c(1,0,0,0)
> contrs <- cbind(AvsB,AvsC,AvsD,AvsAvg)
> rownames(contrs) <- LETTERS[1:4]
> contrs
  AvsB AvsC AvsD AvsAvg
A    1    1    1      1
B   -1    0    0      0
C    0   -1    0      0
D    0    0   -1      0
> solve(contrs)
       A  B  C  D
AvsB   0 -1  0  0
AvsC   0  0 -1  0
AvsD   0  0  0 -1
AvsAvg 1  1  1  1

So the first 3 contrasts would actually just represent the means of groups B, C, D (with signs reversed), and the last contrast would be the sum of all 4 group means.

Here is what is being testing if we remove AvsAvg and allow the intercept:

> contrs2 <- cbind(intercept=1, contrs[,1:3])
> contrs2
  intercept AvsB AvsC AvsD
A         1    1    1    1
B         1   -1    0    0
C         1    0   -1    0
D         1    0    0   -1
> solve(contrs2)
             A     B     C     D
intercept 0.25  0.25  0.25  0.25
AvsB      0.25 -0.75  0.25  0.25
AvsC      0.25  0.25 -0.75  0.25
AvsD      0.25  0.25  0.25 -0.75

This is of course a standard set of "effect codes", in this case testing the means of groups B, C, D against the mean of the group means (not against the mean of group A!), and the intercept being the mean of the group means.

It is often more convenient to specify this inverted matrix directly, and then invert it again to obtain the set of contrasts to actually use. However, it turns out that there is simply no way of asking your 4 questions about group means simultaneously in a single set of codes, because the AvsAvg contrast is necessarily going to be a perfect linear combination of the other 3 contrasts (if we know the answers to AvsB, AvsC, and AvsD, then we can always figure out the answer to AvsAvg).

As a final remark... none of this has anything to do with the fact that you happen to be working with a mixed model. It is basic ANOVA theory (okay, maybe not that basic).

$\endgroup$
0
$\begingroup$

I am not 100% sure about this one, but I think your contrast "AvsAvg" is producing a column in your design matrix that is exactly the same as the column for group A. This design matrix will then be singular and so cannot be inverted to obtain a solution for your fixed effects and contrasts. You get output because the redundant variable "AvsAvg" is automatically dropped. I suspect your contrast "AvsAvg" is meant to be A versus the average, and perhaps this has been coded incorrectly?

$\endgroup$
  • $\begingroup$ Yes, I AvsAVg is A versus average. How can I set this contrast? I thought that with linear models it's possible to set non-orthogonal contrasts. $\endgroup$ – this.is.not.a.nick Jan 22 '13 at 19:13
  • $\begingroup$ I would try $AvsAvg = L = (1,-1/3,-1/3,-1/3)$, so that if $\beta=(A,B,C,D)'$ then $L\beta=A-(B+C+D)/3$ which is the difference of group A versus the average of the other three groups. $\endgroup$ – dandar Jan 22 '13 at 21:21
  • 1
    $\begingroup$ A correction to the above comment: The $L$ vector should sum to $0$, so if you use $L=(3/4,-1/4,-1/4,-1/4)$ then $L\beta=A-(A+B+C+D)/4$ which is group A versus the average of the other four groups. $\endgroup$ – dandar Jan 22 '13 at 21:56
  • $\begingroup$ @dandar Unfortunately that won't work because this new $AvsAvg$ is a perfect linear combination of the other 3 contrasts, specifically, $AvsAvg = (AvsB+AvsC+AvsD)/4$. $\endgroup$ – Jake Westfall May 22 '13 at 22:10

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.