Power analysis from scratch

If pilot data is not available, simr can be used to create lme4 objects from scratch as a starting point. This requires more paramters to be specified by the user. Values for these parameters might come from the literature or the user’s own knowledge and experience.

library(simr)

Covariates and parameters

First set up some covariates with expand.grid.

x <- 1:10
g <- letters[1:3]

X <- expand.grid(x=x, g=g)

Specify some fixed and random parameters.

b <- c(2, -0.1) # fixed intercept and slope
V1 <- 0.5 # random intercept variance
V2 <- matrix(c(0.5,0.05,0.05,0.1), 2) # random intercept and slope variance-covariance matrix
s <- 1 # residual standard deviation

Build a model object

Use the makeLmer or makeGlmer function to build an artificial lme4 object.

model1 <- makeLmer(y ~ x + (1|g), fixef=b, VarCorr=V1, sigma=s, data=X)
print(model1)

## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ x + (1 | g)
##    Data: X
## REML criterion at convergence: 95.0003
## Random effects:
##  Groups   Name        Std.Dev.
##  g        (Intercept) 0.7071  
##  Residual             1.0000  
## Number of obs: 30, groups:  g, 3
## Fixed Effects:
## (Intercept)            x  
##         2.0         -0.1

model2 <- makeGlmer(z ~ x + (x|g), family="poisson", fixef=b, VarCorr=V2, data=X)
print(model2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: z ~ x + (x | g)
##    Data: X
##      AIC      BIC   logLik deviance df.resid 
## 206.3039 213.3099 -98.1520 196.3039       25 
## Random effects:
##  Groups Name        Std.Dev. Corr
##  g      (Intercept) 0.7071       
##         x           0.3162   0.22
## Number of obs: 30, groups:  g, 3
## Fixed Effects:
## (Intercept)            x  
##         2.0         -0.1

Start the power analysis

Now we have “pilot” models, which can be used with simr.

powerSim(model1, nsim=20)

## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')

## Power for predictor 'x', (95% confidence interval):
##       30.00% (11.89, 54.28)
## 
## Test: Kenward Roger (package pbkrtest)
##       Effect size for x is -0.10
## 
## Based on 20 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 30
## 
## Time elapsed: 0 h 0 m 1 s

powerSim(model2, nsim=20)

## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')

## Power for predictor 'x', (95% confidence interval):
##       35.00% (15.39, 59.22)
## 
## Test: z-test
##       Effect size for x is -0.10
## 
## Based on 20 simulations, (0 warnings, 0 errors)
## alpha = 0.05, nrow = 30
## 
## Time elapsed: 0 h 0 m 2 s