Search code examples
rsyntax-errorbayesianmodel-fitting

Error when fitting Bayesian model with the Rethinking R package


I'm trying to fit a very simple model to estimate the prevalence of disease as described here https://www.ecdc.europa.eu/sites/default/files/documents/Methodology-estimating-point-prevalence%20-SARS-CoV-2-infection-pooled-RT-PCR-testing.pdf on page 6, using the Rethinking R package

Here's my code:

quap(alist(
    p ~ dbeta(.3, .3),
    p_test ~ 1 - dbinom(0, s, p), # I tried also p_test <- 1 - dbinom(0, s, p)
    k ~ dbinom(w, p_test)
), data = list(s = 10, k = 30, w = 200))

but I receive the error: Error in pars[[i]] : subscript out of bounds in quap

What am I doing wrong?


Solution

  • Reordering of definitions in alist does the trick.

    quap(
     alist(
        k ~ dbinom(w, p_test),
        p ~ dbeta(.3, .3),
        p_test ~ 1 - dbinom(0, s, p) 
     ), 
     data = list(s = 10, k = 30, w = 200)
    )
    

    This code returns:

    Quadratic approximate posterior distribution
    
    Formula:
    k ~ dbinom(w, p_test)
    p ~ dbeta(0.3, 0.3)
    p_test ~ 1 - dbinom(0, s, p)
    
    Posterior means:
             p 
    0.01575975 
    
    Log-likelihood: -2.55