[R] lme, glmmPQL, multiple random effects

Andrew Beckerman
2005-01-05T16:33:42Z

Hi all -
R2.0.1, OS X

Perhaps while there is some discussion of lme going on.....

I am trying to execute a glmm using glmmPQL from the MASS libray, using
the example data set from McCullagh and Nelder's (1989, p442) table
14.4 (it happens to be the glmm example for GENSTAT as well).  The data
are binary, representing mating success (1,0) for crosses between males
and females from two populations of salamanders.  The idea is to fit a
fixed effect of Cross, and estimate random effects for females and
males. (data available as an *.rda for anyone who wants it...)

Following the advice from various postings from R-help and from Pinhero
and Bates,  I can successfully (I think) code an lme() model, using
pdBlocked, pdIdent and a dummy grouping variable in a grouped data
object.

 > load("~/RbinaryData.rda")
 > dd<-data.frame(RbinaryData,dum=factor(rep(1,120)))
 > names(dd)
[1] "Cross"  "Female" "Male"   "Mate1"
"dum"
 > summary(dd)
  Cross       Female        Male        Mate1        dum
  RR:30   1      : 6   1      : 6   Min.   :0.0000   1:120
  RW:30   2      : 6   2      : 6   1st Qu.:0.0000
  WR:30   3      : 6   3      : 6   Median :1.0000
  WW:30   4      : 6   4      : 6   Mean   :0.5833
          5      : 6   5      : 6   3rd Qu.:1.0000
          6      : 6   6      : 6   Max.   :1.0000
          (Other):84   (Other):84
m1<-lme(Mate1~Cross,data=gd,random=pdBlocked(list(pdIdent(~Female
-1),pdIdent(~Male-1))))

This model executes and estimates separate intercepts for both females
and males.  Good.

However, when I try and run this as a glmmPQL, the following error
ensues (traceback provided below). Am I trying something that is
impossible at the moment, or just doing something wrong?

Cheers
andrew

 >
m2<-glmmPQL(Mate1~Cross,data=gd,random=pdBlocked(list(pdIdent(~Female
-1),pdIdent(~Male-1))),family="binomial")
iteration 1
Error in getGroups.data.frame(dataMix, groups) :
	Invalid formula for groups

 > traceback()
