Perturbing Observed Outcomes Without Altering Missing Values

Why Perturb Only the Recorded Outcomes?

A selection model describes the joint behavior of outcomes and the observation process. For a subject with incomplete longitudinal data, the likelihood includes both the outcome density and the probability of the subject’s observed missingness pattern. A useful local diagnostic can ask: if one or more recorded outcomes were slightly different, how would the fitted parameters respond? The question is about sensitivity to the observed values, not about replacing missing values with guesses.

This distinction matters because the selection probability may depend on outcomes that were not observed. Those outcomes must remain integrated out in the observed-data likelihood. Substituting estimated missing values, or perturbing the full unobserved trajectory as if it had been measured, defines a different calculation.

Definition: For subject $i$, let $Y_i=(Y_{i,\mathrm{obs}},Y_{i,\mathrm{mis}})$ denote the complete longitudinal outcome vector, partitioned according to the observed pattern $R_i$. Let $f_i(y_i;\theta)$ be the outcome density and $g_i(r_i\mid y_i;\psi)$ the selection-model probability of the observed pattern. The observed-data likelihood contribution is $$ L_i(\eta;0)=\int f_i(y_{i,\mathrm{obs}},y_{i,\mathrm{mis}};\theta)\, g_i(r_i\mid y_{i,\mathrm{obs}},y_{i,\mathrm{mis}};\psi)\,dy_{i,\mathrm{mis}}, \qquad \eta=(\theta,\psi). $$ An outcome perturbation replaces only the observed coordinates by $y_{i,\mathrm{obs}}(\omega)=y_{i,\mathrm{obs}}+D_i\omega$. Its contribution is $$ L_i(\eta;\omega)=\int f_i(y_{i,\mathrm{obs}}+D_i\omega,y_{i,\mathrm{mis}};\theta)\, g_i(r_i\mid y_{i,\mathrm{obs}}+D_i\omega,y_{i,\mathrm{mis}};\psi)\,dy_{i,\mathrm{mis}}. $$ The missing coordinates remain variables of integration; they are not fixed, imputed, or shifted.

Here $\omega$ is a vector of small perturbation amounts and $D_i$ specifies which recorded coordinates are shifted and by how much per unit of $\omega$. The unperturbed data correspond to $\omega=0$. In a one-coordinate diagnostic, $D_i\omega$ might add $\omega$ to one observed measurement. The selection probability is evaluated at the perturbed observed values too, because those values enter the selection model. This is not a perturbation of the selection equation’s covariates or coefficients.

Assumptions and Local Parameter Change

Let $\ell(\eta;\omega)=\sum_i\log L_i(\eta;\omega)$ be the observed-data log likelihood. The diagnostic is defined around an unperturbed maximum-likelihood estimate $\hat{\eta}$, satisfying $\ell_\eta(\hat{\eta};0)=0$. For the local derivative to be meaningful, assume the likelihood is twice differentiable near $(\hat{\eta},0)$, its derivatives can be moved through the missing-outcome integrals, and the fitted solution is an interior point with nonsingular observed information. The outcome and selection models must also define valid densities and probabilities, and the observed-data likelihood must use the actual pattern of observed coordinates.

Write $J=-\ell_{\eta\eta}(\hat{\eta};0)$ for the observed information matrix, and $C=\ell_{\eta\omega}(\hat{\eta};0)$ for the matrix of cross-derivatives of the score with respect to the perturbation. Differentiating the score equations gives the local parameter response.

Result: If $J$ is nonsingular, the derivative of the fitted parameter vector with respect to the perturbation is $$ \left.\frac{d\hat{\eta}(\omega)}{d\omega^\mathsf{T}}\right|_{\omega=0} =J^{-1}C. $$ For a scalar target $q(\eta)$, with gradient $a=\nabla_\eta q(\hat{\eta})$, its local response is $a^\mathsf{T}J^{-1}C$. This is a first-order sensitivity, not the exact change for a large perturbation.

The sign follows directly from implicit differentiation. At a perturbed maximum, the score is zero. Differentiating that equation with respect to $\omega$ gives $\ell_{\eta\eta}\,d\hat{\eta}/d\omega+\ell_{\eta\omega}=0$. Since $\ell_{\eta\eta}=-J$, rearrangement yields the displayed expression.

0=\frac{d}{d\omega^\mathsf{T}}\ell_\eta(\hat{\eta}(\omega);\omega) =\ell_{\eta\eta}\frac{d\hat{\eta}}{d\omega^\mathsf{T}}+\ell_{\eta\omega}, \qquad \frac{d\hat{\eta}}{d\omega^\mathsf{T}}=J^{-1}C.

