[R-sig-ME] Proportional odds model w/ random effects

[Ken Knoblauch]

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