Monte Carlo Integration for Non-Gaussian Outcomes

Why integrate over an unobserved outcome?

In a current-outcome selection model, the chance of dropout at a visit may depend on the outcome value at that visit. If dropout occurs before that outcome is recorded, the outcome is unobserved, but it still enters the model for dropout. The observed-data contribution must therefore average the dropout probability over the possible outcome values, weighted by their probabilities under the outcome model. When that outcome model is non-Gaussian, this integral may not have a convenient closed form.

Monte Carlo integration replaces the integral with an average over simulated outcomes. The simulation approximates a model-based expectation; it does not fill in the missing outcome with one guessed value. The introductory material on quadrature for logistic-normal selection models covers deterministic numerical integration in a particular setting. Here the emphasis is simulation-based integration for non-Gaussian outcome distributions and the practical checks needed to use it responsibly.

The observed-data contribution

Let $Y_t$ denote the outcome at visit $t$, and let $H_t$ denote the recorded history available before the current visit. Write $D_t=1$ if dropout occurs before $Y_t$ is observed, and $D_t=0$ if the subject remains observed at that visit. Suppose the outcome model has conditional density or probability mass function $f(y\mid H_t;\theta)$, and the dropout model has conditional probability $\pi_D(y,H_t;\gamma)=P(D_t=1\mid Y_t=y,H_t)$. The parameters $\theta$ and $\gamma$ index the outcome and dropout models.

Definition: When $D_t=1$ and $Y_t$ is unobserved, the current-visit observed-data contribution is $$ L_t(\theta,\gamma)=\int f(y\mid H_t;\theta)\,\pi_D(y,H_t;\gamma)\,dy =E_{\theta}\!\left[\pi_D(Y_t,H_t;\gamma)\mid H_t\right]. $$ For a discrete outcome, the integral is interpreted as a sum over its possible values.

This expression assumes that the specified outcome distribution describes $Y_t$ conditional on the recorded history, and that the dropout probability is correctly specified conditional on both $Y_t$ and that history. It also assumes the timing represented by the model: the current outcome can affect the dropout decision even though it is not observed after dropout. The integration is over the outcome distribution before learning whether dropout occurred, not over a distribution already conditioned on dropout.

For a subject who remains observed, the corresponding current-visit factor, when $Y_t=y_t$ is recorded, is $f(y_t\mid H_t;\theta)\{1-\pi_D(y_t,H_t;\gamma)\}$. Thus the observed outcome density or mass and the response probability both contribute. If dropout occurs and the outcome is missing, use the integrated factor instead; multiplying by a dropout probability evaluated at an invented value of $Y_t$ does not give this contribution.

Monte Carlo approximation and its error

If independent draws $Y_t^{(1)},\ldots,Y_t^{(M)}$ are generated from $f(\,\cdot\mid H_t;\theta)$, then the integral is an expectation under that distribution. Its direct Monte Carlo estimate is the sample average of the dropout probabilities:

\widehat L_t=\frac{1}{M}\sum_{m=1}^{M}\pi_D(Y_t^{(m)},H_t;\gamma).

Each term lies between zero and one, so this estimate is also between zero and one. Under independent draws from the stated outcome model, the law of large numbers implies $\widehat L_t$ converges to $L_t$ as $M$ increases. The estimated Monte Carlo standard error is $s_g/\sqrt{M}$, where $s_g$ is the sample standard deviation of the $M$ simulated dropout probabilities. This standard error describes simulation variability for the specified model and parameters; it does not measure uncertainty in those parameters or error from a misspecified model.

Justification: Set $G_m=\pi_D(Y_t^{(m)},H_t;\gamma)$. The draws are independent and identically distributed with $E(G_m)=L_t$. Therefore $E(\widehat L_t)=L_t$, and $$ \operatorname{Var}(\widehat L_t)=\operatorname{Var}(G_m)/M. $$ Replacing the unknown standard deviation of $G_m$ with its sample standard deviation $s_g$ gives the estimated standard error $s_g/\sqrt{M}$.

The justification depends on drawing from the conditional outcome distribution used in the integral. If the simulation distribution is different, a correction is needed. The normal approximation to Monte Carlo error can also be unreliable with very few draws, especially when nearly all simulated dropout probabilities are close to zero or one. In practice, increasing $M$ and checking whether estimates are stable across independent runs are useful complements to a single reported standard error.

Worked Example: Direct Simulation With a Gamma Outcome

