I'm using the simr package in R to conduct a power analysis for a two-level poisson regression model. I'm trying to plot a power curve to understand the influence of the number of observations per cluster. However, I'm encountering issues with the powerCurve() function.
I'm trying to plot the power curve using
power_curve2 <- powerCurve(model, within="x + Z ", breaks=c(20,50,80,100), nsim = 10)
plot(power_curve2)
This should theoretically plot a power curve with different values for the number of clusters (defined by the breaks). However I get the following Error:
Error in plot.window(...) : need finite 'xlim' values In addition: Warning messages: 1: In min(x) : no non-missing arguments to min; returning Inf 2: In max(x) : no non-missing arguments to max; returning -Inf
Any idea, on how I can plot the power curve?
This is the code setting up the model for the power simulation:
#install.packages("simr")
library(simr)
#### Specification of Input Parameters ####
### Specification of standardized input parameters ###
L1_DE_standardized <- .1 ## standardized L1 direct effect
L2_DE_standardized <- .3 ## standardized L2 direct effect
CLI_E_standardized = .3 ## standardized CLI effect
rand.sl = .09 ## standardized random slope variance component
ICC <- .30 ## standardized intraclass correlation coefficient
cor.i.sl = .00 ## Correlation between random slope and random intercept
alpha.S <- .05 ## significance level
Size.clus <- 100 ## L1 sample size (cluster size)
N.clus <- 30 ## L2 sample size (number of clusters)
### Derivation of a population model for the power analysis ###
## Specification of predictor variables ##
x <- scale(rep(1:Size.clus))
g <- as.factor(1:N.clus)
X <- cbind(expand.grid("x"=x, "g"=g))
X <- data.frame(X, Z=as.numeric(X$g))
X$Z <- scale(X$Z)
## Specification of the outcome variable ##
varL1 <- 1 ## uncond. L1 variance (fixed at 1)
s <- sqrt((varL1)*(1-(L1_DE_standardized^2))) ## cond. L1 variance (Equation 11)
varL2 <- ICC/(1-ICC) ## uncond. L2 variance (Equation 10)
V1 <- varL2*(1-(L2_DE_standardized^2)) ## cond. L2 variance (Equation 11)
## Adjustment of the random slope variance ##
varRS <- rand.sl*varL1 ## uncond. random slope variance (Equation 13)
condRS <- varRS*(1-(CLI_E_standardized^2)) ## cond. Random slope variance (Equation 11)
## Covariance of random intercept and slope (Equation 14) ##
cov.i.sl <- cor.i.sl*sqrt(varL2)*sqrt(varRS)
## Adjustment of fixed effects (Equation 15) ##
L1_DE <- L1_DE_standardized*sqrt(varL1)
L2_DE <- L2_DE_standardized*sqrt(varL2)
CLI_E <- CLI_E_standardized*sqrt(varRS)
#### Implementation of a Power Analysis in a Two-Level Model in SIMR ####
### Vector of fixed effects ###
b <- c(0, L1_DE, L2_DE, CLI_E)
rand_sl.con <- varRS*(1-(CLI_E_standardized^2))
### Random intercept/slope variance-covariance matrix ###
V2 <- matrix(c(V1, cov.i.sl, cov.i.sl, rand_sl.con), 2)
model <- makeGlmer(y ~ x + Z + x:Z + (x|g),
fixef=b, VarCorr=V2,
family="poisson", data=X)
print(model)
### Simulating power for the L1 direct effect ###
sim.ef <- powerSim(model,
fixed("x"),
alpha=alpha.S,
nsim=10)
print(sim.ef)
simdat_Table2 <- cbind(effect="x",
Size.clus,
N.clus,
L1_DE_standardized,
L2_DE_standardized,
CLI_E_standardized,
summary(sim.ef))
power_curve2 <- powerCurve(model, within="x + Z ", breaks=c(20,50,80,100), nsim = 10)
plot(power_curve2)
# trying to plot power curve
I believe the answer lies in the documentation for ?powerCurve. There it says that the breaks argument is only used when along is specified as an argument instead of within. This should give you the desired power curve, as a function of the number of clusters (g is the variable representing cluster identity).
power_curve2 <- powerCurve(model, along = 'g', breaks=c(5,10,15,20), nsim = 10)
plot(power_curve2)