Generalized Additive Models for Propensity Scores

Why Use Additive Smooth Effects?

A logistic propensity model can be too restrictive when it assumes that every continuous covariate has a linear effect on the log odds of treatment. A generalized additive model (GAM) relaxes that restriction by allowing the log odds to change smoothly with each covariate, while retaining an additive structure. This gives the model more flexibility than a linear logistic specification without automatically adding every possible interaction.

The goal remains estimating treatment probabilities that are useful for the intended causal design, not maximizing treatment classification accuracy. Flexible fitted probabilities still need to be examined for overlap and for their performance in balancing measured covariates. The earlier material on functional forms, restricted cubic splines, and treatment-prediction objectives motivates these distinctions.

Definition: Let $D_i$ be a binary treatment indicator and $X_i=(X_{i1},\ldots,X_{ip})$ the measured covariates for unit $i$. The propensity score is $e(X_i)=P(D_i=1\mid X_i)$. A generalized additive logistic propensity model specifies $$ \operatorname{logit}\{e(X_i)\}=\log\left\{\frac{e(X_i)}{1-e(X_i)}\right\} =\beta_0+\sum_{j=1}^{p} f_j(X_{ij}), $$ where each $f_j$ is a smooth function. The fitted propensity score is $\hat e(X_i)=\operatorname{expit}\{\hat\beta_0+\sum_j\hat f_j(X_{ij})\}$, where $\operatorname{expit}(u)=1/(1+\exp(-u))$.

“Additive” describes the predictor: the smooth contributions are summed on the log-odds scale. It does not assert that treatment effects or outcome relationships are additive. Unless interactions are explicitly included, the model also does not represent arbitrary dependence of one covariate’s effect on another covariate.

Basis Functions and Smoothing

A smooth function is represented computationally by a finite basis expansion. For a continuous covariate $X_{ij}$, write $f_j(X_{ij})=\mathbf b_j(X_{ij})^{\mathsf T}\boldsymbol\theta_j$, where $\mathbf b_j$ is a vector of basis functions and $\boldsymbol\theta_j$ is a coefficient vector. Splines are a common basis choice. Identifiability constraints, such as centering a smooth so that its sample contributions sum to zero, separate the smooth terms from the intercept.

A basis with many coefficients can represent a highly wiggly function. A quadratic roughness penalty discourages such complexity. For $n$ independent units, let $\mathbf B$ be the design matrix containing the intercept and all smooth-basis columns, $\boldsymbol\theta$ the combined coefficient vector, and $\boldsymbol\eta=\mathbf B\boldsymbol\theta$. Let $S_j$ be the penalty matrix for smooth $j$, and $\lambda_j\geq 0$ its smoothing parameter. A larger $\lambda_j$ penalizes roughness more strongly; a smaller value allows more flexibility. Each smooth may have its own smoothing parameter.

With $D_i\in\{0,1\}$, the Bernoulli log likelihood and penalized objective are

\ell(\boldsymbol\theta)=\sum_{i=1}^{n}\left[D_i\eta_i-\log\{1+\exp(\eta_i)\}\right], \qquad \ell_{\mathrm{pen}}(\boldsymbol\theta)=\ell(\boldsymbol\theta) -\frac{1}{2}\sum_{j=1}^{p}\lambda_j\boldsymbol\theta_j^{\mathsf T}S_j\boldsymbol\theta_j.

Intercept and any unpenalized basis components have zero penalty blocks. Smoothing parameters are generally selected by a method such as cross-validation or a likelihood-based criterion; their selection is a modeling choice, not evidence that the propensity model is causally adequate. The earlier material on cross-validation covers the general tuning idea. In practice, the basis dimension must also be adequate: a penalty cannot recover features that the chosen basis is unable to represent.

What the Penalized Fit Optimizes

The score equations show how the data and smoothness penalties jointly determine the fitted model. Write $\hat{\mathbf e}=(\hat e_1,\ldots,\hat e_n)^{\mathsf T}$, where $\hat e_i=\operatorname{expit}(\eta_i)$. Differentiating the Bernoulli log likelihood with respect to $\boldsymbol\theta$ gives $\mathbf B^{\mathsf T}(\mathbf D-\hat{\mathbf e})$. Define $P$ as the block-diagonal matrix with zero blocks for unpenalized coefficients and blocks $\lambda_j S_j$ for the smooth coefficients. Then

\nabla_{\boldsymbol\theta}\ell_{\mathrm{pen}} =\mathbf B^{\mathsf T}(\mathbf D-\hat{\mathbf e})-P\boldsymbol\theta.

