1

I estimated a mixed effect model with a nested random effect structure (participants were in different groups) with the lmer command of the lme4 package.
mixed.model <- lmer(ln.v ~ treatment*level+age+income+(1 | group/participant),data=data)
Then I bootstrapped the bootstrap command from the lmeresampler package because of the nested structure. I used the semi-parametric bootstrap.
boot.mixed.model <- bootstrap(model = mixed.model, type = "cgr", fn = extractor, B = 10000, resample=c(data$group,data$participant))
I can obtain bootsrapped confidence intervals via boot.ci (package boot) but in addition I want to report the coefficients' p-values. The output of the bootstrapped model boot.mixed.model provides only the bias and the standard error: 

Bootstrap Statistics :
         original        bias    std. error
t1*   0.658442415 -7.060056e-02  2.34685668
t2*  -0.452128438 -2.755208e-03  0.17041300
…

What is the best way to calculate the p-values based on these values?

user20650
  • 24,654
  • 5
  • 56
  • 91
Julian
  • 240
  • 1
  • 8

1 Answers1

0

I am unaware of the package called lmeresampler, and it seems to have been removed from cran due to compatibility issues (failed cran checks).

Also, the question does not include data and extractor is not defined, so the example is not reproducible. However the output is the same as you would get by using the bootMer function from lme4 so produce and example using the inbuilt function.

Basically this follows the example from the help(bootMer) page, but expanded for the specific problem. If the object returend by the lmeresampler package is similar, it will contain the objects used.

Reproducible example

library(lme4)
data(Dyestuff, package = "lme4")
fm01ML <- lmer(Yield ~ 1|Batch, Dyestuff, REML = FALSE)

Now the bootMer function simply requires a function, which outputs a vector of interesting parameters.

StatFun <- function(merMod){
    pars <- getME(merMod, c("fixef", "theta", "sigma"))
    c(beta = pars$fixef, theta = unname(pars$theta * pars$sigma), sigma = pars$sigma) ### <<== Error corrected
}

We can perform our bootstrapping by using the bootMer, which also contains parametric options in type (i suggest reading the details in the help(bootMer) page for more information)

boo01 <- bootMer(fm01ML, StatFun, nsim = 100, seed = 101)

Now for more precise p-values, I'd advice p-values greater closer to 1000 but for time reasons it might not be feasible in every circumstance.

Regardless the output is stored in a matrix t, which we can use to perform a simple Kolmogorov-supremum test:

H0 <- c(0, 0, 0)
Test <- sweep(abs(boo01$t), 2, H0, "-") <= H0 ###<<=== Error corrected
pVals <- colSums(Test)/nrow(Test)
print(pVals)
#output#
beta.(Intercept)            theta            sigma 
            0.00             0.12             0.00
Community
  • 1
  • 1
Oliver
  • 8,169
  • 3
  • 15
  • 37
  • Thanks for the answer! ```extractor``` was just a function in order to obtain FE coefficients and SDs of the random effects. Unfortunately, I was not able to reproduce your example. Trying to run ```Test``` I get: ```Error in abs(boo01) : non-numeric argument to mathematical function``` I found an interesting article by BMJ (2011; doi: https://doi.org/10.1136/bmj.d2304), the authors illustrate the calculation of p-values from confidence intervals. They suggest this approach: 1) Standard error = (upper limit CI – lower limit CI)/(2×1.96); 2) z = Est/SE; 3) p = exp(−0.717×z − 0.416×z^2). – Julian Oct 14 '19 at 09:54
  • It sounds like you simply copy-pasted the code. `boo01` would be equivalent to `boot.mixed.model` if i understand your question correctly. – Oliver Oct 14 '19 at 12:06
  • 1
    First I tested it with my data with the equivalent ```boot.mixed.model``` and then I tried to reproduce it with ```Dyestuff``` data and the code you provided but in both cases I cannot run ```Test``` because ```abs``` does not accept an ```boot.mixed.model```as input. – Julian Oct 14 '19 at 14:16
  • Ahh my bad. `boo01$t` – Oliver Oct 14 '19 at 15:57
  • I've updated the reproducible example. In my haste i made 2 minor errors, which made it rather hard to reproduce. my bad. – Oliver Oct 14 '19 at 16:00