The derivatives in this result are derivatives of the integrated observed-data likelihood. In general, they are not obtained by treating the missing outcomes as known. In computation, one may evaluate the integral and differentiate it, or use a valid numerical procedure that gives the same observed-likelihood derivatives.

Worked Example: Integrating a Missing Outcome After Shifting an Observed One

Consider a hypothetical subject with one recorded outcome $x$ and a missing binary outcome $Z$. Suppose $P(Z=0)=0.25$, $P(Z=1)=0.75$, and the outcome density for $x$ is $\phi(x;0,1)$, independent of $Z$. The observed pattern is $R=1$, and the selection probability is $\operatorname{expit}(0.5x-1)$ when $Z=0$, and $\operatorname{expit}(0.5x)$ when $Z=1$. The notation $\operatorname{expit}(u)=1/(1+e^{-u})$ denotes the logistic function.

The observed-data contribution integrates over the missing $Z$: $$ L(x)=\phi(x;0,1)\{0.25\operatorname{expit}(0.5x-1)+0.75\operatorname{expit}(0.5x)\}. $$ At the observed value $x=0.4$, the two selection probabilities are $\operatorname{expit}(-0.8)\approx0.3100$ and $\operatorname{expit}(0.2)\approx0.5498$. Thus the integrated selection factor is approximately $0.25(0.3100)+0.75(0.5498)=0.4899$. Since $\phi(0.4;0,1)\approx0.3683$, the likelihood contribution is approximately $0.3683(0.4899)=0.1804$.

Now perturb the recorded outcome to $x^*=0.5$, leaving the missing $Z$ integrated out. The selection probabilities become $\operatorname{expit}(-0.75)\approx0.3208$ and $\operatorname{expit}(0.25)\approx0.5622$. Using the unrounded probabilities, the integrated factor is approximately $0.5018$. Also, $\phi(0.5;0,1)\approx0.3521$, so the perturbed contribution is approximately $0.3521(0.5018)=0.1767$. The decrease is a change in this subject’s observed-data likelihood contribution, not evidence by itself that the fitted parameters or scientific conclusion must change substantially.

The example assumes the stated density and selection probabilities are correctly specified and that the two probabilities for $Z$ sum to one. It also shows why the selection factor must be averaged over the missing outcome rather than evaluated at a chosen value of $Z$.

Reading Parameter Sensitivities

A derivative reports the direction and rate of local change. Its units are the units of the parameter per unit of outcome perturbation. A positive derivative means that a small positive shift in the specified observed measurement moves the estimate upward, to first order; a negative derivative means it moves downward. For a finite perturbation $\Delta\omega$, the approximation is $\hat{\eta}(\Delta\omega)-\hat{\eta}(0)\approx J^{-1}C\Delta\omega$. The approximation should not be treated as an exact refit.

Worked Example: Translating Curvature Into a Contrast Change

Suppose a hypothetical fitted model has parameters $\eta=(\beta,\gamma)^\mathsf{T}$. At the fit, numerical derivatives of the integrated likelihood give $$ J=\begin{pmatrix}40&8\\8&20\end{pmatrix}, \qquad C=\begin{pmatrix}2\\-1\end{pmatrix} $$ for a one-unit upward perturbation of one recorded outcome. These are hypothetical derivative values, not data from a real study. The information matrix is nonsingular: its determinant is $40(20)-8(8)=736>0$, and its leading diagonal entry is positive, so it is positive definite.

The inverse is $J^{-1}=736^{-1}\begin{pmatrix}20&-8\\-8&40\end{pmatrix}$. Therefore, $$ J^{-1}C=\frac{1}{736}\begin{pmatrix}20&-8\\-8&40\end{pmatrix} \begin{pmatrix}2\\-1\end{pmatrix} =\frac{1}{736}\begin{pmatrix}48\\-56\end{pmatrix} \approx\begin{pmatrix}0.06522\\-0.07609\end{pmatrix}. $$

Thus, per one unit of upward shift in that recorded outcome, $\hat{\beta}$ has local change $0.06522$ and $\hat{\gamma}$ has local change $-0.07609$. For a shift of $0.2$ outcome units, the first-order changes are approximately $0.01304$ and $-0.01522$, respectively.

For the target contrast $q(\eta)=\beta-0.5\gamma$, the gradient is $a=(1,-0.5)^\mathsf{T}$. Its sensitivity is $a^\mathsf{T}J^{-1}C=0.06522-0.5(-0.07609)\approx0.1033$ contrast units per outcome unit. A $0.2$-unit shift therefore gives an approximate contrast change of $0.02065$. The interpretation is a local change in the fitted contrast, not a causal effect of changing the recorded measurement.

