The GLIMMIX Procedure

Confidence Bounds Based on Likelihoods

Families of statistical tests can be inverted to produce confidence limits for parameters. The confidence region for parameter theta is the set of values for which the corresponding test fails to reject upper H colon theta equals theta 0. When parameters are estimated by maximum likelihood or a likelihood-based technique, it is natural to consider the likelihood ratio test statistic for H in the test inversion. When there are multiple parameters in the model, however, you need to supply values for these nuisance parameters during the test inversion as well.

In the following, suppose that bold-italic theta is the covariance parameter vector and that one of its elements, theta, is the parameter of interest for which you want to construct a confidence interval. The other elements of bold-italic theta are collected in the nuisance parameter vector bold-italic theta 2. Suppose that ModifyingAbove bold-italic theta With caret is the estimate of bold-italic theta from the overall optimization and that upper L left-parenthesis ModifyingAbove bold-italic theta With caret right-parenthesis is the likelihood evaluated at that estimate. If estimation is based on pseudo-data, then upper L left-parenthesis ModifyingAbove bold-italic theta With caret right-parenthesis is the pseudo-likelihood based on the final pseudo-data. If estimation uses a residual (restricted) likelihood, then L denotes the restricted maximum likelihood and ModifyingAbove bold-italic theta With caret is the REML estimate.

Profile Likelihood Bounds

The likelihood ratio test statistic for testing upper H colon theta equals theta 0 is

2 StartSet log left-brace upper L left-parenthesis ModifyingAbove bold-italic theta With caret right-parenthesis right-brace minus log left-brace upper L left-parenthesis theta 0 comma ModifyingAbove bold-italic theta With caret Subscript 2 Baseline right-parenthesis right-brace EndSet

where ModifyingAbove bold-italic theta With caret Subscript 2 is the likelihood estimate of bold-italic theta 2 under the restriction that theta equals theta 0. To invert this test, a function is defined that returns the maximum likelihood for a fixed value of theta by seeking the maximum over the remaining parameters. This function is termed the profile likelihood (Pawitan 2001, Ch. 3.4),

lamda Subscript p Baseline equals upper L left-parenthesis bold-italic theta 2 vertical-bar theta overTilde right-parenthesis equals sup Underscript bold-italic theta 2 Endscripts upper L left-parenthesis theta overTilde comma bold-italic theta 2 right-parenthesis

In computing lamda Subscript p, theta is fixed at theta overTilde and bold-italic theta 2 is estimated. In mixed models, this step typically requires a separate, iterative optimization to find the estimate of bold-italic theta 2 while theta is held fixed. The left-parenthesis 1 minus alpha right-parenthesis times 100 percent-sign profile likelihood confidence interval for theta is then defined as the set of values for theta overTilde that satisfy

2 StartSet log left-brace upper L left-parenthesis ModifyingAbove bold-italic theta With caret right-parenthesis right-brace minus log left-brace upper L left-parenthesis bold-italic theta 2 vertical-bar theta overTilde right-parenthesis right-brace EndSet less-than-or-equal-to chi Subscript 1 comma left-parenthesis 1 minus alpha right-parenthesis Superscript 2

The GLIMMIX procedure seeks the values theta overTilde Subscript l and theta overTilde Subscript u that mark the endpoints of the set around ModifyingAbove theta With caret that satisfy the inequality. The values left-parenthesis theta overTilde Subscript l Baseline and theta overTilde Subscript u Baseline right-parenthesis are then called the left-parenthesis 1 minus alpha right-parenthesis times 100 percent-sign confidence bounds for theta. Note that the GLIMMIX procedure assumes that the confidence region is not disjoint and relies on the convexity of upper L left-parenthesis ModifyingAbove bold-italic theta With caret right-parenthesis.

