Sampling the dependent variable

27 views
Skip to first unread message

Brian

unread,
Jul 2, 2024, 8:52:28 AMJul 2
to R-inla discussion group
In the spirit of performing the posterior predictive checks, I wish to draw several samples of the dependent variable from the posterior distribution. Is this possible in INLA? 

Finn Lindgren

unread,
Jul 2, 2024, 9:58:41 AMJul 2
to Brian, R-inla discussion group
Hi,
not directly, but you can since it does give you posterior samples of
the linear predictor you can add the observation-level sampling
yourself as a post-processing step.

Finn

On Tue, 2 Jul 2024 at 13:52, Brian <brian...@gmail.com> wrote:
>
> In the spirit of performing the posterior predictive checks, I wish to draw several samples of the dependent variable from the posterior distribution. Is this possible in INLA?
>
> --
> You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/42df8b9f-1b12-4064-91f4-b02b2f7c83cbn%40googlegroups.com.



--
Finn Lindgren
email: finn.l...@gmail.com

Brian

unread,
Jul 2, 2024, 10:08:09 AMJul 2
to R-inla discussion group
Thanks for the response, However I do not understand on 'add the observation-level sampling' part. Assuming the following is my model, and I draw 10 samples, what should follow next?

```
mod <- inla(
  formula,
  family = "binomial",
  Ntrials = numtrials,
  data = inla.stack.data(stk.full),
  control.family = control.family(link = 'logit'),
  control.predictor = list(
    compute = TRUE,
    link = 1,
    A = inla.stack.A(stk.full)
  ),
  control.compute = list(config = TRUE, return.marginals.predictor = TRUE),
  control.inla = list(tolerance = 1e-7),
  num.threads = 7,
  verbose = F
)  

nsamp <- 10
#Posterior sampling
ps <- inla.posterior.sample(nsamp, mod, num.threads = 7, parallel.configs = T)
```

Finn Lindgren

unread,
Jul 2, 2024, 10:14:10 AMJul 2
to R-inla discussion group
Personally,
I'd rewrite it with inlabru, since that eliminates some of the
keeping-manual-track-of indexing. THen I'd do

generate(fit, newdata = ..., formula = ~ rbinom(nrow(.data.), size =
Ntrials, prob = plogis(the + sum + of + model +components)))

That would do the same thing as calling inla.posterior.sample,
extracting/computing the predictor expression for each sample (in a
loop), and then calling rbinom() in a similar way as above.
This simply uses the standard R methods for generating samples.

Finn
> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/8c2b836a-6e80-433b-ada3-50a93daf8a51n%40googlegroups.com.

Håvard Rue

unread,
Jul 2, 2024, 10:15:56 AMJul 2
to Finn Lindgren, R-inla discussion group
something like this

nsamp <- 10
ps <- inla.posterior.sample(nsamp, mod)

ysim <- function(size) {
n <- length(APredictor)
prob <- 1/(1+exp(-APredictor))
return (rbinom(n, prob = prob, size = size))
}
yy <- inla.posterior.sample.eval(
ysim, ps, size = inla.stack.data(stk.full)$numtrials)
Håvard Rue
hr...@r-inla.org

Brian

unread,
Jul 2, 2024, 1:47:07 PMJul 2
to R-inla discussion group
Thanks, this has worked. Since 'yy' is the number of successes, if  I want obtain the probability I can divide the obtained yy with the total number of trials 'numerals'?

Finn Lindgren

unread,
Jul 2, 2024, 1:52:57 PMJul 2
to Brian, R-inla discussion group
To get the simulated _proportion_, yes. Be careful when using the term probability; there are many different probabilities here, and sampling y values isn’t needed to access them! In fact, it’s one of the least efficient methods.
Finn

On 2 Jul 2024, at 18:47, Brian <brian...@gmail.com> wrote:

Thanks, this has worked. Since 'yy' is the number of successes, if  I want obtain the probability I can divide the obtained yy with the total number of trials 'numerals'?

Brian

unread,
Jul 2, 2024, 2:27:13 PMJul 2
to R-inla discussion group
Okay. Thank you. 
Reply all
Reply to author
Forward
0 new messages