How do I carry out a significance test with Tarone's Z-statistic?

Question

Context

In this blog the author suggests using Tarone's Z-statistic to test for overdispersion in a binomial model to determine whether or not it is necessary to use a beta-binomial model instead. In their example they generate some synthetic data from binomial and beta-binomial distributions and then calculate the Z-statistic's for each and plot them, along with a theoretical curve of the null distribution to demonstrate that this metric works.

Question

How do I actually calculate/use this to test for overdispersion? I found the author code difficult to follow and I don't quite understand how I could use this to formally test for over-dispersion.

I have searched around but all I can turn up about Tarone's Z-statistic are the two links I have included.

I am working in R using the lme4 and glmmTMB packages and I would greatly appreciate an answer in this form. I know this question kind of straddles the bounds of CV and stackoverflow, but I considered this a "non-trivial problem" - If the community disagrees I am happy to migrate it!

Update:

I have managed to adapt C.A. Kapourani's code and write a function for calculating Z-statistics for my models (see below), but I still have the problem of how to compare two values. Is it sensible to find the Z-score having a cumulative probability equal to the critical probability (i.e. 0.05) like with other Z-scores? If so, can anyone recommend how I might do this in R?

taronesZStat <- function(model){

  #Extracting model residuals
  res <- residuals(model)

  #Number of residuals
  n <- length(res)

  #Calculating Tarone's Z-statistic
  p_hat <- sum(res)/n

  S <- sum((res - n * p_hat)^2 / (p_hat * (1 - p_hat)))

  Z_score <- (S - sum(n)) / sqrt(2 * sum(n * (n - 1)))

  return(Z_score)
}

Answer 1

This is an addendum to the other answer that is too long for a comment. The original question regarded overdispersion in generalized linear mixed models models (I assume, given OP is using lme4). In this case there are covariates that will induce variation in excess of what would be expected for iid binomial data. This is actually a primary motivation for modelling overdispersion -- in many fields, we are unable to measure all relevant covariates, and the variation that these induce can distort inference.

Without a good reason to do otherwise I would always incorporate overdispersion into a model (in biology -- my field -- there is almost always unexplained variation due to flawed experimental conditions, the inherent complexity of natural systems, etc). If I wanted to test for overdispersion with an arbitrary test statistic and arbitrary model, I would use a parametric bootstrap. The idea is to simulate data under a null model (e.g. the arbitrary model, but without overdispersion) and calculate the test statistic for each simulation. This gives an approximation to the distribution of the test statistic under the null.

Here's an example with lme4:

library(lme4)
set.seed(101)

# "actual" data
intercept <- 0.1; slope <- 0.8
covariate <- runif(1000)
latent_response <- intercept + covariate*slope
observed_response <- rbinom(1000, prob=plogis(latent_response), size=10)
realdat <- data.frame(m = observed_response, n=rep(10,1000), x = covariate, grp = rep(1:100, each = 10))
null <- glmer(cbind(m, n-m) ~ x + (1|grp), data=realdat, family="binomial")

# parametric bootstrap with function from Ben's answer
sims <- replicate(1000,{ fakedat <- simulate(null)[[1]]; Tarone.test(rowSums(fakedat), fakedat[,1])$statistic })
obs <- Tarone.test(realdat$n, realdat$m)$statistic
# approximate one-tailed pval
sum(sims > obs)/length(sims)
hist(sims);abline(v=obs,col="red") #no evidence for overdispersion

# what if covariate was not measured?
null2 <- glmer(cbind(n, n-m) ~ (1|grp), data=realdat, family="binomial")
sims2 <- replicate(1000,{ fakedat <- simulate(null2)[[1]]; Tarone.test(rowSums(fakedat), fakedat[,1])$statistic })
sum(sims2 > obs)/length(sims2)
hist(sims2);abline(v=obs,col="red") #now appears to be evidence for overdispersion

Answer 2

If you would like another explanation of the procedure, you can read the original paper by Tarone (1979), but the blog actually gives a longer and clearer explanation than the original paper. In any case, the function you have programmed does not correspond to the formula for Tarone's Z statistic. I would recommend writing a function that directly calculates the Z statistic from the data, rather than taking an input from a model. To do this, let's start by generating some binomial data.

#Generate example data
N <- c(30, 32, 40, 28, 29, 35, 30, 34, 31, 39);
M <- c( 9, 10, 22, 15,  8, 19, 16, 19, 15, 10);

We can program the Tarone test as a function in R so that it takes the input data directly and outputs all the required elements of a htest object. (This form of function means that the object will contain all the relevant information for the hypothesis test, and it will print in the default mode for a htest object.)

Tarone.test <- function(N, M) {

    #Check validity of inputs
    if(!(all(N == as.integer(N)))) { stop("Error: Number of trials should be integers"); }
    if(min(N) < 1) { stop("Error: Number of trials should be positive"); }
    if(!(all(M == as.integer(M)))) { stop("Error: Count values should be integers"); }
    if(min(M) < 0) { stop("Error: Count values cannot be negative"); }
    if(any(M > N)) { stop("Error: Observed count value exceeds number of trials"); }

    #Set description of test and data
    method      <- "Tarone's Z test";
    data.name   <- paste0(deparse(substitute(M)), " successes from ", 
                          deparse(substitute(N)), " trials");

    #Set null and alternative hypotheses
    null.value  <- 0;
    attr(null.value, "names") <- "dispersion parameter";
    alternative <- "greater";

    #Calculate test statistics
    estimate    <- sum(M)/sum(N);
    attr(estimate, "names") <- "proportion parameter";
    S           <- ifelse(estimate == 1, sum(N),
                          sum((M - N*estimate)^2/(estimate*(1 - estimate))));
    statistic   <- (S - sum(N))/sqrt(2*sum(N*(N-1))); 
    attr(statistic, "names") <- "z";

    #Calculate p-value
    p.value     <- 2*pnorm(-abs(statistic), 0, 1);
    attr(p.value, "names") <- NULL;

    #Create htest object
    TEST        <- list(method = method, data.name = data.name,
                        null.value = null.value, alternative = alternative,
                        estimate = estimate, statistic = statistic, p.value = p.value);
    class(TEST) <- "htest";
    TEST; }

Now let's try applying this function to the example data and see what we get:

#Apply Tarone's test to the example data
TEST <- Tarone.test(N, M);
TEST;

        Tarone's Z test

data:  M successes from N counts
z = 2.5988, p-value = 0.009355
alternative hypothesis: true dispersion parameter is greater than 0
sample estimates:
proportion parameter 
           0.4359756 

In this case we have a low p-value, so we find strong evidence to reject the null hypothesis of a binomial distribution in favour of the alternative of over-dispersion. We can fit the data to the beta-binomial model to estimate the dispersion parameter:

#Fit the example data to a beta-binomial distribution
MODEL <- VGAM::vglm(cbind(DATA$m, DATA$n-DATA$m) ~ 1, betabinomial, trace = FALSE);
Coef(MODEL);

        mu        rho 
0.43507150 0.03355549