Seemingly Unrelated Regression with Tobit regression equations

Modeling
1

I want to code for Seemingly Unrelated Tobit Regression model. I can’t find any example to do that in stan.
Can someone please guide me how can I code it?

Thank you!

2

Hi,
I don’t think “Seemingly Unrelated Tobit Regression” is a popular enough model to be well known to the community. It is usually best to write down mathematical formulas and/or links to references describing the model. I found https://dukespace.lib.duke.edu/dspace/bitstream/handle/10161/1881/Huange_estimation_of_seemingly_unrelated_tobit_regressions.pdf

If I understand the model correctly, it is just censored regression with some correlations - if you are using R, multivariate censored regression with some correlation structures should be possible to implement using brms, so you wouldn’t need to dig deep into Stan code.

Otherwise, you should be able to piece together the necessary Stan code from: User’s guide on linear regression, User’s guide on multivariate outcomes (explicitly mentions seemingly-unrelated regression) and User’s guid on censored data. I am not completely sure how you would handle both correlations in noise and censoring, but I hope this gives you a good starting point.

3

Hi Martin,
Thank you for your reply.
Yes, SUR tobit is a censored regression with tobit. The error of two variables are normally distributed with mean and variance-covariance matrix.

I do am using R. Any help with using brms will be appreciated.

I tried coding the bayesian method myself as mentioned in this paper:

However, I am stuck at drawing from the censored data.
I am not sure what sigma to use for the multiple cases.

Do you think I can use brms to estimate this model?

Many thanks,
Usman

4

So, just to be sure - you are able to use brms to fit both:

  1. Univariate model with a censored outcome
  2. Multivariate model with correlation but without censoring

But you do not know how to use both at the same time? If yes, could you share the code you use for both of the cases so that we have a shared starting point. And - if you can - also share an example of the data you use so that I can run the code myself?

5

Yes, exactly this is the problem and therefore I tried to code the process on my own by trying to draw from the conditional distributions.

For brms i am using this code but it is not working:

bf_1 <- brms::bf(y1 | trunc(lb=0, ub = Inf) ~ x1  + x2)  
bf_2 <- brms::bf(y2 | trunc(lb=0, ub = Inf) ~ x1  + x3)

brm_mv <- brm(
  bf_1 + bf_2 , data = data.frame(data.m),
  chains = 2, cores = 2)

The data i am using is confidential. However, I am not able to run the above function.

6

So maybe you could create some simulated dataset that we could both use?

Could you be more specific? Does this produce an error? If yes, could you share the error message? If no, what is the expected behaviour and what is the actual behavior?

You are fitting a truncated model, not a censored one - so do you want to use a censored or a truncated model? (censored models are fit in brms with the cens addition term).

Also note that you can use triple backticks (```) to format code blocks here. (I’ve edited your previous message for you).

7

My bad. yes i am trying to do censored regression and not truncated. I am creating data and will post soon.

1 Like
8

Here is the data I generated:

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'))

Here is the brms code:

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, 
  data = mydata,
  chains = 2, cores = 2)

and here is the error:

Setting 'rescor' to TRUE by default for this model
Error: The following variables are missing in 'data':
'y2'
In addition: Warning message:
In the future, 'rescor' will be set to FALSE by default for all models. It is thus recommended to explicitely set 'rescor' via 'set_recor' instead of using the default.

However, when i am running the code on th original data i am getting this error:

Setting 'rescor' to TRUE by default for this model
Error: Only 'se', 'weights', 'mi' are supported addition arguments when 'rescor' is estimated.
In addition: Warning message:
In the future, 'rescor' will be set to FALSE by default for all models. It is thus recommended to explicitely set 'rescor' via 'set_recor' instead of using the default.
1 Like
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