[R-sig-ME] Proportional odds model w/ random effects
Ken Knoblauch
ken.knoblauch at inserm.fr
Fri May 14 17:58:13 CEST 2010
Hi,
I've been looking at some data in which observers rated
pairs of images as being of siblings or not. There
were 32 observers and 30 pairs of images, half of
which were of siblings. Treating the ratings as ordinal responses,
I'm interested in fitting the data with a proportional odds model.
I would like to treat both Observer and Image pair as random effects.
I followed the recent discussion on the list on
proportional odds models with mixed-effects.
The clmm function from the ordinal package in the
current version only accepts 1 random effect, as far
as I can tell. I have played with MCMCglmm a bit
but don't feel comfortable with my understanding (yet)
about how best to choose priors. I put (hacked) together
a wrapper (see below) for glmer (or glm if there are no random effect
terms) that builds a model matrix for the fixed effects
part of the model which is combined with the random
effect terms before handing it over to glmer (or glm
if no random effects). It has a lot of imperfections
but it seems to give results for the cutpoints,
predictors and standard errors that are systematically
close to those from polr from MASS (when there are no fixed
effects in the model) or clmm with one random effect and
reasonably otherwise. I say systematically close,
because they always differ in the 2nd decimal place or so.
This may turn out to be an unimportant difference but I worry
because in the past, I've found glm and direct optimization
approaches to be more similar when things are set up right.
The deviances are quite different (usually by about a factor
of 2), but there may be more benign reasons for that (inclusion
of constants, etc.).
I would like to know if this is a reasonable approach or not,
I have played with convergence criteria and number of iterations
and have not found these to explain the differences.
I pasted in my code below for anyone interested and would be happy to
learn of its deficiencies (a few of which I'm already aware).
Thank you.
Ken
polmer <- function(formula, data, lnk = "logit",
cor = FALSE, ...){
# Get fixed and random terms
lTerms <- labels(terms(formula))
rTerms <- sapply(lTerms, function(x){
rT <- grep("|", x, fixed = TRUE)
length(rT) > 0
})
ranTerms <- lTerms[rTerms]
fixTerms <- lTerms[!rTerms]
ranNames <- as.vector(sapply(ranTerms, function(x)
strsplit(x, "\\| ")[[1]][2]
))
NoRanEff <- !(length(ranTerms) > 0)
# Make fixed-effect model matrix
Resp <- ordered(data[[as.character(formula[[2]])]])
Cuts <- seq(0, length(levels(Resp)) - 2)
cumRat <- as.vector(sapply(Cuts,
function(x) Resp <= x))
fRat <- gl(length(Cuts), nrow(data),
nrow(data) * length(Cuts))
X <- model.matrix(~ fRat - 1)
labs <- sapply(Cuts, function(x)
paste(x, "|", x+1, sep = "")
)
colnames(X) <- labs
fxForm <- if(!NoRanEff){
tmp <- formula
tmp[[2]] <- NULL
while( any(xx <- sapply(tmp[[2]],
inherits, what = "(" ))){
RE <- which(xx)
for(ix in RE) tmp[[2]][[ix]] <- NULL
tmp[[2]] <- tmp[[2]][[2]]
}
tmp } else formula
fX <- -model.matrix(fxForm, data)[, -1]
fX <- kronecker(matrix(rep(1, length(Cuts)), nc = 1),
fX)
X <- cbind(X, fX)
colnames(X)[-seq(length(Cuts))] <- fixTerms
p.df <- data.frame(cumRat = cumRat, X = X)
names(p.df) <- c("cumRat", colnames(X))
# Make random-effect model vectors
frm <- if(!NoRanEff){
tmp <- sapply(seq_len(length(ranNames)),
function(x)
rep(data[[ranNames[x]]], 10))
for(ix in seq_len(ncol(tmp)))
assign(ranNames[ix], tmp[, ix])
rxForm <- paste(paste("(", ranTerms, ")",
sep = "", collapse = " + "), " - 1")
as.formula(paste("cumRat ~ . + ", rxForm))
} else as.formula("cumRat ~ . - 1")
res <- if (NoRanEff)
glm(frm, data = p.df,
family = binomial(lnk), ...) else
glmer(frm, data = p.df,
family = binomial(lnk), ...)
if (inherits(res, "mer")) print(res, cor = cor) else
print(res)
invisible(res)
}
--
Ken Knoblauch
Inserm U846
Stem-cell and Brain Research Institute
Department of Integrative Neurosciences
18 avenue du Doyen Lépine
69500 Bron
France
tel: +33 (0)4 72 91 34 77
fax: +33 (0)4 72 91 34 61
portable: +33 (0)6 84 10 64 10
http://www.sbri.fr/members/kenneth-knoblauch.html
More information about the R-sig-mixed-models
mailing list