Multivariate models and mi()

Interfaces brms
1
  • Operating System: Mac
  • brms Version: 2.8.0

Question

Let’s say we have 2 + response variables. One has some missing values and the other(s) are complete. Is there a way to fit a multivariate model where we avoid dropping cases with missingness?

We have examples.

Here’s the typical setup.

dat <-
  tibble(a = rnorm(100),
         b = c(NA, rnorm(99)))

fit1 <-
  brm(data = dat,
      mvbind(a, b) ~ 1)

Before compiling the model, brm() tells us “Rows containing NAs were excluded from the model.” And indeed, the summary confirms we only used 99 out of our 100 cases.

 Family: MV(gaussian, gaussian) 
  Links: mu = identity; sigma = identity
         mu = identity; sigma = identity 
Formula: a ~ 1 
         b ~ 1 
   Data: dat (Number of observations: 99) 
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
a_Intercept    -0.05      0.10    -0.24     0.14       4134 1.00
b_Intercept    -0.13      0.10    -0.31     0.07       4492 1.00

Family Specific Parameters: 
        Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma_a     0.95      0.07     0.82     1.09       4305 1.00
sigma_b     0.98      0.07     0.85     1.15       4603 1.00

Residual Correlations: 
            Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
rescor(a,b)    -0.08      0.10    -0.28     0.12       4583 1.00

If possible, I would like to use all 100.

I did a little fooling around to try to solve the problem on my own. No luck. For example

fit2 <-
  brm(data = dat,
      mvbind(a, b | mi()) ~ 1)

yields a parsing error.

Error in parse(text = x, keep.source = FALSE) : :2:0: unexpected end of input 1: mvbind(a,b~1 ^

This attempt

fit3 <-
  brm(data = dat,
      bf(mvbind(a, b ) | mi() ~ 1) + set_rescor(TRUE))

gives me a whole bunch of error messages

Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can’t start sampling from this initial value.

that finally end with

[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
error occurred during calling the sampler; sampling not done

And I get the same sad results with approach:

bfa <- bf(a ~ 1)
bfb <- bf(b | mi() ~ 1)

fit4 <-
  brm(data = dat,
      mvbf(bfa, bfb, rescor = TRUE))

Is there a nice solution?

1 Like
2

Thanks for reporting these problems. fit4 should now sample as expected and mvbind(a, b | mi()) ~ 1 should return a clearer error message in the latest github version.

3

Thanks @paul.buerkner! I can confirm, the update works like a dream. And yes, the new error message for fit2 is more informative. Is the fit4 method the only way to pull this off or are there other syntax tricks to combine mi() and multivariate models?

4

fit4 is the intended way at least of you want to selectively specify mi() for some but not for others.

2 Likes
5

Got it.

6

This method still works like a champ.

bfa <- bf(a ~ 1)
bfb <- bf(b | mi() ~ 1)

fit4 <-
  brm(data = dat,
      mvbf(bfa, bfb, rescor = TRUE))

However, I get error messages when I try to compute the waic() or the loo(). Both return:

NAs were found in the log-likelihood. Possibly this is because some of your predictors contain NAs. If you use 'mi' terms, try setting 'resp' to those response variables without missing values.Error in validate_ll(x) : NAs not allowed in input.

Is there a fix?

7

The error message has a typo. It should say:

NAs were found in the log-likelihood. Possibly this is because some of your **responses** contain NAs. If you use 'mi' terms, try setting 'resp' to those response variables without missing values.Error in validate_ll(x) : NAs not allowed in input.

In other words, this is currently not ment to work. You may open an issue if you like so that I don’t loose track of it.

8

Will do.

In addition to fixing the typo, do you see a way to get information criteria from models like this?

9

Yes, this is why I said you may want to open an issue. I fixed the typo already ;-)

10

Cool.