Delta-Method Uncertainty for Products of Estimated Probabilities

Why uncertainty in a product needs care

A subject’s inverse probability of treatment weight can involve a product of fitted probabilities from several decision times. Even when each probability is estimated precisely, uncertainty in their product depends on how all the estimates vary together. In particular, probabilities fitted from a shared treatment model, or from parameter estimates that are correlated, are not independent just because they correspond to different times.

The delta method gives a first-order approximation to the uncertainty of a smooth function of an asymptotically normal estimator. Here the function is a product of sequential probabilities, or its inverse. The approximation is useful for understanding how model-parameter uncertainty propagates into a weight; it is not a guarantee that the weight’s sampling distribution is close to normal, especially when the product is very small.

Setup and assumptions

For subject \(i\), let \(u_{it}(\theta)\) be the treatment-model probability assigned to the treatment actually observed at decision time \(t\), conditional on the history used in that model. For \(K\) decision times, define the sequential probability product and its inverse by

q_i(\theta)=\prod_{t=1}^{K}u_{it}(\theta), \qquad w_i(\theta)=\frac{1}{q_i(\theta)}.

The covariate and treatment history for the subject whose weight is being evaluated is held fixed when differentiating with respect to \(\theta\). The parameter \(\theta\) includes all parameters that determine the fitted probabilities. Thus, if a model shares parameters across times, the same parameter appears in multiple factors.

Assumptions: The number of decision times \(K\) is fixed. The true parameter \(\theta_0\) is an interior point of a finite-dimensional parameter space, and \(\hat{\theta}\xrightarrow{p}\theta_0\). For the delta-method limit, assume \(\sqrt{n}(\hat{\theta}-\theta_0)\xrightarrow{d}N(0,V)\), where \(V\) is finite, and that the function being transformed is differentiable at \(\theta_0\). For the product and inverse formulas below, assume each \(u_{it}(\theta_0)>0\) and that the factors are differentiable in a neighborhood of \(\theta_0\). When estimating the variance by plugging in a gradient, also assume that gradient is continuous at \(\theta_0\), and that \(\hat V\xrightarrow{p}V\).

The asymptotic normality assumption is the input to the calculation; its justification depends on the treatment-model estimator and sampling design. For example, an appropriate covariance estimate may come from the model’s information matrix or a sandwich covariance estimator. Whichever is used, \(V\) is the covariance in the limiting distribution of \(\sqrt{n}(\hat{\theta}-\theta_0)\), so the approximate covariance of \(\hat{\theta}\) itself is \(V/n\).

The delta method for the product and its inverse

Theorem: Under the stated asymptotic normality assumption, if \(g\) is differentiable at \(\theta_0\), then \[ \sqrt{n}\{g(\hat{\theta})-g(\theta_0)\} \xrightarrow{d} N\!\left(0,\nabla g(\theta_0)^T V\nabla g(\theta_0)\right). \] If, in addition, \(\nabla g\) is continuous at \(\theta_0\) and \(\hat V\xrightarrow{p}V\), the plug-in variance \[ \widehat{\operatorname{Var}}\{g(\hat{\theta})\} =\frac{1}{n}\nabla g(\hat{\theta})^T\hat V\nabla g(\hat{\theta}) \] is consistent for the first-order asymptotic variance.

For the first statement, differentiability at \(\theta_0\) means that the remainder after the linear approximation is negligible relative to \(\|\hat{\theta}-\theta_0\|\). Asymptotic normality implies \(\sqrt{n}(\hat{\theta}-\theta_0)=O_p(1)\), and therefore \(\hat{\theta}-\theta_0=O_p(n^{-1/2})\). The remainder multiplied by \(\sqrt n\) converges in probability to zero. The linear term then has the stated limiting distribution. The extra continuity assumption in the variance statement ensures \(\nabla g(\hat{\theta})\xrightarrow{p}\nabla g(\theta_0)\); consistency of \(\hat V\) and the continuous mapping theorem give the plug-in result.

For a positive differentiable product, differentiating its logarithm gives a particularly convenient expression:

\nabla q_i(\theta) =q_i(\theta)\sum_{t=1}^{K}\nabla\log u_{it}(\theta), \qquad \nabla w_i(\theta) =-w_i(\theta)\sum_{t=1}^{K}\nabla\log u_{it}(\theta).

The second formula follows by differentiating \(w_i=q_i^{-1}\). Define \(s_i(\theta)=\sum_t\nabla\log u_{it}(\theta)\). The first-order variance approximations are consequently

