|
| 1 | +from typing import Callable, Optional |
| 2 | + |
| 3 | +import numpy as np |
| 4 | +from joblib import Parallel, delayed |
| 5 | +from numpy.typing import ArrayLike |
| 6 | +from sklearn.base import BaseEstimator |
| 7 | +from sklearn.linear_model import LogisticRegression |
| 8 | + |
| 9 | +from pywhy_stats.kernel_utils import _default_regularization |
| 10 | + |
| 11 | + |
| 12 | +def _preprocess_propensity_data( |
| 13 | + group_ind: ArrayLike, |
| 14 | + propensity_model: Optional[BaseEstimator], |
| 15 | + propensity_weights: Optional[ArrayLike], |
| 16 | +): |
| 17 | + if group_ind.ndim != 1: |
| 18 | + raise RuntimeError("group_ind must be a 1d array.") |
| 19 | + if len(np.unique(group_ind)) != 2: |
| 20 | + raise RuntimeError( |
| 21 | + f"There should only be two groups. Found {len(np.unique(group_ind))} groups." |
| 22 | + ) |
| 23 | + if propensity_model is not None and propensity_weights is not None: |
| 24 | + raise ValueError( |
| 25 | + "Both propensity model and propensity estimates are specified. Only one is allowed." |
| 26 | + ) |
| 27 | + if propensity_weights is not None: |
| 28 | + if propensity_weights.shape[0] != len(group_ind): |
| 29 | + raise ValueError( |
| 30 | + f"There are {propensity_weights.shape[0]} pre-defined estimates, while " |
| 31 | + f"there are {len(group_ind)} samples." |
| 32 | + ) |
| 33 | + if propensity_weights.shape[1] != len(np.unique(group_ind.squeeze())): |
| 34 | + raise ValueError( |
| 35 | + f"There are {propensity_weights.shape[1]} group pre-defined estimates, while " |
| 36 | + f"there are {len(np.unique(group_ind))} unique groups." |
| 37 | + ) |
| 38 | + |
| 39 | + |
| 40 | +def _compute_propensity_scores( |
| 41 | + group_ind: ArrayLike, |
| 42 | + propensity_model: Optional[BaseEstimator] = None, |
| 43 | + propensity_weights: Optional[ArrayLike] = None, |
| 44 | + n_jobs: Optional[int] = None, |
| 45 | + random_state: Optional[int] = None, |
| 46 | + **kwargs, |
| 47 | +): |
| 48 | + if propensity_model is None: |
| 49 | + K: ArrayLike = kwargs.get("K") |
| 50 | + |
| 51 | + # compute a default penalty term if using a kernel matrix |
| 52 | + # C is the inverse of the regularization parameter |
| 53 | + if K.shape[0] == K.shape[1]: |
| 54 | + # default regularization is 1 / (2 * K) |
| 55 | + propensity_penalty_ = _default_regularization(K) |
| 56 | + C = 1 / (2 * propensity_penalty_) |
| 57 | + else: |
| 58 | + # defaults to no regularization |
| 59 | + propensity_penalty_ = 0.0 |
| 60 | + C = 1.0 |
| 61 | + |
| 62 | + # default model is logistic regression |
| 63 | + propensity_model_ = LogisticRegression( |
| 64 | + penalty="l2", |
| 65 | + n_jobs=n_jobs, |
| 66 | + warm_start=True, |
| 67 | + solver="lbfgs", |
| 68 | + random_state=random_state, |
| 69 | + C=C, |
| 70 | + ) |
| 71 | + else: |
| 72 | + propensity_model_ = propensity_model |
| 73 | + |
| 74 | + # either use pre-defined propensity weights, or estimate them |
| 75 | + if propensity_weights is None: |
| 76 | + K = kwargs.get("K") |
| 77 | + # fit and then obtain the probabilities of treatment |
| 78 | + # for each sample (i.e. the propensity scores) |
| 79 | + propensity_weights = propensity_model_.fit(K, group_ind.ravel()).predict_proba(K)[:, 1] |
| 80 | + else: |
| 81 | + propensity_weights = propensity_weights[:, 1] |
| 82 | + return propensity_weights |
| 83 | + |
| 84 | + |
| 85 | +def compute_null( |
| 86 | + func: Callable, |
| 87 | + e_hat: ArrayLike, |
| 88 | + X: ArrayLike, |
| 89 | + Y: ArrayLike, |
| 90 | + null_reps: int = 1000, |
| 91 | + n_jobs=None, |
| 92 | + seed=None, |
| 93 | + **kwargs, |
| 94 | +) -> ArrayLike: |
| 95 | + """Estimate null distribution using propensity weights. |
| 96 | +
|
| 97 | + Parameters |
| 98 | + ---------- |
| 99 | + func : Callable |
| 100 | + The function to compute the test statistic. |
| 101 | + e_hat : Array-like of shape (n_samples,) |
| 102 | + The predicted propensity score for ``group_ind == 1``. |
| 103 | + X : Array-Like of shape (n_samples, n_features_x) |
| 104 | + The X (covariates) array. |
| 105 | + Y : Array-Like of shape (n_samples, n_features_y) |
| 106 | + The Y (outcomes) array. |
| 107 | + null_reps : int, optional |
| 108 | + Number of times to sample null, by default 1000. |
| 109 | + n_jobs : int, optional |
| 110 | + Number of jobs to run in parallel, by default None. |
| 111 | + seed : int, optional |
| 112 | + Random generator, or random seed, by default None. |
| 113 | +
|
| 114 | + Returns |
| 115 | + ------- |
| 116 | + null_dist : Array-like of shape (n_samples,) |
| 117 | + The null distribution of test statistics. |
| 118 | + """ |
| 119 | + rng = np.random.default_rng(seed) |
| 120 | + n_samps = X.shape[0] |
| 121 | + |
| 122 | + # compute the test statistic on the conditionally permuted |
| 123 | + # dataset, where each group label is resampled for each sample |
| 124 | + # according to its propensity score |
| 125 | + null_dist = Parallel(n_jobs=n_jobs)( |
| 126 | + [ |
| 127 | + delayed(func)(X, Y, group_ind=rng.binomial(1, e_hat, size=n_samps), **kwargs) |
| 128 | + for _ in range(null_reps) |
| 129 | + ] |
| 130 | + ) |
| 131 | + return np.asarray(null_dist) |
0 commit comments