Correlated error solution

[jebyrnes]

While if you want correlated errors between all endogenous variables, you can set_rescor(TRUE), if it's only a subset, I included a solution from reading several sources at https://jebyrnes.github.io/bayesian_sem/bayesian_sem.html#correlated-error

The example:

library(palmerpenguins)
library(brms)

cor_err_bf <-
  bf(bill_depth_mm ~ species + body_mass_g +
       (1|dl|obs), sigma = 1) +
  bf(bill_length_mm ~ species + body_mass_g +
       (1|dl|obs), sigma = 1) +
  set_rescor(FALSE)

corr_err_brm <- brm(cor_err_bf,
                    data = penguins)

sigma = 1 makes it work, but, is it correct? In comparing it to set_rescor(TRUE), the correlation or errors is different. But, that might be due to the sigma= 1?

cor_err_rescor_bf <-
  bf(bill_depth_mm ~ species + body_mass_g) +
  bf(bill_length_mm ~ species + body_mass_g) +
  set_rescor(TRUE)

corr_err_rescor_brm <- brm(cor_err_bf,
                    data = penguins)

Setting sigma = 0 in the top model causes a failure of convergence. Setting it very small leads to some oddities as well.

Feels like there's something here that needs to be fixed and/or at least get the proper solution to.

[jebyrnes]

paul-buerkner/brms#254 seems like the relevant issue on the brms issue tracker

[jebyrnes]

For reference

