# How to simulate monotonic effects with an ordinal outcome
**URL:** https://discourse.mc-stan.org/t/how-to-simulate-monotonic-effects-with-an-ordinal-outcome/28705
**Category:** brms
**Tags:** brms
**Created:** [September 1, 2022, 1:01am UTC](https://discourse.mc-stan.org/t/how-to-simulate-monotonic-effects-with-an-ordinal-outcome/28705 "2022-09-01T01:01:49Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@shirdekel](https://discourse.mc-stan.org/u/shirdekel)
#### Post date: [September 1, 2022, 1:01am UTC](https://discourse.mc-stan.org/t/how-to-simulate-monotonic-effects-with-an-ordinal-outcome/28705/1 "2022-09-01T01:01:49Z")
I’d like to simulate ordinal data for later analysis with `brms`. At its simplest, I have an ordinal predictor with an ordinal outcome, which I would analyse using something like this:
```r
fit1 <- brm(ordinal_outcome ~ mo(ordinal_predictor), data = d, family = cumulative())
```
The [Estimating Monotonic Effects with brms vignette](https://paul-buerkner.github.io/brms/articles/brms_monotonic.html) uses the following code to simulate an ordinal predictor with a continuous outcome.
```r
income_options <- c("below_20", "20_to_40", "40_to_100", "greater_100")
income <- factor(sample(income_options, 100, TRUE),
levels = income_options, ordered = TRUE)
mean_ls <- c(30, 60, 70, 75)
ls <- mean_ls[income] + rnorm(100, sd = 7)
dat <- data.frame(income, ls)
```
How can I simulate a data set with an ordinal predictor and an ordinal outcome?
- Operating System: macOS 12.5.1
- brms Version: 2.17.0
---
### Author:  [@cgoold](https://discourse.mc-stan.org/u/cgoold)
#### Post date: [September 10, 2022, 7:35am UTC](https://discourse.mc-stan.org/t/how-to-simulate-monotonic-effects-with-an-ordinal-outcome/28705/2 "2022-09-10T07:35:23Z")
To simulate an ordinal outcome, you need to follow the correct data generating assumptions.
An ordinal outcome y is assumed categorically distributed with K category probabilities \mathbf{\pi} = \{\pi\_{1}, \pi\_{2}, ..., \pi\_{K}\}:
y \sim Categorical(\mathbf{\pi})
The category probabilities are generated from an underlying continuous scale split up into K bins using K - 1 thresholds, \mathbf{\theta} = \{\theta\_{1}, ..., \theta\_{K - 1}\}. With a suitable CDF, often the logistic but you could also choose the normal CDF, we can work out the probability in each bin. If we refer to the standard logistic CDF (with mean=0 and scale=1) as F, then:
\pi\_{k} = F(\theta\_{k} - \eta) - F(\theta\_{k-1} - \eta)
where I’m assuming that \theta\_{0} = - \infty and \theta\_{K} = \infty so that F(\theta\_{0} - \eta) = 0 and F(\theta\_{K} - \eta) = 1. The \eta term holds the linear predictor, i.e. your ordinal predictor variable.
In R, you’d have something like (just quickly wrote this, it might contain errors):
```r
# Generate a single ordinal data point
# @K number of ordinal categories
# @theta vector of latent cutpoints
# @eta linear predictor
gen_ordinal_data <- function(K, theta, eta){
if(length(theta) != (K - 1))
stop(paste0("theta must have length ", K - 1, " not ", length(K)))
probs <- c()
for(k in 1:K){
if(k == 1) prev_thresh = -Inf
else prev_thresh = theta[k - 1]
if(k == K) next_thresh = Inf
else next_thresh = theta[k]
probs[k] = plogis(next_thresh - eta) - plogis(prev_thresh - eta)
}
return(sample(K, 1, prob=probs))
}
```
Then, you can plug whatever you want in for `eta` including your ordinal predictor variables.
---
### Author:  [@shirdekel](https://discourse.mc-stan.org/u/shirdekel)
#### Post date: [September 16, 2022, 2:11am UTC](https://discourse.mc-stan.org/t/how-to-simulate-monotonic-effects-with-an-ordinal-outcome/28705/3 "2022-09-16T02:11:44Z")
Thank you so much for that! I really appreciate the detail.
1. How would you explain the reason for taking the difference between \theta and \eta (e.g., within F(\theta\_k - \eta))? Is it just because increasing \eta shifts the distribution of y to the right and decreasing it shifts the distribution to the left?
2. How can I do parameter recovery with this code? That is, how can I set up the simulation so that I can anticipate the monotonic coefficient I will get in the `brms` fit output?
Here is my attempt at using the above code for a simulation and subsequent modelling:
```r
## Generate a single ordinal data point
# @K number of ordinal categories
# @theta vector of latent cutpoints
# @eta linear predictor
gen_ordinal_data <- function(K, theta, eta) {
if (length(theta) != (K - 1)) {
stop(paste0("theta must have length ", K - 1, " not ", length(K)))
}
probs <- c()
for (k in 1:K) {
if (k == 1) {
prev_thresh <- -Inf
} else {
prev_thresh <- theta[k - 1]
}
if (k == K) {
next_thresh <- Inf
} else {
next_thresh <- theta[k]
}
probs[k] <- plogis(next_thresh - eta) - plogis(prev_thresh - eta)
}
return(sample(K, 1, prob = probs))
}
library(purrr)
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(tibble)
library(brms)
#> Loading required package: Rcpp
#> Loading 'brms' package (version 2.17.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
options(mc.cores = parallel::detectCores())
set.seed(1)
K <- 5
theta <- seq(from = -4, to = 4, length.out = 4)
n <- 500
eta <- seq_len(5)
data_sim <-
tibble(
ID = seq_len(n)
) |>
rowwise() |>
mutate(
variable_1 = sample(eta, 1),
variable_2 = gen_ordinal_data(
K = K,
theta = theta,
eta = variable_1
)
)
data_sim
#> # A tibble: 500 × 3
#> # Rowwise:
#> ID variable_1 variable_2
#>
#> 1 1 1 4
#> 2 2 4 5
#> 3 3 1 4
#> 4 4 2 4
#> 5 5 5 5
#> 6 6 3 4
#> 7 7 2 4
#> 8 8 3 4
#> 9 9 3 4
#> 10 10 1 3
#> # … with 490 more rows
fit1 <- brm(variable_2 ~ mo(variable_1), data = data_sim, family = cumulative())
#> Compiling Stan program...
#> Trying to compile a simple C file
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
#> In file included from :1:
#> In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
#> /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
#> namespace Eigen {
#> ^
#> /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
#> namespace Eigen {
#> ^
#> ;
#> In file included from :1:
#> In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
#> /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#> #include
#> ^ ~~~~~~~~
#> 3 errors generated.
#> make: *** [foo.o] Error 1
#> Start sampling
fit1
#> Family: cumulative
#> Links: mu = logit; disc = identity
#> Formula: variable_2 ~ mo(variable_1)
#> Data: data_sim (Number of observations: 500)
#> Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#> total post-warmup draws = 4000
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept[1] -5.16 1.06 -7.72 -3.58 1.00 1602 1062
#> Intercept[2] -2.06 0.26 -2.58 -1.58 1.00 3879 2987
#> Intercept[3] 0.34 0.19 -0.03 0.72 1.00 3307 3324
#> Intercept[4] 2.93 0.25 2.47 3.43 1.00 2356 2708
#> movariable_1 1.03 0.08 0.86 1.20 1.00 2073 2275
#>
#> Simplex Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> movariable_11[1] 0.24 0.06 0.12 0.35 1.00 2286 1677
#> movariable_11[2] 0.20 0.07 0.07 0.34 1.00 2222 1836
#> movariable_11[3] 0.35 0.07 0.21 0.49 1.00 3304 2504
#> movariable_11[4] 0.21 0.07 0.08 0.34 1.00 2892 1425
#>
#> Family Specific Parameters:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> disc 1.00 0.00 1.00 1.00 NA NA NA
#>
#> Draws were sampled 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).
```
Created on 2022-09-16 by the [reprex package](https://reprex.tidyverse.org) (v2.0.1)
> **Session info**
>
> ```r
> sessioninfo::session_info()
> #> ─ Session info ───────────────────────────────────────────────────────────────
> #> setting value
> #> version R version 4.2.1 (2022-06-23)
> #> os macOS Big Sur ... 10.16
> #> system x86_64, darwin17.0
> #> ui X11
> #> language (EN)
> #> collate en_AU.UTF-8
> #> ctype en_AU.UTF-8
> #> tz Australia/Sydney
> #> date 2022-09-16
> #> pandoc 2.19 @ /usr/local/bin/ (via rmarkdown)
> #>
> #> ─ Packages ───────────────────────────────────────────────────────────────────
> #> package * version date (UTC) lib source
> #> abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
> #> assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.2.0)
> #> backports 1.4.1 2021-12-13 [1] CRAN (R 4.2.0)
> #> base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.2.0)
> #> bayesplot 1.9.0 2022-03-10 [1] CRAN (R 4.2.0)
> #> bridgesampling 1.1-2 2021-04-16 [1] CRAN (R 4.2.0)
> #> brms * 2.17.0 2022-04-13 [1] CRAN (R 4.2.0)
> #> Brobdingnag 1.2-7 2022-02-03 [1] CRAN (R 4.2.0)
> #> callr 3.7.1 2022-07-13 [1] CRAN (R 4.2.0)
> #> checkmate 2.1.0 2022-04-21 [1] CRAN (R 4.2.0)
> #> cli 3.3.0 2022-04-25 [1] CRAN (R 4.2.0)
> #> coda 0.19-4 2020-09-30 [1] CRAN (R 4.2.0)
> #> codetools 0.2-18 2020-11-04 [1] CRAN (R 4.2.1)
> #> colorspace 2.0-3 2022-02-21 [1] CRAN (R 4.2.0)
> #> colourpicker 1.1.1 2021-10-04 [1] CRAN (R 4.2.0)
> #> crayon 1.5.1 2022-03-26 [1] CRAN (R 4.2.0)
> #> crosstalk 1.2.0 2021-11-04 [1] CRAN (R 4.2.0)
> #> DBI 1.1.3 2022-06-18 [1] CRAN (R 4.2.0)
> #> digest 0.6.29 2021-12-01 [1] CRAN (R 4.2.0)
> #> distributional 0.3.0 2022-01-05 [1] CRAN (R 4.2.0)
> #> dplyr * 1.0.9 2022-04-28 [1] CRAN (R 4.2.0)
> #> DT 0.24 2022-08-09 [1] CRAN (R 4.2.0)
> #> dygraphs 1.1.1.6 2018-07-11 [1] CRAN (R 4.2.0)
> #> ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.0)
> #> evaluate 0.16 2022-08-09 [1] CRAN (R 4.2.0)
> #> fansi 1.0.3 2022-03-24 [1] CRAN (R 4.2.0)
> #> farver 2.1.1 2022-07-06 [1] CRAN (R 4.2.0)
> #> fastmap 1.1.0 2021-01-25 [1] CRAN (R 4.2.0)
> #> fs 1.5.2 2021-12-08 [1] CRAN (R 4.2.0)
> #> generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.0)
> #> ggplot2 3.3.6 2022-05-03 [1] CRAN (R 4.2.0)
> #> ggridges 0.5.3 2021-01-08 [1] CRAN (R 4.2.0)
> #> glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0)
> #> gridExtra 2.3 2017-09-09 [1] CRAN (R 4.2.0)
> #> gtable 0.3.0 2019-03-25 [1] CRAN (R 4.2.0)
> #> gtools 3.9.3 2022-07-11 [1] CRAN (R 4.2.0)
> #> highr 0.9 2021-04-16 [1] CRAN (R 4.2.0)
> #> htmltools 0.5.3 2022-07-18 [1] CRAN (R 4.2.0)
> #> htmlwidgets 1.5.4 2021-09-08 [1] CRAN (R 4.2.0)
> #> httpuv 1.6.5 2022-01-05 [1] CRAN (R 4.2.0)
> #> igraph 1.3.4 2022-07-19 [1] CRAN (R 4.2.0)
> #> inline 0.3.19 2021-05-31 [1] CRAN (R 4.2.0)
> #> knitr 1.39 2022-04-26 [1] CRAN (R 4.2.0)
> #> later 1.3.0 2021-08-18 [1] CRAN (R 4.2.0)
> #> lattice 0.20-45 2021-09-22 [1] CRAN (R 4.2.1)
> #> lifecycle 1.0.1 2021-09-24 [1] CRAN (R 4.2.0)
> #> loo 2.5.1 2022-03-24 [1] CRAN (R 4.2.0)
> #> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.0)
> #> markdown 1.1 2019-08-07 [1] CRAN (R 4.2.0)
> #> Matrix 1.4-1 2022-03-23 [1] CRAN (R 4.2.1)
> #> matrixStats 0.62.0 2022-04-19 [1] CRAN (R 4.2.0)
> #> mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
> #> miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.0)
> #> munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.0)
> #> mvtnorm 1.1-3 2021-10-08 [1] CRAN (R 4.2.0)
> #> nlme 3.1-157 2022-03-25 [1] CRAN (R 4.2.1)
> #> pillar 1.8.1 2022-08-19 [1] CRAN (R 4.2.0)
> #> pkgbuild 1.3.1 2021-12-20 [1] CRAN (R 4.2.0)
> #> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0)
> #> plyr 1.8.7 2022-03-24 [1] CRAN (R 4.2.0)
> #> posterior 1.3.0 2022-08-15 [1] CRAN (R 4.2.1)
> #> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.0)
> #> processx 3.7.0 2022-07-07 [1] CRAN (R 4.2.0)
> #> promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.0)
> #> ps 1.7.1 2022-06-18 [1] CRAN (R 4.2.0)
> #> purrr * 0.3.4 2020-04-17 [1] CRAN (R 4.2.0)
> #> R.cache 0.16.0 2022-07-21 [1] CRAN (R 4.2.0)
> #> R.methodsS3 1.8.2 2022-06-13 [1] CRAN (R 4.2.0)
> #> R.oo 1.25.0 2022-06-12 [1] CRAN (R 4.2.0)
> #> R.utils 2.12.0 2022-06-28 [1] CRAN (R 4.2.0)
> #> R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0)
> #> Rcpp * 1.0.9 2022-07-08 [1] CRAN (R 4.2.0)
> #> RcppParallel 5.1.5 2022-01-05 [1] CRAN (R 4.2.0)
> #> reprex 2.0.1 2021-08-05 [1] CRAN (R 4.2.0)
> #> reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.2.0)
> #> rlang 1.0.4 2022-07-12 [1] CRAN (R 4.2.0)
> #> rmarkdown 2.14 2022-04-25 [1] CRAN (R 4.2.0)
> #> rstan 2.21.5 2022-04-11 [1] CRAN (R 4.2.0)
> #> rstantools 2.2.0 2022-04-08 [1] CRAN (R 4.2.0)
> #> scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.0)
> #> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0)
> #> shiny 1.7.2 2022-07-19 [1] CRAN (R 4.2.0)
> #> shinyjs 2.1.0 2021-12-23 [1] CRAN (R 4.2.0)
> #> shinystan 2.6.0 2022-03-03 [1] CRAN (R 4.2.0)
> #> shinythemes 1.2.0 2021-01-25 [1] CRAN (R 4.2.0)
> #> StanHeaders 2.21.0-7 2020-12-17 [1] CRAN (R 4.2.0)
> #> stringi 1.7.8 2022-07-11 [1] CRAN (R 4.2.0)
> #> stringr 1.4.0 2019-02-10 [1] CRAN (R 4.2.0)
> #> styler 1.7.0 2022-03-13 [1] CRAN (R 4.2.0)
> #> tensorA 0.36.2 2020-11-19 [1] CRAN (R 4.2.0)
> #> threejs 0.3.3 2020-01-21 [1] CRAN (R 4.2.0)
> #> tibble * 3.1.8 2022-07-22 [1] CRAN (R 4.2.0)
> #> tidyselect 1.1.2 2022-02-21 [1] CRAN (R 4.2.0)
> #> utf8 1.2.2 2021-07-24 [1] CRAN (R 4.2.0)
> #> vctrs 0.4.1 2022-04-13 [1] CRAN (R 4.2.0)
> #> withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.0)
> #> xfun 0.31 2022-05-10 [1] CRAN (R 4.2.0)
> #> xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.0)
> #> xts 0.12.1 2020-09-09 [1] CRAN (R 4.2.0)
> #> yaml 2.3.5 2022-02-21 [1] CRAN (R 4.2.0)
> #> zoo 1.8-10 2022-04-15 [1] CRAN (R 4.2.0)
> #>
> #> [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library
> #>
> #> ──────────────────────────────────────────────────────────────────────────────
>
> ```
---
### Author:  [@cgoold](https://discourse.mc-stan.org/u/cgoold)
#### Post date: [October 4, 2022, 7:04am UTC](https://discourse.mc-stan.org/t/how-to-simulate-monotonic-effects-with-an-ordinal-outcome/28705/4 "2022-10-04T07:04:47Z")
@shirdekel
1. Yes, that’s correct. As the predictor variables increase in size, we subtract because that shifts the thresholds to the left, providing more probability density under the highest categories.
2. If I understand your model, you’d need to fist simulate an ordinal predictor variable, then you’d need to use the `gen_ordinal_data` function to simulate your outcome ordinal variable, where `eta` is a function of the ordinal predictor. If you weren’t using monotonic effects, `eta` would be somthing like `alpha + beta * ordinal_predictor`, producing a vector of linear predictor values to generate the ordinal response variable. But you’d need to change the `beta * ordinal_predictor` variable specification to the structure of the monotonic predictor assumption.