This expression keeps each smooth’s penalty paired with its own smoothing parameter. To see the likelihood part, $\partial\eta_i/\partial\boldsymbol\theta$ is row $i$ of $\mathbf B$, and $\partial\ell_i/\partial\eta_i=D_i-\operatorname{expit}(\eta_i)$; summing those contributions gives $\mathbf B^{\mathsf T}(\mathbf D-\hat{\mathbf e})$. Differentiating the quadratic penalty gives $P\boldsymbol\theta$, which is subtracted because the objective is maximized.

A second derivative gives the penalized Hessian

\nabla^2_{\boldsymbol\theta}\ell_{\mathrm{pen}} =-\mathbf B^{\mathsf T}W\mathbf B-P, \qquad W=\operatorname{diag}\{\hat e_i(1-\hat e_i)\}_{i=1}^{n}.

Both $W$ and the penalty matrices are positive semidefinite, so the Hessian is negative semidefinite. Thus the penalized objective is concave in the coefficients for fixed smoothing parameters. Numerical fitting commonly uses iterative weighted least squares or a related optimizer. Concavity does not guarantee a unique finite estimate in every design: rank deficiencies, unpenalized directions, or separation-like patterns can still cause problems.

Worked Example: Evaluating the Penalized Gradient

Consider four hypothetical units with design-matrix rows $(1,z_{i1},z_{i2})$, where the two smooth-basis columns are $z_1=(-1,-1,1,1)^{\mathsf T}$ and $z_2=(-1,1,-1,1)^{\mathsf T}$. Let $\mathbf D=(1,1,0,0)^{\mathsf T}$, and evaluate the gradient at $\boldsymbol\theta=(0,0.2,0.3)^{\mathsf T}$. Suppose the penalty matrices for the two smooth coefficients are each the scalar 1, with distinct smoothing parameters $\lambda_1=0.2$ and $\lambda_2=0.8$.

The four linear predictors are $(-0.5,0.1,-0.1,0.5)$. Applying the expit function gives fitted probabilities approximately $(0.3775,0.5250,0.4750,0.6225)$. Therefore, $\mathbf D-\hat{\mathbf e}$ is approximately $(0.6225,0.4750,-0.4750,-0.6225)^{\mathsf T}$.

The unpenalized gradient components are the residual sum and the two residual dot products: $0$, $-0.6225-0.4750-0.4750-0.6225=-2.1950$, and $-0.6225+0.4750+0.4750-0.6225=-0.2949$. The penalty gradient components are $0$, $0.2(0.2)=0.0400$, and $0.8(0.3)=0.2400$. Subtracting them yields the penalized gradient $(0,-2.2350,-0.5349)^{\mathsf T}$, up to rounding. The distinct penalty contributions matter: replacing both smoothing parameters with one scalar would not differentiate the stated objective.

A nonzero gradient means these coefficients are not a stationary point for the chosen penalties and data. An optimizer would update them; the calculation alone does not provide a final fitted model.

Worked Probability Comparisons

Worked Example: A Smooth Age Effect Versus a Linear Effect

Suppose, purely for illustration, two fitted models from a hypothetical sample have the following forms. Age is scaled as $t=(\text{age}-50)/10$, and smoker equals 1 for a smoker and 0 otherwise. The GAM predictor is $-2+f(\text{age})+0.6(\text{smoker})$, with fitted smooth values $f(40)=-0.4$, $f(50)=0.1$, and $f(60)=0.8$. The comparison linear logistic predictor is $-2+0.2t+0.6(\text{smoker})$. These stipulated coefficients and smooth values illustrate predictions; they are not claimed to have been estimated from the displayed information.

For a 50-year-old non-smoker, the GAM linear predictor is $-2+0.1=-1.9$, so the fitted probability is $1/(1+\exp(1.9))\approx0.1301$. The linear model gives $-2+0.2(0)=-2$, or $1/(1+\exp(2))\approx0.1192$.

For a 60-year-old smoker, the GAM predictor is $-2+0.8+0.6=-0.6$, giving $1/(1+\exp(0.6))\approx0.3543$. The linear model gives $-2+0.2(1)+0.6=-1.2$, giving $1/(1+\exp(1.2))\approx0.2315$. The difference arises because the GAM’s stipulated age contribution rises more sharply between ages 50 and 60 than the linear contribution does. These are estimated treatment probabilities under the two hypothetical fitted specifications, not treatment effects.

Worked Example: Checking Grouped Calibration and Balance Separately

Imagine a hypothetical analysis with four age bands, each containing 10 units. The observed treatment counts are 1, 3, 7, and 8, so the observed treatment fractions are $0.10,0.30,0.70,0.80$. Suppose a linear logistic model predicts average probabilities $0.18,0.38,0.58,0.78$ in these bands, while a GAM predicts $0.11,0.32,0.72,0.80$.

