The variance of random interception effects is zero

I am running a power analysis using normal LMM in R. I have seven input parameters, two of which I don't need to test (number of years and number of sites). The other 5 parameters are intercept, tilt and standard deviation of the random effects of residual, intercept and tilt.

Given that my response data (year is the only explanatory variable in the model) is related between (-1, +1), interception also falls into this range. However, what I find is that if I run, say, 1000 simulations with a given intercept and slope (which I consider to be constant for over 10 years), then if the random intercept SD falls below a certain value, the simulations, where interception of random effects SD is zero. If I inflate the SD intercept then it seems to simulate correctly (see below where I use the residual Sd = 0.25, SD intercept = 0.10 and SD slope = 0.05 if I increase the SD intercept to 0. 2, it's modeled correctly, or if I drop the residual SD to say 0.05, the variance parameters are correctly modeled).

Is this problem due to my coercion of the range (-1, +1)?

I am including the code for my function and the simulations handling below in case it helps:

Function: data creation:

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))
}

      

        
        
        
      

    

Simulated code processing:

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)

      

        
        
        
      

    

Thanks in advance for any help you can offer.