Calculating grand mean from brms cumulative-link model with sum coded categorical variables

Dear all,

I can’t get my head around this.
I fitted a cumulative-link mixed model with factors situation (levels: Internist, Surgeon), disfigurement (levels: None, Strabism, Acne, Piercing, Tattoo) and face (levels: Man_Average, Woman_Average, Man_Attractive, Woman_Attractive) effect coded (-1, 0, 1) with interactions between all factors. The answer is a likert scale with 11 items (numbers 0-10), hence 10 thresholds.

summary(model)
#> Loading required package: readr
#>  Family: cumulative
#>   Links: mu = logit; disc = identity
#> Formula: answer ~ 1 + situation * face * disfigurement + (1 | participant_id)
#>    Data: original_data (Number of observations: 16036)
#> Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
#>          total post-warmup samples = 10000
#>
#> Group-Level Effects:
#> ~participant_id (Number of levels: 4009)
#>               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
#> sd(Intercept)     1.39      0.03     1.33     1.44       3058 1.00
#>
#> Population-Level Effects:
#>                                 Estimate Est.Error l-95% CI u-95% CI
#> Intercept[1]                       -3.07      0.04    -3.15    -2.99
#> Intercept[2]                       -2.26      0.04    -2.33    -2.19
#> Intercept[3]                       -1.50      0.03    -1.56    -1.44
#> Intercept[4]                       -0.77      0.03    -0.83    -0.71
#> Intercept[5]                       -0.15      0.03    -0.21    -0.09
#> Intercept[6]                        0.57      0.03     0.51     0.63
#> Intercept[7]                        1.22      0.03     1.16     1.28
#> Intercept[8]                        1.98      0.03     1.92     2.05
#> Intercept[9]                        3.02      0.04     2.94     3.10
#> Intercept[10]                       3.88      0.05     3.78     3.97
#> situation1                         -0.23      0.03    -0.28    -0.17
#> face1                              -0.75      0.03    -0.80    -0.70
#> face2                              -0.08      0.03    -0.13    -0.03
#> face3                               0.05      0.03     0.00     0.10
#> disfigurement1                      0.68      0.05     0.58     0.78
#> disfigurement2                     -0.66      0.03    -0.72    -0.60
#> disfigurement3                     -0.54      0.03    -0.60    -0.48
#> disfigurement4                      0.52      0.03     0.46     0.58
#> situation1:face1                    0.10      0.02     0.05     0.15
#> situation1:face2                   -0.01      0.02    -0.06     0.03
#> situation1:face3                   -0.04      0.03    -0.09     0.00
#> situation1:disfigurement1           0.00      0.04    -0.07     0.08
#> situation1:disfigurement2          -0.37      0.03    -0.42    -0.31
#> situation1:disfigurement3           0.35      0.03     0.29     0.40
#> situation1:disfigurement4           0.00      0.02    -0.04     0.05
#> face1:disfigurement1               -0.19      0.05    -0.28    -0.09
#> face2:disfigurement1               -0.01      0.04    -0.09     0.06
#> face3:disfigurement1               -0.17      0.05    -0.27    -0.08
#> face1:disfigurement2                0.35      0.07     0.22     0.49
#> face2:disfigurement2               -0.24      0.07    -0.37    -0.11
#> face3:disfigurement2                0.20      0.07     0.06     0.33
#> face1:disfigurement3                0.13      0.07     0.00     0.26
#> face2:disfigurement3                0.46      0.07     0.33     0.59
#> face3:disfigurement3                0.05      0.06    -0.05     0.17
#> face1:disfigurement4               -0.17      0.07    -0.30    -0.03
#> face2:disfigurement4                0.05      0.06    -0.05     0.18
#> face3:disfigurement4               -0.15      0.07    -0.28    -0.01
#> situation1:face1:disfigurement1    -0.01      0.03    -0.08     0.06
#> situation1:face2:disfigurement1     0.00      0.03    -0.07     0.07
#> situation1:face3:disfigurement1    -0.02      0.04    -0.10     0.05
#> situation1:face1:disfigurement2     0.11      0.06    -0.00     0.23
#> situation1:face2:disfigurement2     0.03      0.05    -0.06     0.13
#> situation1:face3:disfigurement2     0.00      0.04    -0.08     0.09
#> situation1:face1:disfigurement3    -0.12      0.06    -0.25    -0.00
#> situation1:face2:disfigurement3     0.08      0.06    -0.01     0.20
#> situation1:face3:disfigurement3     0.01      0.04    -0.08     0.11
#> situation1:face1:disfigurement4     0.03      0.05    -0.05     0.14
#> situation1:face2:disfigurement4    -0.04      0.05    -0.15     0.05
#> situation1:face3:disfigurement4    -0.01      0.04    -0.11     0.07
#>                                 Eff.Sample Rhat
#> Intercept[1]                          5053 1.00
#> Intercept[2]                          5241 1.00
#> Intercept[3]                          5107 1.00
#> Intercept[4]                          5012 1.00
#> Intercept[5]                          5077 1.00
#> Intercept[6]                          5321 1.00
#> Intercept[7]                          5520 1.00
#> Intercept[8]                          5834 1.00
#> Intercept[9]                          6429 1.00
#> Intercept[10]                         6855 1.00
#> situation1                            4373 1.00
#> face1                                 8374 1.00
#> face2                                 8222 1.00
#> face3                                 6520 1.00
#> disfigurement1                        4736 1.00
#> disfigurement2                        8346 1.00
#> disfigurement3                        7682 1.00
#> disfigurement4                        8303 1.00
#> situation1:face1                      8468 1.00
#> situation1:face2                     10101 1.00
#> situation1:face3                      7254 1.00
#> situation1:disfigurement1             5276 1.00
#> situation1:disfigurement2             9401 1.00
#> situation1:disfigurement3             8789 1.00
#> situation1:disfigurement4             9774 1.00
#> face1:disfigurement1                  9936 1.00
#> face2:disfigurement1                 12024 1.00
#> face3:disfigurement1                  9456 1.00
#> face1:disfigurement2                  4408 1.00
#> face2:disfigurement2                  5068 1.00
#> face3:disfigurement2                  4457 1.00
#> face1:disfigurement3                  4173 1.00
#> face2:disfigurement3                  4164 1.00
#> face3:disfigurement3                  5064 1.00
#> face1:disfigurement4                  4027 1.00
#> face2:disfigurement4                  3905 1.00
#> face3:disfigurement4                  3645 1.00
#> situation1:face1:disfigurement1      11025 1.00
#> situation1:face2:disfigurement1      11474 1.00
#> situation1:face3:disfigurement1       9777 1.00
#> situation1:face1:disfigurement2       3942 1.00
#> situation1:face2:disfigurement2       6539 1.00
#> situation1:face3:disfigurement2       7183 1.00
#> situation1:face1:disfigurement3       3990 1.00
#> situation1:face2:disfigurement3       3861 1.00
#> situation1:face3:disfigurement3       7052 1.00
#> situation1:face1:disfigurement4       5240 1.00
#> situation1:face2:disfigurement4       4725 1.00
#> situation1:face3:disfigurement4       5493 1.00
#>
#> Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
#> is a crude measure of effective sample size, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).

