[R-sig-ME] getting pvalues of fixed effects coefficients in a lmer model with bootMer

[espesser]

  Dear List,
I want to get the pvalues  of the coefficients of a lmer model using 
Bootmer.
I am specially interested in obtaining pvalues for  factor(s) with more 
than 2 levels.

Here is an example:

library(lme4)
data(cake)
str(cake)
'data.frame':   270 obs. of  5 variables:
  $ replicate  : Factor w/ 15 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 
1 1 ...
  $ recipe     : Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 2 2 2 2 ...
  $ temperature: Ord.factor w/ 6 levels "175"<"185"<"195"<..: 1 2 3 4 5 
6 1 2 3 4 ...
  $ angle      : int  42 46 47 39 53 42 39 46 51 49 ...
  $ temp       : num  175 185 195 205 215 225 175 185 195 205 ...

# fitting  the model
   cakenum.lmer_new = lmer(angle ~ recipe + temp + (1|recipe:replicate), 
cake)

  summary(cakenum.lmer_new)
Linear mixed model fit by REML ['merModLmerTest']
Formula: angle ~ recipe + temp + (1 | recipe:replicate)
    Data: cake

REML criterion at convergence: 1693.844

Random effects:
  Groups           Name        Variance Std.Dev.
  recipe:replicate (Intercept) 41.80    6.465
  Residual                     20.71    4.551
Number of obs: 270, groups: recipe:replicate, 45

Fixed effects:
              Estimate Std. Error        df t value Pr(>|t|)
(Intercept)   1.51587    3.67895 257.81000   0.412    0.681
recipeB      -1.47778    2.45625  42.00000  -0.602    0.551
recipeC      -1.52222    2.45625  42.00000  -0.620    0.539
temp          0.15803    0.01622 224.00000   9.746 <2e-16 ***

(please dismiss the SAS df and pvalues )


###   bootstrap confidence intervals

confint(cakenum.lmer_new, method="boot",nsim=N) -> cakenum.lmer_new.confint
Computing bootstrap confidence intervals ...

                                      2.5 %    97.5 %
sd_(Intercept)|recipe:replicate  4.9334498 7.9448858
sigma                            4.1213029 4.9908540
(Intercept)                     -5.5832955 8.6921784
recipeB                         -6.2738450 3.1064975
recipeC                         -6.5293123 3.2290811
temp                             0.1270242 0.1882827


###  bootstrap pvalues

# definition of the private external function

mafun<- function(fit) {
      return(fixef(fit))
  }

N=1000

  resu.boot=bootMer(cakenum.lmer_new, mafun, nsim=N)
  fixeboot = resu.boot$t

  # the code below is derived from pvals.fnc()
  prop <- colSums(fixeboot > 0)/N
  prop
(Intercept)     recipeB     recipeC        temp
       0.631       0.281       0.269       1.000

ans <- 2 * pmax(0.5/N, pmin(prop, 1 - prop))
  ans
[1] 0.738    0.562         0.538         0.001
Inter        recipeB      recipeC      temp        # corresponding fixed 
coeff

which looks plausible (compared to old pvals.fnc() results) and
in line with the preceding  confint results.


a) Is this procedure valid ?
b) Is it valid whatever are the complexity of the model and/or the 
random terms ?

Thanks a lot

R. Espesser

-- 
Robert Espesser
CNRS UMR  7309 - Université Aix-Marseille
5 Avenue Pasteur
13100 AIX-EN-PROVENCE

Tel: +33 (0)413 55 36 26