lmer / bootMer - calculating p-values?

[David Reitter]

Hi all,

The latest versions of the popular 'lme4' package no longer provide an MCMC sampling function to generate p-values and confidence intervals. You may recall that this was problematic with any bot the most basic random effects structures anyway, and lme4 authors point to random effects with low variance as the culprit. If you use languageR to do your MCMC sampling ("pvals.fnc"), you are affected.

One of their suggestions is to use parametric bootstrapping via "bootMer".

I would like to run the code below by you. Questions: shouldn't the "mySumm2" function just return the fixed effects as the statistic of interest? (This was somewhat blindly taken from the bootMer documentation). I then calculate the t value on the fixed effects models and read a p value from the t distribution. I'm surprised there is no function provided - so, are there any caveats?

Thanks,
David



# straight from the bootMer documentation
mySumm2 <- function(.) {
c(beta=fixef(.),sigma=sigma(.),sig01=unlist(VarCorr(.)))
}

# calculate prob for one-tailed t distribution from sample SAMP
# assuming the effect is in the predicted direction
dr.tfun = function (samp) {pt((mean(samp)-0)/(sd(samp)/sqrt(length(samp))), length(samp)-1, lower.tail=(mean(samp)<0))}

# print P values [and more, without labels - fixme]
dr.boot.print = function (boot) { # print prob
for (row in as.data.frame(boot)) {print(dr.tfun(row))}
}

# bootstrap an lmer or glmer model
dr.bootstrap = function (model) {
# bootstrap, fixed nsim
bm = bootMer(model, mySumm2, nsim=20, .progress="txt")
dr.boot.print(bm);
return(bm)
}

[Moreno Coco]

Hi David,

I have also been looking into this issue very recently,
as I am re-analyzing some time-course fixation data.

(another hot-topic that still sparkles tons of controversy,
but nobody really discuss).

I have not found yet, a "solution" to the issue of extracting
the p-values for the estimates of the fixed effect, as there
really is no clear consensus on the usefulness of it, as well as,
there is yet no technically stable solution to this issue.

there is this post:

http://stackoverflow.com/questions/18443127/r-bootstrapped-binary-mixed-model-logistic-regression-using-bootmer-of-the-ne

where Ben Bolker (one of the lme4 developers) gives some
ideas on how to perform bootstrapping, and gives a
hint on how to extract the confidence intervals of
parameter estimates (?confint.merMod).

He also keeps a wiki where these issues are often discussed:

http://glmm.wikidot.com/faq

An interesting discussion by Douglas Bates and other people
about this issue can be found in:

http://rwiki.sciviews.org/doku.php?id=guides:lmer-tests

I hope this helps, and if you do find a solution, please
share it ... I will do the same if I find one :)

Moreno

Post by David Reitter
Hi all,
The latest versions of the popular 'lme4' package no longer provide
an MCMC sampling function to generate p-values and confidence
intervals. You may recall that this was problematic with any bot
the most basic random effects structures anyway, and lme4 authors
point to random effects with low variance as the culprit. If you
use languageR to do your MCMC sampling ("pvals.fnc"), you are
affected.
One of their suggestions is to use parametric bootstrapping via "bootMer".
I would like to run the code below by you. Questions: shouldn't the
"mySumm2" function just return the fixed effects as the statistic of
interest? (This was somewhat blindly taken from the bootMer
documentation). I then calculate the t value on the fixed effects
models and read a p value from the t distribution. I'm surprised
there is no function provided - so, are there any caveats?
Thanks,
David
# straight from the bootMer documentation
mySumm2 <- function(.) {
c(beta=fixef(.),sigma=sigma(.),sig01=unlist(VarCorr(.)))
}
# calculate prob for one-tailed t distribution from sample SAMP
# assuming the effect is in the predicted direction
dr.tfun = function (samp)
{pt((mean(samp)-0)/(sd(samp)/sqrt(length(samp))), length(samp)-1,
lower.tail=(mean(samp)<0))}
# print P values [and more, without labels - fixme]
dr.boot.print = function (boot) { # print prob
for (row in as.data.frame(boot)) {print(dr.tfun(row))}
}
# bootstrap an lmer or glmer model
dr.bootstrap = function (model) {
# bootstrap, fixed nsim
bm = bootMer(model, mySumm2, nsim=20, .progress="txt")
dr.boot.print(bm);
return(bm)
}
--
Dr. David Reitter
Assistant Professor of Information Sciences and Technology
Penn State University
http://www.david-reitter.com

[Jonathan Baron]

I think you omitted the most help page I know, which is the help page
on pvalues in the lme4 package itself:

http://finzi.psych.upenn.edu/library/lme4/html/pvalues.html

(That said, I did not check the stackoverflow page. The url is too
big.)

For me, where I used to use pvals.fnc(), I now use PBmodcomp() in the
pbkrtest package.

Jon

Post by Moreno Coco
Hi David,
I have also been looking into this issue very recently,
as I am re-analyzing some time-course fixation data.
(another hot-topic that still sparkles tons of controversy,
but nobody really discuss).
I have not found yet, a "solution" to the issue of extracting
the p-values for the estimates of the fixed effect, as there
really is no clear consensus on the usefulness of it, as well as,
there is yet no technically stable solution to this issue.
http://stackoverflow.com/questions/18443127/r-bootstrapped-binary-mixed-model-logisti
c-regression-using-bootmer-of-the-ne
where Ben Bolker (one of the lme4 developers) gives some
ideas on how to perform bootstrapping, and gives a
hint on how to extract the confidence intervals of
parameter estimates (?confint.merMod).
http://glmm.wikidot.com/faq
An interesting discussion by Douglas Bates and other people
http://rwiki.sciviews.org/doku.php?id=guides:lmer-tests
I hope this helps, and if you do find a solution, please
share it ... I will do the same if I find one :)
Moreno

Post by David Reitter
Hi all,
The latest versions of the popular 'lme4' package no longer provide
an MCMC sampling function to generate p-values and confidence
intervals. You may recall that this was problematic with any bot
the most basic random effects structures anyway, and lme4 authors
point to random effects with low variance as the culprit. If you
use languageR to do your MCMC sampling ("pvals.fnc"), you are
affected.
One of their suggestions is to use parametric bootstrapping via "bootMer".
I would like to run the code below by you. Questions: shouldn't the
"mySumm2" function just return the fixed effects as the statistic of
interest? (This was somewhat blindly taken from the bootMer
documentation). I then calculate the t value on the fixed effects
models and read a p value from the t distribution. I'm surprised
there is no function provided - so, are there any caveats?
Thanks,
David
# straight from the bootMer documentation
mySumm2 <- function(.) {
c(beta=fixef(.),sigma=sigma(.),sig01=unlist(VarCorr(.)))
}
# calculate prob for one-tailed t distribution from sample SAMP
# assuming the effect is in the predicted direction
dr.tfun = function (samp)
{pt((mean(samp)-0)/(sd(samp)/sqrt(length(samp))), length(samp)-1,
lower.tail=(mean(samp)<0))}
# print P values [and more, without labels - fixme]
dr.boot.print = function (boot) { # print prob
for (row in as.data.frame(boot)) {print(dr.tfun(row))}
}
# bootstrap an lmer or glmer model
dr.bootstrap = function (model) {
# bootstrap, fixed nsim
bm = bootMer(model, mySumm2, nsim=20, .progress="txt")
dr.boot.print(bm);
return(bm)
}
--
Dr. David Reitter
Assistant Professor of Information Sciences and Technology
Penn State University
http://www.david-reitter.com

--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.