Generalized linear model

A simulated example, with the code and its output

The numbers below come from made-up data, made by the code on this page. They show the method. They are not output from this tool and not a finding about your data.

Question: When the outcome is a count (orders per customer), can a Poisson generalized linear model recover the effects we built in, and how often does its 95% interval hold the built-in value?

Inputs: Simulated data. 200 customers per dataset. Three predictors: x1 is a standardized spend score (normal), x2 is a 0/1 new-customer flag (half are 1), x3 is a noise column with no effect. The count outcome is drawn from a Poisson distribution whose mean is exp(0.5 + 0.4*x1 - 0.3*x2 + 0*x3). Built-in coefficients (log scale): intercept 0.5, x1 0.4, x2 -0.3, x3 0.0. So x1 and x2 truly matter and x3 truly does not; those results are expected by construction. Seed 1 gives the single worked fit. Seeds 1000 to 1999 (1000 fresh datasets) give the repeat check.

Method: Generalized linear model (Poisson, log link). A generalized linear model links the average outcome to the predictors through a link function; for counts the log link keeps predictions positive, so a coefficient is a change in the log of the expected count. The script fits it by IRLS (iteratively reweighted least squares): start with starting values for the coefficients, compute the fitted means, build a working response, solve a weighted least squares problem, and repeat until the coefficients stop moving. The standard errors come from the inverse of the Fisher information (a measure of how sharply the data pin down each coefficient). A 95% interval is estimate +/- 1.96 standard errors. The script prints one fit, then repeats on 1000 new datasets and counts (a) how often the interval holds the built-in value and (b) how often the estimate differs from zero at the 5% level (a Wald test).

Result: One fit (seed 1): the estimates 0.583, 0.343, -0.340 and -0.015 sit near the built-in 0.5, 0.4, -0.3 and 0, and every interval holds its built-in value. Across 1000 seeds the intervals hold the built-in value 94.5%, 95.3%, 94.6% and 94.7% of the time, close to the 95% the interval promises; the average estimates match the built-in values. Detection rate by simulation: x1 (built-in 0.4) is flagged in 100% of seeds, x2 (built-in -0.3, a smaller effect with a noisier flag) in 72.3%, and the noise column x3 in 5.3%, which is the false-alarm rate expected at the 5% level. Those x1 and x2 detections and the x3 near-5% are true because they were built in. The spread of the estimate (0.058 to 0.118) is the size of change to treat as noise: one dataset's x3 estimate of -0.015 is inside it and is not a finding. This shows the method on made-up data, not what the tool will find in yours.

Limit: The data are clean: the counts really are Poisson, with no extra spread, no missing values and no outliers, so the intervals hold close to 95% by design. Real counts are often more spread out than Poisson, and then a Poisson fit can look more certain than it should be.

Python script (numpy and scipy only):

import numpy as np
from scipy import stats

BETA = np.array([0.5, 0.4, -0.3, 0.0])   # built in: intercept, x1, x2, x3 (x3 has NO effect)
NAMES = ["intercept", "x1 (spend)", "x2 (new customer)", "x3 (noise)"]
N, SEEDS = 200, 1000

def make_data(seed):
    rng = np.random.default_rng(seed)
    x1 = rng.normal(size=N)                  # standardized spend
    x2 = rng.binomial(1, 0.5, size=N)        # 0/1 flag
    x3 = rng.normal(size=N)                  # unrelated to the outcome
    X = np.column_stack([np.ones(N), x1, x2, x3])
    y = rng.poisson(np.exp(X @ BETA))        # counts, e.g. orders per customer
    return X, y

def poisson_irls(X, y, iters=25):
    """Fit a Poisson GLM (log link) by iteratively reweighted least squares."""
    b = np.zeros(X.shape[1]); b[0] = np.log(y.mean() + 0.1)
    for _ in range(iters):
        mu = np.exp(X @ b)                   # fitted mean
        z = X @ b + (y - mu) / mu            # working response
        XtW = X.T * mu                       # weights are mu for Poisson
        b_new = np.linalg.solve(XtW @ X, XtW @ z)
        if np.max(np.abs(b_new - b)) < 1e-10:
            b = b_new; break
        b = b_new
    mu = np.exp(X @ b)
    cov = np.linalg.inv((X.T * mu) @ X)      # inverse Fisher information
    return b, np.sqrt(np.diag(cov))

zc = stats.norm.ppf(0.975)
print("Simulated data. Built-in coefficients (log scale):", dict(zip(NAMES, BETA.tolist())))
X, y = make_data(1)
b, se = poisson_irls(X, y)
print(f"\nSeed 1, n={N}: estimate, 95% interval, holds built-in value?")
for i, nm in enumerate(NAMES):
    lo, hi = b[i] - zc * se[i], b[i] + zc * se[i]
    print(f"{nm:18s} built-in {BETA[i]:5.2f}  est {b[i]:6.3f}  [{lo:6.3f}, {hi:6.3f}]  {lo <= BETA[i] <= hi}")

hold = np.zeros(4); reject = np.zeros(4); est = np.zeros((SEEDS, 4))
for s in range(SEEDS):
    X, y = make_data(1000 + s)
    b, se = poisson_irls(X, y)
    est[s] = b
    hold += (np.abs(b - BETA) <= zc * se)
    reject += (np.abs(b / se) > zc)          # Wald test against zero
print(f"\nAcross {SEEDS} seeds (seeds 1000..{999 + SEEDS}):")
print(f"{'coefficient':18s} {'mean est':>8s} {'sd of est':>9s} {'interval holds':>15s} {'p<0.05 vs 0':>12s}")
for i, nm in enumerate(NAMES):
    print(f"{nm:18s} {est[:, i].mean():8.3f} {est[:, i].std():9.3f} {hold[i] / SEEDS:15.3f} {reject[i] / SEEDS:12.3f}")

Printed output from running the code above:

Simulated data. Built-in coefficients (log scale): {'intercept': 0.5, 'x1 (spend)': 0.4, 'x2 (new customer)': -0.3, 'x3 (noise)': 0.0}

Seed 1, n=200: estimate, 95% interval, holds built-in value?
intercept          built-in  0.50  est  0.583  [ 0.441,  0.724]  True
x1 (spend)         built-in  0.40  est  0.343  [ 0.223,  0.464]  True
x2 (new customer)  built-in -0.30  est -0.340  [-0.569, -0.112]  True
x3 (noise)         built-in  0.00  est -0.015  [-0.124,  0.094]  True

Across 1000 seeds (seeds 1000..1999):
coefficient        mean est sd of est  interval holds  p<0.05 vs 0
intercept             0.493     0.080           0.945        1.000
x1 (spend)            0.401     0.058           0.953        1.000
x2 (new customer)    -0.301     0.118           0.946        0.723
x3 (noise)            0.003     0.058           0.947        0.053

What's in the report?

What drives an outcome, with the model that fits how it is measured: logistic for yes/no, Poisson for counts, gamma for amounts, linear otherwise.