diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 2421ad6..ec4e173 100644 --- a/.github/workflows/python-app.yml +++ b/.github/workflows/python-app.yml @@ -17,7 +17,7 @@ jobs: fail-fast: false matrix: os: [macos-26] - python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14", "3.15"] steps: - uses: actions/checkout@v7 @@ -28,6 +28,19 @@ jobs: - name: Set up Python ${{ matrix.python-version }} run: uv python install ${{ matrix.python-version }} + # tabmat and glum have no Python 3.15 wheels yet, so they build from source. + # lifelines caps pandas below 3, and pandas 2.x is broken on 3.15, so override + # the cap and stop uv run from re-syncing it back. + - name: Prepare Python 3.15 build + if: matrix.python-version == '3.15' + run: | + brew install libomp xsimd jemalloc + echo "CPPFLAGS=-I$(brew --prefix libomp)/include -I$(brew --prefix xsimd)/include -I$(brew --prefix jemalloc)/include" >> "$GITHUB_ENV" + echo "LDFLAGS=-L$(brew --prefix libomp)/lib -L$(brew --prefix jemalloc)/lib" >> "$GITHUB_ENV" + echo "pandas>=3" > "$RUNNER_TEMP/uv-overrides.txt" + echo "UV_OVERRIDE=$RUNNER_TEMP/uv-overrides.txt" >> "$GITHUB_ENV" + echo "UV_NO_SYNC=1" >> "$GITHUB_ENV" + - name: Install dependencies run: | uv venv --python ${{ matrix.python-version }} diff --git a/CITATION.cff b/CITATION.cff index 28d5beb..9315375 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -18,7 +18,7 @@ authors: affiliation: "CAUSALab, Department of Epidemiology, Harvard T.H. Chan School of Public Health; Department of Biostatistics, Harvard T.H. Chan School of Public Health" orcid: https://orcid.org/0000-0003-1619-8456 title: "pySEQTarget: Sequential Target Trial Emulation" -version: 0.14.2 +version: 0.14.3 url: https://pyseqtarget.readthedocs.io/ repository-code: https://github.com/CausalInference/pySEQTarget repository-artifact: https://pypi.org/project/pySEQTarget/ diff --git a/docs/vignettes/SEQuential_results.png b/docs/vignettes/SEQuential_results.png index 7388757..f996c93 100644 Binary files a/docs/vignettes/SEQuential_results.png and b/docs/vignettes/SEQuential_results.png differ diff --git a/docs/vignettes/end_of_followup.md b/docs/vignettes/end_of_followup.md index b25b977..97af634 100644 --- a/docs/vignettes/end_of_followup.md +++ b/docs/vignettes/end_of_followup.md @@ -51,8 +51,8 @@ model.eof_data.select(["A", "Proportion", "SE", "95% LCI", "95% UCI", ┌─────┬────────────┬──────────┬──────────┬──────────┬────────────┬──────────┐ │ A ┆ Proportion ┆ SE ┆ 95% LCI ┆ 95% UCI ┆ % Censored ┆ Subjects │ ╞═════╪════════════╪══════════╪══════════╪══════════╪════════════╪══════════╡ -│ 0 ┆ 0.019865 ┆ 0.003014 ┆ 0.013956 ┆ 0.025773 ┆ 19.687712 ┆ 266 │ -│ 1 ┆ 0.027256 ┆ 0.002342 ┆ 0.022666 ┆ 0.031847 ┆ 28.46412 ┆ 236 │ +│ 0 ┆ 0.019865 ┆ 0.003014 ┆ 0.013956 ┆ 0.025773 ┆ 17.098809 ┆ 266 │ +│ 1 ┆ 0.027256 ┆ 0.002342 ┆ 0.022666 ┆ 0.031847 ┆ 24.81102 ┆ 236 │ └─────┴────────────┴──────────┴──────────┴──────────┴────────────┴──────────┘ ``` @@ -108,8 +108,8 @@ windowed.eof_data.select(["A", "Proportion", "95% LCI", "95% UCI", │ A ┆ Proportion ┆ 95% LCI ┆ 95% UCI ┆ Trial-period ┆ Trial-period ┆ Trial-period ┆ % Censored │ │ ┆ ┆ ┆ ┆ s (Eligible) ┆ s (Analysed) ┆ s (Censored) ┆ │ ╞═════╪════════════╪══════════╪══════════╪══════════════╪══════════════╪══════════════╪════════════╡ -│ 0 ┆ 0.075704 ┆ 0.058602 ┆ 0.092805 ┆ 2946 ┆ 2523 ┆ 423 ┆ 14.358452 │ -│ 1 ┆ 0.102475 ┆ 0.0864 ┆ 0.11855 ┆ 6257 ┆ 4889 ┆ 1368 ┆ 21.863513 │ +│ 0 ┆ 0.075704 ┆ 0.058602 ┆ 0.092805 ┆ 2894 ┆ 2523 ┆ 371 ┆ 12.819627 │ +│ 1 ┆ 0.102475 ┆ 0.0864 ┆ 0.11855 ┆ 6073 ┆ 4889 ┆ 1184 ┆ 19.49613 │ └─────┴────────────┴──────────┴──────────┴──────────────┴──────────────┴──────────────┴────────────┘ ``` @@ -150,8 +150,8 @@ windowed.diagnostics["nonunique_eof"] │ tx_init_bas ┆ Eligible ┆ At k ┆ In window ┆ Excluded (outside ┆ Excluded (no measurement) │ │ ┆ ┆ ┆ ┆ window) ┆ │ ╞═════════════╪══════════╪══════╪═══════════╪═══════════════════╪═══════════════════════════╡ -│ 0 ┆ 2946 ┆ 2366 ┆ 157 ┆ 423 ┆ 0 │ -│ 1 ┆ 6257 ┆ 4476 ┆ 413 ┆ 1368 ┆ 0 │ +│ 0 ┆ 2894 ┆ 2366 ┆ 157 ┆ 371 ┆ 0 │ +│ 1 ┆ 6073 ┆ 4476 ┆ 413 ┆ 1184 ┆ 0 │ └─────────────┴──────────┴──────┴───────────┴───────────────────┴───────────────────────────┘ ``` @@ -272,8 +272,8 @@ perprotocol.eof_data.select(["A", "Proportion", "95% LCI", "95% UCI", ┌─────┬────────────┬──────────┬──────────┬────────────┬──────────┐ │ A ┆ Proportion ┆ 95% LCI ┆ 95% UCI ┆ % Censored ┆ Subjects │ ╞═════╪════════════╪══════════╪══════════╪════════════╪══════════╡ -│ 0 ┆ 0.045069 ┆ 0.011689 ┆ 0.07845 ┆ 82.959946 ┆ 86 │ -│ 1 ┆ 0.088733 ┆ 0.065572 ┆ 0.111894 ┆ 49.592456 ┆ 211 │ +│ 0 ┆ 0.043282 ┆ 0.012381 ┆ 0.074182 ┆ 82.653766 ┆ 86 │ +│ 1 ┆ 0.088984 ┆ 0.065753 ┆ 0.112215 ┆ 48.065207 ┆ 211 │ └─────┴────────────┴──────────┴──────────┴────────────┴──────────┘ ``` diff --git a/docs/vignettes/exploring_results.md b/docs/vignettes/exploring_results.md index f78ad06..f9fad98 100644 --- a/docs/vignettes/exploring_results.md +++ b/docs/vignettes/exploring_results.md @@ -12,72 +12,40 @@ Let us now go over what the dump to md looks like and explore our output in furt ## Weighting -We begin by exploring the weight models, this gives us general information about the numerator and denominator models, as well as weight statistics before applying any limits. If you recall, we imposed weight bounds at the 99th percentile. This means that in the outcome model our weights will be bound at [0.273721, 423.185]. Note that in real analysis, we would hope the weights are stabilized further. Using the generated data, specifically with excused analysis, often result in larger-than-intended weights. +We begin by exploring the weight models, this gives us general information about the denominator model, as well as weight statistics before applying any limits. With pre-expansion weights and excused conditions the numerator is fixed at 1, so there is no numerator model and the weights are not stabilized; they have a long right tail, which is why we bound them at the 1st and 99th percentiles. This means that in the outcome model our weights will be bound at [0.999954, 19.8889]. -We should also note here that in excused-censoring analysis, our adherance models hold `switch` as the dependent variable. In all non-excused cases, this would normally be your treatment value. - -### Numerator Model - -``` - MNLogit Regression Results -============================================================================== -Dep. Variable: switch No. Observations: 65375 -Model: MNLogit Df Residuals: 65366 -Method: MLE Df Model: 8 -Date: Wed, 10 Dec 2025 Pseudo R-squ.: 0.008332 -Time: 10:18:38 Log-Likelihood: -13986. -converged: True LL-Null: -14103. -Covariance Type: nonrobust LLR p-value: 2.560e-46 -=============================================================================== - switch=1 coef std err z P>|z| [0.025 0.975] -------------------------------------------------------------------------------- -Intercept -0.8797 0.692 -1.271 0.204 -2.236 0.477 -sex[T.1] -0.0461 0.035 -1.325 0.185 -0.114 0.022 -N_bas 0.0026 0.003 0.741 0.459 -0.004 0.009 -L_bas 0.3368 0.032 10.553 0.000 0.274 0.399 -P_bas -0.1864 0.073 -2.556 0.011 -0.329 -0.043 -followup -0.0211 0.006 -3.822 0.000 -0.032 -0.010 -followup_sq 0.0001 0.000 0.637 0.524 -0.000 0.000 -trial -0.0624 0.014 -4.430 0.000 -0.090 -0.035 -trial_sq 0.0004 0.000 2.309 0.021 5.78e-05 0.001 -=============================================================================== -``` +We should also note here that with pre-expansion weights the denominator model has the treatment as its dependent variable. With post-expansion weights (`weight_preexpansion=False`), excused-censoring models hold `switch` as the dependent variable instead. ### Denominator Model ``` - MNLogit Regression Results + Logit Regression Results +============================================================================== +Dep. Variable: tx_init No. Observations: 6137 +Model: Logit Df Residuals: 6130 +Method: MLE Df Model: 6 +Date: Tue, 13 Oct 2026 Pseudo R-squ.: 0.0006500 +Time: 10:31:40 Log-Likelihood: -3017.9 +converged: True LL-Null: -3019.8 +Covariance Type: nonrobust LLR p-value: 0.6867 +============================================================================== + coef std err z P>|z| [0.025 0.975] +------------------------------------------------------------------------------ +Intercept -2.2157 1.214 -1.824 0.068 -4.596 0.165 +sex 0.0331 0.065 0.511 0.609 -0.094 0.160 +N -0.0105 0.006 -1.648 0.099 -0.023 0.002 +L -0.0177 0.052 -0.343 0.732 -0.119 0.083 +P 0.0945 0.128 0.741 0.459 -0.156 0.345 +time 0.0147 0.024 0.609 0.542 -0.033 0.062 +time_sq -1.602e-05 0.000 -0.066 0.947 -0.000 0.000 ============================================================================== -Dep. Variable: switch No. Observations: 65375 -Model: MNLogit Df Residuals: 65363 -Method: MLE Df Model: 11 -Date: Wed, 10 Dec 2025 Pseudo R-squ.: 0.01586 -Time: 10:18:38 Log-Likelihood: -13880. -converged: True LL-Null: -14103. -Covariance Type: nonrobust LLR p-value: 5.374e-89 -=============================================================================== - switch=1 coef std err z P>|z| [0.025 0.975] -------------------------------------------------------------------------------- -Intercept -1.4384 0.715 -2.012 0.044 -2.840 -0.037 -sex[T.1] -0.0447 0.035 -1.281 0.200 -0.113 0.024 -N -0.0195 0.003 -5.655 0.000 -0.026 -0.013 -L 0.3719 0.062 6.025 0.000 0.251 0.493 -P 0.9362 0.139 6.723 0.000 0.663 1.209 -N_bas 0.0023 0.003 0.674 0.501 -0.004 0.009 -L_bas -0.1703 0.092 -1.842 0.066 -0.351 0.011 -P_bas -0.9966 0.139 -7.166 0.000 -1.269 -0.724 -followup 0.0906 0.022 4.164 0.000 0.048 0.133 -followup_sq -0.0007 0.000 -3.678 0.000 -0.001 -0.000 -trial -0.0695 0.014 -4.934 0.000 -0.097 -0.042 -trial_sq 0.0006 0.000 3.788 0.000 0.000 0.001 -=============================================================================== ``` ### Weighting Statistics | weight_min | weight_max | weight_mean | weight_std | weight_p01 | weight_p25 | weight_p50 | weight_p75 | weight_p99 | |-------------:|-------------:|--------------:|-------------:|-------------:|-------------:|-------------:|-------------:|-------------:| -| 3.11308e-08 | 9.08003e+30 | 1.97762e+26 | 3.39822e+28 | 0.260691 | 0.853367 | 1.02192 | 1.28444 | 30119 | +| 0.999923 | 24934.2 | 2.05553 | 33.0393 | 0.999954 | 0.999995 | 1.00002 | 1.26351 | 19.8889 | ## Outcome @@ -88,29 +56,26 @@ After weight information, we begin to gather information about the outcome model ``` Generalized Linear Model Regression Results ============================================================================== -Dep. Variable: outcome No. Observations: 658971 -Model: GLM Df Residuals: 658961 -Model Family: Binomial Df Model: 9 +Dep. Variable: outcome No. Observations: 670413 +Model: GLM Df Residuals: 670406 +Model Family: Binomial Df Model: 6 Link Function: Logit Scale: 1.0000 -Method: IRLS Log-Likelihood: -2844.7 -Date: Wed, 10 Dec 2025 Deviance: 5689.4 -Time: 10:18:38 Pearson chi2: 6.80e+05 -No. Iterations: 11 Pseudo R-squ. (CS): 0.0001638 +Method: IRLS Log-Likelihood: -3096.8 +Date: Tue, 13 Oct 2026 Deviance: 6193.6 +Time: 10:31:40 Pearson chi2: 1.02e+06 +No. Iterations: 13 Pseudo R-squ. (CS): 0.0002477 Covariance Type: nonrobust -==================================================================================== - coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------------- -Intercept -23.1102 2.697 -8.569 0.000 -28.396 -17.824 -tx_init_bas[T.1] -0.2221 0.185 -1.203 0.229 -0.584 0.140 -sex[T.1] -0.5588 0.113 -4.942 0.000 -0.780 -0.337 -followup 0.0060 0.015 0.416 0.678 -0.022 0.035 -followup_sq 0.0001 0.000 0.465 0.642 -0.000 0.001 -trial 0.3377 0.054 6.274 0.000 0.232 0.443 -trial_sq -0.0022 0.001 -4.223 0.000 -0.003 -0.001 -N_bas -0.0007 0.011 -0.066 0.947 -0.022 0.021 -L_bas -0.3595 0.072 -4.962 0.000 -0.501 -0.217 -P_bas 1.6141 0.281 5.752 0.000 1.064 2.164 -==================================================================================== +============================================================================================= + coef std err z P>|z| [0.025 0.975] +--------------------------------------------------------------------------------------------- +Intercept -9.5140 0.371 -25.647 0.000 -10.241 -8.787 +tx_init_bas[T.1] -0.0750 0.343 -0.219 0.827 -0.748 0.598 +followup -0.0560 0.024 -2.346 0.019 -0.103 -0.009 +tx_init_bas[T.1]:followup 0.0581 0.022 2.633 0.008 0.015 0.101 +followup_sq 0.0002 0.000 0.677 0.498 -0.000 0.001 +trial 0.1249 0.018 7.092 0.000 0.090 0.159 +trial_sq -0.0017 0.000 -4.600 0.000 -0.002 -0.001 +============================================================================================= ``` ### Survival @@ -125,21 +90,21 @@ survival_data = my_output.retrieve_data("km_data") #### Risk Differences -| A_x | A_y | Risk Difference | RD 95% LCI | RD 95% UCI | -|------:|------:|------------------:|-------------:|-------------:| -| 0 | 1 | 0.00859802 | -0.169438 | 0.186634 | -| 1 | 0 | -0.00859802 | -0.186634 | 0.169438 | +| Followup | A_x | A_y | Risk Difference | RD 95% LCI | RD 95% UCI | RD SE | +|-----------:|------:|------:|------------------:|-------------:|-------------:|----------:| +| 60 | 0 | 1 | -0.0244359 | -0.126522 | -0.00410246 | 0.0270593 | +| 60 | 1 | 0 | 0.0244359 | 0.00410246 | 0.126522 | 0.0270593 | #### Risk Ratios -| A_x | A_y | Risk Ratio | RR 95% LCI | RR 95% UCI | -|------:|------:|-------------:|-------------:|-------------:| -| 0 | 1 | 1.24069 | 0.0121904 | 126.272 | -| 1 | 0 | 0.806005 | 0.00791939 | 82.032 | +| Followup | A_x | A_y | Risk Ratio | RR 95% LCI | RR 95% UCI | log(RR) SE | +|-----------:|------:|------:|-------------:|-------------:|-------------:|-------------:| +| 60 | 0 | 1 | 0.255543 | 0.114726 | 0.632205 | 0.48761 | +| 60 | 1 | 0 | 3.91324 | 1.58177 | 8.71642 | 0.48761 | -#### Survival Curves +#### Cumulative Risk -![Kaplan-Meier Survival Curves](SEQuential_results.png) +![Cumulative Risk](SEQuential_results.png) ## Diagnostic Tables @@ -147,42 +112,56 @@ After all of our primary results, we are met with a few diagnostic tables. These Because we have an excused-censoring analysis, we are also provided with information about switches from treatment as well as how many of these switches were excused. -### Unique Outcomes +### Unique Outcomes (distinct subjects who had the outcome) | tx_init | outcome | len | |----------:|----------:|------:| -| 0 | 0 | 249 | +| 1 | 0 | 713 | +| 0 | 1 | 5 | | 1 | 1 | 8 | -| 0 | 1 | 4 | -| 1 | 0 | 715 | +| 0 | 0 | 252 | -### Nonunique Outcomes +### Nonunique Outcomes (total outcome events) | tx_init | outcome | len | |----------:|----------:|-------:| -| 0 | 1 | 73 | -| 1 | 0 | 546644 | -| 1 | 1 | 227 | -| 0 | 0 | 117007 | +| 1 | 0 | 553803 | +| 0 | 0 | 121031 | +| 0 | 1 | 70 | +| 1 | 1 | 209 | + +### Unique Followup (distinct subjects contributing follow-up) + +| tx_init_bas | len | +|--------------:|------:| +| 0 | 828 | +| 1 | 913 | + +### Nonunique Followup (person-time intervals) + +| tx_init_bas | len | +|--------------:|-------:| +| 0 | 59827 | +| 1 | 610586 | -### Unique Switches +### Unique Switches (distinct subjects who switched) | tx_init | isExcused | switch | len | |----------:|:------------|---------:|------:| -| 0 | True | 1 | 30 | -| 1 | False | 1 | 47 | -| 0 | False | 1 | 91 | -| 1 | True | 1 | 32 | -| 0 | False | 0 | 132 | -| 1 | False | 0 | 644 | +| 1 | False | 0 | 646 | +| 0 | True | 1 | 29 | +| 1 | True | 1 | 33 | +| 0 | False | 1 | 94 | +| 0 | False | 0 | 134 | +| 1 | False | 1 | 42 | -### Nonunique Switches +### Nonunique Switches (total switch intervals) | tx_init | isExcused | switch | len | |----------:|:------------|---------:|-------:| -| 0 | True | 0 | 22056 | -| 0 | False | 1 | 3724 | -| 1 | False | 1 | 1256 | -| 1 | False | 0 | 527107 | -| 1 | True | 0 | 18508 | -| 0 | False | 0 | 91300 | +| 0 | True | 0 | 23012 | +| 0 | False | 1 | 3394 | +| 1 | False | 1 | 1306 | +| 1 | False | 0 | 533916 | +| 0 | False | 0 | 94695 | +| 1 | True | 0 | 18790 | diff --git a/docs/vignettes/more_advanced_models.md b/docs/vignettes/more_advanced_models.md index d281a45..4231d1b 100644 --- a/docs/vignettes/more_advanced_models.md +++ b/docs/vignettes/more_advanced_models.md @@ -1,7 +1,7 @@ # More Advanced Analysis In getting started, we covered some of the basics for getting up and running on a simple analysis, but there are many options stored within `SEQuential`, or more aptly, many more parameters to play with in {py:class}`~pySEQTarget.SEQopts`. Let's cover a more in-depth analysis. -In this case, let's go over a censoring analysis with excused conditions and stabilized weighting, limiting weights to the 99th percentile, and adjusting for losses-to-followup. Futhermore, we are interested in bootstrapping our results to get a risk estimate with confidence bounds and for ease of computation, we are going to randomly downsample 30% of trials which did not initiate treatment. Because we are downsampling, we are additionally going to turn off the lag condition for our adherance weights. +In this case, let's go over a censoring analysis with excused conditions and inverse probability weighting, limiting weights to the 99th percentile, and adjusting for losses-to-followup. Futhermore, we are interested in bootstrapping our results to get a risk estimate with confidence bounds and for ease of computation, we are going to randomly downsample 30% of trials which did not initiate treatment. If you are coming from the R version, many arguments have been streamlined or inferred - take R's `bootstrap`, and `bootstrap.nboot` - these have been merged such that any `bootstrap_nboot` over 0 automatically starts the bootstrap initiation. @@ -16,14 +16,15 @@ from pySEQTarget.data import load_data data = load_data("SEQdata_LTFU") my_options = SEQopts( bootstrap_nboot = 20, # 20 bootstrap iterations (for demonstration only — use 500+ in practice) + bootstrap_CI_method = "percentile", # percentile intervals stay within [0, 1] cense_colname = "LTFU", # control for losses-to-followup as a censor excused = True, # allow excused treatment swapping excused_colnames = ["excusedZero", "excusedOne"], km_curves = True, # run survival estimates + plot_type = "risk", # plot cumulative risk rather than survival selection_random = True, # randomly sample treatment non-initiators selection_sample = 0.30, # sample 30% of treatment non-initiators weighted = True, # enables the weighting - weight_lag_condition=False, # turn off lag condition when weighting for adherance weight_p99 = True, # bounds weights by the 1st and 99th percentile weight_preexpansion = True # weights are predicted using pre-expansion data ) diff --git a/pySEQTarget/SEQoutput.py b/pySEQTarget/SEQoutput.py index 1d8ab5a..c283271 100644 --- a/pySEQTarget/SEQoutput.py +++ b/pySEQTarget/SEQoutput.py @@ -1,7 +1,7 @@ import tempfile from dataclasses import dataclass from pathlib import Path -from typing import List, Literal, Optional +from typing import Any, Dict, List, Literal, Optional, Union import matplotlib.figure import matplotlib.pyplot as plt @@ -25,8 +25,12 @@ class SEQoutput: :type numerator_models: List[ResultsWrapper] or None :param denominator_models: Denominator models, if applicable, from the weighting process :type denominator_models: List[ResultsWrapper] or None - :param compevent_models: Competing event models, if applicable - :type compevent_models: List[ResultsWrapper] or None + :param outcome_models: Outcome models, one per bootstrap replicate with the + main fit first; for subgroup analyses, a dict of such lists keyed by subgroup + :type outcome_models: List[ResultsWrapper], Dict[Any, List[ResultsWrapper]] or None + :param compevent_models: Competing event models, if applicable, in the same form + as ``outcome_models`` + :type compevent_models: List[ResultsWrapper], Dict[Any, List[ResultsWrapper]] or None :param weight_statistics: Weight statistics once returned back to the expanded dataset :type weight_statistics: dict or None :param hazard: Hazard ratio if applicable @@ -67,8 +71,10 @@ class SEQoutput: method: str = None numerator_models: List[ResultsWrapper] = None denominator_models: List[ResultsWrapper] = None - outcome_models: List[List[ResultsWrapper]] = None - compevent_models: List[List[ResultsWrapper]] = None + outcome_models: Union[List[ResultsWrapper], Dict[Any, List[ResultsWrapper]]] = None + compevent_models: Union[List[ResultsWrapper], Dict[Any, List[ResultsWrapper]]] = ( + None + ) weight_statistics: pl.DataFrame = None hazard: pl.DataFrame = None km_data: pl.DataFrame = None @@ -82,7 +88,7 @@ class SEQoutput: def plot(self) -> None: """ - Displays the Kaplan-Meier graph + Displays the survival, risk or incidence plot """ if self.km_graph is None: raise ValueError( @@ -96,9 +102,10 @@ def summary( type: Optional[ Literal["numerator", "denominator", "outcome", "compevent"] ] = None, - ) -> List: + ) -> Union[List, Dict[Any, List]]: """ - Returns a list of model summaries of either the numerator, denominator, outcome, or competing event models + Returns a list of model summaries of either the numerator, denominator, outcome, or competing event models. + For subgroup analyses, outcome and competing event summaries are a dict of lists keyed by subgroup. :param type: Indicator for which model list you would like returned :type type: str @@ -121,14 +128,19 @@ def summary( if self.options is not None and self.options.offload: loader = Offloader(enabled=True, dir=self.options.offload_dir) - summaries = [] - for model in models: - if model is None: - continue - if loader is not None: - model = loader.load_model(model) - summaries.append(model.summary()) - return summaries + def summarise(group): + summaries = [] + for model in group: + if model is None: + continue + if loader is not None: + model = loader.load_model(model) + summaries.append(model.summary()) + return summaries + + if isinstance(models, dict): + return {val: summarise(group) for val, group in models.items()} + return summarise(models) def retrieve_data( self, diff --git a/pySEQTarget/SEQuential.py b/pySEQTarget/SEQuential.py index 784d37f..f2f578b 100644 --- a/pySEQTarget/SEQuential.py +++ b/pySEQTarget/SEQuential.py @@ -16,6 +16,7 @@ from .expansion import _binder, _diagnostics, _dynamic, _random_selection from .helpers import (Offloader, _bake_model_formulas, _col_string, _format_time, bootstrap_loop) +from .helpers._fix_categories import _cast_fixed_categories from .initialization import (_cense_denominator, _cense_numerator, _denominator, _numerator, _outcome) from .plot import _survival_plot @@ -30,6 +31,16 @@ _DEFAULT_SEED = 0 +def _collect_models(self, kind): + """One model per replicate, main fit first; keyed by subgroup for subgroup analyses.""" + if self.subgroup_colname is None: + return [models[kind] for models in self.outcome_model] + return { + val: [models[val][kind] for models in self.outcome_model if val in models] + for val in self._unique_subgroups + } + + class SEQuential: """ Primary class initializer for SEQuential target trial emulation @@ -41,7 +52,8 @@ class SEQuential: :param treatment_col: Column name specifying treatment per time_col :param outcome_col: Column name specifying outcome per time_col :param time_varying_cols: Time-varying column names as covariates (BMI, Age, etc.) - :param fixed_cols: Fixed column names as covariates (Sex, YOB, etc.) + :param fixed_cols: Fixed column names as covariates (Sex, YOB, etc.). Numeric + columns enter as continuous; supply categorical ones as string or Categorical :param method: Method for analysis ['ITT', 'dose-response', or 'censoring'] :param parameters: Parameters to augment analysis, specified with ``pySEQTarget.SEQopts`` :type parameters: Optional[SEQopts] or None @@ -145,6 +157,7 @@ def expand(self): self.cense_eligible_colname, self.compevent_colname, self.visit_colname, + self.subgroup_colname, *self.weight_eligible_colnames, *self.excused_colnames, ] @@ -306,15 +319,18 @@ def fit(self) -> None: # models on `self` - so we keep the original polars frame for # the downstream steps rather than paying a pl.from_pandas() # round-trip per replicate. - WDT_pd = WDT_pl.to_pandas() - for col in self.fixed_cols: - if col in WDT_pd.columns: - WDT_pd[col] = WDT_pd[col].astype("category") + WDT_pd = _cast_fixed_categories(self.fixed_cols, WDT_pl.to_pandas()) _fit_LTFU(self, WDT_pd) _fit_visit(self, WDT_pd) - _fit_numerator(self, WDT_pd) - _fit_denominator(self, WDT_pd) + # A subject's first row has no previous treatment (tx_lag is a fill) + tx_WDT_pd = ( + WDT_pd[WDT_pd.groupby(self.id_col).cumcount() > 0] + if self.weight_preexpansion + else WDT_pd + ) + _fit_numerator(self, tx_WDT_pd) + _fit_denominator(self, tx_WDT_pd) if self.offload: _offload_weights(self, boot_idx) @@ -355,7 +371,7 @@ def fit(self) -> None: key: (m.params.values, list(m.model.exog_names)) for key, m in sg.items() } - for val, sg in zip(self._unique_subgroups, models_list) + for val, sg in models_list.items() } return models_list @@ -552,9 +568,9 @@ def collect(self) -> SEQoutput: outcome_models = None compevent_models = None elif self.outcome_model is not None: - outcome_models = [model["outcome"] for model in self.outcome_model] + outcome_models = _collect_models(self, "outcome") if self.compevent_colname is not None: - compevent_models = [model["compevent"] for model in self.outcome_model] + compevent_models = _collect_models(self, "compevent") else: compevent_models = None else: diff --git a/pySEQTarget/analysis/_endoffup.py b/pySEQTarget/analysis/_endoffup.py index 2fdf095..bb2e9de 100644 --- a/pySEQTarget/analysis/_endoffup.py +++ b/pySEQTarget/analysis/_endoffup.py @@ -12,6 +12,7 @@ import polars as pl from scipy import stats +from ..weighting._weight_stats import _weight_bounds from ._risk_estimates import _ci_label @@ -85,18 +86,26 @@ def _eof_measure(self, DT): def _eof_period_flags(self, DT): - """One row per trial-period, flagged with what it has available.""" + """One row per trial-period that can reach the window, flagged with what it has.""" k = self.end_of_fup_time w = self.end_of_fup_window valid = pl.col(self.outcome_col).is_not_null() in_window = valid & (pl.col("followup") >= k - w) & (pl.col("followup") <= k + w) + last = DT["period"].max() - return DT.group_by([self.id_col, "trial"] + _eof_group_cols(self, DT)).agg( - [ - (valid & (pl.col("followup") == k)).any().alias("at_k"), - in_window.any().alias("in_window"), - valid.any().alias("measured"), - ] + return ( + DT.group_by([self.id_col, "trial"] + _eof_group_cols(self, DT)) + .agg( + [ + (pl.col("period") - pl.col("followup")).min().alias("_start"), + (valid & (pl.col("followup") == k)).any().alias("at_k"), + in_window.any().alias("in_window"), + valid.any().alias("measured"), + ] + ) + # Trial-periods starting too late to reach the window before the data end + .filter(pl.col("_start") + k - w <= last) + .drop("_start") ) @@ -117,8 +126,9 @@ def _eof_estimate(self): DT = _eof_frame(self) by = _eof_group_cols(self, DT) + lower, upper = _weight_bounds(self) measured = _eof_measure(self, DT).with_columns( - pl.col("weight").clip(lower_bound=self.weight_min, upper_bound=self.weight_max) + pl.col("weight").clip(lower_bound=lower, upper_bound=upper) ) est = measured.group_by(by).agg( diff --git a/pySEQTarget/analysis/_hazard.py b/pySEQTarget/analysis/_hazard.py index 2a0f9c2..a3cbc52 100644 --- a/pySEQTarget/analysis/_hazard.py +++ b/pySEQTarget/analysis/_hazard.py @@ -8,28 +8,30 @@ from ..helpers._predict_model import _safe_predict from ._outcome_fit import _cast_categories +from ._subgroup_fit import _warn_missing_subgroup def _calculate_hazard(self): if self.subgroup_colname is None: - return _calculate_hazard_single(self, self.DT, idx=None, val=None) + return _calculate_hazard_single(self, self.DT, val=None) all_hazards = [] original_DT = self.DT - for i, val in enumerate(self._unique_subgroups): + for val in self._unique_subgroups: + _warn_missing_subgroup(self, val) subgroup_DT = original_DT.filter(pl.col(self.subgroup_colname) == val) - hazard = _calculate_hazard_single(self, subgroup_DT, i, val) + hazard = _calculate_hazard_single(self, subgroup_DT, val) all_hazards.append(hazard) self.DT = original_DT return pl.concat(all_hazards) -def _calculate_hazard_single(self, data, idx=None, val=None): +def _calculate_hazard_single(self, data, val=None): if self.seed is not None: self._rng = np.random.RandomState(self.seed) - full_log_hr = _hazard_handler(self, data, idx, 0, self._rng) + full_log_hr = _hazard_handler(self, data, val, 0, self._rng) if full_log_hr is None or np.isnan(full_log_hr): return _create_hazard_output(None, None, None, val, self) @@ -47,11 +49,11 @@ def _calculate_hazard_single(self, data, idx=None, val=None): # seed (always set since SEQuential pins a default) so each replicate's # RNG — and therefore the result — is identical to the serial path. if getattr(self, "parallel", False) and self.seed is not None: - boot_log_hrs = _parallel_boot_log_hrs(self, data, idx, boot_sample_idx) + boot_log_hrs = _parallel_boot_log_hrs(self, data, val, boot_sample_idx) else: boot_log_hrs = [] for model_pos, sample_idx in enumerate(boot_sample_idx): - boot_log_hr = _one_boot_log_hr(self, data, idx, model_pos, sample_idx) + boot_log_hr = _one_boot_log_hr(self, data, val, model_pos, sample_idx) if boot_log_hr is not None and not np.isnan(boot_log_hr): boot_log_hrs.append(boot_log_hr) @@ -108,7 +110,7 @@ def _truncate_to_first_event(tmp, id_col, event_col): ) -def _one_boot_log_hr(self, data, idx, model_pos, sample_idx): +def _one_boot_log_hr(self, data, val, model_pos, sample_idx): """Build one bootstrap resample of ``data`` and return its log hazard ratio. The RNG is rebuilt from ``seed + sample_idx + 1`` (matching the serial loop @@ -134,7 +136,7 @@ def _one_boot_log_hr(self, data, idx, model_pos, sample_idx): .drop("_count", "_rep") .collect() ) - return _hazard_handler(self, boot_data, idx, model_pos + 1, rng) + return _hazard_handler(self, boot_data, val, model_pos + 1, rng) # Process-pool worker state. Set once per worker process by the initializer so @@ -150,11 +152,11 @@ def _hazard_pool_init(obj, data_ref): _HZ_WORKER_DATA = obj._offloader.load_dataframe(data_ref) -def _hazard_pool_task(idx, model_pos, sample_idx): - return _one_boot_log_hr(_HZ_WORKER_OBJ, _HZ_WORKER_DATA, idx, model_pos, sample_idx) +def _hazard_pool_task(val, model_pos, sample_idx): + return _one_boot_log_hr(_HZ_WORKER_OBJ, _HZ_WORKER_DATA, val, model_pos, sample_idx) -def _parallel_boot_log_hrs(self, data, idx, boot_sample_idx): +def _parallel_boot_log_hrs(self, data, val, boot_sample_idx): """Run the bootstrap hazard simulations over a process pool. The analysis frame is handed to each worker process once (via the offloader @@ -162,7 +164,7 @@ def _parallel_boot_log_hrs(self, data, idx, boot_sample_idx): models. Results are gathered in submission order, matching the serial loop; NaN/None replicates are dropped the same way. """ - data_ref = self._offloader.save_dataframe(data, f"_haz_DT_{idx}") + data_ref = self._offloader.save_dataframe(data, f"_haz_DT_{val}") # Slim copy for the pool: drop the large frames workers reload from data_ref; # keep the fitted models, bootstrap samples, and config. _GlumFit and @@ -179,7 +181,7 @@ def _parallel_boot_log_hrs(self, data, idx, boot_sample_idx): initargs=(slim, data_ref), ) as executor: futures = [ - executor.submit(_hazard_pool_task, idx, model_pos, sample_idx) + executor.submit(_hazard_pool_task, val, model_pos, sample_idx) for model_pos, sample_idx in enumerate(boot_sample_idx) ] for future in futures: @@ -190,7 +192,7 @@ def _parallel_boot_log_hrs(self, data, idx, boot_sample_idx): return boot_log_hrs -def _hazard_handler(self, data, idx, boot_idx, rng): +def _hazard_handler(self, data, val, boot_idx, rng): exclude_cols = [ "followup", f"followup{self.indicator_squared}", @@ -215,8 +217,10 @@ def _hazard_handler(self, data, idx, boot_idx, rng): ) ) - if idx is not None: - model_dict = self.outcome_model[boot_idx][idx] + if val is not None: + model_dict = self.outcome_model[boot_idx].get(val) + if model_dict is None: + return None else: model_dict = self.outcome_model[boot_idx] @@ -226,9 +230,13 @@ def _hazard_handler(self, data, idx, boot_idx, rng): ce_model = self._offloader.load_model(model_dict["compevent"]) all_treatments = [] - for val in self.treatment_level: + for treatment_val in self.treatment_level: tmp = trials.with_columns( - [pl.lit(val).alias(f"{self.treatment_col}{self.indicator_baseline}")] + [ + pl.lit(treatment_val).alias( + f"{self.treatment_col}{self.indicator_baseline}" + ) + ] ) tmp_pd = _cast_categories(self, tmp.to_pandas()) diff --git a/pySEQTarget/analysis/_outcome_fit.py b/pySEQTarget/analysis/_outcome_fit.py index 99cf945..ec64d4b 100644 --- a/pySEQTarget/analysis/_outcome_fit.py +++ b/pySEQTarget/analysis/_outcome_fit.py @@ -1,26 +1,34 @@ import re +import numpy as np # noqa: F401 - formulas are evaluated here, e.g. np.log(x) import pandas as pd import polars as pl import statsmodels.api as sm import statsmodels.formula.api as smf from pandas.api.types import is_numeric_dtype +from ..helpers._check_design import _check_design +from ..helpers._fix_categories import _cast_fixed_categories from ..helpers._spline import _compute_spline_knots, _cr_term +from ..weighting._weight_stats import _weight_bounds def _apply_spline_formula(formula, indicator_squared, spline_knots): inner_knots, lower, upper = spline_knots spline = _cr_term("followup", inner_knots, lower, upper) - formula = re.sub(r"(\w+)\s*\*\s*followup\b", rf"\1*{spline}", formula) - formula = re.sub(r"\bfollowup\s*\*\s*(\w+)", rf"{spline}*\1", formula) + # Mark interactions first, so dropping the bare followup terms below + # cannot reach into the spline, which contains "followup" itself + mark = "_SPLINE_" + formula = re.sub(r"(\w+)\s*\*\s*followup\b", rf"\1*{mark}", formula) + formula = re.sub(r"\bfollowup\s*\*\s*(\w+)", rf"{mark}*\1", formula) formula = re.sub(rf"\bfollowup{re.escape(indicator_squared)}\b", "", formula) formula = re.sub(r"\bfollowup\b", "", formula) formula = re.sub(r"\s+", " ", formula) formula = re.sub(r"\+\s*\+", "+", formula) formula = re.sub(r"^\s*\+\s*|\s*\+\s*$", "", formula).strip() + formula = formula.replace(mark, spline) if formula: return f"{formula} + {spline}" @@ -54,10 +62,7 @@ def _cast_categories(self, df_pd): if squared_col in df_pd.columns: df_pd[squared_col] = df_pd[squared_col].astype("category") - if self.fixed_cols: - for col in self.fixed_cols: - if col in df_pd.columns: - df_pd[col] = df_pd[col].astype("category") + _cast_fixed_categories(self.fixed_cols, df_pd) # Stable factor encoding for categorical time-varying covariates: fix the # level set from the full expanded data (captured on the non-bootstrap @@ -88,10 +93,9 @@ def _outcome_fit( start_params=None, ): if weighted: + lower, upper = _weight_bounds(self) df = df.with_columns( - pl.col(weight_col).clip( - lower_bound=self.weight_min, upper_bound=self.weight_max - ) + pl.col(weight_col).clip(lower_bound=lower, upper_bound=upper) ) if self.method == "censoring": @@ -108,6 +112,8 @@ def _outcome_fit( formula, self.indicator_squared, self._followup_spline_knots ) + if getattr(self, "_current_boot_idx", None) is None: + _check_design(formula, df_pd) full_formula = f"{outcome} ~ {formula}" var_weights = df_pd[weight_col] if weighted else None diff --git a/pySEQTarget/analysis/_subgroup_fit.py b/pySEQTarget/analysis/_subgroup_fit.py index e78e19d..9defdc1 100644 --- a/pySEQTarget/analysis/_subgroup_fit.py +++ b/pySEQTarget/analysis/_subgroup_fit.py @@ -1,13 +1,29 @@ +import warnings + import polars as pl from ._outcome_fit import _outcome_fit +def _warn_missing_subgroup(self, val): + """Warn when bootstrap replicates have no rows in subgroup ``val``.""" + n = sum(val not in models for models in self.outcome_model[1:]) + if n: + warnings.warn( + f"{n} bootstrap replicate(s) had no rows in {self.subgroup_colname} = " + f"{val}; they are left out of its confidence intervals.", + UserWarning, + stacklevel=3, + ) + + def _subgroup_fit(self, start_params=None): + # A resample can miss a subgroup, so key the models by subgroup value subgroups = sorted(self.DT[self.subgroup_colname].unique().to_list()) - self._unique_subgroups = subgroups + if getattr(self, "_current_boot_idx", None) is None: + self._unique_subgroups = subgroups - models_list = [] + models_by_group = {} for val in subgroups: subDT = self.DT.filter(pl.col(self.subgroup_colname) == val) sg_start = (start_params or {}).get(val, {}) or {} @@ -34,5 +50,5 @@ def _subgroup_fit(self, start_params=None): "weight", start_params=sg_start.get("compevent"), ) - models_list.append(models) - return models_list + models_by_group[val] = models + return models_by_group diff --git a/pySEQTarget/analysis/_survival_pred.py b/pySEQTarget/analysis/_survival_pred.py index 81cefa8..2d41db1 100644 --- a/pySEQTarget/analysis/_survival_pred.py +++ b/pySEQTarget/analysis/_survival_pred.py @@ -6,6 +6,7 @@ _model_design_info) from ..helpers._predict_model import _safe_predict from ._outcome_fit import _cast_categories +from ._subgroup_fit import _warn_missing_subgroup def _store_boot_risks(obj, treatment_val, TxDT, boot_cols, is_survival=False): @@ -52,22 +53,47 @@ class data: return np.asarray(dmatrix(design_info, fixed)) -def _cached_predict(model, X_cached, ref_column_names, data): +def _same_state(a, b): + """Deep equality of memorised transform state (dicts, arrays, objects).""" + if isinstance(a, dict): + return ( + isinstance(b, dict) + and a.keys() == b.keys() + and all(_same_state(a[k], b[k]) for k in a) + ) + if isinstance(a, np.ndarray) or isinstance(b, np.ndarray): + return isinstance(a, np.ndarray) and np.array_equal(a, b) + if hasattr(a, "__dict__"): + return type(a) is type(b) and _same_state(vars(a), vars(b)) + return a == b + + +def _transform_state(design_info): + """What stateful transforms such as center() or cr() memorised at fit.""" + return { + f.name(): fi.state.get("transforms", {}) + for f, fi in design_info.factor_infos.items() + if isinstance(fi.state, dict) + } + + +def _cached_predict(model, X_cached, ref_dinfo, data): """ - Predict using a pre-built design matrix when the model's design_info - column structure matches the reference, falling back to patsy via - _safe_predict on mismatch (e.g. a bootstrap resample that dropped a - categorical level). + Predict using a pre-built design matrix when the model's design matches the + reference — same columns and same memorised transform state — falling back + to patsy via _safe_predict otherwise (e.g. a resample that dropped a level). """ dinfo = _model_design_info(model) - if list(dinfo.column_names) == ref_column_names: + if list(dinfo.column_names) == list(ref_dinfo.column_names) and _same_state( + _transform_state(dinfo), _transform_state(ref_dinfo) + ): probs = np.asarray(model.predict(X_cached, transform=False)) if not np.any(np.isnan(probs)): return np.clip(probs, 0, 1) return _safe_predict(model, data) -def _get_outcome_predictions(self, TxDT, idx=None): +def _get_outcome_predictions(self, TxDT, val=None): data = _cast_categories(self, TxDT.to_pandas()) predictions = {"outcome": []} if self.compevent_colname is not None: @@ -77,58 +103,57 @@ def _get_outcome_predictions(self, TxDT, idx=None): # Each bootstrap model that shares the same column structure can then # reuse it, skipping patsy entirely on the predict path. main = self.outcome_model[0] - main_dict = main[idx] if idx is not None else main + main_dict = main[val] if val is not None else main main_outcome = self._offloader.load_model(main_dict["outcome"]) outcome_dinfo = _model_design_info(main_outcome) X_outcome = _build_design_matrix(outcome_dinfo, data) - outcome_cols = list(outcome_dinfo.column_names) - X_compevent = compevent_cols = None + X_compevent = compevent_dinfo = None if self.compevent_colname is not None: main_compevent = self._offloader.load_model(main_dict["compevent"]) compevent_dinfo = _model_design_info(main_compevent) X_compevent = _build_design_matrix(compevent_dinfo, data) - compevent_cols = list(compevent_dinfo.column_names) for boot_model in self.outcome_model: - model_dict = boot_model[idx] if idx is not None else boot_model + if val is not None and val not in boot_model: + continue + model_dict = boot_model[val] if val is not None else boot_model outcome_model = self._offloader.load_model(model_dict["outcome"]) predictions["outcome"].append( - _cached_predict(outcome_model, X_outcome, outcome_cols, data) + _cached_predict(outcome_model, X_outcome, outcome_dinfo, data) ) if self.compevent_colname is not None: compevent_model = self._offloader.load_model(model_dict["compevent"]) predictions["compevent"].append( - _cached_predict(compevent_model, X_compevent, compevent_cols, data) + _cached_predict(compevent_model, X_compevent, compevent_dinfo, data) ) return predictions def _pred_risk(self): - has_subgroups = ( - isinstance(self.outcome_model[0], list) if self.outcome_model else False - ) + has_subgroups = self.subgroup_colname is not None and bool(self.outcome_model) self._boot_risks = {} if not has_subgroups: - return _calculate_risk(self, self.DT, idx=None, val=None) + return _calculate_risk(self, self.DT, val=None) all_risks = [] original_DT = self.DT - for i, val in enumerate(self._unique_subgroups): + for val in self._unique_subgroups: + _warn_missing_subgroup(self, val) subgroup_DT = original_DT.filter(pl.col(self.subgroup_colname) == val) - risk = _calculate_risk(self, subgroup_DT, i, val) + risk = _calculate_risk(self, subgroup_DT, val) all_risks.append(risk) self.DT = original_DT return pl.concat(all_risks) -def _calculate_risk(self, data, idx=None, val=None): +def _calculate_risk(self, data, val=None): a = 1 - self.bootstrap_CI lci = a / 2 uci = 1 - lci @@ -178,7 +203,7 @@ def _calculate_risk(self, data, idx=None, val=None): ] ) - preds = _get_outcome_predictions(self, TxDT, idx=idx) + preds = _get_outcome_predictions(self, TxDT, val=val) # Drop original data columns — only followup and TID needed from here TxDT = TxDT.select(["followup", "TID"]) diff --git a/pySEQTarget/expansion/_binder.py b/pySEQTarget/expansion/_binder.py index 3c58e40..4613c4e 100644 --- a/pySEQTarget/expansion/_binder.py +++ b/pySEQTarget/expansion/_binder.py @@ -22,22 +22,13 @@ def _binder(self, kept_cols): cols = kept_cols.union({self.eligible_col, self.outcome_col, self.treatment_col}) cols = {col for col in cols if col is not None} - regular = { - col - for col in cols - if not (self.indicator_baseline in col or self.indicator_squared in col) - and col not in excluded - } - - baseline = { - col for col in cols if self.indicator_baseline in col and col not in excluded - } - bas_kept = {col.replace(self.indicator_baseline, "") for col in baseline} - - squared = { - col for col in cols if self.indicator_squared in col and col not in excluded - } - sq_kept = {col.replace(self.indicator_squared, "") for col in squared} + # The indicators are suffixes, so e.g. "bmi_baseline" is a regular column + bas_sfx, sq_sfx = self.indicator_baseline, self.indicator_squared + baseline = {c for c in cols if c.endswith(bas_sfx) and c not in excluded} + squared = {c for c in cols if c.endswith(sq_sfx) and c not in excluded} + regular = cols - baseline - squared - excluded + bas_kept = {c[: -len(bas_sfx)] for c in baseline} + sq_kept = {c[: -len(sq_sfx)] for c in squared} kept = list(regular.union(bas_kept).union(sq_kept)) @@ -72,12 +63,11 @@ def _binder(self, kept_cols): if squared: squares = [] - for sq in squared: - col = sq.replace(self.indicator_squared, "") - squares.append((pl.col(col) ** 2).alias(f"{col}{self.indicator_squared}")) + for col in sq_kept: + squares.append((pl.col(col) ** 2).alias(f"{col}{sq_sfx}")) DT = DT.with_columns(squares) - baseline_cols = {bas.replace(self.indicator_baseline, "") for bas in baseline} + baseline_cols = set(bas_kept) needed = {self.eligible_col, self.treatment_col} baseline_cols.update({c for c in needed}) diff --git a/pySEQTarget/expansion/_dynamic.py b/pySEQTarget/expansion/_dynamic.py index 52060a8..ba62d37 100644 --- a/pySEQTarget/expansion/_dynamic.py +++ b/pySEQTarget/expansion/_dynamic.py @@ -33,8 +33,9 @@ def _dynamic(self): for i, val in enumerate(self.treatment_level): colname = self.excused_colnames[i] if colname is not None: + # Column i excuses a switch away from level i, as in the weight models conditions.append( - (pl.col(colname) == 1) & (pl.col(self.treatment_col) == val) + (pl.col(colname) == 1) & (pl.col("tx_lag") == val) ) if conditions: diff --git a/pySEQTarget/helpers/_check_design.py b/pySEQTarget/helpers/_check_design.py new file mode 100644 index 0000000..af6f8ca --- /dev/null +++ b/pySEQTarget/helpers/_check_design.py @@ -0,0 +1,30 @@ +import numpy as np +import patsy + + +def _undefined_terms(formula, data): + """Design-matrix columns with missing or infinite values; ``formula`` may be a DesignInfo.""" + if isinstance(formula, str): + formula = formula.split("~", 1)[-1] + if formula.strip() in ("", "1"): + return [] + # Invalid values are what we are looking for, so don't warn about them + with np.errstate(divide="ignore", invalid="ignore", over="ignore"): + X = patsy.dmatrix( + formula, + data, + NA_action=patsy.NAAction(NA_types=[]), + return_type="dataframe", + ) + return [col for col in X.columns if not np.isfinite(X[col]).all()] + + +def _check_design(formula, data): + """Raise, naming the terms, if a model's design matrix has missing or infinite values.""" + bad = _undefined_terms(formula, data) + if bad: + raise ValueError( + f"Model term(s) {', '.join(bad)} produced missing or infinite values, " + "e.g. log() or sqrt() of a non-positive value; transform the variable " + "so the term is defined for every row." + ) diff --git a/pySEQTarget/helpers/_col_string.py b/pySEQTarget/helpers/_col_string.py index 1f6c07c..5914033 100644 --- a/pySEQTarget/helpers/_col_string.py +++ b/pySEQTarget/helpers/_col_string.py @@ -9,11 +9,14 @@ def _term_columns(expression): except SyntaxError: return set(expression.replace("+", " ").replace("*", " ").split()) - called = { - node.func - for node in ast.walk(tree) - if isinstance(node, ast.Call) and isinstance(node.func, ast.Name) - } + called = set() + for node in ast.walk(tree): + if isinstance(node, ast.Call): + # The function, or the module it comes from, as in np.log(x) + func = node.func + while isinstance(func, ast.Attribute): + func = func.value + called.add(func) return { node.id for node in ast.walk(tree) diff --git a/pySEQTarget/helpers/_fix_categories.py b/pySEQTarget/helpers/_fix_categories.py index ea4fe2d..dc1e590 100644 --- a/pySEQTarget/helpers/_fix_categories.py +++ b/pySEQTarget/helpers/_fix_categories.py @@ -1,4 +1,13 @@ import pandas as pd +from pandas.api.types import is_numeric_dtype + + +def _cast_fixed_categories(fixed_cols, df): + """Cast non-numeric fixed covariates to category; numeric ones stay continuous.""" + for col in fixed_cols: + if col in df.columns and not is_numeric_dtype(df[col]): + df[col] = df[col].astype("category") + return df def _model_design_info(model): diff --git a/pySEQTarget/helpers/_output_files.py b/pySEQTarget/helpers/_output_files.py index 8c85164..c4fbd27 100644 --- a/pySEQTarget/helpers/_output_files.py +++ b/pySEQTarget/helpers/_output_files.py @@ -6,7 +6,7 @@ def _build_md(self, img_path: str = None) -> str: Builds markdown content for SEQuential analysis results. :param self: SEQoutput instance - :param img_path: Path to saved KM graph image (if any) + :param img_path: Path to the saved survival, risk or incidence plot (if any) :return: Markdown string """ @@ -16,14 +16,22 @@ def _model_section(title, kind): # SEQoutput.summary handles absent model lists (e.g. no numerator # models under weighted ITT) and loads offloaded path refs. summaries = self.summary(kind) - if not summaries: - return - lines.append(f"### {title}") - lines.append("") - lines.append("```") - lines.append(str(summaries[0])) - lines.append("```") - lines.append("") + if isinstance(summaries, dict): + groups = summaries.items() + else: + groups = [(None, summaries)] + for val, group in groups: + if not group: + continue + label = title + if val is not None: + label = f"{title} ({self.options.subgroup_colname} = {val})" + lines.append(f"### {label}") + lines.append("") + lines.append("```") + lines.append(str(group[0])) + lines.append("```") + lines.append("") lines.append(f"# SEQuential Analysis: {datetime.date.today()}: {self.method}") lines.append("") @@ -83,9 +91,11 @@ def _model_section(title, kind): lines.append("") if self.km_graph is not None and img_path is not None: - lines.append("#### Survival Curves") + # Model-based curves of the plotted quantity, as titled on the plot + label = f"Cumulative {self.options.plot_type.title()}" + lines.append(f"#### {label}") lines.append("") - lines.append(f"![Kaplan-Meier Survival Curves]({img_path})") + lines.append(f"![{label}]({img_path})") lines.append("") if self.diagnostic_tables: diff --git a/pySEQTarget/helpers/_predict_model.py b/pySEQTarget/helpers/_predict_model.py index 13c2655..29dd962 100644 --- a/pySEQTarget/helpers/_predict_model.py +++ b/pySEQTarget/helpers/_predict_model.py @@ -1,6 +1,16 @@ import numpy as np -from ._fix_categories import _fix_categories_for_predict +from ._check_design import _check_design +from ._fix_categories import (_cast_fixed_categories, + _fix_categories_for_predict, _model_design_info) + + +def _check_predictions(model, data, probs): + """If predictions are NaN because a model term is undefined, say which.""" + if np.any(np.isnan(probs)): + design_info = _model_design_info(model) + if design_info is not None: + _check_design(design_info, data) def _safe_predict(model, data, clip_probs=True): @@ -29,6 +39,7 @@ def _safe_predict(model, data, clip_probs=True): if clip_probs: probs = np.array(probs) + _check_predictions(model, data, probs) if np.any(np.isnan(probs)): raise ValueError( "NaN values in predicted probabilities. This typically indicates " @@ -42,29 +53,26 @@ def _safe_predict(model, data, clip_probs=True): def _prep_predict_frame(self, newdata): - """Convert a polars frame to pandas with fixed_cols cast to category. + """Convert a polars frame to pandas with non-numeric fixed_cols cast to category. Split out from _predict_model so callers predicting with several models on the same rows (e.g. numerator + denominator in _weight_predict) can pay the conversion once and share the frame. """ - newdata = newdata.to_pandas() - for col in self.fixed_cols: - if col in newdata.columns: - newdata[col] = newdata[col].astype("category") - return newdata + return _cast_fixed_categories(self.fixed_cols, newdata.to_pandas()) def _predict_model_pd(model, newdata): """Predict on an already-prepared pandas frame, with category fix retry.""" try: - return np.array(model.predict(newdata)) + probs = np.array(model.predict(newdata)) except Exception as e: - if "mismatching levels" in str(e): - newdata = _fix_categories_for_predict(model, newdata) - return np.array(model.predict(newdata)) - else: + if "mismatching levels" not in str(e): raise + newdata = _fix_categories_for_predict(model, newdata) + probs = np.array(model.predict(newdata)) + _check_predictions(model, newdata, probs) + return probs def _predict_model(self, model, newdata): diff --git a/pySEQTarget/helpers/_spline.py b/pySEQTarget/helpers/_spline.py index 2333e2e..d0a48c0 100644 --- a/pySEQTarget/helpers/_spline.py +++ b/pySEQTarget/helpers/_spline.py @@ -1,6 +1,7 @@ import re import numpy as np +import patsy # The only cr() form whose basis depends on the rows, so the only one to bake. _CR_DF = re.compile( @@ -27,8 +28,10 @@ def _compute_spline_knots(arr, df=3): def _cr_term(var, inner_knots, lower, upper, constraints=None): - """A ``cr()`` term with its knots written out.""" - constraint = f', constraints="{constraints}"' if constraints else "" + """A ``cr()`` term with its knots, and any linear constraint, written out.""" + constraint = ( + f", constraints=[{list(constraints)}]" if constraints is not None else "" + ) return ( f"cr({var}, knots={list(inner_knots)}, " f"lower_bound={lower}, upper_bound={upper}{constraint})" @@ -79,6 +82,13 @@ def _bake_spline_knots(formula, data): f"{len(inner_knots)} interior knot(s) for cr(df={df}); use a " "smaller df, or model this term without a spline." ) + if constraints: + # patsy's centring constraint (the basis column means) also depends + # on the rows it sees, so fix it from the full data too + basis = patsy.dmatrix( + _cr_term("x", inner_knots, lower, upper) + "-1", {"x": arr} + ) + constraints = np.asarray(basis).mean(axis=0).tolist() formula = formula.replace( token, _cr_term(var, inner_knots, lower, upper, constraints) ) diff --git a/pySEQTarget/initialization/_outcome.py b/pySEQTarget/initialization/_outcome.py index cd693f5..de9b702 100644 --- a/pySEQTarget/initialization/_outcome.py +++ b/pySEQTarget/initialization/_outcome.py @@ -14,7 +14,9 @@ def _outcome(self) -> str: if self.time_varying_cols else None ) - fixed = "+".join(self.fixed_cols) if self.fixed_cols else None + # The subgroup column is constant within each subgroup's outcome model + fixed_cols = [c for c in self.fixed_cols if c != self.subgroup_colname] + fixed = "+".join(fixed_cols) if fixed_cols else None trial = ( "+".join(["trial", f"trial{self.indicator_squared}"]) if self.trial_include diff --git a/pySEQTarget/weighting/_weight_fit.py b/pySEQTarget/weighting/_weight_fit.py index 851bff8..a198f99 100644 --- a/pySEQTarget/weighting/_weight_fit.py +++ b/pySEQTarget/weighting/_weight_fit.py @@ -1,7 +1,9 @@ +import numpy as np # noqa: F401 - formulas are evaluated here, e.g. np.log(x) import statsmodels.api as sm import statsmodels.formula.api as smf from ..error._check_separation import _check_separation +from ..helpers._check_design import _check_design from ..helpers._glum_fit import _fit_glum @@ -59,6 +61,8 @@ def _fit_pair( setattr(self, out, None) continue formula = f"{outcome}~{rhs}" + if getattr(self, "_current_boot_idx", None) is None: + _check_design(formula, WDT) if getattr(self, "glm_package", "statsmodels") == "glum": fitted = _fit_glum(formula, WDT) else: @@ -114,6 +118,8 @@ def _fit_numerator(self, WDT): if len(DT_subset[predictor].unique()) < 2: fits.append(None) continue + if getattr(self, "_current_boot_idx", None) is None: + _check_design(formula, DT_subset) # Use logit for binary 0/1 censoring, mnlogit otherwise if is_binary and getattr(self, "glm_package", "statsmodels") == "glum": model_fit = _fit_glum(formula, DT_subset) @@ -125,6 +131,8 @@ def _fit_numerator(self, WDT): model_fit = smf.mnlogit(formula, DT_subset).fit( disp=0, method=self.weight_fit_method ) + # Column order of mnlogit predictions, to look classes up by value + model_fit.seq_classes = sorted(DT_subset[predictor].unique().tolist()) _check_separation(model_fit, label=f"numerator (level {level})") fits.append(model_fit) @@ -155,6 +163,8 @@ def _fit_denominator(self, WDT): if len(DT_subset[predictor].unique()) < 2: fits.append(None) continue + if getattr(self, "_current_boot_idx", None) is None: + _check_design(formula, DT_subset) # Use logit for binary 0/1 censoring, mnlogit otherwise if is_binary and getattr(self, "glm_package", "statsmodels") == "glum": model_fit = _fit_glum(formula, DT_subset) @@ -166,6 +176,8 @@ def _fit_denominator(self, WDT): model_fit = smf.mnlogit(formula, DT_subset).fit( disp=0, method=self.weight_fit_method ) + # Column order of mnlogit predictions, to look classes up by value + model_fit.seq_classes = sorted(DT_subset[predictor].unique().tolist()) _check_separation(model_fit, label=f"denominator (level {level})") fits.append(model_fit) diff --git a/pySEQTarget/weighting/_weight_pred.py b/pySEQTarget/weighting/_weight_pred.py index 56341b5..c6310f4 100644 --- a/pySEQTarget/weighting/_weight_pred.py +++ b/pySEQTarget/weighting/_weight_pred.py @@ -5,21 +5,13 @@ from ..helpers._predict_model import _predict_model_pd, _prep_predict_frame -def _extract_class_probability(p, level_idx, is_binary): - """ - Extract the probability for a specific class from model predictions. - Handles both binary (logit) and multinomial (mnlogit) output formats. - """ +def _extract_class_probability(p, model, level, is_binary): + """P(treatment == level): logit predicts P(A = 1), mnlogit one column per class.""" if is_binary: - # logit returns P(Y=1) directly as 1D array - # For level 0: want P(stay at 0) = 1 - P(Y=1) - # For level 1: want P(switch to 1) = P(Y=1) - return p if level_idx == 1 else (1 - p) - else: - # mnlogit returns [P(Y=0), P(Y=1), ...] as 2D array - if p.ndim == 1: - p = p.reshape(-1, 1) - return p[:, level_idx] + return p if level == 1 else 1 - p + if level not in model.seq_classes: + return np.zeros(len(p)) + return p.reshape(len(p), -1)[:, model.seq_classes.index(level)] def _weight_predict(self, WDT): @@ -69,14 +61,16 @@ def _weight_predict(self, WDT): if denom_model is not None: p = _predict_model_pd(denom_model, subset_pd) - p_class = _extract_class_probability(p, i, is_binary) + p_class = _extract_class_probability( + p, denom_model, level, is_binary + ) pred_denom[lag_mask] = np.where( switched_treatment[lag_mask], 1.0 - p_class, p_class ) if num_model is not None: p = _predict_model_pd(num_model, subset_pd) - p_class = _extract_class_probability(p, i, is_binary) + p_class = _extract_class_probability(p, num_model, level, is_binary) pred_num[lag_mask] = np.where( switched_treatment[lag_mask], 1.0 - p_class, p_class ) @@ -91,84 +85,65 @@ def _weight_predict(self, WDT): else: for i, level in enumerate(self.treatment_level): col = self.excused_colnames[i] - - if col is not None: - denom_model = self._offloader.load_model(self.denominator_model[i]) - denom_mask = ((WDT["tx_lag"] == level) & (WDT[col] != 1)).to_numpy() - - if denom_model is not None and denom_mask.sum() > 0: - pred_denom = np.ones(WDT.height) - subset = WDT.filter(pl.Series(denom_mask)) - p = _predict_model(self, denom_model, subset) - - if p.ndim == 1: - prob_switch = p - else: - prob_switch = p[:, 1] if p.shape[1] > 1 else p.flatten() - - pred_denom[denom_mask] = prob_switch - - WDT = WDT.with_columns( - pl.when(pl.Series(denom_mask)) - .then(pl.Series(pred_denom)) - .otherwise(pl.col("denominator")) - .alias("denominator") - ) - - if i == 0: - flip_mask = ( - (WDT["tx_lag"] == level) - & (WDT[col] == 0) - & (WDT[self.treatment_col] == level) - ).to_numpy() - else: - flip_mask = ( - (WDT["tx_lag"] == level) - & (WDT[col] == 0) - & (WDT[self.treatment_col] != level) - ).to_numpy() - - WDT = WDT.with_columns( - pl.when(pl.Series(flip_mask)) - .then(1.0 - pl.col("denominator")) - .otherwise(pl.col("denominator")) - .alias("denominator") + if col is None: + continue + denom_model = self._offloader.load_model(self.denominator_model[i]) + denom_mask = ((WDT["tx_lag"] == level) & (WDT[col] != 1)).to_numpy() + if denom_model is None or denom_mask.sum() == 0: + continue + subset = WDT.filter(pl.Series(denom_mask)) + p = _predict_model(self, denom_model, subset) + if self.weight_preexpansion: + # Treatment model: probability of the treatment received + p_level = _extract_class_probability( + p, denom_model, level, is_binary ) + stay = (subset[self.treatment_col] == level).to_numpy() + prob = np.where(stay, p_level, 1.0 - p_level) + else: + # Censoring model: probability of remaining uncensored + p_cens = _extract_class_probability(p, denom_model, 1, is_binary) + prob = np.where( + subset["switch"].to_numpy() == 0, 1.0 - p_cens, p_cens + ) + pred_denom = np.ones(WDT.height) + pred_denom[denom_mask] = prob + WDT = WDT.with_columns( + pl.when(pl.Series(denom_mask)) + .then(pl.Series(pred_denom)) + .otherwise(pl.col("denominator")) + .alias("denominator") + ) if self.weight_preexpansion: WDT = WDT.with_columns(pl.lit(1.0).alias("numerator")) else: for i, level in enumerate(self.treatment_level): col = self.excused_colnames[i] + if col is None: + continue + num_model = self._offloader.load_model(self.numerator_model[i]) + num_mask = ( + (WDT[self.treatment_col] == level) & (WDT[col] == 0) + ).to_numpy() + if num_model is None or num_mask.sum() == 0: + continue + subset = WDT.filter(pl.Series(num_mask)) + p = _predict_model(self, num_model, subset) + pred_num = np.ones(WDT.height) + pred_num[num_mask] = _extract_class_probability( + p, num_model, 1, is_binary + ) + WDT = WDT.with_columns( + pl.when(pl.Series(num_mask)) + .then(pl.Series(pred_num)) + .otherwise(pl.col("numerator")) + .alias("numerator") + ) - if col is not None: - num_model = self._offloader.load_model(self.numerator_model[i]) - num_mask = ( - (WDT[self.treatment_col] == level) & (WDT[col] == 0) - ).to_numpy() - - if num_model is not None and num_mask.sum() > 0: - pred_num = np.ones(WDT.height) - subset = WDT.filter(pl.Series(num_mask)) - p = _predict_model(self, num_model, subset) - - if p.ndim == 1: - prob_switch = p - else: - prob_switch = p[:, 1] if p.shape[1] > 1 else p.flatten() - - pred_num[num_mask] = prob_switch - - WDT = WDT.with_columns( - pl.when(pl.Series(num_mask)) - .then(pl.Series(pred_num)) - .otherwise(pl.col("numerator")) - .alias("numerator") - ) - - first_level = self.treatment_level[0] + # As for the denominator, the model predicts P(censored) WDT = WDT.with_columns( - pl.when(pl.col(self.treatment_col) == first_level) + pl.when(pl.col("switch") == 0) .then(1.0 - pl.col("numerator")) .otherwise(pl.col("numerator")) .alias("numerator") @@ -182,8 +157,9 @@ def _weight_predict(self, WDT): cense_denom_model = self._offloader.load_model(self.cense_denominator_model) if cense_num_model is not None and cense_denom_model is not None: WDT_pd = _prep_predict_frame(self, WDT) - p_num = _predict_model_pd(cense_num_model, WDT_pd).flatten() - p_denom = _predict_model_pd(cense_denom_model, WDT_pd).flatten() + # The models predict P(lost); the weight uses P(remaining), as in R + p_num = 1.0 - _predict_model_pd(cense_num_model, WDT_pd).flatten() + p_denom = 1.0 - _predict_model_pd(cense_denom_model, WDT_pd).flatten() WDT = ( WDT.with_columns( [ diff --git a/pySEQTarget/weighting/_weight_stats.py b/pySEQTarget/weighting/_weight_stats.py index 6004884..988a054 100644 --- a/pySEQTarget/weighting/_weight_stats.py +++ b/pySEQTarget/weighting/_weight_stats.py @@ -16,8 +16,12 @@ def _weight_stats(self): ] ) - if self.weight_p99: - self.weight_min = float(stats["weight_p01"][0]) - self.weight_max = float(stats["weight_p99"][0]) - return stats + + +def _weight_bounds(self): + """Truncation bounds for the current fit; with weight_p99, its 1st/99th percentiles.""" + stats = getattr(self, "weight_stats", None) + if self.weight_p99 and stats is not None: + return float(stats["weight_p01"][0]), float(stats["weight_p99"][0]) + return self.weight_min, self.weight_max diff --git a/pyproject.toml b/pyproject.toml index 592513b..2cea093 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "pySEQTarget" -version = "0.14.2" +version = "0.14.3" description = "Sequential Target Trial Emulation" readme = "README.md" license = {text = "MIT"} @@ -19,6 +19,7 @@ classifiers = [ "Programming Language :: Python :: 3.12", "Programming Language :: Python :: 3.13", "Programming Language :: Python :: 3.14", + "Programming Language :: Python :: 3.15", ] authors = [ @@ -53,8 +54,8 @@ output = [ "tabulate" ] gpu = [ - "jax[cuda12]; sys_platform == 'linux'", - "jax; sys_platform != 'linux'", + "jax[cuda12]; sys_platform == 'linux' and python_version < '3.15'", + "jax; sys_platform != 'linux' and python_version < '3.15'", ] dev = [ "pySEQTarget[output,gpu]", @@ -76,7 +77,7 @@ Repository = "https://github.com/CausalInference/pySEQTarget" dev = [ "black", "isort", - "jax", + "jax; python_version < '3.15'", "pytest", "myst-parser", "piccolo_theme", diff --git a/tests/test_cense_weights.py b/tests/test_cense_weights.py new file mode 100644 index 0000000..1e21128 --- /dev/null +++ b/tests/test_cense_weights.py @@ -0,0 +1,34 @@ +"""Censoring weights use the probability of remaining under follow-up.""" + +import numpy as np +import pytest + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data +from pySEQTarget.helpers._predict_model import _predict_model + + +@pytest.mark.filterwarnings("ignore") +def test_ltfu_weight_is_the_ratio_of_probabilities_of_remaining(): + s = SEQuential( + load_data("SEQdata_LTFU"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="ITT", + parameters=SEQopts( + weighted=True, weight_preexpansion=False, cense_colname="LTFU" + ), + ) + s.expand() + s.fit() + # Under ITT the treatment weights are 1, so each weight is the LTFU factor + rows = s.DT.filter(s.DT["followup"] > 0).head(500) + lost_num = _predict_model(s, s.cense_numerator_model, rows) + lost_den = _predict_model(s, s.cense_denominator_model, rows) + expected = (1 - lost_num) / (1 - lost_den) + assert np.allclose(rows["weight"].to_numpy(), expected) diff --git a/tests/test_coefficients.py b/tests/test_coefficients.py index bef2738..2da644f 100644 --- a/tests/test_coefficients.py +++ b/tests/test_coefficients.py @@ -1,3 +1,4 @@ +import pytest from pytest import approx from pySEQTarget import SEQopts, SEQuential @@ -25,11 +26,11 @@ def test_ITT_coefs(): expected = [ -6.828506035553407, 0.18935003090041902, - 0.12717241010542563, 0.033715156987629266, -0.00014691202235029346, 0.044566165558944326, 0.0005787770439053261, + 0.12717241010542563, 0.0032906669395295026, -0.01339242049205771, 0.20072409918428052, @@ -56,14 +57,14 @@ def test_PreE_dose_response_coefs(): s.fit() matrix = s.outcome_model[0]["outcome"].summary2().tables[1]["Coef."].to_list() expected = [ - -4.842735939069144, - 0.14286755770151904, - 0.055221018477671927, - -0.000581657931537684, - -0.008484541900408258, - 0.00021073328759912806, - 0.010537967151467553, - 0.0007772316818101141, + -4.838361606118745, + 0.0555200193349395, + -0.000588987903925482, + -0.008977364567147181, + 0.00022025518488454974, + 0.01046828606393553, + 0.0007774475011663997, + 0.13972199132132426, ] assert matrix == approx(expected, rel=1e-3) @@ -89,13 +90,13 @@ def test_PostE_dose_response_coefs(): matrix = s.outcome_model[0]["outcome"].summary2().tables[1]["Coef."].to_list() expected = [ -6.265901713761531, - 0.14065954021957594, 0.048626017624679704, -0.0004688287307505834, -0.003975906839775267, 0.00016676441745740924, 0.03866279977096397, 0.0005928449623613982, + 0.14065954021957594, 0.0030001459817949844, -0.02106338184559446, 0.14867250693140854, @@ -122,13 +123,13 @@ def test_PreE_censoring_coefs(): s.fit() matrix = s.outcome_model[0]["outcome"].summary2().tables[1]["Coef."].to_list() expected = [ - -4.872373936951975, - 0.48389186624295133, - 0.0477349453301334, - 0.029127276869076173, - 4.784054961154824e-05, - -0.013614654772205668, - 0.0011281734101133744, + -4.866811032893474, + 0.48052384773425627, + 0.029219605971331802, + 4.518185255729768e-05, + -0.013811946311246277, + 0.0011309893093644648, + 0.04763133686352373, ] assert matrix == approx(expected, rel=1e-3) @@ -155,11 +156,11 @@ def test_PostE_censoring_coefs(): expected = [ -9.172266785106519, 0.4707554720116625, - 0.08162617232478116, 0.029021087196430605, 7.8937226861939e-05, 0.06700192286925702, 0.0005834323664644191, + 0.08162617232478116, 0.004870212765388434, 0.013503198983327514, 0.4466573801510379, @@ -191,12 +192,12 @@ def test_PreE_censoring_excused_coefs(): s.fit() matrix = s.outcome_model[0]["outcome"].summary2().tables[1]["Coef."].to_list() expected = [ - -5.028261715903588, - 0.09661040854758277, - -0.029861423750765226, - 0.0014186936955145387, - 0.08365564531281737, - -0.0006220464783614585, + -4.804650815594313, + 0.15964610209528257, + 0.029097252181086327, + -4.654535663714614e-05, + 0.007817954062830172, + 0.0008537393051722129, ] assert matrix == approx(expected, rel=1e-3) @@ -219,23 +220,22 @@ def test_PostE_censoring_excused_coefs(): weight_preexpansion=False, excused=True, excused_colnames=["excusedZero", "excusedOne"], - weight_max=1, ), ) s.expand() s.fit() matrix = s.outcome_model[0]["outcome"].summary2().tables[1]["Coef."].to_list() expected = [ - -7.722441228318476, - 0.25040421685828396, - 0.08370244078073162, - 0.03644249151697272, - -0.00019169394285363785, - 0.053677366381589396, - 0.0005643189202781975, - 0.005250478928581509, - 0.0014679503081325516, - 0.3008769969502361, + -8.311960526025656, + 0.1820521596172769, + 0.03422386974709779, + -0.00013723383158569053, + 0.06646241563319619, + 0.0004999277868749847, + 0.15032828057125242, + 0.0012618850925656703, + 0.0028819214525432277, + 0.3590674018732999, ] assert matrix == approx(expected, rel=1e-3) @@ -261,16 +261,16 @@ def test_PreE_LTFU_ITT(): s.fit() matrix = s.outcome_model[0]["outcome"].summary2().tables[1]["Coef."].to_list() expected = [ - -21.636346991788276, - 0.06813705852786496, - -0.1939555961858531, - 0.02874152772603635, - -0.0005734047013500563, - 0.2854740212699898, - -0.0013729662310668182, - 0.006501915963316852, - -0.4467079969655381, - 1.3870473474960576, + -21.640562053561073, + 0.06852515972631311, + 0.028751019561197418, + -0.0005762181500554517, + 0.2855438174370584, + -0.0013730471122256052, + -0.19004771044056396, + 0.006589451846069475, + -0.44899973509857344, + 1.3875130951733505, ] assert matrix == approx(expected, rel=1e-3) @@ -296,16 +296,16 @@ def test_PostE_LTFU_ITT(): s.fit() matrix = s.outcome_model[0]["outcome"].summary2().tables[1]["Coef."].to_list() expected = [ - -21.847198431385877, - 0.07786703138967718, - -0.15461370944416225, - 0.030140057462437704, - -0.0006287338029348562, - 0.287393206037481, - -0.0013719595115633126, - 0.007295485861066434, - -0.42797049565882755, - 1.4082102322835948, + -21.638981392234413, + 0.06850547874247659, + 0.02874572640721577, + -0.0005760278247524966, + 0.2855200268142661, + -0.0013729693434970476, + -0.1901847051810351, + 0.006586770631265142, + -0.4490558452282004, + 1.3873546219096944, ] assert matrix == approx(expected, rel=1e-3) @@ -331,11 +331,11 @@ def test_ITT_multinomial(): expected = [ -47.505262164163625, 1.76628017234151, - 22.79205044396338, 0.14473536056627245, -0.003725499516376173, 0.2893070991930884, -0.004266608123938117, + 22.79205044396338, 0.05574429164512122, 0.7847862691929901, 1.4703411759229423, @@ -343,6 +343,7 @@ def test_ITT_multinomial(): assert matrix == approx(expected, rel=1e-3) +@pytest.mark.filterwarnings("ignore:overflow encountered in exp:RuntimeWarning") def test_weighted_multinomial(): data = load_data("SEQdata_multitreatment") @@ -362,17 +363,21 @@ def test_weighted_multinomial(): ) s.expand() s.fit() - matrix = s.outcome_model[0]["outcome"].summary2().tables[1]["Coef."].to_list() - expected = [ - -109.99715622379995, - -12.536816769546702, - 9.22013733949143, - -0.6129380297017852, - 0.01597877250531723, - 5.743984176710672, - -0.08478678955657822, - ] - assert matrix == approx(expected, rel=1e-3) + # The outcome model has no events in arm 2, so its fit diverges (and + # overflows); pin the multinomial weights instead + stats = s.weight_stats.row(0, named=True) + expected = { + "weight_min": 0.6247271882268575, + "weight_max": 1.3120226973615323, + "weight_mean": 0.9947273760015201, + "weight_std": 0.05493659250756025, + "weight_p01": 0.8251337794242161, + "weight_p25": 0.972097606970619, + "weight_p50": 0.9968156083037403, + "weight_p75": 1.0164134507964089, + "weight_p99": 1.1626440219480307, + } + assert stats == approx(expected, rel=1e-3) def test_ITT_visit(): @@ -398,11 +403,11 @@ def test_ITT_visit(): expected = [ -21.636346991788276, 0.06813705852786496, - -0.1939555961858531, 0.02874152772603635, -0.0005734047013500563, 0.2854740212699898, -0.0013729662310668182, + -0.1939555961858531, 0.006501915963316852, -0.4467079969655381, 1.3870473474960576, diff --git a/tests/test_end_of_fup.py b/tests/test_end_of_fup.py index 10fd2d2..0c58a5c 100644 --- a/tests/test_end_of_fup.py +++ b/tests/test_end_of_fup.py @@ -435,3 +435,27 @@ def test_end_of_followup_requires_eof_mode_and_fit(): s.expand() with pytest.raises(ValueError, match="fit"): s.end_of_followup() + + +def test_trial_periods_that_cannot_reach_the_window_are_not_eligible(): + k, w = 12, 3 + s = _eof_run(end_of_fup=True, end_of_fup_time=k, end_of_fup_window=w) + starts = s.DT.group_by(["ID", "trial", "tx_init_bas"]).agg( + (pl.col("period") - pl.col("followup")).min().alias("start") + ) + reachable = ( + starts.filter(pl.col("start") + k - w <= s.DT["period"].max()) + .group_by("tx_init_bas") + .len() + .sort("tx_init_bas") + ) + eof = s.eof_data.sort("A") + assert eof["Trial-periods (Eligible)"].to_list() == reachable["len"].to_list() + assert reachable["len"].sum() < starts.height + # Analysed, censored and never measured still partition the eligible total + parts = ( + eof["Trial-periods (Analysed)"] + + eof["Trial-periods (Censored)"] + + eof["Trial-periods (No measurement)"] + ) + assert (parts == eof["Trial-periods (Eligible)"]).all() diff --git a/tests/test_excused_direction.py b/tests/test_excused_direction.py new file mode 100644 index 0000000..cfc00a5 --- /dev/null +++ b/tests/test_excused_direction.py @@ -0,0 +1,31 @@ +"""excused_colnames[i] excuses a switch away from treatment_level[i].""" + +import polars as pl + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data + + +def test_expansion_excuses_switches_away_from_each_level(): + s = SEQuential( + load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="censoring", + parameters=SEQopts( + excused=True, excused_colnames=["excusedZero", "excusedOne"] + ), + ) + s.expand() + excused = s.DT.filter(pl.col("isExcused")) + assert excused.height > 0 + # A switch away from 0 lands on 1 (excusedZero), and away from 1 on 0 (excusedOne) + away = ((pl.col("excusedZero") == 1) & (pl.col("tx_init") == 1)) | ( + (pl.col("excusedOne") == 1) & (pl.col("tx_init") == 0) + ) + assert excused.select(away.all()).item() diff --git a/tests/test_excused_weights.py b/tests/test_excused_weights.py new file mode 100644 index 0000000..99d0e1b --- /dev/null +++ b/tests/test_excused_weights.py @@ -0,0 +1,49 @@ +"""Post-expansion excused weights use P(uncensored) in every arm.""" + +import polars as pl +import pytest + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data + + +@pytest.mark.filterwarnings("ignore") +def test_uncensored_rows_get_the_probability_of_staying_uncensored(): + s = SEQuential( + load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="censoring", + parameters=SEQopts( + weighted=True, + weight_preexpansion=False, + excused=True, + excused_colnames=["excusedZero", "excusedOne"], + ), + ) + s.expand() + s.fit() + # Uncensored rows on their baseline arm, which that arm's models predict + own_col = ( + pl.when(pl.col("tx_init_bas") == 0) + .then(pl.col("excusedZero")) + .otherwise(pl.col("excusedOne")) + ) + rows = s.DT.filter( + (pl.col("switch") == 0) + & (pl.col("followup") > 0) + & (pl.col("tx_init") == pl.col("tx_init_bas")) + & (own_col == 0) + ) + medians = rows.group_by("tx_init_bas").agg( + pl.col("numerator").median(), pl.col("denominator").median() + ) + # Switching is rare, so P(uncensored) is near 1 in both arms + assert medians.height == 2 + assert (medians["numerator"] > 0.5).all() + assert (medians["denominator"] > 0.5).all() diff --git a/tests/test_fixed_covariates.py b/tests/test_fixed_covariates.py new file mode 100644 index 0000000..46a6402 --- /dev/null +++ b/tests/test_fixed_covariates.py @@ -0,0 +1,75 @@ +"""Numeric fixed covariates stay continuous; non-numeric ones are categorical.""" + +import numpy as np +import pandas as pd +import polars as pl + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data +from pySEQTarget.helpers._fix_categories import _cast_fixed_categories + + +def _with_age(data): + rng = np.random.default_rng(1) + ids = data["ID"].unique().sort() + ages = pl.DataFrame({"ID": ids, "age": np.round(rng.uniform(20, 80, len(ids)), 1)}) + return data.join(ages, on="ID") + + +def test_cast_leaves_numeric_fixed_covariates_numeric(): + df = pd.DataFrame( + { + "age": [34.5, 61.2, 47.8], + "sex": [0, 1, 1], + "region": ["north", "south", "north"], + } + ) + out = _cast_fixed_categories(["age", "sex", "region"], df) + assert out["age"].tolist() == [34.5, 61.2, 47.8] + assert out["sex"].tolist() == [0, 1, 1] + assert isinstance(out["region"].dtype, pd.CategoricalDtype) + + +def test_continuous_fixed_covariate_is_a_single_term(): + s = SEQuential( + _with_age(load_data("SEQdata")), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex", "age"], + method="censoring", + parameters=SEQopts(weighted=True, weight_preexpansion=False, km_curves=True), + ) + s.expand() + s.fit() + s.survival() + names = list(s.outcome_model[0]["outcome"].params.index) + assert [n for n in names if n.startswith("age")] == ["age"] + assert "sex" in names + assert "age" in s.denominator_model[0].params.index + + +def test_subgroup_drops_only_the_subgroup_column_from_the_outcome_model(): + s = SEQuential( + _with_age(load_data("SEQdata")), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex", "age"], + method="ITT", + parameters=SEQopts(subgroup_colname="sex", km_curves=True), + ) + terms = s.covariates.split("+") + assert "age" in terms + assert "sex" not in terms + # Still carried through the expansion to split the analysis by + s.expand() + s.fit() + s.survival() + assert set(s.km_data["sex"].unique()) == {0, 1} diff --git a/tests/test_followup_options.py b/tests/test_followup_options.py index 4cd6c65..c33ecc5 100644 --- a/tests/test_followup_options.py +++ b/tests/test_followup_options.py @@ -30,9 +30,9 @@ def test_followup_class(): 0.11375627612408938, 0.14496108664292745, 0.1798424095611678, - 0.09066206802273916, 0.015738693166264354, 0.0009258560187318309, + 0.09066206802273916, 0.011267393559242982, 0.022194521411244304, 0.19115237121222872, @@ -70,9 +70,9 @@ def test_followup_spline(): expected = [ -4.804282252748607, 0.19115933860001255, - 0.12717121164606823, 0.044310717515918724, 0.0005814999431447507, + 0.12717121164606823, 0.0032948355025455216, -0.013371824500839971, 0.19972467861548412, @@ -105,11 +105,40 @@ def test_no_followup(): expected = [ -6.062350570326165, 0.17748844870984498, - 0.11209431124681817, 0.03344595751001804, 0.0005457002039545119, + 0.11209431124681817, 0.0032236473201563585, -0.014463448024337773, 0.20398559747503964, ] assert [round(x, 3) for x in matrix] == [round(x, 3) for x in expected] + + +@pytest.mark.filterwarnings( + "ignore::statsmodels.tools.sm_exceptions.SingularMatrixWarning" +) +@pytest.mark.parametrize("method", ["ITT", "dose-response"]) +def test_followup_spline_with_km_curves_interactions(method): + # km_curves adds treatment- or dose-by-followup interactions, which must + # interact with the spline rather than lose its variable + s = SEQuential( + load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method=method, + parameters=SEQopts( + followup_spline=True, followup_include=False, km_curves=True + ), + ) + s.expand() + s.fit() + s.survival() + names = s.outcome_model[0]["outcome"].params.index + assert any("cr(followup" in n and ":" in n for n in names) + assert s.km_data.height > 0 diff --git a/tests/test_indicator_suffixes.py b/tests/test_indicator_suffixes.py new file mode 100644 index 0000000..282596a --- /dev/null +++ b/tests/test_indicator_suffixes.py @@ -0,0 +1,46 @@ +"""The _bas/_sq indicators are matched only as name suffixes.""" + +import polars as pl + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data + + +def test_indicator_text_mid_name_is_a_regular_column(): + data = load_data("SEQdata").with_columns( + (pl.col("N") + 1).alias("bmi_baseline"), (pl.col("L") + 1).alias("bmi_sqrt") + ) + s = SEQuential( + data, + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex", "bmi_baseline", "bmi_sqrt"], + method="ITT", + ) + s.expand() + assert {"bmi_baseline", "bmi_sqrt", "N_bas"} <= set(s.DT.columns) + s.fit() + names = s.outcome_model[0]["outcome"].params.index + assert {"bmi_baseline", "bmi_sqrt"} <= set(names) + + +def test_suffixed_columns_are_still_derived(): + s = SEQuential( + load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="ITT", + parameters=SEQopts(covariates="tx_init_bas+followup+N_bas+L_sq"), + ) + s.expand() + assert {"N_bas", "L_sq"} <= set(s.DT.columns) + assert (s.DT["L_sq"] == s.DT["L"] ** 2).all() diff --git a/tests/test_output_files.py b/tests/test_output_files.py index 8ff47ac..6ec167e 100644 --- a/tests/test_output_files.py +++ b/tests/test_output_files.py @@ -69,3 +69,59 @@ def test_to_md_with_offloaded_models(tmp_path): assert "Numerator Model" in content assert "Denominator Model" in content assert "Outcome Model" in content + + +def test_subgroup_models_are_collected_and_reported_by_subgroup(tmp_path): + s = SEQuential( + load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="ITT", + parameters=SEQopts(subgroup_colname="sex", bootstrap_nboot=2, seed=42), + ) + s.expand() + s.bootstrap() + s.fit() + out = s.collect() + + # Per subgroup: the main fit, then one model per replicate + assert set(out.outcome_models) == {0, 1} + assert all(len(models) == 3 for models in out.outcome_models.values()) + assert set(out.summary("outcome")) == {0, 1} + + md_file = tmp_path / "report.md" + out.to_md(str(md_file)) + content = md_file.read_text() + assert "Outcome Model (sex = 0)" in content + assert "Outcome Model (sex = 1)" in content + + +def test_report_labels_the_plot_by_what_it_shows(tmp_path): + s = SEQuential( + load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="ITT", + parameters=SEQopts(km_curves=True, plot_type="risk", seed=42), + ) + s.expand() + s.fit() + s.survival() + s.plot() + md_file = tmp_path / "report.md" + s.collect().to_md(str(md_file)) + content = md_file.read_text() + # Model-based curves, so no Kaplan-Meier label + assert "#### Cumulative Risk" in content + assert "![Cumulative Risk](report.png)" in content + assert "Kaplan-Meier" not in content diff --git a/tests/test_parallel.py b/tests/test_parallel.py index a6384f9..f0e8c74 100644 --- a/tests/test_parallel.py +++ b/tests/test_parallel.py @@ -32,11 +32,11 @@ def test_parallel_ITT(): [ -6.828506035553407, 0.18935003090041902, - 0.12717241010542563, 0.033715156987629266, -0.00014691202235029346, 0.044566165558944326, 0.0005787770439053261, + 0.12717241010542563, 0.0032906669395295026, -0.01339242049205771, 0.20072409918428052, diff --git a/tests/test_preexpansion_first_row.py b/tests/test_preexpansion_first_row.py new file mode 100644 index 0000000..36e56ec --- /dev/null +++ b/tests/test_preexpansion_first_row.py @@ -0,0 +1,30 @@ +"""Each subject's first row is left out of the pre-expansion treatment models.""" + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data + + +def _model(): + data = load_data("SEQdata") + s = SEQuential( + data, + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="censoring", + parameters=SEQopts(weighted=True), + ) + s.expand() + s.fit() + return s, data + + +def test_first_rows_are_not_fitted(): + s, data = _model() + n_ids = data["ID"].n_unique() + fitted = sum(int(m.nobs) for m in s.denominator_model) + assert fitted == data.height - n_ids diff --git a/tests/test_stateful_transforms.py b/tests/test_stateful_transforms.py new file mode 100644 index 0000000..e3589e3 --- /dev/null +++ b/tests/test_stateful_transforms.py @@ -0,0 +1,74 @@ +"""Terms that remember their fitting data, e.g. centred splines or center(x).""" + +import numpy as np +import pytest + +import pySEQTarget.analysis._survival_pred as survival_pred +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data + + +def _model(**opts): + s = SEQuential( + load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="censoring", + parameters=SEQopts(seed=42, **opts), + ) + s.expand() + if s.bootstrap_nboot: + s.bootstrap() + s.fit() + return s + + +def _weights(s): + return s.DT.sort(["ID", "trial", "followup"])["weight"].to_numpy() + + +@pytest.mark.filterwarnings("ignore") +def test_glum_weight_spline_survives_offloading(tmp_path): + # Unpickled glum models rebuild their design from two rows, so the + # centring constraint must be baked into the formula, not memorised + opts = dict( + weighted=True, + weight_preexpansion=False, + weight_spline=True, + glm_package="glum", + offload_dir=str(tmp_path), + ) + in_memory = _weights(_model(**opts)) + offloaded = _weights(_model(offload=True, **opts)) + assert np.allclose(in_memory, offloaded) + + +@pytest.mark.filterwarnings("ignore") +@pytest.mark.parametrize("term", ["center(N_bas)", "standardize(L_bas)"]) +def test_bootstrap_survival_uses_each_replicates_own_transform(term, monkeypatch): + opts = dict( + covariates=f"tx_init_bas+followup+trial+{term}", + km_curves=True, + bootstrap_nboot=2, + ) + + def km(s): + s.survival() + return s.km_data.sort(["estimate", "tx_init", "followup"])[ + ["pred", "LCI", "UCI"] + ].to_numpy() + + cached = km(_model(**opts)) + # Without the cached design matrix every model builds its own + monkeypatch.setattr( + survival_pred, + "_cached_predict", + lambda model, X, ref, data: survival_pred._safe_predict(model, data), + ) + uncached = km(_model(**opts)) + assert np.allclose(cached, uncached, equal_nan=True) diff --git a/tests/test_subgroup_bootstrap.py b/tests/test_subgroup_bootstrap.py new file mode 100644 index 0000000..69e8b00 --- /dev/null +++ b/tests/test_subgroup_bootstrap.py @@ -0,0 +1,41 @@ +"""Subgroup bootstrap models are matched by subgroup value, not position.""" + +from collections import Counter + +import polars as pl +import pytest + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data + + +@pytest.mark.parametrize("estimate", ["survival", "hazard"]) +def test_replicate_missing_a_subgroup_is_left_out_of_it(estimate): + data = load_data("SEQdata") + opts = dict(subgroup_colname="sex", bootstrap_nboot=3, seed=42) + opts["km_curves" if estimate == "survival" else "hazard_estimate"] = True + s = SEQuential( + data, + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="ITT", + parameters=SEQopts(**opts), + ) + s.expand() + s.bootstrap() + # The last replicate resamples only sex == 1, so it has no sex == 0 rows + ids = data.filter(pl.col("sex") == 1)["ID"].unique().to_list() + s._boot_samples[-1] = Counter(ids) + s.fit() + + assert s._unique_subgroups == [0, 1] + assert set(s.outcome_model[-1]) == {1} + with pytest.warns(UserWarning, match="1 bootstrap replicate.*sex = 0"): + getattr(s, estimate)() + out = s.km_data if estimate == "survival" else s.hazard_ratio + assert set(out["sex"].unique()) == {0, 1} diff --git a/tests/test_treatment_level_order.py b/tests/test_treatment_level_order.py new file mode 100644 index 0000000..a34417a --- /dev/null +++ b/tests/test_treatment_level_order.py @@ -0,0 +1,55 @@ +"""Treatment weights do not depend on the order of treatment_level.""" + +import numpy as np +import pytest + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data + + +def _weights(dataset, treatment_level, preexpansion, **opts): + s = SEQuential( + load_data(dataset), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="censoring", + parameters=SEQopts( + weighted=True, + weight_preexpansion=preexpansion, + treatment_level=treatment_level, + **opts, + ), + ) + s.expand() + s.fit() + return s.DT.sort(["ID", "trial", "followup"])["weight"].to_numpy() + + +@pytest.mark.filterwarnings("ignore") +@pytest.mark.parametrize("preexpansion", [False, True]) +@pytest.mark.parametrize( + "dataset, levels", + [("SEQdata", [0, 1]), ("SEQdata_multitreatment", [1, 2])], +) +def test_weights_match_with_treatment_level_reversed(dataset, levels, preexpansion): + forward = _weights(dataset, levels, preexpansion) + reverse = _weights(dataset, levels[::-1], preexpansion) + assert np.allclose(forward, reverse) + + +@pytest.mark.filterwarnings("ignore") +@pytest.mark.parametrize("preexpansion", [False, True]) +def test_excused_weights_match_with_treatment_level_reversed(preexpansion): + cols = ["excusedZero", "excusedOne"] + forward = _weights( + "SEQdata", [0, 1], preexpansion, excused=True, excused_colnames=cols + ) + reverse = _weights( + "SEQdata", [1, 0], preexpansion, excused=True, excused_colnames=cols[::-1] + ) + assert np.allclose(forward, reverse) diff --git a/tests/test_undefined_terms.py b/tests/test_undefined_terms.py new file mode 100644 index 0000000..df72f10 --- /dev/null +++ b/tests/test_undefined_terms.py @@ -0,0 +1,86 @@ +"""Transformed terms such as np.log(x), and terms undefined for some rows.""" + +import polars as pl +import pytest + +from pySEQTarget import SEQopts, SEQuential +from pySEQTarget.data import load_data +from pySEQTarget.helpers._col_string import _col_string + +UNDEFINED = r"Model term\(s\) np.log\(x\) produced missing or infinite values" + + +def _data(): + # x is positive for N above its median and negative otherwise + return ( + load_data("SEQdata") + .with_columns( + (pl.col("N") > pl.col("N").median()).cast(pl.Int8).alias("welig"), + (pl.col("N") - pl.col("N").min() + 1).alias("pos"), + ) + .with_columns( + pl.when(pl.col("welig") == 1).then(pl.col("pos")).otherwise(-1.0).alias("x") + ) + ) + + +def _model(**opts): + s = SEQuential( + _data(), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["x", "pos", "N", "L", "P"], + fixed_cols=["sex"], + method="censoring", + parameters=SEQopts(**opts), + ) + s.expand() + return s + + +def _weights(term, **opts): + return dict( + weighted=True, + weight_preexpansion=False, + numerator=f"sex+followup+trial+{term}", + denominator=f"sex+N+L+P+followup+trial+{term}", + **opts, + ) + + +def test_column_reader_skips_the_module_in_a_call(): + assert _col_string(["followup+np.log(pos)+np.sqrt(N)"]) == {"followup", "pos", "N"} + + +@pytest.mark.filterwarnings("ignore") +def test_log_terms_fit_end_to_end(): + s = _model(covariates="tx_init_bas+followup+np.log(pos)", km_curves=True) + s.fit() + s.survival() + assert "np.log(pos)" in s.outcome_model[0]["outcome"].params.index + s = _model(**_weights("np.log(pos)")) + s.fit() + assert "np.log(pos)" in s.denominator_model[0].params.index + + +def test_undefined_outcome_term_is_named(): + s = _model(covariates="tx_init_bas+followup+np.log(x)") + with pytest.raises(ValueError, match=UNDEFINED): + s.fit() + + +def test_undefined_weight_term_is_named(): + s = _model(**_weights("np.log(x)")) + with pytest.raises(ValueError, match=UNDEFINED): + s.fit() + + +@pytest.mark.filterwarnings("ignore") +def test_term_undefined_only_where_predicted_is_named(): + # Fitted on eligible rows, where x > 0, but predicted on every row + s = _model(**_weights("np.log(x)", weight_eligible_colnames=["welig", "welig"])) + with pytest.raises(ValueError, match=UNDEFINED): + s.fit() diff --git a/tests/test_weight_p99.py b/tests/test_weight_p99.py index be856f2..3d55451 100644 --- a/tests/test_weight_p99.py +++ b/tests/test_weight_p99.py @@ -28,8 +28,8 @@ def _coefs(s): def test_weight_p99_truncates_at_p01_p99_percentile_weights(): - # weight_p99=True overrides weight_min/weight_max with the p01/p99 of the - # (untruncated) weights -- these are reported in weight_stats. So it must be + # weight_p99=True truncates at the p01/p99 of the (untruncated) weights -- + # these are reported in weight_stats. So it must be # equivalent to passing those percentile values as explicit bounds, and must # differ from an untruncated weighted fit. p99 = _fit(weight_p99=True) @@ -42,3 +42,26 @@ def test_weight_p99_truncates_at_p01_p99_percentile_weights(): assert np.allclose(_coefs(p99), _coefs(explicit), atol=1e-8) assert not np.allclose(_coefs(p99), _coefs(untruncated), atol=1e-6) + + +def test_weight_p99_leaves_the_supplied_bounds_alone(): + # Each fit, main or replicate, truncates at its own percentiles, but the + # options keep what was supplied, as in SEQTaRget + s = SEQuential( + load_data("SEQdata"), + id_col="ID", + time_col="time", + eligible_col="eligible", + treatment_col="tx_init", + outcome_col="outcome", + time_varying_cols=["N", "L", "P"], + fixed_cols=["sex"], + method="censoring", + parameters=SEQopts(weighted=True, weight_p99=True, bootstrap_nboot=2, seed=1), + ) + s.expand() + s.bootstrap() + s.fit() + assert (s.weight_min, s.weight_max) == (0.0, None) + options = s.collect().options + assert (options.weight_min, options.weight_max) == (0.0, None) diff --git a/tests/test_weight_spline.py b/tests/test_weight_spline.py index 61c1996..ba7976a 100644 --- a/tests/test_weight_spline.py +++ b/tests/test_weight_spline.py @@ -28,7 +28,7 @@ ] BAKED_PLAIN = r"cr\({var}, knots=\[[^]]*\], lower_bound=[-\d.]+, upper_bound=[-\d.]+" -BAKED = BAKED_PLAIN + r', constraints="center"\)' +BAKED = BAKED_PLAIN + r", constraints=\[\[[^]]*\]\]\)" def _model(**opts): @@ -264,7 +264,7 @@ def test_weight_spline_fits_on_either_glm_backend(glm_package): # The basis is centred, which keeps it full rank alongside the model # intercept — glum's solver, unlike statsmodels' pinv, needs that s = _model(weight_preexpansion=False, weight_spline=True, glm_package=glm_package) - assert 'constraints="center"' in s.denominator + assert "constraints=[[" in s.denominator assert s.DT["weight"].to_numpy().min() > 0