Power analysis from scratch

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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.

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')
## 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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.