> tidy(corr_err_brm)
# A tibble: 13 × 9
   response     effect   component group         term  estimate std.error conf.low conf.high
   <chr>        <chr>    <chr>     <chr>         <chr>    <dbl>     <dbl>    <dbl>     <dbl>
 1 billdepthmm  fixed    cond      NA            (Int… 12.8      0.449    11.9      13.7    
 2 billlengthmm fixed    cond      NA            (Int… 24.9      1.06     22.8      27.0    
 3 billdepthmm  fixed    cond      NA            spec…  0.0249   0.152    -0.281     0.321  
 4 billdepthmm  fixed    cond      NA            spec… -5.44     0.204    -5.84     -5.04   
 5 billdepthmm  fixed    cond      NA            body…  0.00151  0.000119  0.00127   0.00174
 6 billlengthmm fixed    cond      NA            spec…  9.91     0.347     9.23     10.6    
 7 billlengthmm fixed    cond      NA            spec…  3.56     0.488     2.61      4.50   
 8 billlengthmm fixed    cond      NA            body…  0.00375  0.000283  0.00319   0.00430
 9 obs          ran_pars cond      billdepthmm_… sd__…  0.242    0.0648    0.117     0.370  
10 obs          ran_pars cond      billlengthmm… sd__…  2.19     0.104     1.99      2.40   
11 obs          ran_pars cond      billdepthmm_… cor_…  0.860    0.119     0.560     0.995  
12 billdepthmm  ran_pars cond      Residual      sd__…  1        0         1         1      
13 billlengthmm ran_pars cond      Residual      sd__…  1        0         1         1      

and

> tidy(corr_err_rescor_brm)
# A tibble: 10 × 9
   response     effect   component group    term       estimate std.error conf.low conf.high
   <chr>        <chr>    <chr>     <chr>    <chr>         <dbl>     <dbl>    <dbl>     <dbl>
 1 billdepthmm  fixed    cond      NA       (Intercep… 12.8      0.391    12.0      13.5    
 2 billlengthmm fixed    cond      NA       (Intercep… 24.9      1.07     22.8      26.9    
 3 billdepthmm  fixed    cond      NA       speciesCh…  0.0259   0.129    -0.227     0.279  
 4 billdepthmm  fixed    cond      NA       speciesGe… -5.44     0.181    -5.78     -5.08   
 5 billdepthmm  fixed    cond      NA       body_mass…  0.00151  0.000104  0.00130   0.00171
 6 billlengthmm fixed    cond      NA       speciesCh…  9.92     0.367     9.21     10.6    
 7 billlengthmm fixed    cond      NA       speciesGe…  3.56     0.495     2.58      4.50   
 8 billlengthmm fixed    cond      NA       body_mass…  0.00375  0.000286  0.00321   0.00432
 9 billdepthmm  ran_pars cond      Residual sd__Obser…  0.883    0.0338    0.821     0.953  
10 billlengthmm ran_pars cond      Residual sd__Obser…  2.42     0.0941    2.24      2.61 

Ach - doesn't show the rescor correlation in the second, so, here's the comparison

OK, from the corr_err model

                                                  Estimate Est.Error l-95% CI u-95% CI Rhat
sd(billdepthmm_Intercept)                             0.24      0.06     0.12     0.37 1.00
sd(billlengthmm_Intercept)                            2.19      0.10     1.99     2.40 1.00
cor(billdepthmm_Intercept,billlengthmm_Intercept)     0.86      0.12     0.56     0.99 1.03
                                                  Bulk_ESS Tail_ESS
sd(billdepthmm_Intercept)                             3611     2463
sd(billlengthmm_Intercept)                            2279     2897
cor(billdepthmm_Intercept,billlengthmm_Intercept)      120      186

Further Distributional Parameters:
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_billdepthmm      1.00      0.00     1.00     1.00   NA       NA       NA
sigma_billlengthmm     1.00      0.00     1.00     1.00   NA       NA       NA

And from the rescor model

Further Distributional Parameters:
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_billdepthmm      0.88      0.03     0.82     0.95 1.00     3920     2939
sigma_billlengthmm     2.42      0.09     2.24     2.61 1.00     4083     2358

Residual Correlations: 
                                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
rescor(billdepthmm,billlengthmm)     0.26      0.05     0.16     0.35 1.00     3789
                                 Tail_ESS
rescor(billdepthmm,billlengthmm)     2867

[jebyrnes]

So, inference about relationships between variables is pretty much unchanged. But, the error structure is..... different. Hrm....

[amynang]

Does fixing sigma come from this post? (ctrl+F "If you wanna get silly" for the relevant part) In that case sigma ~ 0 makes sense because he's working with everything standardized to unit variance.

I think your (1|dl|obs) plus set_rescor(FALSE) cannot be equivalent with set_rescor(TRUE) because in the former case you are fitting an OLRE model.

[jebyrnes]

Might be, but that makes sense. I see that the models are not quite equivalent. I guess the question is, is it possible to make them so? Although the coefficients are the same, which makes me think the relevant problem is being handled. But, yeah, I'm trying to hack my way to arbitrary correlated errors in different parts of a model here - a generalized solution rather than relying on set_rescor().

[jebyrnes]

Of note

cor_err_zero_bf <-
    bf(bill_depth_mm ~ species + body_mass_g +
           (1|dl|obs), sigma = ~0) +
    bf(bill_length_mm ~ species + body_mass_g +
           (1|dl|obs), sigma = ~0) +
    set_rescor(FALSE)

gives the same results as sigma = 1.

In looking at the models, the more I think about this, I think the big difference is that in one model, the 'group' - really observation - level SDs are correlated, but then there is some additional residual. In the second, the residuals at the observation level are correlated, with no additional residual.

If so, it would lead to the correlation in the sigma=1 model to be greater in magnitude than that of the one with just a correlation matrix (as some of the uncorrelated sigma is pulled away). But, in both cases, the beta should be correctly estimated - which is the goal.

Looking at the stan code I think backs this up, although I'm not super-fluent in stan.

The big issue is that you cannot set sigma=0 as the likelihood is a normal_id_glm_lpdf with sigma in it - which cannot be 0.

So, this is a hack to be sure...... but one that should work? Although the error/residuals reported will be wonky. I'm AFK largely, but I wonder if you could posthoc recover the correct error structure.

[jebyrnes]

So, does the hack work? As the model suggested by the hack should produce both response-level residuals and 'obs' level residuals, theoretically, summing the residuals from both and then recalculating sigma and the correlation should produce the same results as obtained from rescor=TRUE

Using corr_err_brm we get the following using tidybayes, tidyr, and dplyr to make things understandable.

b <- residual_draws(corr_err_brm, 
                    newdata = penguins |> 
                      filter(!is.na(bill_length_mm)),
                    re_formula = ~obs)
b1 <- b |>
  pivot_wider(names_from = .category, 
              values_from = .residual) |>
  group_by(.draw) |>
  summarise(sigma_billdepthmm = sd(billdepthmm),
            sigma_billlengthmm = sd(billlengthmm),
            rescor = cor(billlengthmm, billdepthmm))

b1 |>  select(-.draw) |> colMeans()

This yields

 sigma_billdepthmm sigma_billlengthmm             rescor 
         1.3304900          2.6023616          0.1573824 

Which is a bit off

tidy_draws(corr_err_rescor_brm) |>
  select(sigma_billdepthmm,
         sigma_billlengthmm,
         rescor__billdepthmm__billlengthmm) |>
  colMeans()

                sigma_billdepthmm                sigma_billlengthmm 
                        0.8826025                         2.4169913 
rescor__billdepthmm__billlengthmm 
                        0.2546789 

Hrm.......

Closer? But why the difference?

[jebyrnes]

Huh. Making sigma smaller - 0.5 - in the initial model makes them nearly identical.

with sigma = 0.1. Here are some quantiles of the posteriors.

# A tibble: 3 × 4
      q sigma_billdepthmm sigma_billlengthmm rescor
  <dbl>             <dbl>              <dbl>  <dbl>
1  0.25             0.880               2.40  0.251
2  0.5              0.884               2.41  0.256
3  0.75             0.889               2.42  0.261

the set_rescor model

# A tibble: 3 × 4
      q sigma_billdepthmm sigma_billlengthmm
  <dbl>             <dbl>              <dbl>
1  0.25             0.859               2.35
2  0.5              0.881               2.41
3  0.75             0.905               2.48

OK - hack works? But why with 0.5 and not 1...... must be something to do with the scale....