random effect variance increases with the inclusion of predictor What does it mean when the inclusion of a fixed effect (sex) increases variance of the random effect?
dat <- data.frame(ID=rep(1:12), sex=c(0,1,1,1,0,0,1,0,1,1,0,1), region=rep(c("A","B"),each=6), y=c(12,13,15,7,8,4,8,8,9,3,4,5))

m <- lmer(y~1+(1|region),data=dat)
VarCorr(m)
 Groups   Name        Std.Dev.
 region   (Intercept) 2.1820    

m = lmer(y~sex+(1|region),data=dat)
VarCorr(m)
 Groups   Name        Std.Dev.
 region   (Intercept) 2.4544 

I thought including an independent variable could only decrease random effect variance, but I am obviously wrong. Is it a good idea to include sex in the model even though it increases random effect variance?
 A: In this case, the expected behaviour is that the variance of the random intercepts will be unchanged after adding the predictor. In reality, it will change, but the expected change is still zero. If we modify your simulation slightly, we can see that it can decrease slightly:
> dat$region <- rep(c("A","B","C"),each=4)

> lmer(y~1+(1|region),data=dat) %>% VarCorr()
## Groups   Name        Std.Dev.
## region   (Intercept) 3.0721  
## Residual             2.7386  

> lmer(y~sex+(1|region),data=dat) %>% VarCorr()
## Groups   Name        Std.Dev.
## region   (Intercept) 3.0236  
## Residual             2.8577 

If we were to repeat this, with different data, we would find that, on average the difference will be zero.  We can do a simple simulation to show this:
> N <- 100
> n.sim <- 200
> rints_dif <- as.numeric(n.sim)
> for(i in 1:n.sim) {
+   set.seed(i)
+ 
+   X <- rnorm(N) 
+   G <- rbinom(N, 10, 0.5)
+   Y <- G + rnorm(N)
+ 
+   m0 <- lmer(Y ~ 1 + (1|G))
+ 
+   m1 <- lmer(Y ~ 1 + X + (1|G))
+ 
+   rints_dif[i] <- as.data.frame(VarCorr(m0))$sdcor[1] - as.data.frame(VarCorr(m1))$sdcor[1]
+ }
> hist(rints_dif)
> mean(rints_dif)
[1] -0.0002529647


