Alternatives to pvals.fnc to compute confidence intervals for fixed effects?

Lately I keep encountering the same problem and I'm wondering whether other people have been able to get around it. I'm running a mixed effects model using lmer(). My model has by-subject and by-item intercepts and slopes, and random correlation parameters between them. Since the current version of lmer() does not have MCMC sampling implemented, I cannot use pvals.fnc(). I get this message:

Error in pvals.fnc(m, withMCMC = T) :
for models with random correlation parameters


pvals.fnc() is also the function I use to get confidence intervals (HPD95lower and HPD95upper were two columns in the pvals.fnc output). Does anyone know of an alternative way of getting confidence intervals for the fixed effects estimates in the model? Or does using models with random correlations means that we can no longer get CIs from R?

Thanks!

NOTE: I've seen this question asked in other forums in slightly different ways. However, the answers always seem to involve (1) calculating something different as an alternative to the confidence intervals, (2) some complicated solution that is unclear (at least to me) how to implement. I would like to know if there is some alternative way of computing CIs that is both mainstream (so that other researchers can use it) and has a function to do it in R, since I am not a programmer and I feel that trying to create that function myself would be error prone.

• Note that a confidence interval/CI is not the same as a Highest Posterior Density/HPD interval, or a Bayesian Credible Interval. – jona Oct 20 '13 at 13:50
• I would also note that the ?pvalues help page in the the new version of lme4 includes a lot of information on this topic. – Ben Bolker Oct 21 '13 at 18:54
• @BenBolker: Thanks! I hadn't updated lme4 in a while and I didn't know about these useful features. I am trying to use confint now, trying to get profile CIs. Do you know if there is a simple way to have it estimate CIs only for the fixed effects? Otherwise, given my the structure of model, it takes forever to run and it also stops if any of the parameters cannot be estimated (which is likely when you have many, I think) – Sol Oct 22 '13 at 4:41
• I know you have params= to specify the parameters you want to get CIs for, but this requires you knowing which number corresponds to which parameter, and I don't know how to do that – Sol Oct 22 '13 at 4:42
• that's something we wanted to improve. The number of random-effects parameters (which come first) is n_ran <- length(getME(model,"theta")); the number of fixed-effect parameters is n_fix <- length(fixef(model)). Thus you should be able to use params=(n_ran+1):(n_ran+n_fix) to get just the fixed effect parameters profiled. – Ben Bolker Oct 22 '13 at 13:03

Most probably the packages lmerTest and lsmeans provide readily available routines for what you are looking for. Mind you, neither of them uses MCMC methodology. If you want to use something resampling-based, you can use lme4's native bootMer() function to bootstrap your model and get parametric bootstrap estimates (ver. 1.0-4 or newer).

This answer would have deserved comment status at best, but comments are too short and don't really lend themselves to extended mock code. Also, it seems you already got a more sensible answer by @user11852, but I wanted to give a more general answer (though see the comments below!).

If all else fails, one may always (?) obtain CIs from bootstrapping. In your specific case, this may be computationally infeasible, since running the model 1000 or so times may take half a century, but it should be fairly fool proof. I don't know R well, so here is some mock code in fake matlab for the 95% CI for the output generated by some parameter estimation function, such as intercepts in lmer. As a special feature, it bootstraps individual subjects and for each subject generates bootstrapped samples of data points for this subject.

a = data_set
s = (# of bootstrap iterations, e.g. 1000)
n = (# of subjects)

% main loop over bootstrap iterations
for x = 1:s

% bootstrap over subjects

% sample with replacement from your subjects (just collect n indices)
z = random_sample_with_replacement(1:n)

% loop over bootstrap selection of subjects
% to bootstrap sample individual data points within subject
for y = 1:length(z):
% for each subject selected for resampling, draw from their data points with replacement
for w = 1:length(a(z(y)))
within_subj_boot(w) = random_sample_with_replacement(a(z(y)))
end
bootsample(y) = within_subj_boot
end

% perform the respective calculation (e.g., lme4) for the bootstrap sample
% and store the relevant parameter you want a CI for
output(x) = parameter_estimation_function(bootsample)

% check the relevant percentiles of your bootstrapped parameter estimates
CI = percentile(output,[5,95])


Most languages will have some wrapper function for something equivalent to this (but a lot less inefficient), I think for R it might be the Coin package?

• jona, your train of thought is correct but take notice two things: 1. What you described in $non$-parametric bootstrap, not "simple" bootstrap; there is an inherent bias-variance trade-off between the two. 2. You need to be a bit careful how you resample your sample. You may accidentally end up missing a grouping especially if you have a lot of clusters. That is not "the end of the world" and would asymptotically "not happen" but this might mess-up your calculation procedures slightly. – usεr11852 Oct 20 '13 at 14:34
• Regarding 2, I agree - it's just an example that would need to be adapted to the individual model. (You might get by simply sampling from the grouping variables too.) Regarding 1 - I don't understand, can you elaborate? – jona Oct 20 '13 at 14:41
• With any bootstrapping technique the simulations get processed just like the real data. With non-parametric bootstrapping (what you described) you resample your original data. With parametric bootstrapping you simulate a new sample based on the original model you fitted. Non-parametric btsp. makes less assumptions but usually has more variance. Parametric btsp. assumes that the model you fit is "correct", so it makes more assumptions, but it usually has less variance. Param. btsp. also eliminate issues regarding the resampling. (Cont.) – usεr11852 Oct 20 '13 at 15:01
• Given the fact you are making parametric assumptions to start with, when you fit your original model you might as well use them and get a better estimate (ie. if you don't believe the model why bootstrap it anyway). You have the correct idea; I just want to highlight though that there is a trade-off between the non-parametric bootstrap you outlined and the parametric bootstrap that lme4's native bootMer() function offers. – usεr11852 Oct 20 '13 at 15:04
• @jona: Thanks! And thanks for taking the time to explain this as well. I am comfortable with bootstrapping (much more than with mixed effects models anyway) and I usually use the boot package in R. And thanks to user11852, now I know that I can also use bootMer() as well. Thanks both! – Sol Oct 20 '13 at 16:41