Model Checking

Today

So far, we have mostly asked how to estimate parameters and quantify uncertainty about them.

Today we ask: is the model we fit any good?

Recall what we did for parametric bootstrap

In the parametric bootstrap, we replace the unknown truth \(P_{\theta_0}\) with the fitted model \(P_{\hat\theta}\):

\[ Y^{*(b)}\sim P_{\hat\theta}. \]

This tells us how our estimator would behave if the fitted model were generating the data, which gave us an approximation to the sampling distribution of the estimator.

Today, we’ll use the same sorts of simulated data to check whether the model “fits” the data.

What does it mean for a model to fit?

• Different from getting a precise parameter estimate

• If the bootstrap distribution of \(\hat \theta\) is very narrow, then the interpretation is that the estimation procedure would give similar answers across repeated datasets from \(P_{\hat \theta}\). This doesn’t tell us much about whether the fitted model will reproduce aspects of the data we didn’t use for fitting.

• Model “fit” can be quantified in different ways, but essentially the question is: for any characteristic of the data that we can observe, does the fitted model give values that are roughly in line with the observed values?

The two questions:

The basic simulation experiment

Choose a feature of the data, \(T(Y)\), that we want to check.

1. Fit the model: \[ \hat\theta=h(Y_{\mathrm{obs}}). \]

2. Simulate replicated datasets: \[ Y^{\operatorname{rep},(b)}\sim P_{\hat\theta}. \]

3. Calculate the same feature in every dataset: \[ T_{\mathrm{obs}}=T(Y_{\mathrm{obs}}), \qquad T^{\operatorname{rep},(b)}=T(Y^{\operatorname{rep},(b)}). \]

4. Compare \(T_{\mathrm{obs}}\) with the distribution of the replicated values.

What must stay the same?

The replicated data should reproduce the original data-collection experiment:

• the same sample size;
• the same predictor values or experimental design, when those are treated as fixed;
• the same observation schedule;
• the same missing-data or censoring mechanism, if it is part of the analysis;
• the same preprocessing and calculation of \(T\).

Otherwise, differences between the observed and replicated data might come from changing the experiment rather than from a failure of the model.

What makes a useful check?

A useful check should be:

• Observable: We can calculate it from both the real and simulated data.
• Relevant: It describes something scientifically or practically important.
• Capable of failing: It was not automatically forced to agree during fitting.
• Diagnosable: If it fails, the direction of failure helps us understand what the model is missing.

There is no single universally best check. Different checks probe different possible failures.

Why not check the summaries we fitted?

Suppose we chose \(\hat\theta\) specifically to make

\[ s(Y_{\mathrm{obs}}) \approx \overline{s(Y_{\hat \theta})}. \]

Discovering afterward that these summaries agree is not much evidence for the model. We chose the parameters to make them agree.

A check does not need to estimate the parameter well

For estimation, we want a summary whose distribution changes when \(\theta\) changes.

For model checking, we want a feature that represents a consequence of the model that could be wrong; it doesn’t have to change with the paramater we’re estimating.

Section 1: A Normal location-scale model

Start with our very simple model from before:

\[ Y_i\stackrel{\mathrm{iid}}{\sim}N(\mu,\sigma^2). \]

We are measuring one quantity 500 times, and we fit \(\mu\) and \(\sigma\) to the measurements.

The story this time is that we’re observing some quality of a manufacturing process (maybe the location of a pre-drilled hole in a piece of Ikea furniture).

Most of the time, the hole is in the right place, but a small minority of the time the machine/worker messes up and you see a much wider spread in the measurements.

The observed measurements

Almost all of the readings sit in a tight band. A handful are far away.

The fitted model

Suppose we fit the model somehow and obtain \(\hat\mu\) and \(\hat\sigma\).

It does not matter for what follows how we got them.

The fitted generative model is

\[ Y_i^{\operatorname{rep}}\stackrel{\mathrm{iid}}{\sim} N(\hat\mu,\hat\sigma^2), \]

and we can now simulate as many replicate datasets from it as we like.

Simulate and compare

Two things are wrong at once: the simulated bulk is too wide, and the simulated data never produce a reading anywhere near the far ones.

Another way to see this: Q-Q plot

Sort both datasets. Plot the \(k\)-th smallest simulated value against the \(k\)-th smallest observed value, for every \(k\).

