# Puzzling results when modeling heterogeneous variance in the residuals

I’m trying to capture different variances in the residuals through `brms`. Here is the code:

``````require(MASS); library(brms)
ns <- 200
S  <- matrix(c(  1, 0.5,    0,    0,
0.5,   1,    0,    0,
0,   0,    1, 0.25,
0,   0, 0.25,    1), nrow=4, ncol=4)
dat <- mvrnorm(n=ns, mu=rep(0,4), Sigma=S)
dat <- data.frame(f = c(rep(paste0('P',1:ns), 2), rep(paste0('T',1:ns), 2)),
R=c(rep('P',2*ns), rep('T',2*ns)),
y = c(dat))
dat\$P <- as.numeric(dat\$R=='P'); dat\$T <- as.numeric(dat\$R=='T')
options(mc.cores = parallel::detectCores())
fm <- brm(bf(y ~ 1+(0+P|f)+(0+T|f), sigma~0+R), data=dat,
chains = 4, iter=1000))
``````

Note that the heterogeneous variance is modeled through `sigma~0+R` inside the formula `bf`. With the result from the above code, I have the following summary:

``````summary(fm)
``````
``````Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    -0.00      0.04    -0.08     0.07 1.00     1461     1297
sigma_RP     -0.41      0.05    -0.50    -0.31 1.00      934     1318
sigma_RT     -0.13      0.05    -0.23    -0.03 1.00      626      633
``````

With 2 group levels (`P` and `T`) in the factor `R`, we see 2 sigmas of `sigma_RP` and `sigma_RT` in the last two rows in the table above. My question is: why are the 2 sigma values negative?

a very quick answer is that the sigma parameter has a log link. So the sigma estimates you are getting here are on a log scale. If you exponentiate them, you should get values comparable to those you used to generate your data.

1 Like

Thanks for the pointer, @fusaroli ! Now everything goes through the square hole.

1 Like