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))length(sites1))

height=rep(NA,n

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)