Missing binomial data in a meta-analytic epidemiological model

[Jon Fawcett]

Hey Everyone!

Background

I am developing an epidemiological model to estimate the prevalence of a series of eight separate but related mental disorders. The available data represent published studies measuring two or more of those disorders – although the number of disorders measured by a given study varies (e.g., some studies measure only 2 disorders, a few measure all 8). The estimated prevalence for any given disorder is quite rare in the target population (some are as low as 1%) although there is considerable heterogeneity across estimates. The data itself is structured as follows within R:

head(dat)


#   studyid disorderid samplesize count
# 1       1          1        453     6
# 2       1          2        453     4
# 3       1          3        453     1
# 4       1          4        453     4
# 5       1          5        453    12
# 6       1          6        453    90

In this respect, studyid represents the study identifier (e.g., studyid=1 means Study 1); disorderid represents the disorder identifier (e.g., disorderid=1 means Disorder 1); samplesize represents the size of the relevant sample; and, count represents the number of patients found to suffer from a given disorder in a given sample.

Model Summary

The model applies a binomial likelihood function to each value of count with N equal to samplesize and theta equal to the estimated prevalence for each disorder within a given study. To account for between-study variability, I adopt a random-effects/hierarchical approach whereby the prevalence for a given study is modelled as a multivariate normal distribution with unknown means and a covariance matrix representing the combination of the between-study SDs and the between-study correlations amongst the disorders. I have implemented several optimizations (e.g., using a Cholesky factorization) recommended in the manual to improve convergence. The present model converges quickly with acceptable effective N values. Application to simulated data also reveals this model to perform quite well. The Stan code is appended. If needed, I could provide some R code to simulate data suitable to test the model.

Problem

My question has to do with missing data. Despite the model performing well in its current form, it occurs to me that any disorders unmeasured by a given study represent a case of missing data that should be incorporated into the model. I already estimate the logit-transformed prevalence for each disorder within each study (in the study_estimates parameter) – including the missing values; however, I do not model a corresponding count estimate for those missing values. I understand how this would be accomplished for a continuous dependent measure – but I understand this is not so simple for discrete measures. Instead, the recommendation appears to be marginalizing out the missing data. I was hoping someone could advise me on how to do so for a binomial measure – or point me in the direction of an example dealing with a binomial likelihood function. Just to be clear, I have the samplesize, studyid and disorderid for each missing value. 

Of course, any advice on further improvements (or cautions on things I may have done wrong) are also most welcome. This message board (and the excellent documentation for Stan itself) has already proven invaluable!

Cheers,

Jon

[Bob Carpenter]

What exactly is missing? If samplesize is
known, but count is missing, imputing the counts
won't add any information to the model, as they only
come into the model as end observations without
constraints:

for (i in 1:nobs)
count[i] ~ binomial_logit(samplesize[i], study_estimates[studyid[i]][disorderid[i]]);

More technically, if you include them as discrete variables,
then marginalize them out, the remaining inferences aren't
affected.

But you could generate estimates for the missing counts this way
if you knew their samplesize, studyid and disorderid:

count_missing[i]
<- binomial_rng(samplesize[i],
inv_logit(study_estimates[studyid[i]]
[disorderid[i]]));

I had to really strip your program's comments and parens down
to be able to read it:

data {
int<lower=1> nobs;
int<lower=1> nstudies;
int<lower=1> ndisorders;
int<lower=1,upper=nstudies> studyid[nobs];
int<lower=1> samplesize[nobs];
int<lower=0> count[nobs];
int<lower=1,upper=ndisorders> disorderid[nobs];
}
parameters {
vector[ndisorders] prevalence;
vector[ndisorders] zstudy_estimates[nstudies];
cholesky_factor_corr[ndisorders] L_Omega;
vector<lower=0, upper=pi()/2>[ndisorders] ucauchy_tau;
}
transformed parameters {
vector[ndisorders] study_estimates[nstudies];
vector<lower=0>[ndisorders] tau;
for (i in 1:ndisorders)
tau[i] <- .25 * tan(ucauchy_tau[i]);
for (i in 1:nstudies)
study_estimates[i] <- prevalence + tau .* (L_Omega * zstudy_estimates[i]);
}
model {
ucauchy_tau ~ uniform(0, pi()/2);
prevalence ~ normal(-3.5, 1.5);
for (i in 1:nstudies)
zstudy_estimates[i] ~ normal(0, 1);
for (i in 1:nobs)
count[i] ~ binomial_logit(samplesize[i], study_estimates[studyid[i]][disorderid[i]]);
}


You don't need the uniform on ucauchy_tau --- it's implicit
given the constraint, but it'd won't hurt anything.

I think you can also turn zstudy_estimates into a matrix and
vectorize the study_es4timates calculation. That'd then
make it easy to vectorize the zstudy_estimates variable
using

to_vector(zstudy_estimates) ~ normal(0, 1);

But all that's just for efficiency.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <mv_prevalence_model.stan>

[Michael Betancourt]

As Bob notes, imputation of completely missing observations
doesn’t actually add any information to your model and hence
doesn’t improve your inferences at all. If you really wanted to
see what the model thinks are likely values for the missing
observations then you can just sample them in generated
quantities (once you model the missingness as part of the
measurement model).

[Jon Fawcett]

Thank you for the feedback, Bob and Mike!

I had read somewhere online about using multiple imputation to deal with missing dependent measures – however, I must have misread the context. I am just as happy not to worry about it. In any event, at least your recommendations have made my model even faster!

[Bob Carpenter]

The point is just that sometimes even if you can do perfect
imputation, it doesn't affect your inferences for quantities
of interest. It all depends on what happens when you provide
a model for potentially missing data / outcomes and waht it
looks like when you marginalize out missingness.

- Bob

On Jun 20, 2016, at 5:48 PM, Jon Fawcett <jmfa...@gmail.com> wrote:
>
> Thank you for the feedback, Bob!

> I had read somewhere online about using multiple imputation to deal with missing dependent measures – however, I must have misread the context. I am just as happy not to worry about it. In any event, at least your recommendations have made my model even faster!
>