[R] specification for glmmPQL

Andrew R. Criswell r-stats at arcriswell.com
Sun Sep 4 10:48:26 CEST 2005


Hello All,

I have a question regarding how glmmPQL should be specified. Which of 
these two is correct?

summary(fm.3 <- glmmPQL(cbind(response, 100 - response) ~ expt,
                        data = data.1, random = ~ 1 | subject,
                        family = binomial))

summary(fm.4 <- glmmPQL(response ~ expt, data = data.2,
                        random = ~ 1 | subject, family = binomial))

One might think it makes no difference, but it does.

I have an experiment in which 8 individuals were subjected to two types 
of treatment, 100 times per day for 4 consecutive days. The response 
given is binary--yes or no--for each treatment.

I constructed two types of data sets. On Rfile-01.Rdata (attached here) 
are data frames, data.1 and data.2. The information is identical but the 
data are arranged differently between these two data frames. Data frame, 
data.1, groups frequencies by subject, day and treatment. Data frame, 
data.2, is ungrouped.

Consistency of these data frames is substantiated by computing these 
tables:

ftable(xtabs(response ~ expt + subject + day,
             data = data.1))
ftable(xtabs(as.numeric(response) - 1 ~ expt + subject + day,
             data = data.2))

If I ignore the repeated measurement aspect of the data, I get, using 
glm, identical results (but for deviance and df).

summary(fm.1 <- glm(cbind(response, 100 - response) ~ expt,
                    data = data.1, family = binomial))

summary(fm.2 <- glm(response ~ expt, data = data.2,
                    family = binomial))

However, if I estimate these two equations as a mixed model using 
glmPQL, I get completely different results between the two 
specifications, fm.3 and fm.4. Which one is right? The example which 
accompanies help(glmmPQL) suggests fm.4.

summary(fm.3 <- glmmPQL(cbind(response, 100 - response) ~ expt,
                        data = data.1, random = ~ 1 | subject,
                        family = binomial))

summary(fm.4 <- glmmPQL(response ~ expt, data = data.2,
                        random = ~ 1 | subject, family = binomial))

Thank you,
Andrew




More information about the R-help mailing list