8: stop("Invalid formula for groups")
7: getGroups.data.frame(dataMix, groups)
6: getGroups(dataMix, groups)
5: lme.formula(fixed = zz ~ Cross, random = list(numeric(0),
numeric(0)),
        data = list(Mate1 = c(1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0,
        1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0,
        1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0,
        0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1,
        0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0), Cross = c(1, 2, 1,
        2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2,
        1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2,
        1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1,
        3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3,
        4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 4, 3, 4, 3, 4, 3, 4, 3,
        4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4,
        3, 4, 3), Female = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
        3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6,
        6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9,
        9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11,
        12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 14, 14, 14,
        14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16,
        17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 19, 19, 19,
        19, 19, 19, 20, 20, 20, 20, 20, 20), Male = c(1, 14, 5, 11,
        4, 15, 5, 15, 3, 13, 1, 12, 2, 11, 1, 14, 3, 13, 4, 12, 2,
        15, 5, 14, 3, 13, 4, 12, 2, 11, 19, 9, 20, 7, 16, 8, 18,
        8, 19, 9, 17, 6, 16, 6, 17, 10, 20, 9, 20, 7, 18, 6, 19,
        10, 17, 10, 16, 8, 18, 7, 9, 19, 7, 20, 10, 18, 7, 16, 9,
        17, 6, 20, 8, 17, 6, 19, 7, 16, 10, 20, 8, 18, 9, 19, 6,
        18, 10, 16, 8, 17, 15, 2, 13, 4, 12, 1, 14, 1, 15, 2, 11,
        5, 11, 4, 12, 5, 15, 3, 13, 3, 11, 1, 14, 4, 12, 5, 14, 3,
        13, 2), wts = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), zz = c(2.37523506338983,
        2.19314662972470, 2.37523506338983, -2.30685171777016,
2.37523506338983,
        2.19314662972470, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        2.19314662972470, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        -2.30685171777016, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        2.19314662972470, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        -2.30685171777016, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        2.19314662972470, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        2.19314662972470, 2.19314662972470, 2.37523506338983,
2.19314662972470,
        -2.73839300553829, 2.19314662972470, 2.37523506338983,
-2.30685171777016,
        -2.73839300553829, -2.30685171777016, 2.37523506338983,
-2.30685171777016,
        -2.73839300553829, -2.30685171777016, 2.37523506338983,
-2.30685171777016,
        -2.73839300553829, 2.19314662972470, 2.37523506338983,
-2.30685171777016,
        -2.73839300553829, 2.19314662972470, 2.37523506338983,
2.19314662972470,
        2.37523506338983, -2.30685171777016, -2.73839300553829,
2.19314662972470,
        -2.73839300553829, 2.19314662972470, -2.73839300553829,
-2.49392989051313,
        2.27586790674987, 3.09612363937901, 2.27586790674987,
-2.49392989051313,
        2.27586790674987, -2.49392989051313, -2.48603224779624,
-2.49392989051313,
        2.27586790674987, -2.49392989051313, -2.48603224779624,
-2.49392989051313,
        -2.48603224779624, -2.49392989051313, -2.48603224779624,
        -2.49392989051313, 2.27586790674987, -2.49392989051313,
2.27586790674987,
        3.09612363937901, 2.27586790674987, -2.49392989051313,
2.27586790674987,
        -2.49392989051313, 2.27586790674987, -2.49392989051313,
-2.48603224779624,
        -2.49392989051313, -2.48603224779624, -2.48603224779624,
        -2.49392989051313, 2.27586790674987, -2.49392989051313,
-2.48603224779624,
        -2.49392989051313, 2.27586790674987, 3.09612363937901,
2.27586790674987,
        -2.49392989051313, 2.27586790674987, 3.09612363937901,
2.27586790674987,
        -2.49392989051313, 2.27586790674987, -2.49392989051313,
2.27586790674987,
        -2.49392989051313, 2.27586790674987, 3.09612363937901,
2.27586790674987,
        3.09612363937901, 2.27586790674987, -2.49392989051313,
2.27586790674987,
        -2.49392989051313, -2.48603224779624, 3.09612363937901,
2.27586790674987,
        -2.49392989051313), invwt = c(5.11362806892812, 4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214)), method = "ML", weights =
