Skip to contents

Given an already-fitted logistic regression, computes a Huber-White/Eicker sandwich (heteroskedasticity-robust, "HC0") coefficient covariance matrix and, by the delta method, standard errors for the G-computed standardized risk difference and log risk ratio (see gcomp_logistic_point_estimate_cpp for the standardization point estimates this builds inference around). The sandwich covariance is \(\widehat{\mathrm{Var}}(\hat\beta) = B\,M\,B\), with "bread" \(B = (X^\top W X)^{-1}\) (\(W = \mathrm{diag}(\hat\mu_i(1-\hat\mu_i))\), the model-based Fisher information weights) and "meat" \(M = X^\top \mathrm{diag}((y_i - \hat\mu_i)^2) X\) (the empirical score outer-product, making this robust to model misspecification, not just relying on the working Bernoulli variance). The risk-difference standard error is obtained by propagating this sandwich covariance through the standardized risks' gradients with respect to \(\beta\) (\(\nabla_\beta \bar{\mathrm{risk}}_1 - \nabla_\beta \bar{\mathrm{risk}}_0\), each a population-averaged logistic-derivative weighted design-matrix sum, with the treatment column's gradient entry replaced by the sum of the standardized-risk derivative directly since every subject's treatment indicator is held fixed at 1 or 0 in the counterfactual averages); the log risk ratio's standard error is obtained the same way via the gradient of \(\log(\overline{\mathrm{risk}}_1) - \log(\overline{\mathrm{risk}}_0)\), and is only computed (non-NA) when both standardized risks are strictly positive. Aborts with an R error (rather than returning NAs) if mu_hat contains non-finite or boundary (0 or 1) values, if the weighted design crossproduct \(X^\top W X\) is not invertible, if the sandwich covariance comes out non-finite anywhere, or if the treatment coefficient's variance is non-positive.

Usage

gcomp_logistic_post_fit_cpp(X_fit, y, coef_hat, mu_hat, j_treat)

Arguments

X_fit

Numeric matrix of predictors used to fit the model, including an intercept column if the model has one.

y

The observed binary (0/1) response used to fit the model.

coef_hat

Numeric vector of fitted logistic regression coefficients \(\hat\beta\), same length and column order as X_fit.

mu_hat

Numeric vector of fitted probabilities \(\hat\mu_i = \mathrm{logit}^{-1}(x_i^\top\hat\beta)\), one per row of X_fit.

j_treat

1-based column index of the treatment indicator in X_fit.

Value

A list with components vcov (the \(p \times p\) sandwich covariance matrix), std_err and z_vals (per-coefficient standard errors and Wald z-statistics, NA for any coefficient with non-finite or non-positive variance), risk1/risk0 (the standardized mean risks under \(T=1\)/\(T=0\) for everyone), rd (risk1 - risk0) and se_rd (its delta-method standard error), and log_rr/rr/ se_log_rr (the log risk ratio, risk ratio, and the log risk ratio's delta-method standard error — all NA if either standardized risk is not strictly positive).

See also

gcomp_logistic_point_estimate_cpp for the point-estimate computation this function's variances are built around; gcomp_fractional_logit_post_fit_cpp() for the analogous fractional-logit/quasi-binomial post-fit inference.