Seemingly Unrelated Regression with Tobit regression equations

Modeling
9

So the error is probably a minor bug in brms, but if I set rescor excplicitly and call the model as:

bf_1 <- brms::bf(d1 | cens(cen_d1) ~ x1  + x2)  
bf_2 <- brms::bf(d2 | cens(cen_d2) ~ x1  + x3)

brm_mv <- brm(
  bf_1 + bf_2 + set_rescor(FALSE), 
  data = mydata,
  chains = 2, cores = 2)

I get a working model - assuming d1 and d2 are independent. That’s not what you wanted, but this will be a useful baseline to compare against, because the full solution is going to be more involved.

Simple approximate solution

A simple (but imperfect) way to mimic the correlations for censored outcomes is to explicitly add the correlated residuals as latent parameters in the model as a random effect. First we augment the data with an unique_id column that is unique for each row:

mydata <- data.frame(d1 = c(0,0,1,0,3,5,7,0,10,11),
           d2 = c(0,7,2,4,9,0,0,2,5,0),
           x1 = c(0.25,0.45,0.67,0.89,1.2,4.6,8.9,1.1,0.9,1.1),
           x2 = c(1,1,1,0,0,0,1,0,1,0),
           x3 = c(0,1,1,0,0,1,1,0,1,0)
           ) %>% 
  mutate(cen_d1 = ifelse(d1 == 0,'left','none'),
         cen_d2 = ifelse(d2 == 0,'left', 'none'),
         unique_id = 1:n() )

then we change the univariate formula to have a correlated varying intercept for each row:

bf_1 <- brms::bf(d1 | cens(cen_d1) ~ x1  + x2 + (1|p|unique_id))  
bf_2 <- brms::bf(d2 | cens(cen_d2) ~ x1  + x3 + (1|p|unique_id))

However, now the model is not identified as we have the sigma parameter and the sd of the unique_id varying intercept fill almost the same role. What we can do is to fix the sigma of the response to a small constant:

priors <- prior(constant(0.5), class = "sigma", resp = "d1") +
    prior(constant(0.5), class = "sigma", resp = "d2")

brm_mv <- brm(
  bf_1 + bf_2 + set_rescor(FALSE), 
  data = mydata,
  prior =  priors)

The choice of constant (here 0.5) is a bit arbitrary: smaller constants would mean you get closer to the exact SUTR model, but are also more likely to cause problems in sampling.

The problem with this approach is also that it is likely to be slow and has potential to introduce convergence issues.

The advantage is that we are staying within the relative safety of officialy supported and tested brms features.

Precise model - an outline

We could also use the exact version of the model. This would be quite involved, and would require you to check the results much more thoroughly (e.g. by simulating data with known parameters and checking that you can recover those, preferably more as discussed in the Bayesian workflow preprint. So I am giving just an outline - we can dig deeper into this if you are in a position to put extra work there. I’ll start with the math than show an outline how to put that math in action using the custom_family feature of brms and a few hacks.

If only one value in a pair is censored, one can compute the conditional distribution of the second variable as in Multivariate normal distribution - Wikipedia - this transforms the problem into computing the probability of a univariate censored variable, which can be done directly with normal_lcdf (see Stan User’s guide on censoring.

If I understand the literature correctly, the biggest problem becomes when more than one response is censored at the same time. So if your data are only bivariate and contain no (0,0) pairs, the model would become much easier. The reason for this is that if I understand the literature correctly (my main reference is https://arxiv.org/pdf/2012.14331.pdf), then there is no known analytic expression for P(x_1 < 0 \land x_2 < 0) given (x_1, x_2) \sim MVN(\bf{\mu}, \Sigma) which is exactly the value we need to add to the likelihood for the “both censored” case. One could definitely solve the relevant integral numerically, either using integrate_1d in Stan or (possibly more efficiently) implementing the purpose-made method described in the reference above. In any case, this is going to be computationally expensive (but less than the approximate method).

Once we have this in place, we use the same trick as I described at Using a single mixture weight parameter in a multivariate mixture model - #2 by martinmodrak to define a bivariate custom family. We would then use the addition terms vint1 to confer censoring information.

So I think the model is possible to implement, but it would be a bit challenging. We would be happy to help you go through this if you decide to try it, but is likely to be non-trivial amount of effort.

Best of luck with the model!

3 Likes