Separate variances for different levels of a fixed effect


#1

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)


Varying Variances
#2

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


#3

This link is helpful, thanks. Is it possible to do what I’m asking in STAN?


#4

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.


#5

Thanks very much! I might actually like to try coding the model myself in Stan. Is anyone aware of good sample code for this?


#6

You should be able to start with a standard random intercept model, then
change sigma’s definition to vector[n_groups] sigma and use the
likelihood y ~ normal(eta, sigma[group_id]), where eta is replaced with
your linear predictor and group_id is an integer array that indexes which
group each lizard belongs to.


#7

You can also learn writing Stan code with brms for such models by using its make_stancode function.