I am trying to implement one of the methods described in Valen Johnson's A Bayesian Chi-Squared Test for Goodness of Fit. It presents a couple of variants depending on whether the random variable of interest is continuous or discrete but I am specifically interested in a binomial outcome.

The central idea is that the proposed $R^B$ statistic's posterior approaches a $\chi^2_{K-1}$ distribution, where $K$ is the number of discrete values a variable can take. It is defined as

$$R^B(\tilde\theta) = \sum_{k=1}^K \left[{m_k - n p_k(\tilde\theta) \over \sqrt{np_k(\tilde\theta)}}\right]^2$$

where $\tilde\theta$ is a single posterior draw from the parameter vector, $p_k$ is the expected count, calculated over the $n$ observations as

$$p_k(\tilde\theta) = {1 \over n}\sum_{j=1}^n \sum_{y \in \text{bin} k} f_j(y\mid\tilde\theta).$$

The above notation is directly transcribed from the paper, but the notion of bins is irrelevant for my binomial scenario, so a slightly clearer way of denoting this is

$$p_k(\tilde\theta) = {1 \over n}\sum_{j=1}^n f(k -1\mid\tilde\theta_j),$$

since there's only one possible $y$ value at each level $k$ and counting starts at $0$. Also, I shifted the $j$ subindex from the pmf $f$ to the parameter $\tilde\theta$, as the pmf has a fixed functional form, but parameters can be observation-dependent (e.g. the mean in a regression model).

Finally, we have $m_k$, which corresponds to observed counts. Using $I(.)$ to denote the indicator function and $a_k$ for the corresponding quantile, we have

$$m_k(\tilde\theta)=\sum_{j=1}^n I(F(y_j\mid\tilde\theta_j) \in (a_{k-1}, a_k]).$$

For reference, the equations above are numbered $(2)$ through $(5)$ in the paper.

Having implemented this measure for a simple intercept-only logistic regression model in R, the distribution is far from what the paper says it should look like. Here's the code:


# Calculate R^B statistic for single posterior draw
iRB <- function(b, n, data) {

  y <- data*n
  p <- (1/(1+exp(-rep(b, length(y)))))

  pmf <- function(X) dbinom(X, n, p)
  cdf <- function(X) pbinom(X, n, p)

  rbk <- lapply(1:(n+1), function(k) {
      # eq (5)
    pk <- sum(pmf(k-1))
      # eq (4)
    Fy <- cdf(y)
      # eq (2)
    ak <- cdf(k-1.1);aK <- cdf(k-1)
    mk <- sum(ifelse(ak <Fy&Fy<= aK, 1, 0))

    data.frame(pk = pk, mk = mk)
    }) %>% do.call(rbind, .)

  with(rbk,sum(((mk - pk)/sqrt(pk))**2))

# Simulate data
m <- 7
set.seed(1);binomdat <- data.frame(y=rbinom(100, m, 0.5)/m, m  = m)
# Fit intercept-only logistic regression
binomfit <- stan_glm(y ~ 1, family=binomial(), data=binomdat, weights = m)
# Extract posterior
ps <- as.matrix(binomfit)
# Calculate R^B for each posterior draw
chi2b <- sapply(ps, iRB, m, binomdat$y)

# Check results
curve(dchisq(x, m), from = 0, to = 80)
hist(chi2b, probability = T, add = T)

I've already gone over this with a professor and we're both perplexed by the results. Not sure if we're misreading the paper or overlooking an error in the implementation.


1 Answer 1

  • The value/distribution that you have been calculating is the posterior distribution of $R^B$ (based on the posterior distribution of the parameter $\theta$ determined from a single sample/test $y$) but this is not chi-square distributed.

    When $n$ grows the posterior distribution for $\theta$ will approach some degenerate distribution around the true parameter value, and $R^B$ will be as well a degenerate distribution for a fixed value/sample/test $y$. This is like the sharp peak that you observe in your histogram, it is the chi-square statistic for a single sample/draw $y$, where you wiggle a little bit the predicted value $\theta$.

  • The theorem 1 from the article relates to $R^B$ as a frequentist concept. It is the sample distribution of $R^B$ (based on multiple samples/tests $y$) which is chi-square distributed.

    E.g. if you do the test many times and take a single posterior sample $\tilde{\theta}$, for each test, from which you calculate $R^B(\tilde{\theta})$, then you will have $R^B \xrightarrow[\text{}]{\text{d}} \chi^2_{K-1}$.

output image

#do this N times
N <- 500
chi2b <- rep(0,N)

pb <- txtProgressBar(title = "progress bar", min = 0,
                     max = N, style=3)
for (ii in 1:N) {
  binomdat <- data.frame(y=rbinom(100, m, 0.5)/m, m  = m)
  # Fit intercept-only logistic regression  # refresh = 0 suppresses output
  binomfit <- stan_glm(y ~ 1, family=binomial(), data=binomdat, weights = m, refresh= 0)
  # Extract SINGLE posterior
  ps <- as.matrix(binomfit)[1]
  # Calculate R^B for each posterior draw
  chi2b[ii] <- sapply(ps, iRB, m, binomdat$y)
  setTxtProgressBar(pb, ii)

# Check results
curve(dchisq(x, m), from = 0, to = 80)
hist(chi2b, probability = T, add = T, breaks = seq(0,80,1))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.