Why Interval Estimation Needs Special Care
A frailty variance describes unobserved heterogeneity: clusters with larger frailty tend to have higher event rates, conditional on measured covariates and the specified event model. An estimate alone does not show how precisely that heterogeneity has been measured. An interval is especially important when the estimated variance is small, because zero is a genuine boundary value, not an ordinary point on an unrestricted real line.
This tutorial develops profile-likelihood and bootstrap intervals for a frailty variance. The likelihood must be based on the assumed frailty model and observation process. The methods do not remove uncertainty caused by a poor model choice; they quantify uncertainty conditional on the model and assumptions used.
For example, a shared-frailty event model may specify the conditional event intensity for member \(j\) of cluster \(i\) as \(\lambda_{ij}(t\mid Z_i)=Z_i\lambda_0(t)\exp(x_{ij}^{T}\beta)\). The marginal likelihood integrates over the frailty distribution. Its exact form depends on the event process, censoring, covariates, baseline hazard, and chosen frailty family. Here the distributional family is treated as fixed; uncertainty intervals for \(\theta\) do not by themselves assess sensitivity to that choice.
The usual working assumptions include independent clusters, a correctly specified conditional event and censoring model, noninformative censoring given the modeled information, and a frailty distribution indexed by \(\theta\) with mean one. For a nonparametric cluster bootstrap, the sampled clusters must also reasonably represent the target population. If clusters are dependent, or censoring is informative in a way the model does not handle, the procedures below do not automatically give valid intervals.
Profile-Likelihood Confidence Sets
Write the full parameter as \((\theta,\eta)\), where \(\eta\) contains all nuisance parameters, such as regression coefficients and baseline-hazard parameters. Let \(\ell(\theta,\eta)\) be the marginal log likelihood and define the profile log likelihood by maximizing over the nuisance parameters separately for each fixed candidate variance.
The likelihood-ratio statistic comparing a candidate \(\theta\) with the best-fitting model is
For an interior true value, under standard likelihood regularity conditions, \(R(\theta)\) is approximately distributed as a chi-square random variable with one degree of freedom. A familiar approximate 95% profile-likelihood confidence set is therefore
The set is obtained by finding the candidate values whose profile likelihood is not too far below its maximum. If the profile is smooth and unimodal, the endpoints solve \(R(\theta)=3.841\), except that the lower endpoint is zero if the set includes zero. It is better to define the interval as a set than to assume it always has two finite endpoints: a very flat profile can produce a wide interval or an unbounded upper endpoint.
The justification for the likelihood-ratio cutoff is the local quadratic approximation to the log likelihood near an interior maximum. After nuisance parameters are accounted for, the loss in log likelihood for a one-dimensional parameter is approximately one half of a squared standard-normal displacement. Multiplying that loss by two gives an approximate \(\chi^2_1\) statistic. This argument relies on an interior parameter value and regularity conditions; it does not apply unchanged at \(\theta=0\).
Worked Example: Finding a Profile-Likelihood Interval
Suppose a hypothetical shared-frailty analysis gives a maximum profile log likelihood of \(-82.40\) at \(\hat{\theta}=0.12\). The fitted models at several fixed values of \(\theta\) have the following profile log likelihoods. The variance is dimensionless under the chosen frailty scale.
| Candidate \(\theta\) | \(\ell_p(\theta)\) | \(R(\theta)=2[-82.40-\ell_p(\theta)]\) |
|---|---|---|
| 0.00 | -83.55 | 2.30 |
| 0.05 | -82.62 | 0.44 |
| 0.12 | -82.40 | 0.00 |
| 0.25 | -82.86 | 0.92 |
| 0.40 | -83.50 | 2.20 |
| 0.50 | -84.15 | 3.50 |
| 0.55 | -84.55 | 4.30 |
For instance, at \(\theta=0.50\), \(R(0.50)=2\{ -82.40-(-84.15)\}=2(1.75)=3.50\). At \(\theta=0\), the statistic is \(2\{ -82.40-(-83.55)\}=2.30\), below 3.841, so the 95% set includes zero and its lower endpoint is zero.
The upper endpoint lies between 0.50, where \(R=3.50\), and 0.55, where \(R=4.30\). Linear interpolation of the likelihood-ratio values gives \(0.50+0.05(3.841-3.50)/(4.30-3.50)=0.521\), approximately. Thus the conventional profile interval is about \([0,0.521]\). The interpolation is a numerical approximation; a final analysis should evaluate the profile more finely and solve for the crossing directly.
Under this model, the data are compatible with no frailty variation and also with a range of positive variances. The interval is not a 95% probability statement about the fixed parameter; it is a procedure whose repeated-sampling coverage is approximately 95% under its calibration assumptions.
Boundary Estimates and Calibration
Because \(\theta\) cannot be negative, the usual regularity conditions for the chi-square likelihood-ratio result fail at zero. In the standard single-variance boundary setting, the likelihood-ratio statistic for testing zero variance has an asymptotic mixture distribution: probability one half at zero and probability one half from \(\chi^2_1\). This is the boundary issue encountered in testing zero frailty variance.
The mixture affects calibration at the boundary. Its 95th percentile is the 90th percentile of \(\chi^2_1\), about 2.706, because half the probability is already concentrated at zero. Consequently, using 3.841 for an interval endpoint when the fitted maximum is at zero is often conservative at the boundary in this idealized asymptotic setting: the mixture distribution is about 97.5% below 3.841. This does not mean 2.706 is a universally correct cutoff for every endpoint. For positive candidate values, the parameter being tested is generally interior, and model complexity or small-sample behavior can also affect calibration.
A boundary-adjusted or candidate-specific parametric bootstrap can be more informative. For each candidate value, fit the model with that value fixed, simulate data under that fitted model, refit both constrained and unconstrained models, and use the resulting likelihood-ratio distribution to decide whether the candidate is retained. This directly approximates the sampling distribution under the candidate model. It can account for a boundary and finite-sample features better than a single asymptotic cutoff, provided the model and fitting algorithm are adequate.
Worked Example: A Boundary-Aware Bootstrap Profile Set
In a hypothetical small study, the unconstrained fit has \(\hat{\theta}=0\). For each candidate variance, the analyst simulates 2,000 data sets from the fitted model with that candidate fixed, refits both models, and obtains the 95th percentile of the simulated likelihood-ratio statistic. The observed profile statistic and the bootstrap cutoff near the upper endpoint are:
| Candidate \(\theta\) | Observed \(R(\theta)\) | Bootstrap 95th-percentile cutoff | Decision |
|---|---|---|---|
| 0.00 | 0.00 | 2.70 | Retain |
| 0.20 | 2.20 | 3.10 | Retain |
| 0.24 | 2.82 | 3.05 | Retain |
| 0.25 | 3.00 | 2.96 | Reject |
The candidate is retained when \(R(\theta)\) does not exceed its bootstrap cutoff. Thus 0.24 is retained because \(2.82<3.05\), while 0.25 is rejected because \(3.00>2.96\). A finer search would locate the crossing between them; based on this grid, the upper endpoint is approximately 0.25 and the confidence set begins at zero. The specific cutoffs are hypothetical bootstrap results, not universal constants. In particular, the cutoff at zero need not equal the asymptotic mixture value exactly.
The conclusion is that this sample does not resolve a positive variance from zero, and the upper endpoint indicates the approximate range of larger frailty variances not excluded by this model and calibration. It does not show that the true variance is zero.
Bootstrap Intervals and Practical Choices
A simpler bootstrap interval applies a sampling procedure to the variance estimator itself. For independent clusters, a nonparametric bootstrap samples whole clusters with replacement, preserving each selected cluster’s event histories and covariates. It then refits the complete model in each bootstrap sample. Resampling individual event records separately would break the within-cluster dependence that the shared frailty represents and is generally inappropriate.
A percentile interval uses the empirical quantiles of the bootstrap estimates. A parametric bootstrap instead simulates complete data sets from the fitted model, including the cluster structure and censoring mechanism, then refits the model. The parametric version depends more directly on correct specification; the cluster bootstrap depends on having enough independent clusters to represent the population. Either procedure can be unstable when many refits fail or estimates pile up at zero.
Worked Example: Reading a Percentile Bootstrap Interval
Suppose a hypothetical analysis has 80 independent clusters and gives \(\hat{\theta}=0.08\). The analyst resamples the 80 clusters with replacement 2,000 times and successfully refits the model each time. In the sorted bootstrap estimates, the empirical 2.5th percentile is 0 and the 97.5th percentile is 0.34. Using the common order-statistic convention \(k\approx p(B+1)\), where \(B=2{,}000\), the percentile ranks are \(0.025(2001)=50.0\) and \(0.975(2001)=1951.0\). The corresponding estimates give the interval \([0,0.34]\).
This interval reflects the observed skewness and boundary pile-up more naturally than a symmetric interval centered at 0.08. It remains an approximate interval: quantile estimates have Monte Carlo error, and with only 80 independent clusters the empirical cluster distribution may poorly represent rare cluster types. The bootstrap interval does not establish that zero is the true variance; it indicates that the percentile procedure includes zero for these data.
A Wald interval, \(\hat{\theta}\pm 1.96\,SE\), is usually a poor default here. It is symmetric, can extend below the allowable value zero, and relies on an approximately normal estimator. Truncating a negative lower endpoint to zero does not repair a strongly skewed or boundary-concentrated sampling distribution. Profile likelihood uses the shape of the likelihood, while bootstrap intervals approximate a sampling distribution; each can be useful, but neither compensates for too few independent clusters or a misspecified frailty model.
Diagnostics, Limitations, and Common Mistakes
A profile interval is only as reliable as the fitted likelihood at each candidate value. Check convergence across the entire profile, not only at the unrestricted maximum. Repeating fits from different starting values, evaluating a sufficiently fine grid, and inspecting the profile curve can reveal numerical problems. For integrated frailty likelihoods, verify that quadrature accuracy is adequate; otherwise small numerical errors can shift the likelihood-ratio crossing.
For bootstrap methods, report the number of attempted and successful refits, the fraction of estimates equal to zero, and the interval method used. A large failure fraction or a very discrete bootstrap distribution is a warning, not a reason to silently discard inconvenient fits. More clusters, improved optimization, or a carefully specified parametric bootstrap may be needed. The effective information about a variance can be low even when the total number of event records is large, because the independent sampling units for shared frailty are clusters.
- Common mistake: treating the variance as an unrestricted parameter. A normal Wald interval may include negative values. The correction is to respect \(\theta\ge 0\) and use a method suited to the boundary and likelihood shape.
- Common mistake: fixing nuisance estimates during profiling. This generally lowers the candidate log likelihood relative to re-optimizing nuisance parameters and inflates the apparent likelihood-ratio cost. Refit nuisance parameters at every candidate value.
- Common mistake: interpreting inclusion of zero as proof of no heterogeneity. A wide interval including zero can reflect weak information. Report the interval and its uncertainty rather than declaring the variance absent.
- Common mistake: resampling individual records for shared frailty. That destroys the cluster dependence structure. Resample independent clusters, or simulate complete clustered data under a parametric model.
The main alternatives answer related but distinct questions. Profile likelihood summarizes the relative support in the fitted likelihood; a cluster bootstrap approximates repeated-sampling behavior by resampling observed clusters; a parametric bootstrap approximates it under a fully specified model. Boundary-aware bootstrap calibration can be helpful when the estimate is zero or nearly zero, while conventional profile likelihood is often a useful compact summary when the profile is well behaved. Sensitivity to the assumed frailty distribution and uncertainty in model structure remain separate concerns.
Check Your Understanding
Answer each question, then compare your reasoning with the brief solutions.
- Why must nuisance parameters be re-optimized at every candidate frailty variance when constructing a profile-likelihood interval?
- A profile-likelihood statistic is 3.2 at a candidate value, using the conventional 95% cutoff 3.841. Is that value included in the approximate confidence set?
- Why can a symmetric Wald interval be misleading when the estimated frailty variance is zero?
- For shared frailty across clusters, what should be resampled in a nonparametric bootstrap, and why?
Solutions. (1) The profile likelihood gives the best fit among nuisance values at each candidate variance; fixing them generally makes that candidate appear more costly than after re-optimization. (2) Yes, because \(3.2<3.841\), subject to the approximation and assumptions of that cutoff. (3) The variance is restricted to be nonnegative and its estimator may be skewed or concentrated at the boundary, so a symmetric normal approximation can extend below zero and misrepresent uncertainty. (4) Resample whole independent clusters so that within-cluster event dependence and histories remain intact.