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.


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!


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)))


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:


# "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

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);

        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 

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);

        mu        rho 
0.43507150 0.03355549
  • 1
    $\begingroup$ Good answer ... but this procedure doesn't involve covariates. If there are covariates (as in the typical binomial glm). With even a slightly strong signal, these will induce overdispersion (as covariates are ignored and the test is used on the raw data). Why not use a parametric bootstrap to estimate the distribution of the test statistic? $\endgroup$ – Nate Pope May 28 '19 at 0:45
  • $\begingroup$ For this answer I am just going by the simpler application in the blog that was linked by the OP, which is the application applied in Tarone (1979). No doubt there are other methods for other types of question. $\endgroup$ – Ben May 28 '19 at 0:50
  • $\begingroup$ @NatePope how would I incorporate parametric boostrapping to estimate the distribution of the test stats? This is more what I was after but hadn't realised that covariates could induce over-dispersion before your comment. $\endgroup$ – André.B May 28 '19 at 0:53
  • $\begingroup$ @ReinstateMonica I've pulled your function in to the stackoverflow package github.com/nfultz/stackoverflow/blob/master/R/Tarone.R - I've cited you as 'ReinstateMonica' but feel free to message me if you'd prefer to be listed in a different way. Thanks for providing a clear, reusable function in your answer! $\endgroup$ – Neal Fultz Dec 22 '19 at 19:04
  • $\begingroup$ @Neal Fultz: Thanks for that --- I hope the code is useful to others. My username is temporarily set to that as part of a show of solidarity with another SE user, but my actual name is Ben O'Neill. If you could cite my code by my actual name, that would be preferable. Thanks heaps. $\endgroup$ – Ben Dec 23 '19 at 1:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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