Answer

How do you test whether event counts vary by month?

Inspired by a question on Cross Validated · October 19, 2026

chi-squaretime serieshypothesis testing

The short answer

Suppose you have several years of monthly counts of some event (hospital admissions, accidents, births) and want to know whether some months have more events than others. The simplest test pools the years, adds up each calendar month, and compares those totals with what you would expect if the daily rate were the same all year. Because months differ in length, the expected share of each month is its number of days divided by the total number of days, with leap-year Februaries counted correctly.

That chi-square goodness-of-fit test is a legitimate test of "no monthly variation". But it asks a very general question: are the 12 months different in any way? If what you suspect is a seasonal wave, with a high season and a low season half a year apart, a test aimed at that shape is more powerful and gives you an effect size (how big the swing is and when it peaks) instead of just a p value.

What the goodness-of-fit test assumes

• Events occur independently, so counts behave like Poisson counts. If events cluster (one incident causing several, or unmeasured year-to-year shocks), the counts vary more than the test expects and its p values are too small. This is called overdispersion.
• The months are unordered categories. The test does not know that June sits next to July, so a smooth rise and fall across the year is spread over 11 degrees of freedom. A modest seasonal wave can easily come out non-significant.
• A long-term trend is handled only roughly. Pooling works because every calendar month contains the same years, so a trend mostly cancels out. A steep trend still makes December slightly larger than January of the same year, which can leak into the monthly comparison.
• Expected counts are large. With several years of data each month's expected total is usually in the hundreds, so the small-count issues covered in chi-square tests with expected counts below 5 rarely arise here.

A better model: Poisson regression with a seasonal cycle

Fit the monthly counts (not the pooled totals) with a Poisson regression. Include log(days) as an offset, so the model describes the daily rate and month length is handled automatically. Add a factor for year (or a smooth function of time) to absorb the trend. Then compare models:

1. Month as a factor (11 df). This is the regression version of the goodness-of-fit test and gives almost the same answer, but it already adjusts for the trend.
2. One cosine and one sine term (2 df). Together they describe one smooth wave per year whose peak can fall in any month. Their combined size gives the amplitude and their ratio the timing of the peak. The idea is the same as in classic seasonality tests such as Edwards' test, and it is far more powerful when the pattern really is one wave a year.
3. A second cosine-sine pair (period six months) if you expect two peaks a year. Compare it with the one-wave model in the same way.
4. Check the dispersion. Divide the Pearson chi-square of the fitted model by its residual degrees of freedom. Values near 1 support the Poisson assumption; clearly above 1, refit with a quasi-Poisson or negative binomial model, which widens the standard errors.
5. Check the residuals over time. If neighbouring months are still correlated after the trend and season are modelled (look at the autocorrelation function of the residuals), standard errors that assume independence are too small.

Tools such as STL decomposition are useful for looking at the data, but they are not tests: with a fixed seasonal window, STL will draw a seasonal component even from pure noise. Seasonal ARIMA models are built mainly for forecasting; for the question "is there seasonality, and how big is it?", the regression approach answers directly. For trend and interruption questions in time series more broadly, see interrupted time series vs event study vs difference-in-differences.

See it in R and Python

The code simulates 12 years (2013 to 2024) of monthly counts with a daily rate of about 2 that rises 3% a year, plus a mild seasonal wave peaking in July, about 6% above the trend at its peak and 6% below at its trough. It then runs the pooled goodness-of-fit test, the month-factor Poisson regression and the cosine-sine regression.

set.seed(144745)
# Monthly event counts for 12 years (2013-2024): an upward trend plus a mild
# seasonal cycle peaking in July (about 6% above and below the trend)
starts <- seq(as.Date("2013-01-01"), as.Date("2025-01-01"), by = "month")
d <- data.frame(year  = as.numeric(format(starts[-145], "%Y")),
                month = as.numeric(format(starts[-145], "%m")),
                days  = as.numeric(diff(starts)))
