Parametric Bootstrap for Multivariate Likelihoods

Why Bootstrap a Multivariate Likelihood?

Likelihood-ratio tests often rely on a large-sample reference distribution, such as a chi-square distribution. That approximation can be inaccurate when the sample is modest, the model is complicated, or nuisance parameters strongly affect the statistic’s finite-sample distribution. The parametric bootstrap offers a targeted alternative: fit the model under the null, simulate data from that fitted model, and repeat the same fitting and testing procedure on every simulated dataset.

The refitting step is essential. A simulated statistic must reflect the estimation that would occur in a new dataset, not merely evaluate a statistic with the original data’s fitted parameters held fixed. The method does not make a model assumption-free: its calibration is conditional on the chosen parametric model being a reasonable description of the data.

Definition: Let \(Y_1,\ldots,Y_n\) be independent observations from a multivariate model with density or mass function \(f_\theta\), where \(\theta=(\psi,\lambda)\). Here \(\psi\) is the parameter of interest and \(\lambda\) is a nuisance parameter. To test \(H_0:\psi=\psi_0\), a parametric bootstrap fits the model under \(H_0\), simulates new datasets from that fitted null model, and recomputes the chosen test statistic after refitting each dataset.

The Bootstrap Likelihood-Ratio Procedure

Write \(\ell(\theta)\) for the log-likelihood. Let \(\hat\theta\) maximize it over the full parameter space, and let \(\hat\theta_0=(\psi_0,\hat\lambda_0)\) maximize it subject to the null hypothesis. The likelihood-ratio statistic is

T_{\mathrm{obs}} =2\{\ell(\hat\theta)-\ell(\hat\theta_0)\}.

A large value is evidence against the null. For bootstrap calibration, simulate each dataset from \(f_{\hat\theta_0}\), not from the unrestricted fit \(f_{\hat\theta}\). For replicate \(b\), refit both the unrestricted and restricted models and calculate

T_b^* =2\{\ell_b^*(\hat\theta_b^*)-\ell_b^*(\hat\theta_{0,b}^*)\}.

The bootstrap \(p\)-value is commonly estimated by

\hat p_{\mathrm{boot}} =\frac{1+\sum_{b=1}^B I(T_b^*\ge T_{\mathrm{obs}})}{B+1}.

The added one prevents a reported value of zero when the number \(B\) of replicates is finite. It does not make a plug-in parametric bootstrap an exact test: the data and bootstrap samples are not generally exchangeable, because the null nuisance parameter was estimated from the observed data.

Fit under the null.
Find \(\hat\lambda_0\), keeping \(\psi=\psi_0\). This restricted fit supplies the bootstrap data-generating parameter.

Simulate complete datasets.
Generate \(Y_1^*,\ldots,Y_n^*\) from \(f_{(\psi_0,\hat\lambda_0)}\), preserving the sample size and the model’s multivariate structure.

Repeat the analysis.
For each dataset, refit the restricted and unrestricted models and calculate the same statistic used on the observed data.

Calibrate the result.
Compare the observed statistic with the simulated null distribution, or use its upper quantile as a critical value.

Worked Examples

Worked Example: Testing a Multivariate Mean With Known Covariance

Suppose \(Y_i\sim N_2(\mu,\Sigma)\), independently, with known \(\Sigma=\begin{pmatrix}1&0\\0&4\end{pmatrix}\). Test \(H_0:\mu=(0,0)^\mathsf{T}\) using \(n=4\) observations whose sample mean is \(\bar Y=(1,2)^\mathsf{T}\).

For known covariance, the unrestricted estimate is \(\hat\mu=\bar Y\), and under the null the mean is fixed at zero. The likelihood-ratio statistic reduces to

T_{\mathrm{obs}} =n(\bar Y-\mu_0)^\mathsf{T}\Sigma^{-1}(\bar Y-\mu_0).

Here \(\Sigma^{-1}=\operatorname{diag}(1,1/4)\), so

T_{\mathrm{obs}} =4\left\{1^2+\frac{2^2}{4}\right\} =4(2)=8.

The null bootstrap simulates four independent \(N_2((0,0)^\mathsf{T},\Sigma)\) observations per replicate. It then calculates the replicate sample mean and applies the same quadratic form. In this example the statistic has exactly a \(\chi^2_2\) distribution under the null, so simulation is not needed for accuracy; it illustrates the procedure and provides a check on an implementation. A simulated mean such as \((0.5,-1)^\mathsf{T}\) gives \(T^*=4\{0.5^2+(-1)^2/4\}=2\). A parametric bootstrap would compare 8 with the distribution of these replicate values.

Worked Example: Testing Zero Correlation With Unknown Marginal Parameters

