2
$\begingroup$

I am trying to create a contrast in R for a rectangular matrix.

pos_c1 = c(1/2,1/2,-1/6,-1/6,-1/6,-1/6,-1/6,-1/6)
posmat.temp = rbind(constant = 1/8,pos_c1)

in order to solve the matrix, I have read somewhere to use the function 'pseudoinverse'

posmat = pseudoinverse(posmat.temp)
posmat = posmat[, -1]
posmat

Then I have associated my matrix with the data:

contrasts(datasheet.complete$LastPosition) = posmat

but when I make my mixed-effect analysis with lmer(), I do not get the result for the contrast, but for each condition separately. To give an idea, with a square-matrix that calls the function solve(), in the output I get a p_value that tells me how much my contrast hypothesis is confirmed by the data analysis.

cl_c1 = c(0,0,-1,1)
cl_c2 = c(0,-1,0,1)
cl_c3 = c(-1,0,0,1)
closuremat.temp = rbind(constant = 1/4,cl_c1,cl_c2,cl_c3)
closuremat = solve(closuremat.temp)
closuremat = closuremat[, -1]
closuremat

contrasts(datasheet.complete$Closure) = closuremat

model = lmer(Score~Closure*ExpertiseType*LastPosition+(1|Participant)+(1|Item), data = datasheet.complete, REML = TRUE)
summary(model)

-> Closurecl_c1   ....   ....   p_value

I would expect therefore in the output a parameter called LastPositionpos_c1, but instead I get:

LastPosition1   ....   .....   p_value
Lastposition2   ....   .....   p_value
Lastposition3   ....   .....   p_value
....

Do you know what is the right function to use to solve a matrix which is not squared, and to obtain the result for the contrast hypothesis? Or can you suggest me how to move?


Ok, if I write the following:

conditions <- data.frame(Pos=factor(c(2,4,1,2,5,6,7,2,2,2,5,6,3,3,3,8,5,3,4,2,1,4,3,3,2,6,1,8,3,7,5,7,8,3,6,6,1,6,3,7)))
pos_c1 = c(1/2,1/2,-1/6,-1/6,-1/6,-1/6,-1/6,-1/6)
posmat.temp = cbind(pos_c1)
posmat = pseudoinverse(posmat.temp)
posmat = posmat[, -1, drop = FALSE]
posmat

I get:

     [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]
[1,] 0.75 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25

Looking at the other squared contrasts that work fine, it looks to me that this contrast should present 8 rows (conditions) and 1 column (the contrast hypothesis). Am I correct?

Following this fact, I get an error, which makes sense:

contrasts(conditions$Pos) = posmat
Error in `contrasts<-`(`*tmp*`, value = c(0.75, -0.25, -0.25, -0.25, -0.25,  : 
  wrong number of contrast matrix rows

Ok, here is an example of the issue:

library(lme4)
library(lmerTest)
library(corpcor)

database <- data.frame(
Clos=factor(c(4,4,1,4,4,3,2,1,2,1,2,2,4,3,1,2,1,4,1,3,2,2,4,4,4,4,2,1,4,2,2,1,4,2,4,2,1,4,4,3)),
Pos=factor(c(2,4,1,2,5,6,7,2,2,2,5,6,3,3,3,8,5,3,4,2,1,4,3,3,2,6,1,8,3,7,5,7,8,3,6,6,1,6,3,7)),
RF=c(8,6,2,9,7,1,7,6,3,4,6,4,5,2,5,5,3,4,1,3,1,2,3,1,2,2,3,1,8,5,2,2,7,1,9,4,5,6,4,2),
Score=c(4,3,3,5,4,3,2,4,5,2,2,3,3,4,4,4,3,2,3,3,5,4,3,4,4,2,3,4,3,4,1,2,2,2,3,4,5,3,1,2)
)
clos_c1 = c(0,0,-1,1)
clos_c2 = c(0,-1,0,1)
clos_c3 = c(-1,0,0,1)
closmat.temp = rbind(constant = 1/4,clos_c1,clos_c2,clos_c3)
closmat = solve(closmat.temp)
closmat = closmat[, -1]
closmat

pos_c1 = c(1/2,1/2,-1/6,-1/6,-1/6,-1/6,-1/6,-1/6)
posmat.temp = rbind(pos_c1)
posmat = pseudoinverse(posmat.temp)
posmat

contrasts(database$Clos) = closmat
contrasts(database$Pos) = posmat

model = lmer(Score~Clos+Pos+(1|RF), data = database, REML = TRUE)
summary(model)

The result of the summary(model) is

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  3.16682    0.23250  8.79795  13.621 3.26e-07 ***
Closclos_c1 -0.04582    0.61284 28.45101  -0.075   0.9409    
Closclos_c2 -0.25304    0.45563 28.59711  -0.555   0.5830    
Closclos_c3 -0.07088    0.49369 28.72373  -0.144   0.8868    
Pos1         0.99056    0.39208 26.60481   2.526   0.0178 *  
Pos2         0.03166    0.45087 28.87030   0.070   0.9445    
Pos3         0.45739    0.56255 26.22149   0.813   0.4235    
Pos4        -0.38971    0.51039 27.77335  -0.764   0.4516    
Pos5         0.19800    0.48213 28.97190   0.411   0.6843    
Pos6        -0.46234    0.52282 26.99040  -0.884   0.3843    
Pos7         0.28973    0.56143 25.70950   0.516   0.6102

As you can see, I have the p_value for the 3 contrast hypothesis named Closclos_c1, _c2 and _c3. I would therefore expect to have the p_value also for the rectangular contrast hypothesis called pos_c1, but instead I have the p_value for 7 out of 8 conditions.

What I would expect instead would be something like:

Closclos_c1
Closclos_c2
Closclos_c3
Pospos_c1
$\endgroup$
  • $\begingroup$ Probably posmat[, -1, drop = FALSE] will do. However, you don't need to rbind the constant when using the pseudoinverse. So posmat <- pseudoinverse(rbind(pos_c1)) or equivalently MASS::ginv(rbind(pos_c1)) should be all you need. $\endgroup$ – statmerkur May 18 '18 at 11:22
  • $\begingroup$ Also, you can name the contrast via colnames(posmat) <- ".pos_c1". $\endgroup$ – statmerkur May 18 '18 at 11:32
  • $\begingroup$ Thank you @statmerkur, now I have another problem. When setting the contrast for my parameter, I get the error "Error in cm[, 1L:nc] <- value : replacement has length zero" $\endgroup$ – Luca Danieli May 18 '18 at 12:31
  • $\begingroup$ I can have a look into it if you provide a reproducible example. $\endgroup$ – statmerkur May 18 '18 at 12:55
  • $\begingroup$ @statmerkur, you are right. I guess that there was a problem that something was "remembered", and got a false feedback. Now that I had to restart R, I get a different message, First of all, I understand that my 'posmat' (8 conditions) presents 7 columns (that I guess should be 8 rows). Secndly, the result of the contrast is "Error in contrasts<-(*tmp*, value = c(0.75, -0.25, -0.25, -0.25, -0.25, : wrong number of contrast matrix rows". I update now this problem in the main text $\endgroup$ – Luca Danieli May 19 '18 at 18:18
3
$\begingroup$

When computing the generalized inverse (aka pseudoinverse) of a matrix in order to specify custom contrasts in R one does not need to rbind (not to be confused with cbind) a constant to that matrix.

The following code shows how the desired matrix can be computed. Note that I use MASS::ginv() here but corpcor::pseudoinverse() works equally well.

conditions <- data.frame(Pos=factor(c(2,4,1,2,5,6,7,2,2,2,5,6,3,3,3,8,5,3,4,2,1,4,3,3,2,6,1,8,3,7,5,7,8,3,6,6,1,6,3,7)))
pos_c1 <- c(1/2,1/2,-1/6,-1/6,-1/6,-1/6,-1/6,-1/6)

posmat.temp1 <- rbind(pos_c1)
posmat1 <- MASS::ginv(posmat.temp1)
posmat1
      [,1]
[1,]  0.75
[2,]  0.75
[3,] -0.25
[4,] -0.25
[5,] -0.25
[6,] -0.25
[7,] -0.25
[8,] -0.25

If you really want to rbind the constant to the matrix and delete it afterwards, this can be achieved with the following code. drop = FALSE is important in this case as it will prevent R from coercing the matrix to a vector.

posmat.temp2 <- rbind(constant = 1/8, pos_c1)
posmat2 <- MASS::ginv(posmat.temp2)
posmat2 <- posmat2[, -1, drop = FALSE]

all.equal(posmat1, posmat2)
[1] TRUE

It is easy to name the contrast with:

colnames(posmat1) <- ".pos_c1"

A nice feature of R comes in handy now, it automatically adds some orthogonal contrasts to fill in the missing columns of the contrast matrix:

contrasts(conditions$Pos) <- posmat1
contrasts(conditions$Pos)
     .pos_c1                                                                  
1    0.75 -0.2886751 -0.2886751 -0.2886751 -0.2886751 -0.2886751 -0.2886751
2    0.75  0.2886751  0.2886751  0.2886751  0.2886751  0.2886751  0.2886751
3   -0.25  0.8333333 -0.1666667 -0.1666667 -0.1666667 -0.1666667 -0.1666667
4   -0.25 -0.1666667  0.8333333 -0.1666667 -0.1666667 -0.1666667 -0.1666667
5   -0.25 -0.1666667 -0.1666667  0.8333333 -0.1666667 -0.1666667 -0.1666667
6   -0.25 -0.1666667 -0.1666667 -0.1666667  0.8333333 -0.1666667 -0.1666667
7   -0.25 -0.1666667 -0.1666667 -0.1666667 -0.1666667  0.8333333 -0.1666667
8   -0.25 -0.1666667 -0.1666667 -0.1666667 -0.1666667 -0.1666667  0.8333333

Regarding your last example:

You do get the p-value for the contrast pos_c1 but corpcor::pseudoinverse() and MASS::ginv(), unlike solve(), drop the column/row names. Therefore you have to manually set the desired name for your contrast.
If you add colnames(posmat) <- "pos_c1" before you assign the matrix to the factor you get the same fixed effects table, but now "Pos1" has the expected name "Pospos_c1".

[...]
pos_c1 = c(1/2,1/2,-1/6,-1/6,-1/6,-1/6,-1/6,-1/6)
posmat.temp = rbind(pos_c1)
posmat = pseudoinverse(posmat.temp)
colnames(posmat) <- "pos_c1"
contrasts(database$Clos) = closmat
contrasts(database$Pos) = posmat
model = lmer(Score~Clos+Pos+(1|RF), data = database, REML = TRUE)
summary(model)
[...]
Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  3.16682    0.23250  8.79800  13.621 3.26e-07 ***
Closclos_c1 -0.04582    0.61284 28.45100  -0.075   0.9409    
Closclos_c2 -0.25304    0.45563 28.59700  -0.555   0.5830    
Closclos_c3 -0.07088    0.49369 28.72400  -0.144   0.8868    
Pospos_c1    0.99056    0.39208 26.60500   2.526   0.0178 *  
Pos          0.03166    0.45087 28.87000   0.070   0.9445    
Pos          0.45739    0.56255 26.22100   0.813   0.4235    
Pos         -0.38971    0.51039 27.77300  -0.764   0.4516    
Pos          0.19800    0.48213 28.97200   0.411   0.6843    
Pos         -0.46234    0.52282 26.99000  -0.884   0.3843    
Pos          0.28973    0.56143 25.71000   0.516   0.6102 
[...]
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ that's great! I have solved one issue. I am going towards my goal. Unfortunately I am still unable to get the p_value for the contrast. Any advice? $\endgroup$ – Luca Danieli May 20 '18 at 10:16
  • $\begingroup$ @Luca Danieli: Please, again, provide a reproducible example. Are you using lmerTest? $\endgroup$ – statmerkur May 20 '18 at 17:24
  • $\begingroup$ ok, I have added a complete reproducible example above. Yes, I am using lmerTest $\endgroup$ – Luca Danieli May 21 '18 at 16:14
  • $\begingroup$ @Luca Danieli: see my edit. Hope this clears things up. $\endgroup$ – statmerkur May 21 '18 at 17:06
  • $\begingroup$ Yes! This is exactly the situation I had in the beginning, but as I did not understand that the name is fundamental to read the output, I did not give it much importance. I get indeed an error about 'duplicate rows': 'Error in row.names<-.data.frame(*tmp*, value = value) : duplicate 'row.names' are not allowed'. How's that you don't? $\endgroup$ – Luca Danieli May 21 '18 at 17:31

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.