Created on 2019-06-18 by the reprex package (v0.2.1)

Now I would like to calculate the grand mean but I don’t get the idea how it is done. Could someone please explain how I can achieve this?

Thank you very much in advance!

Kindest regards,
Pascal

Session info
devtools::session_info()
#> ─ Session info ──────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 3.5.1 (2018-07-02)
#>  os       CentOS Linux 7 (Core)
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Zurich
#>  date     2019-06-18
#>
#> ─ Packages ──────────────────────────────────────────────────────────────
#>  package        * version date       lib
#>  abind            1.4-5   2016-07-21 [2]
#>  assertthat       0.2.1   2019-03-21 [1]
#>  backports        1.1.4   2019-04-10 [1]
#>  base64enc        0.1-3   2015-07-28 [2]
#>  bayesplot        1.7.0   2019-05-23 [1]
#>  bridgesampling   0.6-0   2018-10-21 [1]
#>  brms             2.9.0   2019-06-07 [1]
#>  Brobdingnag      1.2-6   2018-08-13 [1]
#>  callr            3.2.0   2019-03-15 [1]
#>  cli              1.1.0   2019-03-19 [1]
#>  coda             0.19-2  2018-10-08 [2]
#>  colorspace       1.4-1   2019-03-18 [1]
#>  colourpicker     1.0     2017-09-27 [1]
#>  crayon           1.3.4   2017-09-16 [2]
#>  crosstalk        1.0.0   2016-12-21 [2]
#>  desc             1.2.0   2018-05-01 [2]
#>  devtools         2.0.2   2019-04-08 [1]
#>  digest           0.6.19  2019-05-20 [1]
#>  dplyr            0.8.1   2019-05-14 [1]
#>  DT               0.6     2019-05-09 [1]
#>  dygraphs         1.1.1.6 2018-07-11 [1]
#>  evaluate         0.12    2018-10-09 [2]
#>  fs               1.3.0   2019-05-02 [1]
#>  ggplot2          3.1.1   2019-04-07 [1]
#>  ggridges         0.5.1   2018-09-27 [2]
#>  glue             1.3.1   2019-03-12 [1]
#>  gridExtra        2.3     2017-09-09 [2]
#>  gtable           0.3.0   2019-03-25 [1]
#>  gtools           3.8.1   2018-06-26 [2]
#>  hms              0.4.2   2018-03-10 [2]
#>  htmltools        0.3.6   2017-04-28 [2]
#>  htmlwidgets      1.3     2018-09-30 [2]
#>  httpuv           1.5.1   2019-04-05 [1]
#>  igraph           1.2.4.1 2019-04-22 [1]
#>  inline           0.3.15  2018-05-18 [2]
#>  knitr            1.20    2018-02-20 [2]
#>  later            0.8.0   2019-02-11 [1]
#>  lattice          0.20-35 2017-03-25 [2]
#>  lazyeval         0.2.2   2019-03-15 [1]
#>  loo              2.1.0   2019-03-13 [1]
#>  magrittr         1.5     2014-11-22 [2]
#>  markdown         0.9     2018-12-07 [1]
#>  Matrix           1.2-14  2018-04-13 [2]
#>  matrixStats      0.54.0  2018-07-23 [2]
#>  memoise          1.1.0   2017-04-21 [2]
#>  mime             0.6     2018-10-05 [2]
#>  miniUI           0.1.1.1 2018-05-18 [2]
#>  munsell          0.5.0   2018-06-12 [2]
#>  mvtnorm          1.0-10  2019-03-05 [1]
#>  nlme             3.1-137 2018-04-07 [2]
#>  pillar           1.4.1   2019-05-28 [1]
#>  pkgbuild         1.0.3   2019-03-20 [1]
#>  pkgconfig        2.0.2   2018-08-16 [2]
#>  pkgload          1.0.1   2018-10-11 [2]
#>  plyr             1.8.4   2016-06-08 [2]
#>  prettyunits      1.0.2   2015-07-13 [2]
#>  processx         3.3.1   2019-05-08 [1]
#>  promises         1.0.1   2018-04-13 [2]
#>  ps               1.3.0   2018-12-21 [1]
#>  purrr            0.3.2   2019-03-15 [1]
#>  R6               2.4.0   2019-02-14 [1]
#>  Rcpp             1.0.1   2019-03-17 [1]
#>  readr          * 1.1.1   2017-05-16 [2]
#>  remotes          2.0.4   2019-04-10 [2]
#>  reshape2         1.4.3   2017-12-11 [2]
#>  rlang            0.3.4   2019-04-07 [1]
#>  rmarkdown        1.10    2018-06-11 [2]
#>  rprojroot        1.3-2   2018-01-03 [2]
#>  rsconnect        0.8.13  2019-01-10 [1]
#>  rstan            2.18.2  2018-11-07 [1]
#>  rstantools       1.5.1   2018-08-22 [1]
#>  scales           1.0.0   2018-08-09 [2]
#>  sessioninfo      1.1.1   2018-11-05 [1]
#>  shiny            1.3.2   2019-04-22 [1]
#>  shinyjs          1.0     2018-01-08 [1]
#>  shinystan        2.5.0   2018-05-01 [1]
#>  shinythemes      1.1.2   2018-11-06 [1]
#>  StanHeaders      2.18.1  2019-01-28 [1]
#>  stringi          1.4.3   2019-03-12 [1]
#>  stringr          1.4.0   2019-02-10 [1]
#>  testthat         2.0.1   2018-10-13 [2]
#>  threejs          0.3.1   2017-08-13 [1]
#>  tibble           2.1.3   2019-06-06 [1]
#>  tidyselect       0.2.5   2018-10-11 [1]
#>  usethis          1.5.0   2019-04-07 [1]
#>  withr            2.1.2   2018-03-15 [2]
#>  xtable           1.8-4   2019-04-21 [1]
#>  xts              0.11-2  2018-11-05 [1]
#>  yaml             2.2.0   2018-07-25 [2]
#>  zoo              1.8-6   2019-05-28 [1]
#>  source

