Puzzling results when modeling heterogeneous variance in the residuals

Modeling
,
1

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?

2

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
3

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

1 Like