Enforcing conditional monotonicity

Interfaces brms
1

In Bürkner, Paul-Christian, and Emmanuel Charpentier. “Monotonic Effects: A Principled Approach for Including Ordinal Predictors in Regression Models,” it is stated that conditional monotonicity can be achieved when including interaction effects by using a cell means parameterization. However, it seems to me that that isn’t the case because although there is no reference category for main effects, there is still an implicit reference category for the interaction effects.

See below for a simple simulation demonstrating this. Am I doing something wrong here? If not, is there a way to enforce conditional monotonicity other than fixing all simplex parameters to the same value for main and interaction effects? Lastly, is there a way in brms to fix simplex parameters to the same value for main and interaction effects?

set.seed(123)
# Monotonic in group1
group1a <- rnorm(10, mean=10, sd=1)
group1b <- rnorm(10, mean=10 +.4*2, sd=1)
group1c <- rnorm(10, mean=10 + 2, sd=1)

# NOT monotonic in group2
group2a <- rnorm(10, mean=8, sd=1)
group2b <- rnorm(10, mean=8 - 1.2*2, sd=1)
group2c <- rnorm(10, mean=8 - 2, sd=1)

sim_data <- data.frame(y=c(group1a, group1b, group1c, group2a, group2b, group2c),
                       X1=rep(c("Group1", "Group2"), each=30),
                       X2=ordered(rep(rep(c("LevelA", "LevelB", "LevelC"), each=10), 2)))

# Even though X2 is not conditionally monotonic, we want to model it as if it is
fit_cellmeans <- brm(bf(y ~ 0 + X1 * mo(X2), cmc=TRUE),
                     data=sim_data,
                     chains=4,
                     cores=4)

fit_cellmeans

Population-Level Effects: 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
X1Group1         10.14      0.29     9.57    10.71       1687 1.00
X1Group2          8.30      0.30     7.71     8.87       2839 1.00
moX2              1.51      0.41     0.72     2.30       1736 1.00
moX2:X1Group2    -3.71      0.54    -4.79    -2.66       1740 1.00

Simplex Parameters: 
                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
moX21[1]              0.48      0.19     0.09     0.85       2478 1.00
moX21[2]              0.52      0.19     0.15     0.91       2478 1.00
moX2:X1Group21[1]     0.89      0.08     0.70     1.00       3368 1.00
moX2:X1Group21[2]     0.11      0.08     0.00     0.30       3368 1.00

sim_data %>% 
  bind_cols(., as.data.frame(fitted(fit_cellmeans))) %>% 
  distinct(X1, X2, .keep_all = T) %>% 
  select(-y)

# X2 is monotonic in Group1, but not in Group2
      X1     X2  Estimate Est.Error      Q2.5     Q97.5
1 Group1 LevelA 10.142398 0.2926905  9.569837 10.706621
2 Group1 LevelB 10.880892 0.2522873 10.394014 11.387516
3 Group1 LevelC 11.649422 0.2748625 11.126107 12.198392
4 Group2 LevelA  8.298298 0.3002720  7.710111  8.872855
5 Group2 LevelB  5.737879 0.2682472  5.225042  6.273904
6 Group2 LevelC  6.097649 0.2689416  5.551706  6.637781
  • Operating System: macOS High Sierra 10.13.1
  • brms Version: brms_2.8.0
2

How does y ~ 0 + X1:mo(X2) work out?

3

Conditional monotonicity is achieved, but no main effects are estimated and as a result, we are not able to get estimates for observations in the lowest level of X2:

fit2_cellmeans <- brm(bf(y ~ 0 + X1 : mo(X2), cmc=TRUE),
                     data=sim_data,
                     chains=4,
                     cores=4,
                     seed=123)

sim_data %>% 
  bind_cols(., as.data.frame(fitted(fit2_cellmeans))) %>% 
  distinct(X1, X2, .keep_all = T) %>% 
  select(-y)

      X1     X2  Estimate Est.Error     Q2.5     Q97.5
1 Group1 LevelA  0.000000  0.000000 0.000000  0.000000
2 Group1 LevelB 10.108732  1.454284 7.079374 12.794374
3 Group1 LevelC 12.259702  1.512303 9.432396 15.374217
4 Group2 LevelA  0.000000  0.000000 0.000000  0.000000
5 Group2 LevelB  4.619737  1.488213 1.508115  7.463217
6 Group2 LevelC  6.704395  1.556781 3.676564  9.821315

Also, FWIW the model above fitted with cmc=T and cmc=F give the same results with the same seed. I think this makes sense given there are no main effects.

4

What happens if we went for 0 + X1 + X1:mo(X2)?

5

That worked, and it makes sense that now, we’re estimating monotonic effects separately for Group1 and Group2. Thanks!

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: y ~ 0 + X1 + X1:mo(X2) 
   Data: sim_data (Number of observations: 60) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
X1Group1         10.13      0.30     9.53    10.72       2826 1.00
X1Group2          8.27      0.30     7.67     8.85       3138 1.00
X1Group1:moX2     1.42      0.43     0.61     2.28       2672 1.00
X1Group2:moX2    -2.44      0.38    -3.17    -1.68       3137 1.00

Simplex Parameters: 
                  Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
X1Group1:moX21[1]     0.59      0.22     0.11     0.97       3460 1.00
X1Group1:moX21[2]     0.41      0.22     0.03     0.89       3460 1.00
X1Group2:moX21[1]     0.92      0.07     0.74     1.00       4311 1.00
X1Group2:moX21[2]     0.08      0.07     0.00     0.26       4311 1.00