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

Question

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) : 
MCMC sampling is not implemented in recent versions of lme4
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.

Answer 1

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).

Answer 2

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?