[R] Bradley Terry model and glmmPQL

[Simon Blomberg]

Dear R-ers,

I am having trouble understanding why I am getting an error using glmmPQL (library MASS).

I am getting the following error:

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

The long story:

I have data from an experiment on pairwise comparisons between 3 treatments (a, b, c). So a typical run of an experiment involves a rater choosing between a versus b, b versus c, or a versus c. The response is either 0 (not chosen) or 1 (chosen) for each treatment. I am interested in using the Bradley-Terry model to analyse this data. The Bradley-Terry model is a reparameterization of an ordinary logistic regression, using dummy variables to represent the treatment contrasts (see e.g. Agresti 1996 or Agresti 1990). It allows you to rank the treatments in accordance to preference. So my data looks something like this:

    rater age  a  b  c success failure
1       a   1  1 -1  0       1       0
2       b   1 -1  0  1       0       1
3       c   1  0  1 -1       1       0
4       d   1  1 -1  0       1       0
5       e   1 -1  0  1       0       1
6       f   1  0  1 -1       0       1
7       g   1  1 -1  0       0       1
8       h   1 -1  0  1       0       1
9       i   1  0  1 -1       1       0
10      j   1  1 -1  0       1       0
11      a   1 -1  0  1       0       1
12      b   1  0  1 -1       0       1
13      c   1  1 -1  0       0       1
14      d   1 -1  0  1       1       0
15      e   1  0  1 -1       1       0
16      f   1  1 -1  0       1       0
17      g   1 -1  0  1       1       0
18      h   1  0  1 -1       0       1
19      i   1  1 -1  0       1       0
20      j   1 -1  0  1       1       0
...

This is simulated data, but it corresponds to what my real data look like. Note: There are 10 raters (a to j), they repeat the experiment at two ages (1 and 2). Each rater rates all 3 treatment combinations at least once (and in my simulations up to 10 times). So for case 1, this corresponds to a trial between treatments a and b. a was chosen, which corresponds to a success. Similarly, for case 2 which is a trial between c and a, a was chosen, which corresponds to a failure. There are no ties.

The Bradley-Terry model can be fit using glm:

fit <- glm(cbind(success, failure) ~ c + b + nb -1, family=binomial, data=dat)

which works fine. I can also include the age factor, which also seems to work ok:

fit2 <- glm(cbind(success, failure) ~ c + b + nb + as.factor(age) -1, family=binomial, data=dat)

Now, since each rater performs ratings on each of the 3 treatment combinations, I was interested in including rater as a random factor. My naive method was to use glmmPQL from library MASS:

fit3 <- glmmPQL(cbind(success, failure) ~ c + b + nb -1, random = ~1|rater, family=binomial, data=dat)

However, I get the following error:

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

Can someone interpret this error for me, and tell me where I have gone wrong? Is there an alternative approach? Or am I thinking about this problem in the wrong way?

Thanks in advance,

Simon.

Simon Blomberg
Depression & Anxiety Consumer Research Unit
Centre for Mental Health Research
Australian National University
http://www.anu.edu.au/cmhr/
Simon.Blomberg at anu.edu.au  +61 (2) 6125 3379