代码之家 › 专栏 › 技术社区 › Linus

在R中使用SIMR绘制每个簇具有不同观测次数的功率曲线时出错

power-analysis mixed-models r

1

Linus · 技术社区 · 11 月前

我在用 simr R中的包,用于对两级泊松回归模型进行幂分析。我试图绘制一条功率曲线,以了解每个集群的观测次数的影响。然而,我遇到了一些问题 powerCurve() 功能。

我正试图用以下公式绘制功率曲线

power_curve2 <- powerCurve(model, within="x + Z ", breaks=c(20,50,80,100), nsim = 10)
plot(power_curve2)

理论上,这应该绘制一条具有不同簇数值的功率曲线(由 breaks ). 但是我得到以下错误:

plot.window(…)错误:需要有限的“xlim”值此外:警告信息: 1:在min(x)中:min没有不缺失的参数;返回Inf 2:在max(x)中:max没有不缺失的参数;返回-Inf

你知道如何绘制功率曲线吗?

这是设置电源模拟模型的代码:

#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

1 回复 | 直到 11 月前

1

qdread 11 月前

我认为答案就在文件中 ?powerCurve 那里写着 breaks 参数仅在以下情况下使用 along 被指定为参数,而不是 within 。这应该会给你提供所需的功率曲线,作为集群数量的函数( g 是表示集群身份的变量)。

power_curve2 <- powerCurve(model,  along = 'g', breaks=c(5,10,15,20), nsim = 10)
plot(power_curve2)