rate <- 2 * exp(0.03 * (d$year - 2013) + 0.06 * cos(2 * pi * (d$month - 7) / 12))
d$count <- rpois(nrow(d), d$days * rate)

# 1. Chi-square goodness of fit on pooled monthly totals, expected shares set by days
tot  <- tapply(d$count, d$month, sum)
days <- tapply(d$days, d$month, sum)
gof  <- chisq.test(tot, p = days / sum(days))
round(c(n = sum(tot), chisq = unname(gof$statistic), df = unname(gof$parameter), p = gof$p.value), 4)
round(rbind(observed = tot, expected = gof$expected), 0)

# 2. Poisson regression: year effects + month effects, log(days) as offset
lrt <- function(small, big) {
  dev <- deviance(small) - deviance(big); df <- df.residual(small) - df.residual(big)
  round(c(LR = dev, df = df, p = pchisq(dev, df, lower.tail = FALSE)), 4)
}
m0 <- glm(count ~ factor(year) + offset(log(days)), poisson, data = d)
m1 <- glm(count ~ factor(year) + factor(month) + offset(log(days)), poisson, data = d)
lrt(m0, m1)
# Overdispersion check: Pearson chi-square / residual df (near 1 if Poisson fits)
round(sum(residuals(m1, "pearson")^2) / df.residual(m1), 3)

# 3. One smooth yearly cycle (cosine + sine) uses 2 df instead of 11
d$c1 <- cos(2 * pi * d$month / 12); d$s1 <- sin(2 * pi * d$month / 12)
m2 <- glm(count ~ factor(year) + c1 + s1 + offset(log(days)), poisson, data = d)
lrt(m0, m2)
b <- coef(m2)[c("c1", "s1")]
amp <- sqrt(sum(b^2))
round(c(amplitude = amp, peak_to_trough = exp(2 * amp),
        peak_month = unname(atan2(b[2], b[1]) * 12 / (2 * pi)) %% 12), 3)

import numpy as np
from scipy import stats

rng = np.random.default_rng(144745)
# Monthly event counts for 12 years (2013-2024): an upward trend plus a mild
# seasonal cycle peaking in July (about 6% above and below the trend)
year = np.repeat(np.arange(2013, 2025), 12)
month = np.tile(np.arange(1, 13), 12)
starts = np.arange("2013-01", "2025-02", dtype="datetime64[M]").astype("datetime64[D]")
days = np.diff(starts).astype(int)
rate = 2 * np.exp(0.03 * (year - 2013) + 0.06 * np.cos(2 * np.pi * (month - 7) / 12))
count = rng.poisson(days * rate)

# 1. Chi-square goodness of fit on pooled monthly totals, expected shares set by days
tot = np.bincount(month, weights=count)[1:]
mdays = np.bincount(month, weights=days)[1:]
expected = tot.sum() * mdays / mdays.sum()
chi2, p = stats.chisquare(tot, expected)
print({"n": int(tot.sum()), "chisq": round(chi2, 4), "df": 11, "p": float(f"{p:.3g}")})
print(np.round(np.vstack([tot, expected])).astype(int))

# 2. Poisson regression (fitted by IRLS): year effects + month effects, log(days) offset
def poisson_glm(X, y, offset):
    beta = np.zeros(X.shape[1])
    beta[0] = np.log(y.sum() / np.exp(offset).sum())
    for _ in range(100):
        mu = np.exp(X @ beta + offset)
        new = np.linalg.solve(X.T @ (mu[:, None] * X), X.T @ (mu * (X @ beta) + y - mu))
        if np.max(np.abs(new - beta)) < 1e-10:
            break
        beta = new
    mu = np.exp(X @ beta + offset)
    return beta, mu, 2 * np.sum(y * np.log(y / mu) - (y - mu))

def lrt(dev_small, dev_big, df):
    lr = dev_small - dev_big
    return {"LR": round(lr, 4), "df": df, "p": float(f"{stats.chi2.sf(lr, df):.3g}")}

