[R] specification for glmmPQL

[Andrew R. Criswell]

Hello Dr. Bates and group,

I understand, the attached data file did not accompany my original 
message. I have listed below the code used to create that file.

data.1 <- data.frame(subject  = factor(rep(c("one", "two", "three", "four",
                                             "five", "six", "seven", 
"eight"),
                                           each = 4),
                                       levels = c("one", "two", "three",
                                                  "four", "five", "six",
                                                  "seven", "eight")),
                     day      = factor(rep(c("one", "two", "three", "four"),
                                           times = 8),
                                       levels = c("one", "two", "three",
                                                  "four")),
                     expt     = rep(c("control", "treatment"), each = 16),
                     response = c(58, 63, 57, 54, 63, 59, 61, 53, 52, 62,
                                  46, 55, 59, 63, 58, 59, 62, 59, 64, 53,
                                  63, 75, 62, 64, 53, 58, 62, 53, 64, 72,
                                  65, 74))

mtrx.1 <- matrix(apply(data.1[, -4], 2, function(x)
                 rep(x, 100 - data.1$response)), ncol = 3, byrow = F)
mtrx.2 <- matrix(apply(data.1[, -4], 2, function(x)
                 rep(x, data.1$response)), ncol = 3, byrow = F)
                
data.2 <- data.frame(subject  = factor(c(mtrx.1[,1], mtrx.2[,1]),
                                       levels = c("one", "two", "three",
                                                  "four", "five", "six",
                                                  "seven", "eight")),
                     day      = factor(c(mtrx.1[,2], mtrx.2[,2]),
                                       levels = c("one", "two", "three",
                                                  "four")),
                     expt     = factor(c(mtrx.1[,3], mtrx.2[,3]),
                                       levels = c("control", "treatment")),
                     response = factor(c(rep("yes", nrow(mtrx.1)),
                                         rep("no", nrow(mtrx.2))),
                                       levels = c("yes", "no")))

#-------------------------------------------------------------------------------#


Douglas Bates wrote:

>On 9/4/05, Andrew R. Criswell <r-stats at arcriswell.com> wrote:
>
>>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.
>>
>
>I don't think your attached .Rdata file made it through the various
>filters on the lists or on my receiving email.  Could you send me a
>copy in a separate email message?
>
>
>>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
>>
>>
>>
>>
>>______________________________________________
>>R-help at stat.math.ethz.ch mailing list
>>https://stat.ethz.ch/mailman/listinfo/r-help
>>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
>>
>>
>>
>
>
>