(Possibly more a question than an issue) Incompatibility of prior_sample and bayes_factor

[EmmanuelCharpentier]

Long justification so the question comes first :

Is it possible to somehow massage a brm fit obtained with sample_prior=TRUE to make it acceptable to bayes_factor (i. e. to bridgesampling, I presume...), or are there fundamental obstacles implying a refit ?

Justification : we may want to be able to do both global model comparisons via bayes_factor AND undirected point hypothesis testing in the "fuller" model, as illustrated in the following example :

Hypothetical situation : analysis of the difference between the means of three groups.

Already known or postulated :

• homoscedasticity,
• normality of within-group distributions.

Accepted "frequentist" analysis plan (at least in medical circles) :

• one-way anova searching for an overall between-groups effect, then, if "significant",
• analysis of intergroup differences by means of ``bilateral'' (i. e. undirected) NHST.

In order to appease a reviewer who is both ignorant and wary of Bayesian analysis, we want to mimic this ``reasoning''. Let's simulate :

> library(brms)
> Sim1 <- data.frame(G=factor(c(rep("A",15), rep("B", 15), rep("C", 15))),
+                    X=c(rnorm(15,0,1), rnorm(15,1,1), rnorm(15,2,1)))
> Full <- brm(X~G, data=Sim1, family=gaussian,
+             prior=prior(student_t(3,0,10), class="b"),
+             save_all_pars=TRUE, sample_prior=TRUE,
+             silent=TRUE, refresh=0)
Compiling the C++ model
Start sampling

Gradient evaluation took 1.6e-05 seconds

[ Snip... ]

Note that Stan still won't shut the fuck up. As far as I know, the problem is at least three years old, has been repeatedly mentioned on the Stan lists/groups/fora and that a solution "in the next version" has been repeatedly promised...

               0.083437 seconds (Total)

> Null <- update(Full, .~.-G, silent=TRUE, refresh=0)
The desired updates require recompiling the model

[ Re-Snip... ]

> summary(Full)
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: X ~ G 
   Data: Sim1 (Number of observations: 45) 
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
Intercept    -0.47      0.26    -0.98     0.05       3347 1.00
GB            1.27      0.37     0.51     2.00       3508 1.00
GC            2.50      0.38     1.75     3.25       3471 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     1.03      0.12     0.83     1.30       4000 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Looks good. Lets do the overall comparison, with information criteria :

> waic(Full,Null)
              WAIC   SE
Full        132.38 8.40
Null        161.88 8.57
Full - Null -29.49 8.63
> loo(Full, Null)
             LOOIC   SE
Full        132.46 8.42
Null        161.90 8.58
Full - Null -29.44 8.64