It is not always possible to find values theta overTilde Subscript l and theta overTilde Subscript u that satisfy the inequalities. For example, when the parameter space is (0 comma normal infinity right-parenthesis and

2 StartSet log left-brace upper L left-parenthesis ModifyingAbove bold-italic theta With caret right-parenthesis right-brace minus log left-brace upper L left-parenthesis bold-italic theta 2 vertical-bar 0 right-parenthesis right-brace EndSet greater-than chi Subscript 1 comma left-parenthesis 1 minus alpha right-parenthesis Superscript 2

a lower bound cannot be found at the desired confidence level. The GLIMMIX procedure reports the right-tail probabilities that are achieved by the underlying likelihood ratio statistic separately for lower and upper bounds.

Effect of Scale Parameter

When a scale parameter phi is eliminated from the optimization by profiling from the likelihood, some parameters might be expressed as ratios with phi in the optimization. This is the case, for example, in variance component models. The profile likelihood confidence bounds are reported on the scale of the parameter in the overall optimization. In case parameters are expressed as ratios with phi or functions of phi, the column RatioEstimate is added to the "Covariance Parameter Estimates" table. If parameters are expressed as ratios with phi and you want confidence bounds for the unscaled parameter, you can prevent profiling of phi from the optimization with the NOPROFILE option in the PROC GLIMMIX statement, or choose estimated likelihood confidence bounds with the TYPE=ELR suboption of the CL option in the COVTEST statement. Note that the NOPROFILE option is automatically in effect with METHOD=LAPLACE and METHOD=QUAD.

Estimated Likelihood Bounds

Computing profile likelihood ratio confidence bounds can be computationally expensive, because of the need to repeatedly estimate bold-italic theta 2 in a constrained optimization. A computationally simpler method to construct confidence bounds from likelihood-based quantities is to use the estimated likelihood (Pawitan 2001, Ch. 10.7) instead of the profile likelihood. An estimated likelihood technique replaces the nuisance parameters in the test inversion with some other estimate. If you choose the TYPE=ELR suboption of the CL option in the COVTEST statement, the GLIMMIX procedure holds the nuisance parameters fixed at the likelihood estimates. The estimated likelihood statistic for inversion is then

lamda Subscript e Baseline equals upper L left-parenthesis theta overTilde comma ModifyingAbove bold-italic theta With caret Subscript 2 Baseline right-parenthesis

where ModifyingAbove bold-italic theta With caret Subscript 2 are the elements of ModifyingAbove bold-italic theta With caret that correspond to the nuisance parameters. As the values of theta overTilde are varied, no reestimation of bold-italic theta 2 takes place. Although computationally more economical, estimated likelihood intervals do not take into account the variability associated with the nuisance parameters. Their coverage can be satisfactory if the parameter of interest is not (or only weakly) correlated with the nuisance parameters. Estimated likelihood ratio intervals can fall short of the nominal coverage otherwise.

Figure 11 depicts profile and estimated likelihood ratio intervals for the parameter sigma in a two-parameter compound-symmetric model, bold-italic theta equals left-bracket sigma comma phi right-bracket prime, in which the correlation between the covariance parameters is small. The elliptical shape traces the set of values for which the likelihood ratio test rejects the hypothesis of equality with the solution. The interior of the ellipse is the "acceptance" region of the test. The solid and dashed lines depict the PLR and ELR confidence limits for sigma, respectively. Note that both confidence limits intersect the ellipse and that the ELR interval passes through the REML estimate of phi. The PLR bounds are found as those points intersecting the ellipse, where phi equals the constrained REML estimate.

Figure 11: PLR and ELR Intervals, Small Correlation between Parameters

PLR and ELR Intervals, Small Correlation between Parameters


The major axes of the ellipse in Figure 11 are nearly aligned with the major axes of the coordinate system. As a consequence, the line connecting the PLR bounds passes close to the REML estimate in the full model. As a result, ELR bounds will be similar to PLR bounds. Figure 12 displays a different scenario, a two-parameter AR(1) covariance structure with a more substantial correlation between the AR(1) parameter (rho) and the residual variance (phi).

Figure 12: PLR and ELR Intervals, Large Correlation between Parameters

PLR and ELR Intervals, Large Correlation between Parameters


The correlation between the parameters yields an acceptance region whose major axes are not aligned with the axes of the coordinate system. The ELR bound for rho passes through the REML estimate of phi from the full model and is much shorter than the PLR interval. The PLR interval aligns with the major axis of the acceptance region; it is the preferred confidence interval.

Last updated: March 08, 2022