Please also provide the following information in addition to your question:

  • Operating System: CentOS 7.0
  • brms Version: 2.9.0

I believe one approach would be to use the function fitted() along with some post-processing.

Specifically:

  1. Define all possible combinations of the predictor variables. For example, if you had two predictors (x1 and x2) with two levels each, that would be:
conditions <- dplyr::tribble(
  ~x1, ~x2
  'x1_level1', 'x2_level1',
  'x1_level1', 'x2_level2',
  'x1_level2', 'x2_level1',
  'x1_level2', 'x2_level2'
)
  1. Use the function fitted() to obtain the mean predictions.
model_predictions <- cbind(
  conditions,
  fitted(model, newdata = conditions)
)

model_predictions$mu <- with(model_predictions,
  `Estimate.P(Y = 1)` * 1 +
  `Estimate.P(Y = 2)` * 2 +
  `Estimate.P(Y = 3)` * 3 # etc.
)
  1. Compute the grand mean.
mean(model_predictions$mu)

Hi @aljaz,

Thank you very much for your answer! I think I asked the wrong question, because the answer is not what I expected :-D

But first another question. I don’t see what you are doing here, could you please elaborate?

Now my proper question:
I want to get something like a “sum coefficient” for every combination of predictors, so I can report “For the average looking man with strabismus being a surgeon, the probability increases/decreases with x.xxx (y.yyy, z.zzz) compared with no disfigurement”. Is there a way to achieve this with proper credible intervals?

Thanks and regards,
Pascal

Hi Pascal,

But first another question. I don’t see what you are doing here, could you please elaborate?

Computing the expected value for each of the four combinations of the two predictors.

I want to get something like a “sum coefficient” for every combination of predictors, so I can report “For the average looking man with strabismus being a surgeon, the probability increases/decreases with x.xxx (y.yyy, z.zzz) compared with no disfigurement”. Is there a way to achieve this with proper credible intervals?

Using the function fitted() you can obtain non-summarized predictions by setting the summary argument to FALSE. Just like I did above, you’d define the various combinations of predictors that you’re interested in and pass them to fitted() as new data. Then you could compute the differences between these different combinations of predictors and summarize them with credible intervals. But I’m not sure what you want to summarize. You’re talking about a probability … of what? Giving a certain response?