For \(n=20\) independent bivariate normal observations, test whether the population correlation is zero. Both means and both variances are unknown nuisance parameters. Suppose the observed sample correlation is \(r=0.40\).

Under the unrestricted model, the fitted covariance matrix may have a nonzero off-diagonal entry. Under the null it must be diagonal, while the two means and variances are still estimated. For this bivariate normal comparison, profiling the means and variances gives

T_{\mathrm{obs}}=-n\log(1-r^2) =-20\log(1-0.40^2) =-20\log(0.84) \approx 3.487.

For each bootstrap replicate, estimate the two means and variances under zero correlation. Generate 20 independent pairs from the resulting fitted diagonal-covariance normal distribution. Then compute the replicate correlation \(r_b^*\), refit both models, and calculate \(T_b^*=-20\log\{1-(r_b^*)^2\}\). For example, if a replicate has \(r_b^*=0.20\), then

T_b^*=-20\log(0.96)\approx 0.816.

The bootstrap distribution reflects the estimated marginal scales and the finite-sample fitting procedure. It is not obtained by simulating from the unrestricted fit, which would build the observed nonzero correlation into the simulated null datasets.

Worked Example: Inverting Bootstrap Tests for a Mean Confidence Set

Suppose \(Y_i\sim N_2(\mu,\Sigma)\) with both the mean and covariance unknown, and \(n=10\). Consider testing a candidate mean \(\mu_0\). Use the maximum-likelihood covariance estimate with divisor \(n\), and suppose the observed sample mean is \(\bar Y\), the unrestricted covariance estimate is \(S=I_2\), and \(\bar Y-\mu_0=(0.2,0)^\mathsf{T}\).

With the mean fixed at \(\mu_0\), the residual covariance estimate is \(S_0=S+n(\bar Y-\mu_0)(\bar Y-\mu_0)^\mathsf{T}\). Thus

S_0= \begin{pmatrix}1&0\\0&1\end{pmatrix} +10\begin{pmatrix}0.04&0\\0&0\end{pmatrix} =\begin{pmatrix}1.4&0\\0&1\end{pmatrix}.

For this model the profiled likelihood-ratio statistic is \(T_{\mathrm{obs}}=n\log\{|S_0|/|S|\}\). Since \(|S_0|=1.4\) and \(|S|=1\),

T_{\mathrm{obs}}=10\log(1.4)\approx 3.365.

To test this candidate mean, simulate datasets from \(N_2(\mu_0,S_0)\), the fitted model under the candidate null. For each dataset, recompute its sample mean and covariance, refit the covariance under the candidate mean constraint, and calculate the likelihood-ratio statistic. If, for illustration, five simulated statistics in increasing order are \(0.8,1.7,2.6,3.1,4.2\), the empirical upper-tail fraction at 3.365 is \(1/5\); with the plus-one convention the estimated \(p\)-value is \(2/6\). In practice many more replicates are needed for a stable tail estimate.

Repeat this test for candidate values of \(\mu_0\). The confidence set consists of candidate means not rejected at the chosen level. The restricted covariance fit must be recalculated for each candidate; using the original candidate’s fitted covariance for all tests would not invert the stated procedure.

Score Statistics and Nuisance Parameters

The same resampling logic applies to a score test. Let \(U_\psi\) and \(U_\lambda\) be the score vectors, and partition the Fisher information into blocks \(I_{\psi\psi}, I_{\psi\lambda}, I_{\lambda\lambda}\). When \(I_{\lambda\lambda}\) is nonsingular, the efficient score for \(\psi\), accounting for estimation of \(\lambda\), is

U_{\psi\cdot\lambda} =U_\psi-I_{\psi\lambda}I_{\lambda\lambda}^{-1}U_\lambda, \qquad I_{\psi\cdot\lambda} =I_{\psi\psi}-I_{\psi\lambda}I_{\lambda\lambda}^{-1}I_{\lambda\psi}.

A score statistic uses \(U_{\psi\cdot\lambda}^{\mathsf T}I_{\psi\cdot\lambda}^{-1}U_{\psi\cdot\lambda}\), evaluated at the restricted fit. In each bootstrap replicate, estimate nuisance parameters under the null and recompute the score and information there. Holding the observed-data information fixed generally fails to reproduce the fitting variability the bootstrap is intended to capture.

Why the Calibration Works Asymptotically

Theorem: Suppose the observations are independent and identically distributed from a regular, identifiable, finite-dimensional parametric model. Assume the true parameter is an interior point, the Fisher information is nonsingular, the log-likelihood is sufficiently smooth for a uniform local quadratic expansion, and the unrestricted and restricted maximum-likelihood estimators are consistent. If \(H_0\) is a smooth restriction of codimension \(q\), then under \(H_0\), \(T_{\mathrm{obs}}\xrightarrow{d}\chi^2_q\). If bootstrap samples are generated from the restricted estimate and both models are refitted, the conditional bootstrap distribution also converges in probability to \(\chi^2_q\).