Suppose that, conditional on the recorded history, the unobserved current outcome has a Gamma distribution with shape 2 and scale 2. The dropout model is logistic: $ \pi_D(y)=\operatorname{logit}^{-1}(-2+0.5y) $, where $\operatorname{logit}^{-1}(u)=1/(1+e^{-u})$. The outcome distribution and logistic dropout model are assumed to be the specified conditional models. Simulate eight outcomes from the Gamma distribution and suppose the resulting dropout probabilities, calculated from the displayed outcomes before rounding, are 0.2000, 0.3000, 0.3500, 0.3114, 0.8000, 0.7000, 0.6500, and 0.6886.

The observed-data contribution is the expectation of the dropout probability, so its Monte Carlo estimate is

\widehat L_t=\frac{0.2000+0.3000+0.3500+0.3114+0.8000+0.7000+0.6500+0.6886}{8} =\frac{4.0000}{8}=0.5000.

The eight probabilities have sample standard deviation approximately 0.2318. The estimated Monte Carlo standard error is therefore $0.2318/\sqrt{8}=0.0820$. The estimate says that the integrated dropout contribution is 0.5000 for these parameters and history, with substantial simulation variability from using only eight draws. It is not evidence that the true dropout probability has a 95% chance of being within a particular range, and the small simulation is not presented as a reliable final computation.

The calculation uses the Gamma distribution to generate outcomes and then evaluates the logistic hazard at each draw. It does not assume that the outcome itself is normally distributed. With more draws, the Monte Carlo error would generally decrease at a rate proportional to $1/\sqrt{M}$, though the model-based estimate itself could still be inaccurate if either conditional model is unsuitable.

When direct simulation is inefficient

Direct simulation is convenient when the outcome model is easy to sample from and the simulated dropout probabilities adequately represent the integral. It can be inefficient if the integral is driven by outcome values that are rare under the outcome distribution. Importance sampling addresses this by drawing from a proposal density $q(y)$ that places more draws in relevant regions. Provided $q(y)>0$ wherever $f(y\mid H_t;\theta)\pi_D(y,H_t;\gamma)$ is nonzero, the same integral can be written as

L_t=\int q(y)\frac{f(y\mid H_t;\theta)}{q(y)}\pi_D(y,H_t;\gamma)\,dy, \qquad \widehat L_t^{\,IS}=\frac{1}{M}\sum_{m=1}^{M} \frac{f(Y^{(m)}\mid H_t;\theta)}{q(Y^{(m)})}\pi_D(Y^{(m)},H_t;\gamma), \quad Y^{(m)}\sim q.

The weights correct for the fact that the proposal does not draw from the target outcome distribution. A proposal that puts too little probability where the integrand is large can produce a highly variable estimate, even if the estimator is mathematically valid. The sample average above is the ordinary importance-sampling estimator; dividing by the sum of weights instead produces a self-normalized estimator, which is a different finite-sample procedure.

Worked Example: Importance Sampling for a Difficult Integral

For one history, suppose the target outcome density is $f$, but direct draws are inconvenient. Draw four independent values from a proposal $q$. At those values, suppose the density ratios $f(Y^{(m)}\mid H_t)/q(Y^{(m)})$ are 0.5, 1.0, 1.5, and 2.0, and the corresponding dropout probabilities are 0.10, 0.30, 0.60, and 0.80. Assume the proposal covers the support of the target integrand, as required.

The weighted terms are $0.5(0.10)=0.05$, $1.0(0.30)=0.30$, $1.5(0.60)=0.90$, and $2.0(0.80)=1.60$. Thus

\widehat L_t^{\,IS}=\frac{0.05+0.30+0.90+1.60}{4} =\frac{2.85}{4}=0.7125.

This estimate is mathematically formed with the required density-ratio correction. But the four weighted terms vary widely: one term, 1.60, accounts for more than half their sum. That is a warning that four draws provide weak evidence of numerical stability. A practical analysis would use many more draws, inspect the weights and weighted terms, and assess sensitivity to the proposal distribution. Importance sampling can reduce error when the proposal is well chosen; it does not guarantee a precise answer merely because the weights have been included.

Putting an integrated factor into a longitudinal contribution

For monotone dropout, a subject’s observed-data contribution combines factors for visits where outcomes are recorded with an integrated factor at the visit of dropout. The outcome distributions and response probabilities at earlier visits are conditional on the history appropriate to those visits. The following example isolates this bookkeeping; it does not prescribe a complete trial model.