• Each point is one pair of matched order statistics (= (smallest in simulated, smallest in observed), or (second-smallest in simulated, second-smallest in observed), (k-th smallest in simulated, k-th smallest in observed)).
• The dashed line is where simulated \(=\) observed. If the two datasets came from the same distribution, the points scatter around it.
• Points below the line mean the simulated values are smaller than the observed in that region, points above the line mean that the simulated values are larger than the observed in that region.
• At the upper tail (far right of the plot), simulated is substantially lower than observed.
• At the lower tail (far left of the plot), simulated is substantially higher than observed.
• Just below the median (\(\approx\) 97-100), simulated is below observed, just above the median (\(\approx\) 100-103), simulated is above observed.

One simulation against the observed data

The bulk crosses the line steeply, and the extreme observed readings have no simulated counterpart.

Is that just one unlucky simulation?

A one-number summary

The plots show the discrepancy.

In this case, it’s clear by eye that there’s a huge discrepancy between observed and simulated, but if we want to be a bit more formal (or let a computer automate it), a statistic lets us measure it.

How far from the middle is the most extreme reading?

\[ T(Y) = \max_i\left|Y_i-\operatorname{median}(Y)\right|. \]

Simulate the reference distribution

For \(b=1,\ldots,B\):

1. Simulate a dataset of the original size: \[ Y^{\operatorname{rep},(b)} \sim N(\hat\mu,\hat\sigma^2). \]

2. Calculate \[ T^{\operatorname{rep},(b)} =T(Y^{\operatorname{rep},(b)}). \]

Then compare \(T(Y_{\mathrm{obs}})\) with the values the fitted model produces.

The check in code

outlier_distance <- function(y)
    max(abs(y - median(y)))

T_obs <- outlier_distance(y_obs)
T_rep <- replicate(2000,
    outlier_distance(rnorm(
        length(y_obs), mu_hat, sigma_hat)))

def outlier_distance(y):
    return np.max(np.abs(y - np.median(y)))


T_obs = outlier_distance(y_obs)
T_rep = np.array([
    outlier_distance(rng.normal(
        mu_hat, sigma_hat, len(y_obs)))
    for _ in range(2000)])

Everything specific to this check lives in outlier_distance(). The surrounding simulate-and-compare loop is the same in all three examples.

The observed value is far outside anything the fitted model produces.

What did this check tell us?

The failure is not evidence that \(\mu\) or \(\sigma\) was estimated imprecisely.

Changing \(\mu\) moves the distribution and changing \(\sigma\) stretches it. Neither lets a Normal model produce a tight bulk and an occasional reading thirteen units away from the middle.

The fitting method was not the important part

We could have obtained \((\hat\mu,\hat\sigma)\) using:

• closed-form maximum likelihood;
• numerical maximum likelihood;
• method of moments;
• simulated method of moments.

Once the fitted generative model is \(N(\hat\mu,\hat\sigma^2)\), the same simulation-based check is available.

Section 2: Linear regression

Consider the ordinary linear regression model

\[ Y_i\mid x_i \sim N(\beta_0+\beta_1x_i,\sigma^2), \]

independently across observations.

This is a generative model for the responses, conditional on the predictor values.

After fitting \((\hat\beta_0,\hat\beta_1,\hat\sigma)\), we can simulate complete alternative response datasets from it.

Reaction velocity and substrate concentration

The Puromycin data in R record reaction velocity at several substrate concentrations.

We will use one treatment group and start with the linear model

\[ \operatorname{velocity}_i = \beta_0+ \beta_1\operatorname{concentration}_i+ \varepsilon_i, \qquad \varepsilon_i\stackrel{\mathrm{iid}}{\sim}N(0,\sigma^2). \]

Assume we have somehow obtained fitted values for \(\beta_0\), \(\beta_1\), \(\sigma\).

The observed data

Velocity rises with concentration, steeply at first and then much less so.

What experiment should we imitate?

The concentrations were part of the experimental design, so we hold them fixed.

For each replicate dataset, simulate

\[ Y_i^{\operatorname{rep},(b)} \sim N\left( \hat\beta_0+\hat\beta_1x_i, \hat\sigma^2 \right) \]

at exactly the observed values \(x_1,\ldots,x_n\).

We are asking what other response datasets could have occurred under the fitted line in this same experiment.

Simulate and compare

Many simulations at once

Each grey shape is the distribution of simulated velocities at that concentration.

Averaging the runs at each concentration

Each concentration was run twice and averaged, and the red dots are the average of the observed data at each concentration.

One potential discrepancy

Every concentration was run exactly twice, and so we can measure the amount of noise using the absolute difference between the two observations at a given concentration.

\[ T(Y)=\sum_{\text{concentrations }c} \left|Y_{c,1}-Y_{c,2}\right|. \]

We can compare observed and simulated on this statistic.

The pairs check in code

pair_spread <- function(y)
    sum(abs(tapply(y, puro$conc, diff)))