dummies = lambda x: (x[:, None] == np.unique(x)[None, 1:]).astype(float)
off = np.log(days)
X0 = np.column_stack([np.ones(len(year)), dummies(year)])
X1 = np.column_stack([X0, dummies(month)])
_, _, dev0 = poisson_glm(X0, count, off)
_, mu1, dev1 = poisson_glm(X1, count, off)
print(lrt(dev0, dev1, 11))
# Overdispersion check: Pearson chi-square / residual df (near 1 if Poisson fits)
print(round(np.sum((count - mu1) ** 2 / mu1) / (len(count) - X1.shape[1]), 3))

# 3. One smooth yearly cycle (cosine + sine) uses 2 df instead of 11
X2 = np.column_stack([X0, np.cos(2 * np.pi * month / 12), np.sin(2 * np.pi * month / 12)])
b2, _, dev2 = poisson_glm(X2, count, off)
print(lrt(dev0, dev2, 2))
c1, s1 = b2[-2:]
amp = np.hypot(c1, s1)
print({"amplitude": round(amp, 3), "peak_to_trough": round(np.exp(2 * amp), 3),
       "peak_month": round((np.arctan2(s1, c1) * 12 / (2 * np.pi)) % 12, 3)})

The figures below come from one seeded run of the R code. Python's random numbers differ from R's, so its simulated counts and results differ (for example a goodness-of-fit chi-square of 43.75 on 11 df and an estimated peak at month 7.07), but they show the same pattern: all three tests reject, the 2-df test gives the smallest p value, and the dispersion ratio is close to 1.

• Pooled goodness of fit. Across 10,547 events, the monthly totals differed from the day-weighted expectation, χ²(11) = 26.68, p = .005. August had the largest excess (991 observed vs 895 expected) and December the largest shortfall (812 vs 895).
• Month as a factor, adjusting for year. The likelihood-ratio test gave 26.63 on 11 df, p = .005: practically the same as the pooled test, because the trend cancels out when every month contains the same years.
• Dispersion. The Pearson chi-square divided by its residual df was 1.035, so the Poisson assumption looks reasonable for these simulated data. Real event data are often more dispersed than this.
• One seasonal wave (2 df). The likelihood-ratio test gave 16.52 on 2 df, p < .001. That p value is much smaller than the 11-df tests' even though the statistic is smaller, because the evidence is concentrated in 2 degrees of freedom instead of spread over 11.
• Size and timing. The estimated amplitude was 0.056 on the log scale, so the peak rate is about 1.12 times the trough rate (true value 1.13). The peak was estimated at month 6.53, between June and July (true peak: July).

Not sure which test fits your counts? Try the test chooser or browse the DASS blog.

How to report a chi-square goodness-of-fit test in APA style (7th edition)

Give the degrees of freedom and the sample size in parentheses after χ², then the statistic and the exact p value, and say what the expected counts were based on. Using the example above:

"Monthly event counts differed from those expected if the daily rate were constant throughout the year, with expected counts proportional to the number of days in each month, χ²(11, N = 10,547) = 26.68, p = .005."

If you also fit the seasonal regression, report it as a model comparison with its effect size: "Adding a single annual cycle to a Poisson model with year effects improved fit, χ²(2) = 16.52, p < .001; the estimated peak fell between June and July, with the peak-month rate 1.12 times the trough-month rate." Put the observed and expected counts by month in a table or figure, and state in the method section how month length and leap years were handled.

Related tools and guides

• Chi-square power and sample size calculator
• APA 7 formatter for chi-square tests
• Check a reported p-value
• Which statistical test should I use?
• Chi-square test with expected counts below 5
• Is the chi-square test one-tailed or two-tailed?
• Poisson regression (Wikipedia)

More answered questions

• Chi-square test with expected counts below 5: what should you do?
• Is the chi-square test one-tailed or two-tailed?
• How do you compare medians when the groups differ in shape and spread?

Written by AskStats with AI assistance. This is general information, not advice for your specific data or study. When the results matter, check your approach with a qualified statistician.