varFixed(~invwt))
4: lme(fixed = zz ~ Cross, random = list(numeric(0), numeric(0)),
        data = list(Mate1 = c(1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0,
        1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0,
        1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0,
        0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1,
        0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0), Cross = c(1, 2, 1,
        2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2,
        1, 2, 1, 2, 1, 2, 1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2,
        1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1,
        3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3,
        4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 4, 3, 4, 3, 4, 3, 4, 3,
        4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4, 3, 4,
        3, 4, 3), Female = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
        3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6,
        6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9,
        9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11,
        12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 14, 14, 14,
        14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 16,
        17, 17, 17, 17, 17, 17, 18, 18, 18, 18, 18, 18, 19, 19, 19,
        19, 19, 19, 20, 20, 20, 20, 20, 20), Male = c(1, 14, 5, 11,
        4, 15, 5, 15, 3, 13, 1, 12, 2, 11, 1, 14, 3, 13, 4, 12, 2,
        15, 5, 14, 3, 13, 4, 12, 2, 11, 19, 9, 20, 7, 16, 8, 18,
        8, 19, 9, 17, 6, 16, 6, 17, 10, 20, 9, 20, 7, 18, 6, 19,
        10, 17, 10, 16, 8, 18, 7, 9, 19, 7, 20, 10, 18, 7, 16, 9,
        17, 6, 20, 8, 17, 6, 19, 7, 16, 10, 20, 8, 18, 9, 19, 6,
        18, 10, 16, 8, 17, 15, 2, 13, 4, 12, 1, 14, 1, 15, 2, 11,
        5, 11, 4, 12, 5, 15, 3, 13, 3, 11, 1, 14, 4, 12, 5, 14, 3,
        13, 2), wts = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), zz = c(2.37523506338983,
        2.19314662972470, 2.37523506338983, -2.30685171777016,
2.37523506338983,
        2.19314662972470, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        2.19314662972470, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        -2.30685171777016, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        2.19314662972470, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        -2.30685171777016, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        2.19314662972470, 2.37523506338983, 2.19314662972470,
2.37523506338983,
        2.19314662972470, 2.19314662972470, 2.37523506338983,
2.19314662972470,
        -2.73839300553829, 2.19314662972470, 2.37523506338983,
-2.30685171777016,
        -2.73839300553829, -2.30685171777016, 2.37523506338983,
-2.30685171777016,
        -2.73839300553829, -2.30685171777016, 2.37523506338983,
-2.30685171777016,
        -2.73839300553829, 2.19314662972470, 2.37523506338983,
-2.30685171777016,
        -2.73839300553829, 2.19314662972470, 2.37523506338983,
2.19314662972470,
        2.37523506338983, -2.30685171777016, -2.73839300553829,
2.19314662972470,
        -2.73839300553829, 2.19314662972470, -2.73839300553829,
-2.49392989051313,
        2.27586790674987, 3.09612363937901, 2.27586790674987,
-2.49392989051313,
        2.27586790674987, -2.49392989051313, -2.48603224779624,
-2.49392989051313,
        2.27586790674987, -2.49392989051313, -2.48603224779624,
-2.49392989051313,
        -2.48603224779624, -2.49392989051313, -2.48603224779624,
        -2.49392989051313, 2.27586790674987, -2.49392989051313,
2.27586790674987,
        3.09612363937901, 2.27586790674987, -2.49392989051313,
2.27586790674987,
        -2.49392989051313, 2.27586790674987, -2.49392989051313,
-2.48603224779624,
        -2.49392989051313, -2.48603224779624, -2.48603224779624,
        -2.49392989051313, 2.27586790674987, -2.49392989051313,
-2.48603224779624,
        -2.49392989051313, 2.27586790674987, 3.09612363937901,
2.27586790674987,
        -2.49392989051313, 2.27586790674987, 3.09612363937901,
2.27586790674987,
        -2.49392989051313, 2.27586790674987, -2.49392989051313,
2.27586790674987,
        -2.49392989051313, 2.27586790674987, 3.09612363937901,
2.27586790674987,
        3.09612363937901, 2.27586790674987, -2.49392989051313,
2.27586790674987,
        -2.49392989051313, -2.48603224779624, 3.09612363937901,
2.27586790674987,
        -2.49392989051313), invwt = c(5.11362806892812, 4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        5.11362806892812, 4.49999834749486, 5.11362806892812,
4.49999834749486,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 4.49999834749486,
5.11362806892812,
        4.49999834749486, 5.11362806892812, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        5.59005352989214, 4.76190015454611, 5.59005352989214,
4.76190015454611,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214, 4.76190015454611,
5.59005352989214,
        4.76190015454611, 5.59005352989214)), method = "ML", weights =
varFixed(~invwt))
3: eval(expr, envir, enclos)
2: eval(mcall)
1: glmmPQL(Mate1 ~ Cross, data = gd, random =
pdBlocked(list(pdIdent(~Female -
        1), pdIdent(~Male - 1))), family = "binomial")

------------------------------------------------------------------------
---------
Dr. Andrew Beckerman
Department of Animal and Plant Sciences, University of Sheffield,
Alfred Denny Building, Western Bank, Sheffield S10 2TN, UK
ph +44 (0)114 222 0026; fx +44 (0)114 222 0002
http://www.shef.ac.uk/beckslab
------------------------------------------------------------------------
----------