\widehat{\operatorname{Var}}\{q_i(\hat{\theta})\} =\frac{q_i(\hat{\theta})^2}{n}\, s_i(\hat{\theta})^T\hat V s_i(\hat{\theta}), \qquad \widehat{\operatorname{Var}}\{w_i(\hat{\theta})\} =\frac{w_i(\hat{\theta})^2}{n}\, s_i(\hat{\theta})^T\hat V s_i(\hat{\theta}).

These are first-order approximations, not exact finite-sample variances. The same quadratic form appears for both functions, but the scale differs: \(w_i=1/q_i\). The full matrix \(\hat V\), including off-diagonal entries, is essential.

Worked calculations

Worked Example: Two Fitted Probabilities With Correlated Parameters

Suppose two observed treatments have fitted probabilities \(u_1(\theta)=\operatorname{expit}(\theta_1)\) and \(u_2(\theta)=\operatorname{expit}(\theta_2)\), with fitted values \(0.8\) and \(0.5\). At the parameter value of interest, let the asymptotic covariance matrix be \(V=\begin{pmatrix}0.04&0.01\\0.01&0.09\end{pmatrix}\). Find the first-order variances of the product and its inverse.

The product is \(q=0.8(0.5)=0.4\), and \(w=1/0.4=2.5\). For an observed treatment coded \(1\), the derivative of the log fitted logistic probability with respect to its linear predictor is \(1-p\). Thus \(s=(0.2,0.5)^T\), and \(\nabla q=q s=(0.08,0.20)^T\). The product variance is

\frac{1}{n}\nabla q^T V\nabla q =\frac{1}{n}\left\{ 0.08^2(0.04)+2(0.08)(0.20)(0.01)+0.20^2(0.09) \right\} =\frac{0.004176}{n}.

For the inverse, \(\nabla w=-w s=(-0.5,-1.25)^T\), so

\frac{1}{n}\nabla w^T V\nabla w =\frac{1}{n}\left\{ 0.5^2(0.04)+2(0.5)(1.25)(0.01)+1.25^2(0.09) \right\} =\frac{0.163125}{n}.

The middle term in each calculation is the contribution of parameter covariance. It is positive here because both the off-diagonal covariance and the corresponding gradient components have the same sign.

Worked Example: A Shared Parameter Across Decision Times

Suppose the observed treatment is \(1\) at both times, and the two model probabilities are \(u_1(\theta)=\operatorname{expit}(\theta)\) and \(u_2(\theta)=\operatorname{expit}(2\theta)\). At \(\theta_0=0\), both probabilities are \(0.5\). Suppose \(\sqrt n(\hat\theta-\theta_0)\xrightarrow{d}N(0,0.16)\). Calculate the approximate variances of the fitted product and inverse.

Here \(q=0.5(0.5)=0.25\) and \(w=4\). The log-product derivative is the sum of the two log-probability derivatives. At \(\theta_0=0\), these are \((1-0.5)(1)=0.5\) and \((1-0.5)(2)=1\), respectively. Therefore \(s=1.5\), \(\nabla q=0.25(1.5)=0.375\), and \(\nabla w=-4(1.5)=-6\). The variances are

\widehat{\operatorname{Var}}(q(\hat\theta)) \approx \frac{0.375^2(0.16)}{n} =\frac{0.0225}{n}, \qquad \widehat{\operatorname{Var}}(w(\hat\theta)) \approx \frac{6^2(0.16)}{n} =\frac{5.76}{n}.

There is only one estimated parameter, so both probability factors necessarily move together as \(\theta\) changes. Treating their uncertainty as independent would discard this shared-parameter dependence and would not describe the estimator used here.

Worked Example: Direct Probability Parameters

For a two-time product, suppose the parameters are the probabilities themselves, \(\theta=(p_1,p_2)^T\). At the fitted value \((0.7,0.6)\), let the asymptotic covariance matrix be \(V=\begin{pmatrix}0.0004&0.0001\\0.0001&0.0009\end{pmatrix}\). Find the product variance and the inverse-weight variance.

The product is \(q=p_1p_2=0.42\), and \(\nabla q=(p_2,p_1)^T=(0.6,0.7)^T\). Hence

\widehat{\operatorname{Var}}(q(\hat\theta)) \approx \frac{1}{n}\left\{ 0.6^2(0.0004)+2(0.6)(0.7)(0.0001)+0.7^2(0.0009) \right\} =\frac{0.000669}{n}.

Since \(w=1/0.42\), its gradient is \(-\nabla q/q^2\), and its variance is the product variance multiplied by \(1/q^4\). As \(q^4=0.42^4=0.03111696\),

