攔截的隨機效應方差爲零

我正在R中使用正常的LMM進行功耗分析。我有七個輸入參數，其中兩個我不需要測試（年數和站點數量）。其他5個參數是殘差，截距和斜率的截距，斜率和隨機效應標準差。由於我的回答數據（年份是模型中唯一的解釋變量）約束在（-1，+1）之間，所以截距也落在此範圍內。然而，我發現如果我運行了1000次模擬，並給定了一個給定的截距和斜率（我認爲這是10年以來的常數），那麼如果隨機效應截距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)

非常感謝您提供任何幫助。

來源

2014-10-31 user3628889