specification for glmmPQL
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