\widehat{\operatorname{Var}}(w(\hat\theta)) \approx \frac{0.000669}{0.03111696\,n} \approx \frac{0.02150}{n}.

This calculation also checks the inverse derivative: differentiating \(1/(p_1p_2)\) gives components \(-1/(p_1^2p_2)\) and \(-1/(p_1p_2^2)\), equal to \(-p_2/q^2\) and \(-p_1/q^2\).

Logistic factors and practical diagnostics

For a binary treatment at time \(t\), write \(A_{it}\in\{0,1\}\), let \(p_{it}(\theta)=\operatorname{expit}\{x_{it}^T\theta\}\), and let the probability of the observed value be \(u_{it}(\theta)=p_{it}(\theta)^{A_{it}}\{1-p_{it}(\theta)\}^{1-A_{it}}\). Differentiating its log gives

\nabla\log u_{it}(\theta) =(A_{it}-p_{it}(\theta))x_{it}.

Thus the derivative of the log product is the sum of these score-like contributions across times. If parameters are divided into time-specific blocks, each contribution has zeros in the other blocks; the covariance matrix can still have off-diagonal blocks. If parameters are shared, a contribution can affect more than one time’s probability.

Useful checks include verifying that every probability factor corresponds to the observed treatment and its intended history, checking that all factors are positive, and comparing the analytic gradient with a numerical derivative at the fitted parameter. Inspect the covariance matrix and the individual contributions to \(s_i^T\hat V s_i\). A large estimated standard error, an extreme inverse weight, or strong sensitivity to small changes in the fitted probabilities warns that the local linear approximation may be fragile.

Small products: As \(q_i\) approaches zero, the inverse \(w_i\) becomes large, and its absolute uncertainty can increase sharply. A normal approximation for the inverse may then be poor even if the parameter estimator is approximately normal.
Boundary probabilities: If a factor is zero, the inverse is undefined; probabilities close to zero can also make numerical derivatives and uncertainty estimates unstable.
Model and data uncertainty: This calculation propagates uncertainty in \(\hat\theta\) for the specified subject and fitted model. It does not, by itself, account for model selection, misspecification, or additional sources of uncertainty not included in \(V\).
Changing number of times: The formulas use fixed \(K\). If the number of decision times grows with sample size, this fixed-dimensional argument alone does not justify the approximation.

The main alternative for assessing sampling uncertainty is bootstrap inference, which resamples the estimation procedure rather than relying only on a local linear approximation. It can reflect nonlinear behavior more directly, but its validity depends on the resampling design and it can be unstable when fitted probabilities are near zero. The delta method is transparent and computationally direct when the parameter covariance and gradient are available; neither approach repairs a poorly supported treatment model.

A common mistake

Common mistake: Compute a separate variance for each fitted probability, then add the variances as if the factors were independent. This is wrong whenever their estimation errors are correlated or they share model parameters. For example, with two factors \(u_1,u_2\), the product gradient has one component for each parameter block, and the variance includes the cross term \(2g_1g_2\operatorname{Cov}(\hat\theta_1,\hat\theta_2)\). The correction is to differentiate the whole product and use the full joint covariance matrix. Independence may be imposed only when it is justified by the estimator’s actual covariance structure.

Check Your Understanding

Answer each question, then compare with the solutions below.

For \(q(\theta)=\prod_t u_t(\theta)\), state its gradient using the gradients of the log factors.
Why is differentiability at \(\theta_0\) enough for the delta-method limit but not, by itself, enough for the plug-in gradient variance estimate to be consistent?
Two estimated log-probability factors have gradient contributions \(0.2\) and \(0.5\), and their parameter covariance matrix is the one in the first worked example. What is the quadratic form \(s^TVs\)?
If a positive product is multiplied by a factor of \(1/2\) while the log-product gradient remains fixed, how does the first-order variance of its inverse change?

Solutions. (1) \(\nabla q=q\sum_t\nabla\log u_t\), by differentiating \(\log q=\sum_t\log u_t\). (2) Differentiability at a point need not make the gradient continuous nearby, so it does not ensure \(\nabla g(\hat\theta)\xrightarrow{p}\nabla g(\theta_0)\); gradient continuity, together with consistency of \(\hat V\), supplies the needed plug-in convergence. (3) \(s^TVs=0.2^2(0.04)+2(0.2)(0.5)(0.01)+0.5^2(0.09)=0.0261\). (4) The inverse variance is \(w^2s^TVs/n\). Halving the product doubles \(w\), so the variance is multiplied by \(2^2=4\).