Simply extracting random effects variances

[strengejacke]

Close #717

[strengejacke]

ok, two things:

1. looks like there are some differences in rounding between the old and the new approach to extract random effects variances.
2. looks like the order of parameters changed in a few places. this means, either fixing this in the code, or - if the new sorting is OK - fix the related tests.

[mattansb]

I don't think the order of terms should matter too much...

[strengejacke]

I fixed the ordering-issue now, and actually improved it. Order might matter, because in your poly-examples, parameters were ordered alphabetically (.C, .L, .Q), which should be .L, .Q and .C (and then .4 and .5). This was previously a wrong behaviour and should now also be fixed.

This PR still needs some testing, especially for profiled-CIs, to see if parameters and SE/CI still match. Feel free to play around with some models using model_parameters(model, ci_method = "profile") and see if the CIs still match those from ci_method = "wald" (they're not identical, but I think we can see if the CIs are "matching" and in the correct order).

[strengejacke]

And we need to check if this also work for other packages, like GLMMadaptive, ordinal:clmm, nlme etc.

[strengejacke]

@mattansb @DominiqueMakowski @bwiernik

If you look at

library(rstanarm)
library(lme4)
data(sleepstudy)
data(cake)
set.seed(123)

m1 <- suppressMessages(stan_lmer(angle ~ temperature + (temperature | recipe) + (temperature | replicate), data = cake))
model_parameters(m1, effects = "all")

you see that correlations are not bounded between [-1, 1]. This is something to be fixed in bayestestR, I think, but do you have any idea how to transform the posterior in order to have a proper "point estimate" here?

[mattansb]

you see that correlations are not bounded between [-1, 1]. This is something to be fixed in bayestestR, I think, but do you have any idea how to transform the posterior in order to have a proper "point estimate" here?

I think the values from m1 are covariances and variances, not correlations and standard deviation. So you would have to:

library(rstanarm)
library(lme4)
data(cake)
set.seed(123)

m1 <- stan_lmer(angle ~ temp + (temp | replicate),
                data = cake,
                chains = 1)

bayestestR::describe_posterior(m1, effects = "all")
#> ...
#> 
#> # Random effects SD/Cor: replicate
#> 
#> Parameter          |   Median |         95% CI |     pd |          ROPE | % in ROPE |  Rhat |    ESS
#> ----------------------------------------------------------------------------------------------------
#> (Intercept)        | 2.90e-03 | [ 0.00,  0.70] |   100% | [-0.82, 0.82] |      100% | 0.999 | 612.00
#> temp ~ (Intercept) | 1.86e-04 | [-0.01,  0.01] | 59.20% | [-0.82, 0.82] |      100% | 0.999 | 395.00
#> temp               | 1.12e-03 | [ 0.00,  0.00] |   100% | [-0.82, 0.82] |      100% | 1.000 | 325.00
#> 
#> ...


posteriors <- as.data.frame(m1)

rho <- with(posteriors, {
  `Sigma[replicate:temp,(Intercept)]` / 
    (sqrt(`Sigma[replicate:temp,temp]`) * sqrt(`Sigma[replicate:(Intercept),(Intercept)]`))
})

# Cor
bayestestR::describe_posterior(rho)
#> Summary of Posterior Distribution
#> 
#> Parameter | Median |        95% CI |     pd |          ROPE | % in ROPE
#> -----------------------------------------------------------------------
#> Posterior |   0.19 | [-0.92, 0.95] | 59.20% | [-0.10, 0.10] |     9.05%


# Slopes
a <- with(posteriors, sqrt(`Sigma[replicate:temp,temp]`))

bayestestR::describe_posterior(a)
#> Summary of Posterior Distribution
#> 
#> Parameter | Median |       95% CI |   pd |          ROPE | % in ROPE
#> --------------------------------------------------------------------
#> Posterior |   0.03 | [0.02, 0.05] | 100% | [-0.10, 0.10] |      100%


# Compare to:
lme4::lmer(angle ~ temp + (temp | replicate),
           data = cake) |> 
  lme4::VarCorr()
#>  Groups    Name        Std.Dev. Corr 
#>  replicate (Intercept) 2.604527      
#>            temp        0.027073 0.134
#>  Residual              4.825927

[mattansb]

Is this maybe an issue in insight?

[strengejacke]

Is this maybe an issue in insight?

I'm not sure, because - if I recall right - we just summarize the posteriors using mean/median etc., no matter whether fixed or random. insight just returns (roughly speaking) as.data.frame(model).

[bwiernik]

How does VarCorr() line up the appropriate parameters?

[strengejacke]

How does VarCorr() line up the appropriate parameters?

What do you mean? How does it extract all parameters, or sorting?

[bwiernik]

Re:

I'm not sure about var1/var2, when we have only one posterior of correlation values (on the "wrong" scale). We also need to grep column names to extract random effects... It's a slightly annoying issue, I must admit. ;-)

VarCorr() lines up the variance and covariance parameters and uses these to compute correlations. Can we copy that logic?

[strengejacke]

@IndrajeetPatil Do you know how I correctly include https://github.com/glmmTMB/glmmTMB/tree/ci_tweaks/glmmTMB in the "Remotes" section of the DESCRIPTION file?

[IndrajeetPatil]

Sorry, I couldn't figure this one out.

The problem is that you need to specify a subdirectory on a given branch.

remotes::install_github("glmmTMB/glmmTMB@ci_tweaks", subdir = "glmmTMB")

I am not sure what URL this resolves to, and that's exactly what you will need to include in Remotes.

[IndrajeetPatil]

Let's see if StackOverflow can help: https://stackoverflow.com/questions/72671996/how-to-include-subdirectory-in-remotes-field-for-r-package-description

[codecov-commenter]

Codecov Report

Merging #720 (2b3b353) into main (b1eaa53) will decrease coverage by 0.32%.
The diff coverage is 60.81%.

@@            Coverage Diff             @@
##             main     #720      +/-   ##
==========================================
- Coverage   53.23%   52.91%   -0.33%     
==========================================
  Files         183      183              
  Lines       12104    12027      -77     
==========================================
- Hits         6444     6364      -80     
- Misses       5660     5663       +3     

Impacted Files	Coverage Δ
R/1_model_parameters.R	87.96% <ø> (ø)
R/utils_format.R	56.64% <0.00%> (-0.40%)	⬇️
R/utils_model_parameters.R	91.53% <ø> (ø)
R/extract_random_variances.R	65.50% <60.49%> (-6.01%)	⬇️
R/extract_parameters.R	77.30% <100.00%> (+0.23%)	⬆️
R/format.R	56.67% <100.00%> (+0.28%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update b1eaa53...2b3b353. Read the comment docs.

[IndrajeetPatil]

That worked!

[strengejacke]

Nice, thanks! :-)