Worked Example: Comparing Two Observed Coordinates

For the same hypothetical information matrix, suppose the cross-derivative vector for perturbing a second recorded measurement is $C_2=(-1,3)^\mathsf{T}$. This is another assumed derivative from the integrated likelihood. Then $$ J^{-1}C_2=\frac{1}{736}\begin{pmatrix}20&-8\\-8&40\end{pmatrix} \begin{pmatrix}-1\\3\end{pmatrix} =\frac{1}{736}\begin{pmatrix}-44\\128\end{pmatrix} \approx\begin{pmatrix}-0.05978\\0.1739\end{pmatrix}. $$

For $q=\beta-0.5\gamma$, the response is $-0.05978-0.5(0.1739)\approx-0.1467$ contrast units per outcome unit. A $0.1$-unit upward shift in the first recorded coordinate gives an approximate contrast change of $0.01033$, while the same-sized shift in the second gives approximately $-0.01467$. These opposite signed changes illustrate that sensitivity depends on which observed coordinate is perturbed and on how it enters both the outcome and selection components. The comparison is most directly interpretable when the outcomes share a common measurement scale.

Diagnostics, Alternatives, and Limits

The local calculation is useful when it identifies which recorded outcomes have the greatest effect on a parameter or contrast near the fitted model. It should be checked against numerical refits. For a small $\Delta\omega$, compare the exact refitted difference with $J^{-1}C\Delta\omega$; close agreement supports the local approximation in that direction. Large discrepancies can indicate that the perturbation is too large, the likelihood is strongly nonlinear, the optimizer has moved to a different solution, or numerical integration and differentiation are inaccurate.

Other warning signs include a nearly singular $J$, large changes under tiny perturbations, unstable derivatives across integration tolerances, and sensitivity that changes sharply when the perturbation size is halved. These can signal weak identification or numerical instability rather than an influential outcome alone. Verify that the perturbed observed values remain in the model’s support, that the same missingness pattern is being evaluated, and that the missing-outcome integral is computed accurately at each perturbation.

This diagnostic differs from several alternatives. Case deletion removes a subject’s entire likelihood contribution; an observed-outcome perturbation retains the subject and changes selected recorded values. Weighted likelihood perturbations change a subject’s contribution weight rather than its outcome values. A complete-data analysis that plugs in a missing outcome does not preserve the observed-data likelihood. Each answers a different sensitivity question, so none should be substituted for another without stating the change being investigated.

A perturbation does not establish that the data are erroneous, that a subject is atypical, or that a model is misspecified. Its size depends on the measurement units and on the chosen perturbation direction. Moreover, it is local: for a substantial change, refit the model at the altered values and report the actual difference. Conclusions also remain conditional on the selected outcome and selection models. In an MNAR analysis, the likelihood’s dependence on missing outcomes makes the selection specification important; an outcome-perturbation diagnostic does not resolve uncertainty about that specification.

Common pitfall: “I will shift the observed outcome and hold the missing outcome at its estimated value, because that keeps the calculation simple.” This replaces the required integral with a plug-in likelihood and can change both the outcome and selection contributions. Instead, shift only the observed coordinate inside the integrand and integrate over the missing coordinate under the specified model at each perturbation.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

When an observed outcome is shifted in this diagnostic, what happens to the missing outcomes in the subject’s likelihood contribution?
Suppose $J$ is nonsingular and $C$ is the score cross-derivative matrix. What is the local derivative of the fitted parameters with respect to the perturbation?
For $q(\eta)=\beta-0.5\gamma$, the local changes per outcome unit are $d\hat{\beta}/d\omega=0.04$ and $d\hat{\gamma}/d\omega=-0.02$. Find the local change in $q$ for a $0.5$-unit upward shift.
Name two diagnostics that can indicate the local approximation or its numerical calculation is unreliable.

Solutions. (1) Missing outcomes remain variables of integration; they are neither fixed nor perturbed. (2) The derivative is $J^{-1}C$, obtained by differentiating the zero-score equation. (3) The contrast derivative is $0.04-0.5(-0.02)=0.05$ contrast units per outcome unit. Multiplying by $0.5$ gives an approximate change of $0.025$ contrast units. (4) Examples include substantial disagreement between the derivative approximation and a small-perturbation refit; unstable derivatives when numerical integration tolerances change; a nearly singular information matrix; or sharply different sensitivities when the perturbation size is halved.

Perturbing Observed Outcomes Without Altering Missing Values

What You'll Learn