R で通常の LMM を使用して電力分析を実行しています。7 つの入力パラメーターがあり、そのうちの 2 つ (年数とサイト数) をテストする必要はありません。他の 5 つのパラメーターは、切片、勾配、および残差、切片、勾配の変量効果標準偏差です。
私の応答データ (モデルの唯一の説明変数は年) が (-1, +1) の間でバインドされていることを考えると、切片もこの範囲に収まります。ただし、私が見つけたのは、たとえば、特定の切片と勾配 (10 年間で一定として扱っている) で 1000 回のシミュレーションを実行すると、ランダム効果の切片 SD が特定の値を下回ると、多くのランダム効果が SD をインターセプトするシミュレーションはゼロです。インターセプト SD を膨らませると、これは正しくシミュレートされるように見えます (残差 Sd=0.25、インターセプト SD = 0.10、勾配 SD = 0.05 を使用する以下を参照してください。インターセプト SD を 0.2 に増やすと、これは正しくシミュレートされます。または、残差 SD を 0.05 に落とすと、分散パラメーターが正しくシミュレートされます)。
この問題は、範囲を強制的に (-1, +1) にしたことが原因ですか?
これが役立つ場合は、関数のコードと以下のシミュレーションの処理を含めます。
機能: データの生成:
normaldata <- function (J, K, beta0, beta1, sigma_resid,
sigma_beta0, sigma_beta1){
year <- rep(rep(0:J),K) # 0:J replicated K times
site <- rep (1:K, each=(J+1)) # 1:K sites, repeated J years
mu.beta0_true <- beta0
mu.beta1_true <- beta1
# random effects variance parameters:
sigma_resid_true <- sigma_resid
sigma_beta0_true <- sigma_beta0
sigma_beta1_true <- sigma_beta1
# site-level parameters:
beta0_true <<- rnorm(K, mu.beta0_true, sigma_beta0_true)
beta1_true <<- rnorm(K, mu.beta1_true, sigma_beta1_true)
# data
y <<- rnorm(n = (J+1)*K, mean = beta0_true[site] + beta1_true[site]*(year),
sd = sigma_resid_true)
# NOT SURE WHETHER TO IMPOSE THE LIMITS HERE OR LATER IN CODE:
y[y < -1] <- -1 # Absolute minimum
y[y > 1] <- 1 # Absolute maximum
return(data.frame(y, year, site))
}
シミュレートされたコードの処理:
vc1 <- as.data.frame(VarCorr(lme.power))
vc2 <- as.numeric(attributes(VarCorr(lme.power)$site)$stddev)
n.sims = 1000
sigma.resid <- rep(0, n.sims)
sigma.intercept <- rep(0, n.sims)
sigma.slope <- rep(0,n.sims)
intercept <- rep(0,n.sims)
slope <- rep(0,n.sims)
signif <- rep(0,n.sims)
for (s in 1:n.sims){
y.data <- normaldata(10,200, 0.30, ((0-0.30)/10), 0.25, 0.1, 0.05)
lme.power <- lmer(y ~ year + (1+year | site), data=y.data)
summary(lme.power)
theta.hat <- fixef(lme.power)[["year"]]
theta.se <- se.fixef(lme.power)[["year"]]
signif[s] <- ((theta.hat + 1.96*theta.se) < 0) |
((theta.hat - 1.96*theta.se) > 0) # returns TRUE or FALSE
signif[s]
betas <- fixef(lme.power)
intercept[s] <- betas[1]
slope[s] <- betas[2]
vc1 <- as.data.frame(VarCorr(lme.power))
vc2 <- as.numeric(attributes(VarCorr(lme.power)$site)$stddev)
sigma.resid[s] <- vc1[4,5]
sigma.intercept[s] <- vc2[1]
sigma.slope[s] <- vc2[2]
cat(paste(s, " ")); flush.console()
}
power <- mean (signif) # proportion of TRUE
power
summary(sigma.resid)
summary(sigma.intercept)
summary(sigma.slope)
summary(intercept)
summary(slope)
あなたが提供できる助けを前もって感謝します。