The key calculation is the local quadratic likelihood expansion. Write \(U_n\) for the score at the true parameter and \(I\) for the per-observation information. For parameter displacements of order \(n^{-1/2}\), regularity gives

\ell_n(\theta_0+h/\sqrt n)-\ell_n(\theta_0) =h^\mathsf{T}\frac{U_n}{\sqrt n} -\frac12h^\mathsf{T}Ih+o_p(1).

The central limit theorem gives \(U_n/\sqrt n\xrightarrow{d}N(0,I)\). Maximizing this quadratic over the full local parameter space and over the null subspace leaves a difference equal to the squared length of the score’s component in the \(q\) excluded directions, after scaling by information. That component has \(q\) independent standard normal coordinates in the limit. Its squared length therefore has a \(\chi^2_q\) distribution, proving the first assertion.

For the bootstrap, consistency of the restricted estimator gives \(\hat\theta_0\xrightarrow{p}\theta_0\). Conditional on the observed data, the bootstrap observations are independent draws from \(f_{\hat\theta_0}\). The same local expansion and central limit argument apply to this sequence of fitted parameters; the assumed uniformity of the expansion near \(\theta_0\) makes the conditional remainder negligible in probability. The bootstrap score has conditional limiting covariance \(I\), and refitting again projects away the \(q\) null restrictions. Thus the conditional likelihood-ratio distribution converges to \(\chi^2_q\). \(\square\)

Consequence: Under these conditions, a bootstrap critical value \(c_{n,1-\alpha}^*\) for the upper \(1-\alpha\) quantile satisfies \(c_{n,1-\alpha}^*\xrightarrow{p}c_{1-\alpha}\), where \(c_{1-\alpha}\) is the \(\chi^2_q\) quantile. Since the limiting distribution is continuous, the test that rejects when \(T_{\mathrm{obs}}>c_{n,1-\alpha}^*\) has asymptotic rejection probability \(\alpha\) under the null.

To see the size conclusion, fix \(\varepsilon>0\). With probability tending to one, the random bootstrap cutoff lies between \(c_{1-\alpha}-\varepsilon\) and \(c_{1-\alpha}+\varepsilon\). The rejection probability is therefore bounded by the probabilities that \(T_{\mathrm{obs}}\) exceeds these two fixed values, up to a term tending to zero. Taking limits and then letting \(\varepsilon\) decrease to zero gives \(P(T_{\mathrm{obs}}>c_{1-\alpha})=\alpha\), by continuity of the chi-square distribution.

Common Pitfalls and Practical Limits

Simulating from the unrestricted fit. This does not generate data under the null. Simulate from the restricted fit, including its nuisance-parameter estimates.
Not refitting each replicate. Reusing the observed-data estimates changes the statistic and typically understates the variability caused by estimation.
Assuming the plug-in bootstrap is exact. Estimating nuisance parameters makes this a plug-in calibration, usually justified asymptotically under regularity conditions rather than by exact exchangeability.
Ignoring model failure. A parametric bootstrap reproduces features of the assumed model, not unmodeled dependence, outliers, or misspecification. The independent-observation setup here must not be used when the sampling design has dependence that the model omits.
Using ordinary chi-square theory at a boundary. If a variance component is zero under the null, or another parameter lies on the boundary, the regular theorem above does not apply. The bootstrap must reproduce the constrained fitting and boundary behavior; even then, its accuracy requires care.
Using too few replicates for a tail probability. With \(B\) replicates, the estimated tail probability changes in increments of \(1/(B+1)\). Small \(B\) makes small \(p\)-values and extreme quantiles especially imprecise.

The central idea is simple but exacting: define the statistic, fit the null model, simulate from that fitted null, and repeat the full computation. The method is especially useful when multivariate nuisance parameters make a standard reference distribution doubtful, provided the fitted model and the bootstrap’s regularity assumptions are credible.

Check Your Understanding

Work through these questions, paying attention to what is estimated and what is simulated.

For a likelihood-ratio test of \(H_0:\psi=\psi_0\), why are bootstrap datasets generated from the restricted fit rather than the unrestricted fit?
In each bootstrap replicate, which two models must be fitted to calculate a likelihood-ratio statistic?
Explain why estimating nuisance parameters under the null usually makes a parametric bootstrap approximate rather than exactly valid in finite samples.
For a score test with nuisance parameters, why must the score and information be recomputed in each bootstrap replicate?
Give one reason the regular chi-square likelihood-ratio approximation may fail, and explain what feature the bootstrap procedure must reproduce in that case.