... and Bayes factor (better known outside statisticians' circles) :

> bayes_factor(Full,Null)
Erreur : Models including prior samples are not usable in method 'bridge_sampler'.

Ouch ! We need to refit.

The problem is that we will assess the overall homogeneity on one fit (which needs prior_sample=FALSE) then test point hypotheses on another fit (which needs prior_sample=TRUE). This is not only a loss of time, but a weak point, open to criticism of Bayesian methods.

Refit :

> CastratedFull <- update(Full, sample_prior=FALSE, silent=TRUE, refresh=0)

[Re-Snip... ]

> CastratedNull <- update(CastratedFull, .~.-G, silent=TRUE, refresh=0)

[Re-Snip... ]

> bayes_factor(CastratedFull, CastratedNull, log=TRUE)
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
Iteration: 5
Iteration: 1
Iteration: 2
Iteration: 3
Iteration: 4
The estimated log Bayes factor in favor of x1 over x2 is equal to: 9.30157

This indicates strong support in favor of the full model, but the numerical values have no obvious relatin with the information criteria.

Note that these information criteria remain consistent with what they were on the model fitted with sample_prior=TRUE :

> waic(CastratedFull, CastratedNull, log=TRUE)
                                WAIC   SE
CastratedFull                 132.32 8.47
CastratedNull                 162.08 8.49
CastratedFull - CastratedNull -29.76 8.61
> loo(CastratedFull, CastratedNull, log=TRUE)
                               LOOIC   SE
CastratedFull                 132.40 8.49
CastratedNull                 162.11 8.50
CastratedFull - CastratedNull -29.71 8.62

Now that we have reassured our reviewer that we had strong evidence of between-groups heterogeneity, we can show him that we may have some evidence about group means differences :

> hypothesis(Full,c("GB=0","GC=0","(GC-GB)=0"))
Hypothesis Tests for class b:
     Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1      (GB) = 0     1.27      0.37     0.51     2.00       0.09      0.08    *
2      (GC) = 0     2.50      0.38     1.75     3.25       0.00      0.00    *
3 ((GC-GB)) = 0     1.23      0.38     0.48     1.95       0.41      0.29    *
---
'*': The expected value under the hypothesis lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

One notes that the apparent `strength of evidence'' is much weaker for the differences between the "median" group and the others than for the difference between the extremes. But I think there is no point in reinventing Newman-Keuls, nor even Tukey'sHSD` : both are nonsense in a Bayesian perspective (as is the NHST stricto sensu, BTW...).

I think that this scenario has some validity, and that not having to refit the same model to switch from bayes_factor to hypothesis would be useful.

What do you think ?

[paul-buerkner]

The reason, this is not allowed at the same time is because of the implementation of sample_prior = "yes". To sample from truncated priors in Stan, we can't use the generated quantitites block, but have to use the MCMC sampler itself. However, sampling from priors this way interferes with the log-posterior and thus with the marginal likelihood computed by the bridgesampling package. That's the only reason both kinds of bayes factors don't work together. As soon as truncated distributions are supported in the generated quantitites block, I can support both BFs for the same fit.

[EmmanuelCharpentier]

Provided you know some analytic form for (a numerical approximation of) your inverse CDF, wouldn't something like this suffice ?

Another "brute force solution is to generate from the whole distribution and reject values out of the accepted region (rejection sampling). This, of course, is costly if the probability of the acceptable region is small. The following example shows how to sample (r2) from a normal truncated to (-mu, mu), where mu is the current value of the mean of a sample passed as data. Running this sample program shows the effect of truncature on r2 (compare plots of r1 and r2, compare their ranges and sds).

HTH,

// GenTrunc.stan : sampling from a truncated distribution as an aside
// from something else ? 

data {
  int       nObs;
  vector[nObs]  X;
}

parameters {
  real              mu;
  real<lower=0>         sigma;
}

model {
  target += student_t_lpdf(mu | 3, 0, 100);
  target += student_t_lpdf(sigma | 3, 0, 100);
  target += normal_lpdf(X | mu, sigma);
}

generated quantities {
  real r1;
  real r2;
  r1=normal_rng(0,1);
  r2=normal_rng(0,1);
  while ((r2< -fabs(mu))||(r2>fabs(mu))) {
    r2=normal_rng(0,1);
  }
}

[paul-buerkner]

We don't have an analytic form of the inverse CDF for all possibe distributions in Stan that we may want to truncate. Thus, this is not a general solution that can be done in brms I think. Even if possible it would require too much special case coding that is not worth the effort of maintaining it.

The brute force options may be worth trying out, though. Let me think about it.

[paul-buerkner]

In the dev version of brms, bayes_factor should now be usable even if sample_prior = "yes".

[EmmanuelCharpentier]

26 hours from suggestion to implementation : Wow... Thanks a lot ! FYI : my (pseudo-) example runs as expected. --Emmanuel Charpentier Le mardi 31 juillet 2018 à 14:56 -0700, Paul-Christian Bürkner a écrit :

In the dev version of brms, bayes_factor should now be usable even if sample_prior = "yes".

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

{"api_version":"1.0","publisher":{"api_key":"05dde50f1d1a384dd78767c5 5493e4bb","name":"GitHub"},"entity":{"external_key":"github/paul- buerkner/brms","title":"paul-buerkner/brms","subtitle":"GitHub repository","main_image_url":"https://assets-cdn.github.com/images/em ail/message_cards/header.png","avatar_image_url":"https://assets- cdn.github.com/images/email/message_cards/avatar.png","action":{"name ":"Open in GitHub","url":"https://github.com/paul-buerkner/brms"}},"u pdates":{"snippets":[{"icon":"PERSON","message":"@paul-buerkner in

485: In the dev version of brms, bayes_factor should now be usable

even if sample_prior = \"yes\"."}],"action":{"name":"View Issue","url":"https://github.com/paul-buerkner/brms/issues/485#issuec omment-409380893"}}} [ { "@context": "http://schema.org", "@type": "EmailMessage", "potentialAction": { "@type": "ViewAction", "target": "https://github.com/paul-buerkner/brms/issues/485#issuecomm ent-409380893", "url": "https://github.com/paul-buerkner/brms/issues/485#issuecomment -409380893", "name": "View Issue" }, "description": "View this Issue on GitHub", "publisher": { "@type": "Organization", "name": "GitHub", "url": "https://github.com" } }, { "@type": "MessageCard", "@context": "http://schema.org/extensions", "hideOriginalBody": "false", "originator": "AF6C5A86-E920-430C-9C59-A73278B5EFEB", "title": "Re: [paul-buerkner/brms] (Possibly more a question than an issue) Incompatibility of prior_sample and bayesfactor (#485)", "sections": [ { "text": "", "activityTitle": "Paul-Christian Bürkner", "activityImage": "https://assets-cdn.github.com/images/email/message cards/avatar.png", "activitySubtitle": "@paul-buerkner", "facts": [

] } ], "potentialAction": [ { "name": "Add a comment", "@type": "ActionCard", "inputs": [ { "isMultiLine": true, "@type": "TextInput", "id": "IssueComment", "isRequired": false } ], "actions": [ { "name": "Comment", "@type": "HttpPOST", "target": "https://api.github.com", "body": "{\n\"commandName\": \"IssueComment\",\n\"repositoryFullName\": \"paul- buerkner/brms\",\n\"issueId\": 485,\n\"IssueComment\": \"{{IssueComment.value}}\"\n}" } ] }, { "targets": [ { "os": "default", "uri": "https://github.com/paul-buerkner/brms/issues/485#issuecomment -409380893" } ], "@type": "OpenUri", "name": "View on GitHub" }, { "name": "Unsubscribe", "@type": "HttpPOST", "target": "https://api.github.com", "body": "{\n\"commandName\": \"MuteNotification\",\n\"threadId\": 362442438\n}" } ], "themeColor": "26292E" } ]