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/
[EMAIL PROTECTED] +61 (2) 6125 3379
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