T_obs <- pair_spread(puro$rate)
T_rep <- replicate(2000,
    pair_spread(rnorm(nrow(puro),
        fitted(fit), summary(fit)$sigma)))

def pair_spread(y):
    return sum(abs(np.diff(y[conc == c])[0])
               for c in cs)


T_obs = pair_spread(rate)
T_rep = np.array([
    pair_spread(rng.normal(fitted, sigma_lm))
    for _ in range(2000)])

Two different ways to see one failure

The violin plots and this statistic detect the same defect in different ways:

• the violins show the true mean function bending where the line is straight;
• the pairs show the noise estimate inflated to cover that bend.

In this case, it’s a bit easier to check on the discrepancy in the noise.

How could we improve the model?

The straight-line fit does not explain the data well.

We need some sort of a curve.

A scientifically motivated revision might use a Michaelis-Menten curve,

\[ E(Y\mid x) = \frac{V_{\max}x}{K+x}, \]

rather than simply adding terms until the residual plot looks better.

Familiar diagnostics fit this framework

Several classical regression diagnostics are model checks:

• residuals versus predictors check the conditional mean and variance;
• a Normal Q-Q plot checks the error distribution;
• residual autocorrelation checks independence;
• influential-point diagnostics check whether conclusions depend heavily on a few observations.

Simulation makes the reference comparison explicit and extends the same logic to models where standard diagnostic theory is unavailable.

Section 3: A single-server queue

Return to the queue from Lecture 3.

In an \(M/M/1\) queue:

\[ A_i\stackrel{\mathrm{iid}}{\sim}\operatorname{Exponential}(\lambda), \]

\[ S_i\stackrel{\mathrm{iid}}{\sim}\operatorname{Exponential}(\mu), \]

where \(A_i\) is the time between arrivals, \(S_i\) is a service time, and there is one first-come, first-served server.

We have data and a fitted model

We observe one queue: the arrival times and the waiting times of a long run of customers.

Suppose we have fitted \(\lambda\) and \(\mu\) somehow, giving \(\hat\lambda\) and \(\hat\mu\). As before, how we fitted them does not matter here.

We can again use the fitted model to generate new queues and check whether the simulated queues look like the observed queues.

The observed waiting times

Simulate and compare

Simulate one complete queue from the fitted model (same number of customers, same fitted rates) and put its waiting times beside the observed ones.

The complete check

For \(b=1,\ldots,B\):

1. Simulate interarrival times from \(\operatorname{Exponential}(\hat\lambda)\).
2. Simulate service times from \(\operatorname{Exponential}(\hat\mu)\).
3. Pass them through the queueing simulator.
4. Calculate the chosen summary (in this case mean waiting time across all customers in the queue).

The queue check in code

mean_wait <- function(q) mean(q$waiting)

T_obs <- mean_wait(queue_obs)
T_rep <- replicate(2000,
    mean_wait(simulate_queue(
        n_customers, lambda_hat, mu_hat)))

def mean_wait(q):
    return q["waiting"].mean()


T_obs = mean_wait(queue_obs)
T_rep = np.array([
    mean_wait(simulate_queue(
        n_customers, lambda_hat, mu_hat))
    for _ in range(2000)])

simulate_queue() is the Lecture 3 simulator. Only the choice of mean_wait() is new.

The fitted model expects a mean wait of about 0.74. The observed queue averages 1.57.

What did this check tell us?

• Mainly that the model is insufficient to explain waiting times

• You would need to know a bit more about the model and play around with other possibilities to see how to fix it.

• In this case, it could be fixed by modifying the service-time distribution (the data were actually generated under a model where arrivals were from an exponential distribution but service times were from a distribution with heavier tails).

What does a failure mean?

If the observed feature is unusual under the fitted model, several explanations are possible:

• the structural story is missing an important mechanism;
• a distributional assumption is wrong;
• the observation or sampling model is wrong;
• an important covariate has been omitted;
• the fitting or simulation code contains an error;
• the observed dataset is an unusual realization of an adequate model.

A useful check narrows these possibilities, but rarely identifies a unique replacement model by itself.

What does passing a check mean?

If the observed feature looks ordinary relative to the simulations, then the model has reproduced that particular feature.

It does not follow that:

• the model is true;
• every other feature of the data is reproduced;
• the scientific mechanism is correct.

What we saw

• Uncertainty quantification and model checking answer different questions.
• A generative model makes predictions beyond the quantities used to fit it.
• We check those predictions by simulating complete datasets and calculating the same features in observed and replicated data.
• A failed check is most valuable when it points toward a specific weakness in the model and a plausible revision.