As a simple descriptive grouped-calibration summary, calculate the mean absolute difference between each band’s observed fraction and its predicted average. For the linear model this is $(|0.10-0.18|+|0.30-0.38|+|0.70-0.58|+|0.80-0.78|)/4=(0.08+0.08+0.12+0.02)/4=0.075$. For the GAM it is $(0.01+0.02+0.02+0)/4=0.0125$. On these four grouped summaries, the GAM is closer to the observed treatment fractions.

This calculation does not establish covariate balance: it checks grouped agreement between predicted and observed treatment frequencies. Balance must be assessed directly for the measured covariates after applying the intended matching or weighting design. Nor does this in-sample grouped comparison establish performance on new data; evaluation should respect the fitting and tuning procedure.

Assumptions and Diagnostics

The response model assumes a binary treatment recorded consistently and a suitable independent-unit likelihood conditional on the included covariates. The additive specification assumes that the log odds can be represented adequately by the chosen sum of smooth functions, unless relevant interactions are included. The basis, centering constraints, and penalty must be implemented correctly, and the smoothing selection procedure should be appropriate for the data. These are modeling assumptions; they are not guarantees of the causal conditions.

For a causal interpretation of propensity-score adjustment, the usual requirements include consistency, conditional exchangeability given measured covariates, and positivity for the target population. Model flexibility cannot repair unmeasured confounding or structural absence of treated or control units at some covariate values. In particular, a GAM can produce numerically stable probabilities in regions where the data provide little support; such smooth extrapolation should not be mistaken for empirical overlap.

Useful diagnostics include plotting each estimated smooth with its uncertainty or effective degrees of freedom, checking sensitivity to basis dimension and smoothing choices, and examining the distribution of fitted scores by treatment group. Very wiggly fits, sharp local changes driven by few observations, probabilities extremely near zero or one, and materially different predictions under modest tuning changes are warning signs. After fitting, assess covariate balance using the design’s relevant balance diagnostics. A visually plausible smooth or good treatment classification does not imply balance.

Specify the additive structure.
Choose measured covariates and decide which continuous effects need smooth terms, which can remain linear, and whether substantive interactions are needed.

Fit with controlled flexibility.
Choose adequate basis dimensions and select each smooth’s penalty according to a stated tuning procedure. Inspect convergence and fitted probabilities.

Evaluate the causal design, not only prediction.
Check overlap and covariate balance under the intended use of the scores. Revise the model or target population if diagnostics show a problem, and report limitations.

How GAMs Compare and Where They Fall Short

A linear logistic propensity model is simpler to communicate and can be stable when its functional form is adequate, but a linear log-odds assumption may miss curvature. A GAM permits flexible univariate shapes while preserving an additive predictor. A fixed spline logistic model also represents nonlinear effects, but its basis and degrees of freedom are typically specified more directly; a GAM’s smoothing penalty offers a systematic way to regulate roughness. These approaches do not remove the need for covariate-balance diagnostics.

The additive restriction is also a limitation: strong effect modification among covariates may not be captured without interactions. Flexible smooths can become unstable in sparse regions, especially with limited data or many covariates, and estimated scores can still be extreme. The additive form may be easier to interpret than a highly interactive model, but it is not automatically the correct form. Model choice should be guided by design knowledge, adequate overlap, and observed balance rather than by classification performance alone.

Common Mistake

Common mistake: “The GAM predicts treatment more accurately, so it must be the better propensity model and will balance the covariates.” Better treatment prediction is not the target by itself, and neither a flexible fit nor a good prediction score proves balance. Compare candidate specifications using relevant balance and overlap diagnostics, while retaining the causal assumptions and target population in view.

Check Your Understanding

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

In a GAM propensity model, are the smooth functions added on the probability scale or the log-odds scale?
What happens to the penalty’s influence on a smooth term as its smoothing parameter increases?
A fitted GAM has excellent grouped calibration but leaves a large standardized difference in a covariate after weighting. Does calibration establish that the design is adequate?
Why can a GAM fail to represent a relationship that depends jointly on two covariates?

Solutions. (1) The smooth functions are added on the log-odds scale, and the expit transformation converts the total predictor to a probability. (2) A larger smoothing parameter penalizes roughness more strongly, usually producing a less variable smooth, subject to the basis and penalty structure. (3) No. Grouped calibration compares predicted and observed treatment frequencies; the remaining covariate imbalance is a separate diagnostic problem. (4) An additive predictor sums functions of individual covariates and does not generally include joint dependence; an interaction or another suitable term must be specified to represent it.