Hi，I am confusing about how to use the time and history data. It seems that we should make a model with data of questionnaires k-1 and k-2, but I couldn't find such code implementation. I found you used data of time k only while fitting you model.

#176 exposure_model #207 outcome_model #274 add_covariate_model

@pzivich, When using Singlecrossfit TMLE for a continuous outcome with sm.Gaussian GLM class. I have encountered the following error:

xxx/lib/python3.7/site-packages/zepid/causal/doublyrobust/crossfit.py:1663: RuntimeWarning: invalid value encountered in log
  log = sm.GLM(ys, np.column_stack((h1ws, h0ws)), offset=np.log(probability_to_odds(py_os)),
xxx/lib/python3.7/site-packages/zepid/causal/doublyrobust/crossfit.py:1669: RuntimeWarning: invalid value encountered in log
  ystar0 = np.append(ystar0, logistic.cdf(np.log(probability_to_odds(py_ns)) - epsilon[1] / pa0s))

Here is how I defined the estimator for superleaner as well as the parameter input.

link_i = sm.genmod.families.links.identity()
SL_glm = GLMSL(family = sm.families.family.Gaussian(link=link_i))
GLMSL(family = sm.families.family.Binomial())

sctmle = SingleCrossfitTMLE(dataset = df, exposure='treatment', outcome='paid_amt', continuous_bound = 0.01)
sctmle.exposure_model('gender_cd_F + prospective_risk + age_nbr', GLMSL(family = sm.families.family.Binomial()), bound=[0.01, 0.99])
sctmle.outcome_model('gender_cd_F + prospective_risk + age_nbr', SL_glm)
sctmle.fit(n_splits = 2, n_partitions=3, random_state=12345, method = 'median')
sctmle.summary()

If I uses any other ML estimates such as Lasso, GBM, RandomForest from Sklearn for outcome model estimator, it will work fine. The error only related to use of GLMSL family.

Could you share any idea of the reason of this error and how I can fix this issue? Much appreciated!

One lofty goal is to implement the G-formula. Would need to code two versions; time-fixed and time-varying. The Chapter by Robins & Hernan is good reference. I have code that implements the g-formula using pandas. It is reasonably fast.

TODO: generalize to a class, allow input models then predict, need to determine how to allow users to input custom treatment regimes (all/none/natural course are easy to do), compare results (https://www.ncbi.nlm.nih.gov/pubmed/25140837)

Time-fixed version will be relatively easy to write up

Time-varying will need the ability to specify a large amount of models and specify the order in which the models are fit.

Note; I am also considering reorganizing in v0.2.0 that IPW/g-formula/doubly robust will all be contained within a folder caused causal, rather than adding to the current ipw folder

In the future, I think a zepid.causal.generalize branch would be a useful addition. This branch would contain some generalizability and transportability tools. Specifically, the g-transport formula, inverse probability of sampling weights, inverse odds of sampling weights, and doubly robust generalizability.

Generally, I think I can repurpose a fair bit of the existing code. I need to consider how best to handle the distinction between generalizability (sample from target population) and transportability (sample not from target population). I am imagining that the user would pass in two data sets, the data to estimate the model on, and the other data set to generalize to.

As far as I know, these resources are largely lacking in all other commonly used epi softwares. Plus this is becoming an increasingly hot topic in epi (and I think it will catch on more widely once people recognize you can go from your biased sample to a full population under an additional set of assumptions)

Resources:

• g-transport and IPSW estimators: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5466356/

• inverse odds of sampling: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5860052/

• Doubly Robust estimator: https://arxiv.org/pdf/1902.06080.pdf

• General introduction: https://arxiv.org/pdf/1805.00550.pdf

Notes: Some record of Issa's presentation at CIRG. This is the more difficult doubly robust estimator. It is when only a subset of the cohort has some measure needed for transportability. Rather than throwing out the individual who don't have X2 measures, you can use the process in the arXiv paper. For the nested-nested population, the robust estimator has three parts. b_a(X_1, S) is harder to estimate but you can use the following algorithm

1. model Y as a function of X={X1, X2} among S=1 and A=a

2. Predict \hat{Y} among those with D=1

3. Model \hat{Y} as X1, S in D=1

4. Predict \hat{Y*} in D={0, 1}

Also hard to estimate Pr(S=1|X) becuase X2 only observed for subset. Can use m-estimator to obtain. Can do this by a weighted regression with 1 for D=1 & S=1 and 1/C(X1, S) for D=1 & S=0. This is a little less obvious to me but seems doable

TMLE is not guaranteed to attain nominal coverage when used with machine learning. A simulation paper showing major problems is: https://arxiv.org/abs/1711.07137 As a result, I don't feel like TMLE can continue to be supported with machine learning, especially since it implies the confidence intervals are way too narrow (sometimes resulting in 0% coverage). I know this is a divergence from R's tmleverse, but I would rather enforce the best practice/standards than allow incorrect use of methods

Due to this issue, I will be dropping support for TMLE with machine learning. In place of this, I plan on adding CrossfitTMLE which will support machine learning approaches. The crossfitting will result in valid confidence intervals / inference.

Tentative plan:

• In v0.8.0, TMLE will throw a warning when using the custom_model argument.

• Once the Crossfit-AIPW and Crossfit-TMLE are available (v0.9.0), TMLE will lose that functionality. If users want to use TMLE with machine learning, they will need to use a prior version

Add SNM to the zepid.causal branch. After this addition, all of Robin's g-methods will be implemented in zEpid.

SNM are discussed in the Causal Inference book (https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/) and The Chapter. SAS code for search-based and closed-form solvers is available at the site. Ideally will have both implemented. Will start with time-fixed estimator

Update for version 0.8.0. Below are listed the planned additions

• [x] ~Update how the weights argument works in applicable causal models (IPTW, AIPW, g-formula)~ #106 No longer using this approach

Inverse Probability of Treatment Weights

• [x] Changing entire structure #102

• [x] Figure out why new structure is giving slightly different values for positivity calculations...

• [x] Add g-bounds option to truncate weights

• [x] Update tests for new structure

• [x] Weight argument behaves slightly different (diagnostics now available for either IPTW alone or with other weights)

• [x] New summary function for results

• [x] ~Allowing for continuous treatments in estimation of MSM~ ...for later

• [x] ~Plots available for binary or continuous treatments~ ...for later

Inverse Probability of Censoring Weights

• [x] ~Correction for pooled logistic weight calculation with late-entry occurring~ Raise ValueError if late-entry is detected. The user will need to do some additional work

• [x] Create better documentation for when late-entry occurs for this model

G-formula

• [x] Add diagnostics (prediction accuracy of model)

• [x] Add run_diagnostics()

Augmented IPTW

• [x] Add g-bounds

• [x] Add diagnostics for weights and outcome model

• [x] Add run_diagnostics()

TMLE

• [x] New warning when using machine learning algorithms to estimate nuisance functions #109

• [x] Add diagnostics for weights and outcome model

• [x] Add run_diagnostics()

S-value

• [x] Add calculator for S-values, a (potentially) more informative measure than p-values #107

ReadTheDocs Documentation

• [x] Add S-value

• [x] Update IPTW

• [x] Make sure run_diagnostics() and bound are sufficiently explained

Previously my code might look like:

rossi_with_splines[['age0', 'age1']] = spline(rossi_with_splines, var='age', term=2, restricted=True)
cph = CoxPHFitter().fit(rossi_with_splines.drop('age', axis=1), 'week', 'arrest')

# this part is nasty
df =rossi_with_splines.drop_duplicates('age').sort_values('age')[['age', 'age0', 'age1']].set_index('age')
(df * cph.hazards_[['age0', 'age1']].values).sum(1).plot()

vs

spline_transform, _ = create_spline_transform(df['age'], term=2, restricted=True)
rossi_with_splines[['age0', 'age1']] = spline_transform(rossi_with_splines['age'])

cph = CoxPHFitter().fit(rossi_with_splines.drop('age', axis=1), 'week', 'arrest')

ages_to_consider = np.arange(20, 45))
y = spline_transform(ages_to_consider).dot(cph.hazards_[['age0', 'age1']].values)
plot(ages_to_consider, y)

Version 0.5.0

Features to be implemented:

• [x] Replace AIPW with the more specific AIPTW #57

• [x] Add support for monotone IPMW #55

• [ ] ~~Add support for nonmonotone IPMW #55~~ As I have read further into this, it gets a little complicated (even for the unconditional scenario). Will save for later implementation

• [ ] Add support for continuous treatments in TimeFixedGFormula #49

• [ ] ~~Add stratify option to measures #56~~

• [x] TMLE with continuous outcomes #39

• [x] TMLE with missing data #39 (only applies to missing outcome data)

• [ ] ~~Add support for stochastic interventions into TMLE #52~~ Above two changes to TMLE will take precedence. Stochastic treatments are to be added later

• [ ] ~~Add support for permutation weighting (TBD depending on complexity)~~ Will open a new branch for this project. No idea on how long implementation may take

• [x] Incorporate random error in MonteCarloRR

Maintenance

• [x] Add Python 3.7 support

• [x] Check to see if matplotlib 3 breaks anything. LGTM via test_graphics_manual.py

• [x] Magic-g warning updates for g-formula #63

Later addition, but since statsmodels 0.9.0 has GLIMMIX, I would like to add something to deal with interference for the causal branch. I don't have any part of this worked out, so I will need to take some time to really learn what is happening in these papers

References: https://www.ncbi.nlm.nih.gov/pubmed/21068053 https://onlinelibrary.wiley.com/doi/abs/10.1111/biom.12184 https://github.com/bsaul/inferference

Branch plan:

---causal
      |
      -interference

Verification: inferference the R package has some datasets that I can compare results with

Other: Will need to update requirements to need statsmodels 0.9.0

As noted in #73 and #77 there are some further optional enhancements I can add to MonteCarloGFormula

Items to add:

• [x] Censoring model

• [ ] Competing risk model

Testing:

• [x] Test censoring model works as intended (compare to Keil 2014)

• [ ] Test competing risks. May be easiest to simulate up a quick data set to compare. Don't have anything on hand

The updates to Monte-Carlo g-formula will be added to a future update (haven't decided which version they will make it into)

Optional:

• [x] Reduce memory burden of unneeded replicants

I sometimes run into a MemoryError when replicating Keil et al 2014 with many resamples. A potential way out of this is to "throw away" the observations that are not the final observation for that individual. Can add option low_memory=True to throw out those unnecessary observations. User could return the full simulated dataframe with False.

Using the latest version of pip 22.2.2 I am unable to install the most recent zEpid 0.9.0 release on python 3.7.0

pip install -Iv zepid==0.9.0

ERROR: Could not find a version that satisfies the requirement zepid==0.9.0 (from versions: 0.1.0, 0.1.1, 0.1.2, 0.1.3, 0.1.5, 0.1.6, 0.2.0, 0.2.1, 0.3.0, 0.3.1, 0.3.2, 0.4.0, 0.4.1, 0.4.2, 0.4.3, 0.5.0, 0.5.1, 0.5.2, 0.6.0, 0.6.1, 0.7.0, 0.7.1, 0.7.2, 0.8.0, 0.8.1, 0.8.2) ERROR: No matching distribution found for zepid==0.9.0

I had corresponded with @pzivich over email and am posting our communication here for the benefit of other users.

JD.

Is it possible to program saving figures of directed acyclic graphs (DAGs) using zEpid? E.g. using the M-bias DAG code in the docs, typing plt.savefig('dag.png') only saves a blank PNG. To save it to disk, I'd need to plot the figure then manually click-and-point on the pop-up to save it.

PZ.

Unfortunately, saving the DAGs draw isn't too easy. In the background, I use NetworkX to organize and plot the diagram. NetworkX uses matplotlib, but it doesn't return the matplotlib axes object. So while you can tweak parts of the graph in various ways, NetworkX doesn't allow you to directly access the drawn part of the image. Normally, this isn't a problem but when it gets wrapped up in a class object that returns the matplotlib axes (which is what DirectedAcyclicGraph. draw_dag(...) does) it can lead to the issues you noted.

Currently, the best work-around is to generate the image by hand. Below is some code that should do the trick to match what is output by DirectedAcyclicGraph

import networkx as nx
import matplotlib.pyplot as plt
from zepid.causal.causalgraph import DirectedAcyclicGraph

dag = DirectedAcyclicGraph(exposure='X', outcome="Y")
dag.add_arrows((('X', 'Y'),
                ('U1', 'X'), ('U1', 'B'),
                ('U2', 'B'), ('U2', 'Y')
               ))

fig = plt.figure(figsize=(6, 5))
ax = plt.subplot(1, 1, 1)
positions = nx.spectral_layout(dag.dag)
nx.draw_networkx(dag.dag, positions, node_color="#d3d3d3", node_size=1000, edge_color='black',
                 linewidths=1.0, width=1.5, arrowsize=15, ax=ax, font_size=12)
plt.axis('off')
plt.savefig("filename.png", format='png', dpi=300)
plt.close()

Thanks Paul for the advice!

For longer term, it seems useful to build this or something similar into zEpid graphics to programatically save (complex) DAGs in Python for publication. Possibly using position values from DAGs generated in dagitty, which is handy to quickly graph and analyse complex DAGs. Just a thought.

Cheers

Currently AIPTW only returns RD and RR. TMLE returns those and OR as well. I should add support for OR with AIPTW (even though I am not a huge fan of OR when we have nicer estimands)

I should also add support for all / none, and things like ATT and ATU for TMLE and AIPTW both. Basically I need to look up the influence curve info in the TMLE book(s)

Currently you need to set the np.random.seed outside of the function for reproducibility (which isn't good). I should use a similar RandomState approach that the cross-fit estimators use

I wrote most of the ReadTheDocs documentation 2-3 years ago now. It is dated (and my understanding has expanded), so I should go back and review everything after the v0.9.0 release

Here are some things to consider

• Use a different split than time-fixed and time-varying exposures
• Add a futures section (rather than having embedded in documents)
• Update the LIPTW / SIPTW info (once done)
• Replace Chat Gitter button with GitHub Discussions
• Add SuperLearner page to docs