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
>>
>>
>>
>
>
>