Worked Example: A Dropout Factor After an Observed Visit

Suppose a hypothetical subject has an observed first-visit outcome with conditional probability mass 0.25. Given that outcome and the recorded history, the probability of remaining observed at that visit is 0.90. At the next visit, dropout occurs before the current outcome is recorded. A Monte Carlo calculation using 10,000 draws from the specified non-Gaussian outcome model gives an estimated integrated dropout factor of 0.4000. The outcome model and dropout probabilities are assumed correctly specified for this illustration.

If the simulated dropout indicators or probabilities have sample standard deviation 0.4900, the estimated Monte Carlo standard error of the integrated factor is $0.4900/\sqrt{10{,}000}=0.0049$. The subject’s observed-data contribution is the product of the first-visit observed outcome mass, the first-visit continuation probability, and the next-visit integrated dropout factor:

$$ 0.25(0.90)(0.4000)=0.0900. $$

Thus the estimated contribution is 0.0900 under the stated model and history. Replacing the next-visit factor by a dropout probability at one fabricated current outcome would change the likelihood contribution and would not implement the required integration. The Monte Carlo standard error applies to the estimated 0.4000 factor; it is not a standard error for the whole fitted model or for a treatment effect.

Diagnostics, alternatives, and limitations

Before interpreting a numerical estimate, check that the simulation matches the mathematical integral. Draws for direct Monte Carlo must come from the outcome distribution conditional on the observed history, and each draw must be evaluated using the dropout equation for that same history. For importance sampling, verify support coverage and examine both density ratios and weighted terms. A few very large weights, estimates that change sharply with the random seed, or materially different estimates across batches signal that the computation may be unstable.

The main alternatives have different strengths:

Approach	Useful when	Important limitation
Direct Monte Carlo	The outcome model is easy to simulate from and the integral is an expectation under that model.	Rare but influential outcome regions may be missed; results vary with the simulated draws.
Importance sampling	A proposal can place more draws where the weighted integrand matters.	Unstable or extreme weights can make the estimate noisy.
Quadrature	The integral is low-dimensional and a suitable deterministic rule is available.	Performance depends on the integration rule and dimension; non-Gaussian shapes may require careful adaptation.
Exact summation	The outcome has a manageable finite or countable support that can be summed accurately.	The number of relevant values can be too large for practical summation.

These methods approximate the same model-based contribution when correctly implemented, but they do not repair a poor outcome or dropout model. Monte Carlo error is numerical error conditional on the specified distributions and parameters; it is distinct from uncertainty in estimated parameters and from uncertainty about assumptions that cannot be learned from the observed data alone. A finite simulation can also make the logarithm of an estimated contribution numerically troublesome if an estimate is zero. More draws or a better integration strategy may help, but any numerical remedy should preserve the intended likelihood contribution.

A common mistake is to simulate one plausible missing outcome, evaluate its dropout probability, and treat that single value as the observed-data contribution. That replaces an expectation with one random draw. The correction is to average the dropout probabilities across draws from the outcome model, or to use another valid integration method. Another mistake is to report a tiny Monte Carlo standard error as proof that the model is appropriate: it only indicates precision of the numerical approximation under the stated model.

Key takeaway: When an unobserved current outcome affects dropout, its observed-data contribution averages the dropout probability over the conditional outcome distribution. Monte Carlo integration estimates that average; adequate sampling and diagnostics are essential, and the result remains conditional on the specified models.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

Why is the dropout probability evaluated at one imputed current outcome not generally the observed-data contribution?
If $M$ independent direct draws give a sample standard deviation of dropout probabilities equal to 0.20, what is the estimated Monte Carlo standard error when $M=400$?
For importance sampling, what support condition must the proposal distribution satisfy?
A direct Monte Carlo estimate changes markedly across independent runs. Name one diagnostic conclusion and one reasonable response.

Solutions. (1) The required contribution is the expectation of the dropout probability over the conditional outcome distribution, whereas one imputed value gives only one random evaluation. (2) $0.20/\sqrt{400}=0.20/20=0.0100$. (3) The proposal must have positive density wherever the target integrand $f(y\mid H_t)\pi_D(y,H_t)$ is nonzero. (4) The variation suggests appreciable Monte Carlo error or an unstable integrand; increase the number of draws and compare batch estimates, or consider a suitable proposal or quadrature method.