Posted by on Wed 23 Jan 2008 at 6:25 AM:
It looks like pvals.fnc may not work with lmer2 objects. I get this error when passing an lmer2 object:
Error in pvals.fnc(bySubj.lmer2) :
argument should be an lmer model object
Posted by Dale Barr on Thu 24 Jan 2008 at 12:03 AM:
Yes, that's correct, I've noticed that too. The pvals.fnc require an lmer object, and the lmer2 function gives you a lmer2 object, which won't work for whatever reason.
So if you want to get the p-values as computed by MCMC, you can do it using the following two functions (I found the first function floating around somewhere on the web, can't quite remember where). I confess to not understanding what it does, but seems to give the right results (checked against lmer). I wrote the second function myself (and it could probably stand some improvement, but it does the job).
mcmcpvalue <- function(samp) {
std <- backsolve(chol(var(samp)),
cbind(0, t(samp)) - colMeans(samp),
transpose = TRUE)
sqdist <- colSums(std*std)
sum(sqdist[-1] > sqdist[1])/nrow(samp)
}
pvalues <- function(mylmer) {
fes <- names(fixef(mylmer))
mx <- data.frame(Parameter=fes, pvalue=round(0,6))
mcmcres <- mcmcsamp(mylmer, n=50000)
for (i in 1:length(fes)) {
mx[i,"pvalue"] <- round(mcmcpvalue(as.matrix(mcmcres[,i])),6)
}
return(mx)
}
Then just pass your lmer2 object (e.g., "subj.lmer2") to the pvalues function like this:
my.p.values <- pvalues(subj.lmer2)
The number of samples for computing p in the example above is 50000, but you can change it by altering the line where it calls the mcmcsamp function.
Posted by on Mon 04 Feb 2008 at 10:33 AM:
In equation (6) of your JML (in press) paper, what is the purpose of adding .5 to numerator and denominator? Is it to avoid cases in which the ratio inside the parentheses = 0 ? The log of 0 is impossible.
Thanks
Posted by Dale Barr on Tue 05 Feb 2008 at 5:31 PM:
Yes, that is correct-- the .5 is to prevent the value from going to negative or positive infinity (when Y=0 or when Y=N).
The fact that it includes a correction factor is why it is called a quasi-logit (or empirical logit).
I found this equation originally in Agresti's textbook (Categorical Data Analysis) as well as in McCullagh & Nelder 1989, but the idea goes back much further (see for instance, Gart & Zweifel 1967, a paper that Florian Jaeger brought to my attention).
Gart & Zweifel (1967). On the bias of various estimators of the logit and its variance with application to quantal bioassay. Biometrica, 54, 181-187.