Meta-analysis models: Binomial vs Normal?

[Burak Gunhan]

Dear Stan experts,

I am using Stan to fit a fixed effect meta-analysis model. I took Stan code for the non-hierarchical model from Chapter 11.2 of the Stan manual. And I also want to use binomial non-hierarchical model. I attached my two Stan codes (NM.stan and BM.stan).

To compare those two models, I simulated data based on the binomial non-hierarchical model (I attached the Stan code, BM_sim.stan). Then I analyzed the simulated data using normal non-hierarchical model and binomial non-hierarchical model. And then I compared them using loo package based on WAIC results. WAIC says that normal model is preferable, but I dont understand, how come true model is not preferable?

Maybe it is not a Stan-related question but a simple statistical question, but how come normal model performs better than binomial model when I simulate data from binomial model? At what point I am doing wrong?

I also simulated 1000 dataset. As a result, bias that I obtained from Normal model is smaller than binomial model.

Here is my Sessioninfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.1 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=de_DE.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=de_DE.UTF-8  
 [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=de_DE.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C           
[11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C      

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base    

other attached packages:
[1] loo_0.1.6          rstan_2.11.1       StanHeaders_2.11.0 ggplot2_2.1.0    

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.6         matrixStats_0.50.2  codetools_0.2-14    lattice_0.20-33     digest_0.6.10       grid_3.3.1        
 [7] plyr_1.8.4          gtable_0.2.0        stats4_3.3.1        scales_0.4.0        INLA_0.0-1470222140 sp_1.2-3          
[13] Matrix_1.2-6        labeling_0.3        splines_3.3.1       tools_3.3.1         munsell_0.4.3       parallel_3.3.1    
[19] rsconnect_0.4.3     inline_0.3.14       colorspace_1.2-6    gridExtra_2.2.1   

Thank you very much for your help,

Best,

Burak

[Bob Carpenter]

I don't understand the transforms you're doing here or why the
normal model has only one parameter whereas the binomial model has
N + 1 parameters. Or why you model rc and rt in the binomial
model and not in the normal model. I may be missing something,
but the models don't look comparable. Also, you have an inverse logit
transform in the binomial model and I don't see anything similar in
the normal model (but may be missing it).

I'm not an expert on WAIC, but it's not the gold standard in model
comparison. If the models are correct, you might want to try more
robust predictive evaluations (like predictive p-values or simple
cross-validation or held-out evaluation).

- 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.
> <BM.stan><NM.stan><Meta_analysis_stan.R><BM_sim.stan>

[Burak Gunhan]

Hi Bob,

As you described in the Chaper 11.2 of the Stan manual, meta-analyses based on binomial data can be modelled using normal model (using log-odds ratio transformation) or using binomial model (using binomial structure of data directly). I took the binomial stan code from one of your BUGS-translation: https://github.com/stan-dev/example-models/blob/master/bugs_examples/vol1/blocker/blocker.stan.

I am interested in sparse datasets. And since normal modeling approach using approximation, I want to show that binomial modelling give better results compared to normal modeling. This is done by here  http://www.ncbi.nlm.nih.gov/pubmed/20827667

I dont think it is becuse of WAIC, but maybe because they are not comparable models. But still I want to compare them?

[Bob Carpenter]

I didn't write that BUGS example. Not sure I wrote that part of
the manual either --- it's a group effort. Also, I'm not the
WAIC or model comparison expert, so hopefully someone else will
jump in.

- Bob

[Jonah Gabry]

I didn't have a chance to go through everything, so I'm not sure if this is the only issue, but I think there's a mistake in the generated quantities block of the binomial model. You have this:

parameters {

  real d;

  vector[N] mu;

}
...

generated quantities {
  vector[N] log_lik;

  for (n in 1:N)
    log_lik[n] = binomial_logit_log(rt[n], nt[n], mu + d);
}

but mu is declared as an N-vector in the parameters block. So I think you want mu[n] + d in the call to binomial_logit_log. I'm surprised this didn't result in an error when you ran it. Bob, is there a reason why there wouldn't be a runtime error if N > 1?

Jonah

[Bob Carpenter]

> On Aug 20, 2016, at 12:33 AM, Jonah Gabry <jga...@gmail.com> wrote:
>
> I didn't have a chance to go through everything, so I'm not sure if this is the only issue, but I think there's a mistake in the generated quantities block of the binomial model. You have this:
>
> parameters {
> real d;
> vector[N] mu;
> }
> ...
> generated quantities {
> vector[N] log_lik;
>
> for (n in 1:N)
> log_lik[n] = binomial_logit_log(rt[n], nt[n], mu + d);
> }
>
> but mu is declared as an N-vector in the parameters block. So I think you want mu[n] + d in the call to binomial_logit_log.

Good catch.

> I'm surprised this didn't result in an error when you ran it. Bob, is there a reason why there wouldn't be a runtime error if N > 1?

Yes, vectorization. We will take any of the
arguments as vectors/arrays or scalars as long as
all the arguments either have the same size container
or are scalars. I don't know how to warn users about
this, but this isn't the only time it's come up. Reminds
me of the some of the "clever" things R does that can lead
to errors, like broadcasting up (repeating entries) to fill
arrays.

- Bob

[Jonah Sol Gabry]

On Friday, August 19, 2016, Bob Carpenter <ca...@alias-i.com> wrote:

>  I'm surprised this didn't result in an error when you ran it. Bob, is there a reason why there wouldn't be a runtime error if N > 1?

Yes, vectorization.  We will take any of the
arguments as vectors/arrays or scalars as long as
all the arguments either have the same size container
or are scalars.  I don't know how to warn users about
this, but this isn't the only time it's come up.  Reminds
me of the some of the "clever" things R does that can lead
to errors, like broadcasting up (repeating entries) to fill
arrays.

Ah yes, thanks. I'm just not used to seeing the first argument get recycled I think. 

Jonah

[Bob Carpenter]

We could potentially remove that feature, though I'd rather
wait for Stan 3, because it'll break backward compatibility.

I added a note to flag this in "pedantic mode".

- Bob