How can I include random effects from a brms car model in the model_parameters output?

[mcewenkhundi]

Apologies if this is not the right place for this question.

library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.15.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').
#> 
#> Attaching package: 'brms'
#> The following object is masked from 'package:stats':
#> 
#>     ar
library(parameters)


  # generate some spatial data
  east <- north <- 1:10
  Grid <- expand.grid(east, north)
  K <- nrow(Grid)

  # set up distance and neighbourhood matrices
  distance <- as.matrix(dist(Grid))
  W <- array(0, c(K, K))
  W[distance == 1] <- 1

  # generate the covariates and response data
  x1 <- rnorm(K)
  x2 <- rnorm(K)
  theta <- rnorm(K, sd = 0.05)
  phi <- rmulti_normal(
    1, mu = rep(0, K), Sigma = 0.4 * exp(-0.1 * distance)
  )
  eta <- x1 + x2 + phi
  prob <- exp(eta) / (1 + exp(eta))
  size <- rep(50, K)
  y <- rbinom(n = K, size = size, prob = prob)
  dat <- data.frame(y, size, x1, x2)

  # fit a CAR model
  fit <- brm(y | trials(size) ~ x1 + x2 + car(W),
             data = dat, data2 = list(W = W),
             family = binomial())
#> Warning: Using CAR terms without a grouping factor is deprecated. Please use
#> argument 'gr' even if each observation represents its own location.
#> Compiling Stan program...
#> Start sampling
#> 
#> SAMPLING FOR MODEL 'ac2c262b3bd98cbce581999034ed54e0' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 2.415 seconds (Warm-up)
#> Chain 1:                1.721 seconds (Sampling)
#> Chain 1:                4.136 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'ac2c262b3bd98cbce581999034ed54e0' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 0 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 2.661 seconds (Warm-up)
#> Chain 2:                3.881 seconds (Sampling)
#> Chain 2:                6.542 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'ac2c262b3bd98cbce581999034ed54e0' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 0 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 2.605 seconds (Warm-up)
#> Chain 3:                1.73 seconds (Sampling)
#> Chain 3:                4.335 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'ac2c262b3bd98cbce581999034ed54e0' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 0 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 2.308 seconds (Warm-up)
#> Chain 4:                3.121 seconds (Sampling)
#> Chain 4:                5.429 seconds (Total)
#> Chain 4:
#> Warning: The largest R-hat is 1.06, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#tail-ess

  summary(fit)
#> Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
#> careful when analysing the results! We recommend running more iterations and/or
#> setting stronger priors.
#>  Family: binomial 
#>   Links: mu = logit 
#> Formula: y | trials(size) ~ x1 + x2 + car(W) 
#>    Data: dat (Number of observations: 100) 
#> Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup samples = 4000
#> 
#> Correlation Structures:
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> car       0.97      0.03     0.89     1.00 1.01      159      255
#> sdcar     0.60      0.08     0.46     0.77 1.01      977     2034
#> 
#> Population-Level Effects: 
#>           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept    -0.32      0.38    -1.35     0.55 1.06       46       32
#> x1            0.97      0.06     0.86     1.08 1.00     2124     2499
#> x2            1.02      0.06     0.90     1.13 1.00     2338     2906
#> 
#> Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

  #How can I get the random effects to be part of the results of the mode_parameters ?

  model_parameters(fit, centrality = "mean", effects = "all", component = "all", ci = 0.95, exponentiate = TRUE)
#> Warning: Following potential variables could not be found in the data: W
#> # Fixed effects
#> 
#> Parameter   | Mean |       95% CI |     pd | % in ROPE |  Rhat |     ESS
#> ------------------------------------------------------------------------
#> (Intercept) | 0.72 | [0.37, 1.89] | 87.45% |    23.05% | 1.075 |   35.00
#> x1          | 2.64 | [2.37, 2.93] |   100% |        0% | 1.000 | 2117.00
#> x2          | 2.76 | [2.47, 3.10] |   100% |        0% | 1.000 | 2315.00
#> 
#> Using highest density intervals as credible intervals.

^{Created on 2021-08-14 by the reprex package (v2.0.0)}

[mattansb]

Hmm... This seems to be an issue to with {insight} not returning the cor parameters.

> insight::get_parameters(fit, effects = "all") |>
+   head()
  b_Intercept    b_x1     b_x2
1    0.643152 1.08399 0.973802
2    0.619776 1.01089 1.027060
3    0.712314 1.09694 1.033150
4    0.659402 1.08305 1.040860
5    0.680401 1.11275 1.038320
6    0.597129 1.02805 0.991711
> as.data.frame(fit) |> colnames()
  [1] "b_Intercept" "b_x1"        "b_x2"        "car"         "sdcar"       "Intercept"   "rcar[1]"     "rcar[2]"     "rcar[3]"    
 [10] "rcar[4]"     "rcar[5]"     "rcar[6]"     "rcar[7]"     "rcar[8]"     "rcar[9]"     "rcar[10]"    "rcar[11]"    "rcar[12]"   
 [19] "rcar[13]"    "rcar[14]"    "rcar[15]"    "rcar[16]"    "rcar[17]"    "rcar[18]"    "rcar[19]"    "rcar[20]"    "rcar[21]"   
 [28] "rcar[22]"    "rcar[23]"    "rcar[24]"    "rcar[25]"    "rcar[26]"    "rcar[27]"    "rcar[28]"    "rcar[29]"    "rcar[30]"   
 [37] "rcar[31]"    "rcar[32]"    "rcar[33]"    "rcar[34]"    "rcar[35]"    "rcar[36]"    "rcar[37]"    "rcar[38]"    "rcar[39]"   
 [46] "rcar[40]"    "rcar[41]"    "rcar[42]"    "rcar[43]"    "rcar[44]"    "rcar[45]"    "rcar[46]"    "rcar[47]"    "rcar[48]"   
 [55] "rcar[49]"    "rcar[50]"    "rcar[51]"    "rcar[52]"    "rcar[53]"    "rcar[54]"    "rcar[55]"    "rcar[56]"    "rcar[57]"   
 [64] "rcar[58]"    "rcar[59]"    "rcar[60]"    "rcar[61]"    "rcar[62]"    "rcar[63]"    "rcar[64]"    "rcar[65]"    "rcar[66]"   
 [73] "rcar[67]"    "rcar[68]"    "rcar[69]"    "rcar[70]"    "rcar[71]"    "rcar[72]"    "rcar[73]"    "rcar[74]"    "rcar[75]"   
 [82] "rcar[76]"    "rcar[77]"    "rcar[78]"    "rcar[79]"    "rcar[80]"    "rcar[81]"    "rcar[82]"    "rcar[83]"    "rcar[84]"   
 [91] "rcar[85]"    "rcar[86]"    "rcar[87]"    "rcar[88]"    "rcar[89]"    "rcar[90]"    "rcar[91]"    "rcar[92]"    "rcar[93]"   
[100] "rcar[94]"    "rcar[95]"    "rcar[96]"    "rcar[97]"    "rcar[98]"    "rcar[99]"    "rcar[100]"   "lp__"      

[strengejacke]

How do we get that estimate? The car consists of 100 columns of posterior draws, and there's just one single estimate.

[mattansb]

The 4th column is the parameter (car). The others (rcar) are (I think) components of the total parameter, kind of like random effects? Really should learn Stan at some point...

[strengejacke]

should work now after updating insight and parameters.

[mcewenkhundi]

Yes I can confirm, it is now working. Thank you so much.