# Separate variances for different levels of a fixed effect

Is there a way to get a glm to evaluate separate variances for different levels of a fixed effect? I am prepping a model for an experiment involving lizards on six islands. My model attempts to predict lizard height based on an experimental treatment (fixed effect with two levels) with a random intercept term to account for differences between islands. I would like for the model to evaluate whether the variance in lizard height is different between treatment levels (ideally in addition to evaluating shifts in mean lizard height). Below are some simulated data and the model. I recognize that my current simulations do not actually introduce different variances for different treatment levels.

#CHOOSE PARAMETERS FOR FAKE DATA
b.0=200 #whole model intercept - represents perch height on control islands
b.1=-100 #fixed treatment affect
g.0=0 #intercept for group level predictor (island)
sigma.y.true=25 #within group (island) SD
sigma.a.true=25 #between group (island) SD
n=50 #lizards caught per island
nsims=1 #number of simulations

#SIMULATE GROUP LEVEL (ISLAND) INTERCEPTS
sites1=1:6
z=data.frame(Site=sites1, Int=rep(NA,length(sites1)))
for(i in 1:length(z\$Site)){
z\$Int[i]=rnorm(1,g.0,sigma.a.true )
}

#CREATE DATA FROM TO STORE SIMULATED DATA
sites2=as.factor(c(rep(1,n),rep(2,n),rep(3,n),rep(5,n),rep(5,n),rep(6,n)))
treat=c(rep(0,(3n)),rep(1,(3n))) #control=0 and removal=1
IslandInt=rep(NA,nlength(sites1))
height=rep(NA,n
length(sites1))
sims=data.frame(Site=sites2,Treatment=treat,IslandInt=IslandInt,Height=height)

#LOOP IN PREVIOUSLY MODELED INSLAND INTERCEPTS
for(j in 1:(length(sims\$Site))){
sims\$IslandInt[j]=z\$Int[z\$Site==as.character(sims\$Site[j])]
}

#SIMULATE FAKE DATA
for(j in 1:(length(sims\$Site))){
sims\$Height[j]=rnorm(1,sims\$IslandInt[j]+b.0+b.1*sims\$Treatment[j],sigma.y.true)
}

#RUN MODEL
fit = stan_glmer(Height ~ Treatment + (1|Site), data=sims)
print(fit)

That isn’t really supported by `[g]lmer` in the lme4 package, so it isn’t really supported by `stan_[g]lmer`. @bbolker has an example
https://rpubs.com/bbolker/factorvar
using lme4. You could squash the covariance to a near zero value with a sufficiently large value of `regularization` in `prior_covariance = decov(regularization = 10)`.

If I’m understanding your question correctly, you want to see if the residual standard deviation/variance is different between the two groups? You should be able to fit these sorts of models with the `brms` package (which also uses Stan as a backend). Take a look at the first example in this vignette.
change sigma’s definition to `vector[n_groups] sigma` and use the
likelihood `y ~ normal(eta, sigma[group_id])`, where eta is replaced with