Nalin Ranjan

Princeton University

Zechao Shang
Work was done mainly when the author was at the University of

Snowflake Inc.

Aaron J. Elmore
Sanjay Krishnan

University of Chicago

We propose , a prototype of a versioned database for efficient collaborative data management.
supports offline collaboration, where users work independently without real-time correspondence. The core of is a critical step of offline collaboration: reconciling divergent branches made by simultaneous data manipulation. We formalize the concept of , a condition under which branches may be reconciled without human intervention, and propose an efficient framework for determining whether two branches are and identifying particular records for manual reconciliation.

#1     #1#1#1     #1     #1     #1     #1#1          #1#1

[Aleen and ClarkAleen and Clark2009]
authorpersonFarhana Aleen and personNathan
Clark. year2009.
Commutativity Analysis for Software
Parallelization: Letting Program Transformations See the Big Picture. In
booktitleProceedings of the 14th International Conference on
Architectural Support for Programming Languages and Operating Systems
(Washington, DC, USA) (seriesASPLOS XIV).
publisherACM, addressNew York, NY, USA,

[Alvaro, Conway, Hellerstein, and MaierAlvaro
et al2017]
authorpersonPeter Alvaro, personNeil
Conway, personJoseph M. Hellerstein, and
personDavid Maier. year2017.
Blazes: Coordination Analysis and Placement for
Distributed Programs.
journalACM Trans. Database Syst.
volume42, number4, Article articleno23
(dateOct. year2017), numpages31 pages.

[Alvaro, Conway, Hellerstein, and
MarczakAlvaro et al2011]
authorpersonPeter Alvaro, personNeil
Conway, personJoseph M. Hellerstein, and
personWilliam R. Marczak. year2011.
Consistency Analysis in Bloom: a CALM and
Collected Approach. In booktitleCIDR 2011, Fifth Biennial
Conference on Innovative Data Systems Research, Asilomar, CA, USA, January
9-12, 2011, Online Proceedings. publisherwww.cidrdb.org,
addressAsilomar, California, pages249–260.

[Attic Labs.Attic Labs.2021]
authorpersonAttic Labs.

authorpersonPeter Bailis.
titleCoordination Avoidance in Distributed
thesistypePh.D. Dissertation. schoolUniversity
of California, Berkeley, USA.

[Bailis, Fekete, Franklin, Ghodsi, Hellerstein,
and StoicaBailis et al2014]
authorpersonPeter Bailis, personAlan
Fekete, personMichael J. Franklin, personAli Ghodsi,
personJoseph M. Hellerstein, and personIon Stoica.
Coordination Avoidance in Database Systems.
journalProc. VLDB Endow. volume8,
number3 (dateNov. year2014),

[Balegas, Duarte, Ferreira, Rodrigues, and
PreguiçaBalegas et al2018]
authorpersonValter Balegas, personSérgio
Duarte, personCarla Ferreira, personRodrigo
Rodrigues, and personNuno Preguiça.
IPA: Invariant-preserving Applications for Weakly
Consistent Replicated Databases.
journalProc. VLDB Endow. volume12,
number4 (dateDec. year2018),

[Bernstein, Lewis, and LuBernstein
et al2000]
authorpersonA. J. Bernstein, personP. M.
Lewis, and personS. Lu. year2000.
Semantic conditions for correctness at different
isolation levels. In booktitleProceedings of 16th
International Conference on Data Engineering (Cat. No.00CB37073).
publisherIEEE Computer Society, addressWashington, DC,
USA, pages57–66.

[Bernstein, Fekete, Guo, Ramakrishnan, and
TammaBernstein et al2006]
authorpersonPhilip A. Bernstein, personAlan
Fekete, personHongfei Guo, personRaghu Ramakrishnan,
and personPradeep Tamma. year2006.
Relaxed-currency Serializability for Middle-tier
Caching and Replication. In booktitleProceedings of the
2006 ACM SIGMOD International Conference on Management of Data (Chicago,
IL, USA) (seriesSIGMOD '06). publisherACM,
addressNew York, NY, USA, pages599–610.

[Bhardwaj, Bhattacherjee, Chavan, Deshpande,
Elmore, Madden, and ParameswaranBhardwaj et al2015]
authorpersonAnant P. Bhardwaj, personSouvik
Bhattacherjee, personAmit Chavan, personAmol
Deshpande, personAaron J. Elmore, personSamuel
Madden, and personAditya G. Parameswaran.
DataHub: Collaborative Data Science & Dataset
Version Management at Scale. In booktitleCIDR 2015,
Seventh Biennial Conference on Innovative Data Systems Research, Asilomar,
CA, USA, January 4-7, 2015, Online Proceedings.

[Brutschy, Dimitrov, Müller, and
VechevBrutschy et al2017]
authorpersonLucas Brutschy, personDimitar
Dimitrov, personPeter Müller, and personMartin
Vechev. year2017.
Serializability for Eventual Consistency:
Criterion, Analysis, and Applications. In
booktitleProceedings of the 44th ACM SIGPLAN Symposium on
Principles of Programming Languages (Paris, France)
(seriesPOPL 2017). publisherACM,
addressNew York, NY, USA, pages458–472.

[Burckhardt, Baldassin, and LeijenBurckhardt
et al2010]
authorpersonSebastian Burckhardt,
personAlexandro Baldassin, and personDaan Leijen.
Concurrent Programming with Revisions and Isolation
Types. In booktitleProceedings of the ACM International
Conference on Object Oriented Programming Systems Languages and
Applications (Reno/Tahoe, Nevada, USA) (seriesOOPSLA
'10). publisherACM, addressNew York, NY, USA,

[Burckhardt, Gotsman, Yang, and
ZawirskiBurckhardt et al2014]
authorpersonSebastian Burckhardt,
personAlexey Gotsman, personHongseok Yang, and
personMarek Zawirski. year2014.
Replicated Data Types: Specification, Verification,
Optimality. In booktitleProceedings of the 41st ACM
SIGPLAN-SIGACT Symposium on Principles of Programming Languages (San Diego,
California, USA) (seriesPOPL '14).
publisherACM, addressNew York, NY, USA,

[Chairunnanda, Daudjee, and
ÖzsuChairunnanda et al2014]
authorpersonPrima Chairunnanda,
personKhuzaima Daudjee, and personM. Tamer
Özsu. year2014.
ConfluxDB: Multi-master Replication for Partitioned
Snapshot Isolation Databases.
journalProc. VLDB Endow. volume7,
number11 (dateJuly year2014),

[Clements, Kaashoek, Zeldovich, Morris, and
KohlerClements et al2015]
authorpersonAustin T. Clements,
personM. Frans Kaashoek, personNickolai Zeldovich,
personRobert T. Morris, and personEddie Kohler.
The Scalable Commutativity Rule: Designing Scalable
Software for Multicore Processors.
journalACM Trans. Comput. Syst.
volume32, number4, Article articleno10
(dateJan. year2015), numpages47 pages.

[Crooks, Pu, Alvisi, and ClementCrooks
et al2017]
authorpersonNatacha Crooks, personYouer Pu,
personLorenzo Alvisi, and personAllen Clement.
Seeing is Believing: A Client-Centric Specification
of Database Isolation. In booktitleProceedings of the ACM
Symposium on Principles of Distributed Computing (Washington, DC, USA)
(seriesPODC '17). publisherAssociation for
Computing Machinery, addressNew York, NY, USA,

[Crooks, Pu, Estrada, Gupta, Alvisi, and
ClementCrooks et al2016]
authorpersonNatacha Crooks, personYouer Pu,
personNancy Estrada, personTrinabh Gupta,
personLorenzo Alvisi, and personAllen Clement.
TARDiS: A Branch-and-Merge Approach To Weak
Consistency. In booktitleProceedings of the 2016
International Conference on Management of Data (San Francisco, California,
USA) (seriesSIGMOD '16). publisherACM,
addressNew York, NY, USA, pages1615–1628.

titleIntroduction to Delta Lake.
Accessed: 2020-09-22.

[DoltHub, Inc.DoltHub, Inc.2021]
authorpersonDoltHub, Inc.

authorpersonPaul Fitzpatrick.

[Gotsman, Yang, Ferreira, Najafzadeh, and
ShapiroGotsman et al2016]
authorpersonAlexey Gotsman, personHongseok
Yang, personCarla Ferreira, personMahsa Najafzadeh,
and personMarc Shapiro. year2016.
'Cause I'm Strong Enough: Reasoning About
Consistency Choices in Distributed Systems. In
booktitleProceedings of the 43rd Annual ACM SIGPLAN-SIGACT
Symposium on Principles of Programming Languages (St. Petersburg, FL, USA)
(seriesPOPL '16). publisherACM,
addressNew York, NY, USA, pages371–384.

[Gramoli and GuerraouiGramoli and
authorpersonVincent Gramoli and
personRachid Guerraoui. year2014.
Democratizing Transactional Programming.
journalCommun. ACM volume57,
number1 (dateJan. year2014),

[Guo, Larson, and RamakrishnanGuo
et al2005]
authorpersonHongfei Guo, personPer-Ake
Larson, and personRaghu Ramakrishnan.
Caching with "Good Enough" Currency, Consistency,
and Completeness. In booktitleProceedings of the 31st
International Conference on Very Large Data Bases
(seriesVLDB '05). publisherVLDB Endowment,
addressTrondheim, Norway, pages457–468.

[Guo, Larson, Ramakrishnan, and GoldsteinGuo
et al2004a]
authorpersonHongfei Guo, personPer-Ake
Larson, personRaghu Ramakrishnan, and
personJonathan Goldstein. year2004a.
Relaxed Currency and Consistency: How to Say "Good
Enough" in SQL. In booktitleProceedings of the 2004 ACM
SIGMOD International Conference on Management of Data (Paris, France)
(seriesSIGMOD '04). publisherACM,
addressNew York, NY, USA, pages815–826.

[Guo, Larson, Ramakrishnan, and GoldsteinGuo
et al2004b]
authorpersonHongfei Guo, personPer-Ake
Larson, personRaghu Ramakrishnan, and
personJonathan Goldstein. year2004b.
Support for Relaxed Currency and Consistency
Constraints in MTCache. In booktitleProceedings of the 2004
ACM SIGMOD International Conference on Management of Data (Paris, France)
(seriesSIGMOD '04). publisherACM,
addressNew York, NY, USA, pages937–938.

authorpersonJoseph M. Hellerstein.
The Declarative Imperative: Experiences and
Conjectures in Distributed Logic.
journalSIGMOD Rec. volume39,
number1 (dateSept. year2010),

[Herlihy and KoskinenHerlihy and
authorpersonMaurice Herlihy and personEric
Koskinen. year2008.
Transactional Boosting: A Methodology for
Highly-concurrent Transactional Objects. In
booktitleProceedings of the 13th ACM SIGPLAN Symposium on
Principles and Practice of Parallel Programming (Salt Lake City, UT, USA)
(seriesPPoPP '08). publisherACM,
addressNew York, NY, USA, pages207–216.

[Huang, Xu, Liu, Elmore, and
ParameswaranHuang et al2017]
authorpersonSilu Huang, personLiqi Xu,
personJialin Liu, personAaron J Elmore, and
personAditya Parameswaran. year2017.
Orpheus DB: bolt-on versioning for relational
journalProceedings of the VLDB Endowment
volume10, number10 (year2017),

[Iterative, Inc.Iterative, Inc.2021]
authorpersonIterative, Inc.
titleDava Version Control.

[Jordan, Banerjee, and BatmanJordan
et al1981]
authorpersonJ. R. Jordan, personJ.
Banerjee, and personR. B. Batman.
Precision Locks. In
booktitleProceedings of the 1981 ACM SIGMOD International
Conference on Management of Data (Ann Arbor, Michigan)
(seriesSIGMOD '81). publisherAssociation for
Computing Machinery, addressNew York, NY, USA,

[Kraska, Hentschel, Alonso, and
KossmannKraska et al2009]
authorpersonTim Kraska, personMartin
Hentschel, personGustavo Alonso, and personDonald
Kossmann. year2009.
Consistency Rationing in the Cloud: Pay Only when
It Matters.
journalProc. VLDB Endow. volume2,
number1 (dateAug. year2009),

[Kulkarni, Nguyen, Prountzos, Sui, and
PingaliKulkarni et al2011]
authorpersonMilind Kulkarni, personDonald
Nguyen, personDimitrios Prountzos, personXin Sui,
and personKeshav Pingali. year2011.
Exploiting the Commutativity Lattice. In
booktitleProceedings of the 32Nd ACM SIGPLAN Conference on
Programming Language Design and Implementation (San Jose, California, USA)
(seriesPLDI '11). publisherACM,
addressNew York, NY, USA, pages542–555.

[Larson, Goldstein, and ZhouLarson
et al2004]
authorpersonP. Ake. Larson, personJ.
Goldstein, and personJ. Zhou.
MTCache: transparent mid-tier database caching in
SQL server. In booktitleProceedings. 20th International
Conference on Data Engineering. publisherIEEE Computer Society,
addressWashington, DC, USA, pages177–188.

[Li, Leitão, Clement, Preguiça,
Rodrigues, and VafeiadisLi et al2014]
authorpersonCheng Li, personJoão
Leitão, personAllen Clement, personNuno
Preguiça, personRodrigo Rodrigues, and
personViktor Vafeiadis. year2014.
Automating the Choice of Consistency Levels in
Replicated Systems. In booktitleProceedings of the 2014
USENIX Conference on USENIX Annual Technical Conference (Philadelphia, PA)
(seriesUSENIX ATC'14). publisherUSENIX
Association, addressBerkeley, CA, USA,

[Li, Porto, Clement, Gehrke, Preguiça, and
RodriguesLi et al2012]
authorpersonCheng Li, personDaniel Porto,
personAllen Clement, personJohannes Gehrke,
personNuno Preguiça, and personRodrigo
Rodrigues. year2012.
Making Geo-replicated Systems Fast As Possible,
Consistent when Necessary. In booktitleProceedings of the
10th USENIX Conference on Operating Systems Design and Implementation
(Hollywood, CA, USA) (seriesOSDI'12).
publisherUSENIX Association, addressBerkeley, CA,
USA, pages265–278.

[Li, Preguiça, and RodriguesLi
et al2018]
authorpersonCheng Li, personNuno Preguiça, and personRodrigo Rodrigues.
Fine-grained consistency for geo-replicated
systems. In booktitle2018 USENIX Annual Technical
Conference (USENIX ATC 18). publisherUSENIX Association,
addressBoston, MA, pages359–372.

[Maddox, Goehring, Elmore, Madden,
Parameswaran, and DeshpandeMaddox et al2016]
authorpersonMichael Maddox, personDavid
Goehring, personAaron J Elmore, personSamuel Madden,
personAditya Parameswaran, and personAmol
Deshpande. year2016.
Decibel: The relational dataset branching system.
journalProceedings of the VLDB Endowment
volume9, number9 (year2016),

[Narula, Cutler, Kohler, and MorrisNarula
et al2014]
authorpersonNeha Narula, personCody Cutler,
personEddie Kohler, and personRobert Morris.
Phase Reconciliation for Contended In-memory
Transactions. In booktitleProceedings of the 11th USENIX
Conference on Operating Systems Design and Implementation (Broomfield, CO)
(seriesOSDI'14). publisherUSENIX Association,
addressBerkeley, CA, USA, pages511–524.

[Olston, Loo, and WidomOlston
et al2001]
authorpersonChris Olston, personBoon Thau
Loo, and personJennifer Widom.
Adaptive Precision Setting for Cached Approximate
Values. In booktitleProceedings of the 2001 ACM SIGMOD
International Conference on Management of Data (Santa Barbara, California,
USA) (seriesSIGMOD '01). publisherACM,
addressNew York, NY, USA, pages355–366.

[Olston and WidomOlston and Widom2000]
authorpersonChris Olston and
personJennifer Widom. year2000.
Offering a Precision-Performance Tradeoff for
Aggregation Queries over Replicated Data. In
booktitleProceedings of the 26th International Conference on
Very Large Data Bases (seriesVLDB '00).
publisherMorgan Kaufmann Publishers Inc., addressSan
Francisco, CA, USA, pages144–155.

[Pachyderm Inc.Pachyderm Inc.2021]
authorpersonPachyderm Inc.

authorpersonCalton Pu.
Relaxing the Limitations of Serializable
Transactions in Distributed Systems. In
booktitleProceedings of the 5th Workshop on ACM SIGOPS
European Workshop: Models and Paradigms for Distributed Systems Structuring
(Mont Saint-Michel, France) (seriesEW 5).
publisherACM, addressNew York, NY, USA,

[Rinard and DinizRinard and Diniz1997]
authorpersonMartin C. Rinard and
personPedro C. Diniz. year1997.
Commutativity Analysis: A New Analysis Technique
for Parallelizing Compilers.
journalACM Trans. Program. Lang. Syst.
volume19, number6 (dateNov.
year1997), pages942–991.

[Shang and YuShang and Yu2017]
authorpersonZechao Shang and
personJeffrey Xu Yu. year2017.
My Weak Consistency is Strong. In
booktitleCIDR 2017, 8th Biennial Conference on Innovative
Data Systems Research, Chaminade, CA, USA, January 8-11, 2017, Online
Proceedings. publisherwww.cidrdb.org.

[Shang, Yu, and ElmoreShang
et al2018]
authorpersonZechao Shang, personJeffrey Xu
Yu, and personAaron J. Elmore.
RushMon: Real-time Isolation Anomalies Monitoring.
In booktitleProceedings of the 2018 International Conference
on Management of Data (Houston, TX, USA) (seriesSIGMOD
'18). publisherACM, addressNew York, NY, USA,

[Shapiro, Preguiça, Baquero, and
ZawirskiShapiro et al2011]
authorpersonMarc Shapiro, personNuno
Preguiça, personCarlos Baquero, and personMarek
Zawirski. year2011.
Conflict-free Replicated Data Types. In
booktitle13th International Conference on Stabilization,
Safety, and Security of Distributed Systems (seriesSSS
2011). publisherSpringer LNCS volume 6976,

[Shasha, Simon, and ValduriezShasha
et al1992]
authorpersonDennis Shasha, personEric
Simon, and personPatrick Valduriez.
Simple Rational Guidance for Chopping Up
journalSIGMOD Rec. volume21,
number2 (dateJune year1992),

[Sivaramakrishnan, Kaki, and
JagannathanSivaramakrishnan et al2015]
authorpersonKC Sivaramakrishnan,
personGowtham Kaki, and personSuresh Jagannathan.
Declarative Programming over Eventually Consistent
Data Stores. In booktitleProceedings of the 36th ACM
SIGPLAN Conference on Programming Language Design and Implementation
(Portland, OR, USA) (seriesPLDI '15).
publisherACM, addressNew York, NY, USA,

[Terry, Theimer, Petersen, Demers, Spreitzer,
and HauserTerry et al1995]
authorpersonDouglas B Terry, personMarvin M
Theimer, personKarin Petersen, personAlan J Demers,
personMike J Spreitzer, and personCarl H Hauser.
Managing update conflicts in Bayou, a weakly
connected replicated storage system.
journalACM SIGOPS Operating Systems Review
volume29, number5 (year1995),

[Tripp, Manevich, Field, and SagivTripp
et al2012]
authorpersonOmer Tripp, personRoman
Manevich, personJohn Field, and personMooly
Sagiv. year2012.
JANUS: Exploiting Parallelism via Hindsight. In
booktitleProceedings of the 33rd ACM SIGPLAN Conference on
Programming Language Design and Implementation (Beijing, China)
(seriesPLDI '12). publisherACM,
addressNew York, NY, USA, pages145–156.

authorpersonTodd L. Veldhuizen.
Transaction Repair: Full Serializability Without
journalCoRR  volumeabs/1403.5645

[von Koch, Manilov, Vasiladiotis, Cole, and
Frankevon Koch et al2018]
authorpersonTobias J. K. Edler von Koch,
personStanislav Manilov, personChristos Vasiladiotis,
personMurray Cole, and personBjörn Franke.
Towards a Compiler Analysis for Parallel
Algorithmic Skeletons. In booktitleProceedings of the 27th
International Conference on Compiler Construction (Vienna, Austria)
(seriesCC 2018). publisherACM,
addressNew York, NY, USA, pages174–184.

[Wang, Dinh, Lin, Xie, Zhang, Cai, Chen, Ooi,
and RuanWang et al2018]
authorpersonSheng Wang, personTien Tuan Anh
Dinh, personQian Lin, personZhongle Xie,
personMeihui Zhang, personQingchao Cai,
personGang Chen, personBeng Chin Ooi, and
personPingcheng Ruan. year2018.
Forkbase: An efficient storage engine for
blockchain and forkable applications.
journalProceedings of the VLDB Endowment
volume11, number10 (year2018),

[Wang, Mu, Cui, Yi, Chen, and LiWang
et al2016]
authorpersonZhaoguo Wang, personShuai Mu,
personYang Cui, personHan Yi, personHaibo
Chen, and personJinyang Li. year2016.
Scaling Multicore Databases via Constrained
Parallel Execution. In booktitleProceedings of the 2016
International Conference on Management of Data (San Francisco, California,
USA) (seriesSIGMOD '16). publisherACM,
addressNew York, NY, USA, pages1643–1658.

authorpersonW. E. Weihl.
Commutativity-based concurrency control for
abstract data types.
journalIEEE Trans. Comput. volume37,
number12 (dateDec year1988),

authorpersonW. E. Weihl.
Local Atomicity Properties: Modular Concurrency
Control for Abstract Data Types.
journalACM Trans. Program. Lang. Syst.
volume11, number2 (dateApril
year1989), pages249–282.

[Whittaker and HellersteinWhittaker and
authorpersonMichael Whittaker and
personJoseph M. Hellerstein. year2018.
Interactive Checks for Coordination Avoidance.
journalProc. VLDB Endow. volume12,
number1 (dateSept. year2018),

[Wu, Yu, and PuWu et al1992]
authorpersonKun-Lung Wu, personPhilip S.
Yu, and personCalton Pu. year1992.
Divergence Control for Epsilon-Serializability. In
booktitleProceedings of the Eighth International Conference
on Data Engineering. publisherIEEE Computer Society,
addressWashington, DC, USA, pages506–515.

[Wu, Yu, and PuWu et al1997]
authorpersonKun-Lung Wu, personPhilip S.
Yu, and personCalton Pu. year1997.
Divergence Control Algorithms for Epsilon
journalIEEE Trans. on Knowl. and Data Eng.
volume9, number2 (dateMarch
year1997), pages262–274.

[Wu, Chan, and TanWu et al2016]
authorpersonYingjun Wu, personChee-Yong
Chan, and personKian-Lee Tan.
Transaction Healing: Scaling Optimistic Concurrency
Control on Multicores. In booktitleProceedings of the 2016
International Conference on Management of Data (San Francisco, California,
USA) (seriesSIGMOD '16). publisherACM,
addressNew York, NY, USA, pages1689–1704.

[Xie, Su, Littley, Alvisi, Kapritsos, and
WangXie et al2015]
authorpersonChao Xie, personChunzhi Su,
personCody Littley, personLorenzo Alvisi,
personManos Kapritsos, and personYang Wang.
High-performance ACID via Modular Concurrency
Control. In booktitleProceedings of the 25th Symposium on
Operating Systems Principles (Monterey, California)
(seriesSOSP '15). publisherACM,
addressNew York, NY, USA, pages279–294.

[Zellag and KemmeZellag and
authorpersonK. Zellag and personB.
Kemme. year2011.
Real-time quantification and classification of
consistency anomalies in multi-tier architectures. In
booktitle2011 IEEE 27th International Conference on Data
Engineering. publisherIEEE Computer Society,
addressWashington, DC, USA, pages613–624.

[Zellag and KemmeZellag and Kemme2012]
authorpersonKamal Zellag and personBettina
Kemme. year2012.
How Consistent is Your Cloud Application?. In
booktitleProceedings of the Third ACM Symposium on Cloud
Computing (San Jose, California) (seriesSoCC '12).
publisherACM, addressNew York, NY, USA, Article
articleno6, numpages14 pages.

[Zellag and KemmeZellag and Kemme2014]
authorpersonKamal Zellag and personBettina
Kemme. year2014.
Consistency Anomalies in Multi-tier Architectures:
Automatic Detection and Prevention.
journalThe VLDB Journal volume23,
number1 (dateFeb. year2014),

[Zhang, Power, Zhou, Sovran, Aguilera, and
LiZhang et al2013]
authorpersonYang Zhang, personRussell
Power, personSiyuan Zhou, personYair Sovran,
personMarcos K. Aguilera, and personJinyang Li.
Transaction Chains: Achieving Serializability with
Low Latency in Geo-distributed Storage Systems. In
booktitleProceedings of the Twenty-Fourth ACM Symposium on
Operating Systems Principles (Farminton, Pennsylvania)
(seriesSOSP '13). publisherACM,
addressNew York, NY, USA, pages276–291.


Causal Inference with Orthogonalized Regression:
Taming the Phantom
Stephen Bates1, Edward Kennedy2, Robert Tibshirani3, Valérie Ventura4 and
Larry Wasserman5

  * 1

Departments of EECS and Statistics, University of California Berkeley;
<EMAIL_ADDRESS>
  * 2

Department of Statistics and Data Science, Carnegie Mellon University;
<EMAIL_ADDRESS>
  * 3

Departments of Biomedical Data Science and Statistics, Stanford University;
<EMAIL_ADDRESS>
  * 4

Department of Statistics and Data Science, Carnegie Mellon University;
<EMAIL_ADDRESS>
  * 5

Departments of Statistics and Data Science and of Machine Learning, Carnegie
Mellon University<EMAIL_ADDRESS>
###### Abstract

Standard regression methods can lead to inconsistent estimates of causal
effects when there are time-varying treatment effects and time-varying
covariates. This is due to certain non-confounding latent variables that
create colliders in the causal graph. These latent variables, which we call
phantoms, do not harm the identifiability of the causal effect, but they
render naive regression estimates inconsistent. Motivated by this, we ask: how
can we modify regression methods so that they hold up even in the presence of
phantoms? We develop an estimator for this setting based on regression
modeling (linear, log-linear, probit and Cox regression), proving that it is
consistent for the causal effect of interest. In particular, the estimator is
a regression model fit with a simple adjustment for collinearity, making it
easy to understand and implement with standard regression software. From a
causal point of view, the proposed estimator is an instance of the parametric
g-formula. Importantly, we show that our estimator is immune to the _null
paradox_ that plagues most other parametric g-formula methods.

##  1 Introduction

Regression adjustment is likely the most commonly-encountered causal inference
method: the analyst regresses an outcome variable against treatment variables
and observed covariates and then interprets the coefficient of the treatment
variables as causal estimates. But when analyzing longitudinal data with time-
varying treatments and covariates, regression methods can lead to inconsistent
estimates of causal effects (Robins,, 1986). This is due to latent variables
that affect both the outcome and covariates but are _not_ unmeasured
confounders. For example, suppose we wish to understand the impact of a
medication on survival, and also have weight and blood pressure as covariates.
Then, the (latent) overall health of a subject affects both their outcome and
covariates, but not necessarily their choice to take medication. We call these
unmeasured variables phantoms, distinguishing them from unmeasured
confounders. Phantoms do not disrupt the identifiability of causal effects,
but they cause problems for standard regression estimates, biasing them even
in the large-sample limit. From the point of view of causal graphs, this is
because phantoms create colliders—variables with two arrows pointing into them
as in Figure 2. Methods have been developed to avoid this problem, such as
marginal structural models (Robins et al.,, 2000), structural nested models
(Robins,, 2000), and g-formula methods (Naimi et al.,, 2016). Taking a
different approach, we study this problem through the lens of collinearity in
regression models. We explain how phantoms invalidate naive estimates and
derive a collinearity-adjusted regression estimate that is valid even in the
presence of phantoms.

###  1.1 Motivating examples

Figure 1: (Left) Linear regression in a Gaussian model. We plot the regression
coefficient of $A_{1}$ from naive linear regression and the proposed method.
(Right) Cox regression with surival outcomes. We report on the regression
coefficient of $A_{1}$ for naive Cox regression and the proposed method. In
both panels, there is no true causal effect (solid horizontal line).

To see the shortcomings of standard regression modeling, we first consider a
simple linear model with Gaussian errors. The outcome is a Gaussian random
variable $Y$. We have three observed explanatory variables: the treatment
level at time $1$, $A_{1}$, an observed covariate $X_{2}$ that is affected by
$A_{1}$, and the treatment level at time $2$, $A_{2}$, which is affected by
$X_{2}$. There is no unmeasured confounding, but $X_{2}$ is correlated with
$Y$. In the language of this paper, there is a phantom variable—see Section 2.
Lastly, there is no treatment effect from $A_{1}$ or $A_{2}$ to $Y$. See the
right side of Figure 2 for a visualization of these dependencies.

With $Y$ Gaussian and mean a linear function of $(A_{1},X_{2},A_{2})$, we
generate $n=100,000$ data points from the model and then fit a standard linear
regression of $Y$ onto $(A_{1},X_{2},A_{2})$. The left side of Figure 1 shows
the estimated regression coefficient of $A_{1}$ versus the correlation $\rho$
between $A_{1}$ and $X_{2}$. The estimate deviates farther from zero as $\rho$
increases, which is misleading because there is no causal effect of $A_{1}$ on
$Y$. By contrast, orthogonalized regression, a method that will be described
in this work, yields the correct estimate up to simulation error. We will
revisit this simple setting again soon in Section 2.2, where we derive
formulae underlying this simulation result.

The behavior we see in the above is widespread. To demonstrate this, we carry
out a second experiment in a parallel setup, but with survival outcomes. In
this, $Y$ is now a survival time, but the setting is otherwise unchanged. In
particular, $Y$ is correlated with $X_{2}$, but there is again no true causal
effect. We run a Cox regression and report the coefficient for $A_{1}$ in the
right side of Figure 1, which reveals the and same trend; as the effect from
$A_{1}$ to $X_{2}$ increases, the regression coefficient for $A_{1}$ becomes
more misleading.

###  1.2 Related work

Causal inference with longitudinal data was first put on firm formal
foundations in Robins, (1986). It has subsequently blossomed into a rich body
of work; see Hernán and Robins, (2020) for a recent overview. Within this line
of work, marginal structural models (MSMs) (Robins et al.,, 2000; Hernán et
al.,, 2001) are commonly used to precisely define the causal quantity of
interest, and are readily fit using IPW estimators. An MSM directly models the
effect of the treatment $A_{t}$ on the outcome $Y_{t}$, without having to
specify a conditional model for $Y_{t}$ given the treatment and covariates
$(A_{t},X_{t})$). Finding optimal estimators for causal quantities is also of
interest, and much effort has gone into deriving semi- and non-parametric
efficient and doubly robust estimators for this setting; see, e.g., Bang and
Robins, (2005); van der Laan and Rose, (2018); Lendle et al., (2017);
Schomaker et al., (2019); Luedtke et al., (2017). In contrast, here we do not
concentrate on MSMs but rather focus on modelling the outcome $Y_{t}$
parametrically and estimating causal effects using regression methods. This
approach does result in estimates of the marginal effect, but they arise from
a conditional model that includes both the observed covariates and treatment
status rather than a marginal model only in terms of the treatment status. Our
motivations for taking this route is that regression methods are widely
available and familiar to data analysts.

In our setting, the causal effect is given by the g-formula (Robins, (1986);
see (3)) and can be estimated by plugging in estimates for the unknown
quantities. When plugging in parametric estimates, we will call such methods
_parametric g-formula_ methods; see Naimi et al., (2016) for a recent overview
and McGrath et al., (2020) for a software library implementing these methods.
Parametric g-formula methods have been scrutinized for their vulnerability to
the _null paradox_ (Robins and Wasserman,, 1997): they can give inconsistent
estimates of causal effects, especially when there is no true effect. McGrath
et al., (2022) evaluate the degree of possible bias and provide
recommendations, and Evans and Didelez, (2021) derive a set of models that
avoid the null paradox. Here, we will also introduce a set of parametric
g-formula methods that avoid the null paradox, although, unlike Evans and
Didelez, (2021), we are motivated by regression modeling and introduce a
distinct set of models. In principle, one can plugin nonparametric estimates
which then avoid the null paradox. But, as the number of time points
increases, we face an explosion of dimensionality and so, once again, the
estimate is inconsistent due to the curse of dimensionality.

In the present work, we consider several families of regression models: linear
regression, log-linear regression, probit regression, and Cox regression. We
will focus on the latter to highlight our method’s applicability to
challenging, nonlinear regression tasks. In more detail, we pay special
attention to the Cox proportional hazards model with time-varying covariates
(Tian et al.,, 2005; Wang et al.,, 2017). Correctly analyzing time-varying
treatment effects with these models does require the additional causal
machinery mentioned above: marginal structural Cox models have been developed
to allow for flexible specification (Xiao et al.,, 2014) and efficient
estimation (Zheng et al.,, 2016; Dáaz,, 2019) of causal effects. Similarly,
parametric g-formula methods for survival data have been developed previously
in Keil et al., (2014); Young and Tchetgen Tchetgen, (2014); Wen et al.,
(2021).

Our approach in the linear model case is closely related to the rich body of
work on identifiability and estimation in linear structural equation models.
In particular, in this setting our estimator is the usual maximum likelihood
estimator (MLE) for Gaussian linear models. However, we use a new
parameterization that enforces certain orthogonality conditions to yield
simple computations. Leveraging orthogonality to compute the MLE in linear
models goes back classically to Frisch and Waugh, (1933), and ideas similar to
those in our work are implicit in the computational schemes in Drton et al.,
(2009) and Nandy et al., (2017); see also Shpitser et al., (2018). Recently,
Guo and Perković, (2021) showed that such estimators are not only optimal for
the Gaussian case, but are also optimal among estimators that use only the
empirical covariance matrix for (non-Gaussian) linear models. More broadly,
ideas of orthogonality are ever-present in econometrics and causal inference;
see, e.g., Lancaster, (2002); Kennedy, (2016); Chernozhukov et al., (2018).
However, the role of orthogonality in this work is distinct from these
previous instances. Here, we use orthogonality to break certain treatment to
covariate edges in the causal graph, thereby removing causal paths that would
otherwise invalidate the causal effect estimates derived from regression
modeling while avoiding the null paradox.

###  1.3 Our contribution

We propose a new estimator, _orthogonalized regression_ , to obtain causal
estimates via regression modeling for data with time-varying treatments and
covariates. We develop the case of linear regression, log-linear regression,
probit regression and Cox regression, and we expect that further cases can
also be handled. The method is computationally and conceptual simple; it
requires only a sequence of standard regression fits and the operations are
easily understood in terms of collinearity, a familiar topic in regression
modeling. Secondly, we prove that our proposed method is a parametric
g-formula method that overcomes the null paradox; it gives valid estimates
even under the null.

Our presentation proceeds as follows. In Section 2 we introduce our notation,
explain the difficulty with regression modeling in our setting, and introduce
our estimator. In Section 3, we show that the estimator is indeed estimating
the causal effect for linear, log-linear, probit, and Cox regression. In
Section 4, we evaluate the proposed method in simple simulated examples, and
then in Section 5 we demonstrate it’s use in a realistic setting.

##  2 What is the problem?

###  2.1 Notation

Suppose we observe data on $n$ individuals over time. The time ordered data on
the $i^{\rm th}$ subject are

$\left(X_{1}^{(i)},A_{1}^{(i)},\right),\ldots,\left(X_{t}^{(i)},A_{t}^{(i)}\right),\ldots,\left(X_{T}^{(i)},A_{T}^{(i)}\right),Y^{(i)}.$

Here, $X_{t}^{(i)}$ is a vector of covariates $A_{t}^{(i)}$ is treatment or
exposure, and $Y^{(i)}$ is the outcome of interest. The vector $X_{t}^{(i)}$
can include past values of the outcome. The reader should keep in mind the
case where $X_{t}^{(i)}$ corresponds to the same covariate measured at
different time points as the representative example of our setup. But this is
not a requirement: $X_{t}^{(i)}$ could also have different dimensions and
represent different covariates for different $t$. In what follows, we omit the
superscript $(i)$ when considering a generic subject and employ the notation
$\overline{X}_{t}=(X_{1},\ldots,X_{t})$, $\overline{X}=(X_{1},\ldots,X_{T})$,
$\overline{A}_{t}=(A_{1},\ldots,A_{t})$ and
$\overline{A}=(A_{1},\ldots,A_{T})$.

We are interested in the effect of treatment on $Y$. Let $Y^{\overline{a}}$ be
the value that $Y$ would have if the treatment variable
$\overline{A}=(A_{1},\ldots,A_{T})$ was set to a particular value
$\overline{a}=(a_{1},\ldots,a_{T})$. We are interested in

$\psi(\overline{a})=\mbox{$\mathbb{E}$}[Y^{\overline{a}}].$ (1)

This is the average value of $Y$ if we were to intervene to set
$\overline{A}=\overline{a}$ for everyone in the population—it answers the
question “if we changed $\overline{a}$ in the same way for every unit, what
would happen to the population?” Our modeling choices will imply that $\psi$
is parameterized by a finite vector, denoted by $\theta\in\mathbb{R}^{d}$. A
basic question of interest is whether or not there is any treatment effect,
that is, whether or not $\psi$ is a constant. A more precise question is
exactly how the different coordinates of $\overline{a}$ affect the outcome.
This corresponds to estimating the whole function $\psi$.

###  2.2 A linear regression example

To understand the shortcomings of standard regression modeling, we return to
the motivating example on the left of Figure 1, described by the DAG in the
right side of Figure 2. To facilitate explicit formulae for the effects seen
in Figure 2, we further assume
$(A_{1},X_{2},A_{2},Y,U)\sim\mathcal{N}(0,\Sigma)$, with covariance matrix
$\Sigma$ determined by size of the effects $A_{1}\to X_{2}$, $X_{2}\to A_{2}$,
$U\to X_{2}$ and $U\to Y$. We standardize all variables to have variance 1,
and refer to these correlations as $\rho_{A\to X},\rho_{X\to A},\rho_{U\to X}$
and $\rho_{U\to Y}$, respectively. It suffices to vary $\rho_{U\to X}$ and set
$\rho_{U\to Y}=1$, which reduces the model to three parameters. The implied
covariance matrix for the observed quantities $(A_{1},X_{2},A_{2},Y)$ is a
function of these three parameters:

$\begin{bmatrix}1&\rho_{A\to X}&\rho_{A\to X}\rho_{X\to A}&0\\\ \rho_{A\to
X}&1&\rho_{X\to A}&\rho_{U\to X}\\\ \rho_{A\to X}\rho_{X\to A}&\rho_{X\to
A}&1&\rho_{X\to A}\rho_{U\to X}\\\ 0&\rho_{U\to X}&\rho_{X\to A}\rho_{U\to
X}&1\end{bmatrix}.$

We assume that the correlations are such that this is a valid covariance
matrix. A calculation yields the population regression coefficients of $Y$ on
$(A_{1},X_{2},A_{2})$:

$\beta=\left(\frac{-\rho_{A\to X}\rho_{U\to X}}{1-\rho_{A\to
X}^{2}},\frac{\rho_{U\to X}}{1-\rho_{A\to X}^{2}},0\right)$ (2)

from which we confirm that the regression coefficient for $A_{1}$ is
misleading: it is nonzero, even though there is no treatment effect. This is
the result of the collider bias from including $X_{2}$ in the regression
model. On the other hand, the regression coefficient for $A_{2}$ is correct in
this case; including $X_{2}$ in the regression model correctly accounts for
the confounding effect of $X_{2}$ on $A_{2}$. If we instead removed $X_{2}$
from the regression model, we would avoid the collider bias, but fail to
account for the confounding, so the coefficient for $A_{2}$ would be biased.
As previewed in our experiment, we will soon introduce an alternative to this
regression model that correctly identifies the causal effects.

Notice that (2) provides a hint about when the regression model will fail
badly. As $|\rho_{A\to X}|$ increases, the regression coefficient for $A_{1}$
deviates farther from zero. The magnitude of the collider induced bias varies
as with impact of the unobserved variable times the magnitude of the effect
from $A_{1}$ to $X_{2}$, $\rho_{A\to X}$. Thus, the collider bias can be of
magnitude comparable to the effect of confounding. This is an appreciable
effect that cannot typically be ignored. In summary, the regression model has
worse bias when there is a strong effect from the time-varying treatment to
the time-varying covariate.

###  2.3 Phantoms variables: why naive modeling fails

Why do straightforward regression models fail with time-varying treatments?
The heart of the issue is what we call _phantom_ variables: unmeasured
variables that link measured variables with the response. For example,
consider the examples from Robins and Wasserman, (1997) with measured
variables $(A_{1},X_{2},A_{2},Y)$ and graphical structures in Figure 2. Here,
the unmeasured $U$ is a phantom, linking the observed variables $X_{2}$ with
the response $Y$. An example of a phantom is undelying health status of a
subject. Then $U$ is likely associated with the covariates $X_{2}$ and outcome
$Y$. Note that, even in randomized studies, there will likely be many
phantoms. Phantoms are _not_ confounders. Unlike the case with unmeasured
confounders, even with unmeasured phantoms the causal effect is identified and
can be estimated using appropriate methods.

$A_{1}$$X_{2}$$A_{2}$$Y$$U$

$A_{1}$$X_{2}$$A_{2}$$Y$$U$

Figure 2: Example from Robins and Wasserman (1997). Treatments
$(A_{1},A_{2})$, covariate $X_{2}$, outcome $Y$ and phantom $U$. (Left) The
treatment affects the outcome $Y$. There are no unmeasured confounders.
(Right) There is no treatment effect as there are no arrows from $A_{1}$ to
$Y$ or from $A_{2}$ to $Y$. The covariate $X_{2}$ is a collider, meaning that
two arrowheads meet at $X_{2}$. This implies that $Y$ and $(A_{1},A_{2})$ are
dependent conditional on $X_{2}$.

However, phantom variables invalidate a naive regression analysis in the
following way. The variable $X_{2}$ in Figure 2 is an observed confounder that
must be included in the model to get a consistent estimate of the treatment
effect. However, by (naively) including $X_{2}$, the analyst is conditioning
on a post-treatment variable, which creates _collider bias_ and leads to
invalid estimates of, and inference for, the causal effect.

For readers familiar with directed graphs, we see from Figure 2 that $X_{2}$
is a collider and this is why $Y$ is conditionally dependent on
$(A_{1},A_{2})$ given $X_{2}$, even when there is no treatment effect. For
those not familiar with directed graphs, recall that we computed the
conditional correlation algebraically in the previous section and saw that $Y$
and $A_{1}$ were not conditionally independent.

###  2.4 Causal identifiability and the g-formula

Naive regression modeling fails not because estimation is impossible in this
setting; the causal effect is indeed identifiable, and so it can be estimated
with the right methods. Robins, (1986) showed that, if there are no unmeasured
confounders (and subject to two technical conditions described in Appendix B)
the causal effect $\psi(\overline{a})$ in (1) is given by

$\psi(\overline{a})=\int\cdots\int\mu(\overline{a},\overline{x})\
\prod_{t=1}^{T}p(x_{t}|\overline{x}_{t-1},\overline{a}_{t-1})dx_{t}$ (3)

where

$\mu(\overline{a},\overline{x})=\mbox{$\mathbb{E}$}[Y|\overline{X}=\overline{x},\overline{A}=\overline{a}].$

The intuition for the formula is this as follows. Take the joint density

$p(x_{1})p(a_{1}|x_{1})p(x_{2}|x_{1},a_{1})p(a_{2}|x_{1},x_{2},a_{1})\cdots.$

Now replace each $p(a_{t}|{\rm past})$ with a point mass at $a_{t}$ (as if we
intervened to set the treatment to be $a_{t}$ for all units) and take the mean
of $Y$ from the resulting distribution.

Using the g-formula as the starting point for estimation has led to the
development of _g-formula methods_. Nonparametric approaches estimate the
density of $p(x_{s}\mid\overline{x}_{s-1},\overline{a}_{s})$, but this has
high statistical complexity: many samples are needed for reliable density
estimation. On the other hand, parametric models are known to be vulnerable to
the _null paradox_ : naively plugging in parametric models can result in
inconsistent estimates even under the null (Robins and Wasserman,, 1997). In
this work, we show how to use regression models that fit well with the
g-formula and avoid the null paradox. This enables the use of familiar
regression tools to correctly estimate causal effects.

###  2.5 Proposed method: orthogonalized regression

To adapt regression models for use with time-varying treatments and
covariates, we propose _orthogonalized regression_. In particular, we proposed
modifications to linear regression, Poisson regression, Cox regression and so
on, that correctly handle time-varying treatments and covariates. The key idea
is that before model fitting, we replace each $X_{t}$ with an _orthogonalized_
version $\widetilde{X}_{t}$, where $\widetilde{X}_{t}$ is the residual of
$X_{t}$ when we regress out all past variables
$(X_{1},A_{1},\dots,X_{t-1},A_{t-1})$, as outlined in Algorithm 1. To obtain
causal estimates, we then regress $Y$ onto
$(\widetilde{X}_{1},A_{1},\dots,\widetilde{X}_{T},A_{T})$ and extract the
fitted coefficients for $A_{1},\dots,A_{T}$.

The motivation is that because naive regression modeling fails due to collider
bias, if we can remove the collider bias we can obtain reliable estimates of
the causal effect. To remove the colider bias, we remove the edges from
$(\overline{X}_{t-1},\overline{A}_{t-1})$ to $X_{t}$ by residualizing out the
past information. For example, in Figure 3, the orthogonalization of $X_{2}$
removes the edge from $A_{1}$ to $X_{2}$, so that $X_{2}$ is no longer a
collider. In the language of regression modeling, the confounding is due to
collinearity between treatment variables and subsequent covariates, so
orthogonalized regression removes this collinearity before fitting the final
regression model.

The remainder of this work will make precise the causal guarantees of this
method in various settings.

Algorithm 1 Orthogonalized regression

Input: data $\left(\overline{X}^{(i)},\overline{A}^{(i)},Y^{(i)}\right)\text{
for subjects }i=1,\dots,n$, regression fitting subroutines
$\textsc{regression}^{(X)}$ and $\textsc{regression}^{(Y)}$

for t = 1,…,T do

$\widehat{g}_{t}=\textsc{regression}^{(X)}\left(X_{t}\sim(\overline{X}_{t-1},\overline{A}_{t-1})\right)$
$\triangleright$ fit regression model with observed data

$\widetilde{X}_{t}^{(i)}=X_{t}^{(i)}-\widehat{g}_{t}\left(\overline{X}_{t-1}^{(i)},\overline{A}_{t-1}^{(i)}\right)$
for $i=1,\dots,n$ $\triangleright$ compute residuals

end for

$\widehat{f}=\textsc{regression}^{(Y)}\left(Y\sim(\widetilde{X}_{1},A_{1},\dots,\widetilde{X}_{T},A_{T})\right)$
$\triangleright$ orthogonalized regression fit

Output:  $\widehat{f}$ $\triangleright$ $A_{t}$ coefficients in the
orthogonalized model are causal estimates

##  3 The causal estimand of orthogonalized regression

The central result of this work is that regression coefficients from Algorithm
1 are valid estimates of causal effects. In this section, we make this precise
for a variety of cases: linear models, log-linear models, probit regression
and Cox proportional hazard model. In each case, given certain conditions on
the time-varying covariates $\overline{X}$, we show that the regression
coefficients from Algorithm 1 correspond to the plug-in estimator from the
g-formula in (3). This implies that the estimates are consistent for a causal
parameter. We will now make explicit for each case.

###  3.1 Warm-up: orthogonalization in the Gaussian linear model

As an initial simple example, we return to the linear structural causal model
from Section 2.2, treating this case in detail. However, we will now only
assume that the variables are jointly Gaussian with graph given in the left
side of Figure 2. Our goal is to estimate the total causal effect (1) when
setting the treatment $\overline{A}$ to some value $\overline{a}$. Since our
main goal is to show that our estimate from Algorithm 1 is consistent for the
true causal effect, for now we will work in the population (infinite data)
setting.

We seek to estimate the causal effect of $(A_{1},A_{2})$ on $Y$, which is
given by the g-formula in (3) with

$\mu(a_{1},x_{2},a_{2})=\mbox{$\mathbb{E}$}[Y|A_{1}=a_{1},X_{2}=x_{2},A_{2}=a_{2}]=\beta^{\mathrm{(naive)}}_{1}a_{1}+\lambda^{\mathrm{(naive)}}_{2}x_{2}+\beta^{\mathrm{(naive)}}_{2}a_{2},$

where

$\left(\beta^{\mathrm{(naive)}}_{1},\lambda^{\mathrm{(naive)}}_{2},\beta^{\mathrm{(naive)}}_{2}\right)=\operatorname*{argmin}_{(\beta_{1},\lambda_{2},\beta_{2})}\mbox{$\mathbb{E}$}\left(Y-\beta_{1}A_{1}-\lambda_{2}X_{2}-\beta_{2}A_{2}\right)^{2}$
(4)

are the population regression coefficients of $Y$ onto $(A_{1},X_{2},A_{2})$.
Plugging this into the g-formula in (3) and simplifying yields

$\psi(a_{1},a_{2})=\underbrace{\left(\beta^{\mathrm{(naive)}}_{1}+\frac{\textnormal{Cov}(A_{1},X_{2})}{\textnormal{Var}(A_{1})}\lambda^{\mathrm{(naive)}}_{2}\right)}_{\theta_{1}}a_{1}+\underbrace{\beta^{\mathrm{(naive)}}_{2}}_{\theta_{2}}a_{2}.$

This calculation shows that $\psi(a_{1},a_{2})$ is linear in $(a_{1},a_{2})$,
with coefficients $(\theta_{1},\theta_{2})$ as the causal effects of $A_{1}$
and $A_{2}$ on $Y$, respectively. Because
$\theta_{1}\neq\beta^{\mathrm{(naive)}}_{1}$, this calculation also shows that
the naive regression coefficient of $A_{1}$ is not consistent for the causal
effect of $A_{1}$. (In this example, naive regression does yield a consistent
estimate of the causal effect of $A_{2}$; this is because there is no collider
path from $A_{2}$ to $Y$.)

Next, we show that orthogonalized regression in Algorithm 1 yields consistent
estimates of $(\theta_{1},\theta_{2})$. The method consists of regressing $Y$
on $A_{1}$, $A_{2}$ and the _orthogonalized_ version of $X_{2}$:

$\widetilde{X}_{2}:=X_{2}-\frac{\textnormal{Cov}(A_{1},X_{2})}{\textnormal{Var}(A_{1})}A_{1},$
(5)

yielding regression coefficient estimates for $(A_{1},A_{2})$:

$({\beta}_{1}^{\mathrm{(ortho)}},{\beta}^{\mathrm{(ortho)}}_{2})=\operatorname*{argmin}_{\beta_{1},\beta_{2}}\min_{\lambda}\mbox{$\mathbb{E}$}\left(Y-\beta_{1}A_{1}-\beta_{2}A_{2}-\lambda\widetilde{X}_{2}\right)^{2}.$
(6)

The intuition is that it is safe to include $\widetilde{X}_{2}$ in the
regression model because it does not include any collider bias:
$\widetilde{X}_{2}$ is independent of $A_{1}$. Our key result is that this is
a consistent estimator for the causal effect:

######  Proposition 1 (Orthogonalized regression for Gaussian linear models)

In the setting above, $\psi$ is linear in $(a_{1},a_{2})$:

$\psi(a_{1},a_{2})=\theta_{1}a_{1}+\theta_{2}a_{2}.$

Moreover, the orthogonalized regression coefficients correspond to the causal
parameters:

$({\beta}_{1}^{\mathrm{(ortho)}},{\beta}_{2}^{\mathrm{(ortho)}})=(\theta_{1},\theta_{2}).$

We include the proof of the preceding proposition because it is
straightforward and illustrates the key idea.

Proof. Starting with (6) and the definition of $\widetilde{X}_{2}$, we have

$({\beta}_{1}^{\mathrm{(ortho)}},{\beta}^{\mathrm{(ortho)}}_{2})=\operatorname*{argmin}_{\beta_{1},\beta_{2}}\min_{\lambda}\mbox{$\mathbb{E}$}\left(Y-\beta_{1}A_{1}-\beta_{2}A_{2}-\left[X_{2}-\frac{\textnormal{Var}(A_{1})}{\textnormal{Cov}(A_{1},X_{2})}A_{1}\right]\lambda\right)^{2}.$
(7)

The optimization problems in (7) and (4) are reparameterizations of each other
so they take the same minimum value. We will thus obtain a set of parameters
for (7) that achieves the optimal value in (4), which proves that they are
also minimizers of (7).

First set
${\beta}^{\mathrm{(ortho)}}_{2}=\beta^{\mathrm{(naive)}}_{3}=\theta_{2}$ and
rewrite

$\displaystyle{\beta}_{1}^{\mathrm{(ortho)}}=\operatorname*{argmin}_{\beta_{1}}\min_{\lambda}\mbox{$\mathbb{E}$}\left(Y-\left[\beta_{1}-\frac{\textnormal{Var}(A_{1})}{\textnormal{Cov}(A_{1},X_{2})}\lambda\right]A_{1}-\theta_{2}A_{2}-X_{2}\lambda\right)^{2},$

from which we see that the minimizer for $\lambda$ above is
$\lambda^{\mathrm{(naive)}}_{2}$, and thus

${\beta}^{\mathrm{(ortho)}}_{1}-\frac{\textnormal{Var}(A_{1})}{\textnormal{Cov}(A_{1},X_{2})}\lambda^{\mathrm{(naive)}}_{2}=\beta^{\mathrm{(naive)}}_{1}.$

$\Box$

###  3.2 Linear model without Gaussianity

We next generalize the estimator from the previous to the case without the
Gaussian assumption, and consider an arbitrary number of time steps $T$.

Let $\widetilde{X}=(\widetilde{X}_{1},\dots,\widetilde{X}_{T})$ be the
orthogonalized covariates where

$\widetilde{X}_{t}:=X_{t}-\mbox{$\mathbb{E}$}[X_{t}\mid\overline{X}_{t-1},\overline{A}_{t-1}].$
(8)

From here, we assume that the mean of $Y$ is linear in the treatment and
orthogonalized covariates:

$\mbox{$\mathbb{E}$}[Y\mid\overline{X},\overline{A}]=\alpha+\beta^{\top}\overline{A}+\lambda^{\top}\widetilde{X}.$
(9)

We emphasize that we are not assuming anything about the distribution of
$(\overline{X},\overline{A})$ except that the conditional mean of $X_{t}$
given the past exists.

In the setting above, it is again the case that a regression fit of $Y$ onto
$(\widetilde{X},\overline{A})$ corresponds to a g-formula plug-in estimator.
Specifically, we consider the orthogonalized regression fit

$\beta^{\mathrm{(ortho)}}=\operatorname*{argmin}_{\beta}\min_{\alpha,\lambda}\mbox{$\mathbb{E}$}\left(Y-\alpha-\beta^{\top}\overline{A}-\lambda^{\top}\widetilde{X}\right)^{2}.$
(10)

These are causal estimates, as stated formally next.

######  Proposition 2 (Orthogonalized regression for linear models)

In the setting above, $\psi(\overline{a})$ is linear:

$\psi(\overline{a})=\theta_{0}+\theta^{\top}a,$

and the orthogonalized regression coefficients in (10) are the causal
parameters: $\beta^{\mathrm{(ortho)}}=\theta$.

Notice again that this is a population level result. In practice, we do not
observe $\widetilde{X}$, so Algorithm 1 instead estimates the distribution of
$X_{t}\mid\overline{X}_{t-1},\overline{A}_{t-1}$ and uses this to obtain an
approximation of $\widetilde{X}_{t}$.

###### Remark 3

The result above for linear models is known to some researchers (personal
correspondence). Nonetheless, we are not aware of any prior writing similar to
our exposition.

###  3.3 Modeling with orthogonalized variables

We pause here to discuss the interpretation of the outcome model assumed in
(9), where the mean of $Y$ is assumed to be linear in the variables
$\widetilde{X}_{t}$ and $A_{t}$ for $t=1,\dots,T$. How does this relate to the
familiar assumption that $Y$ is linear in $X_{t}$ and $A_{t}$? In fact, these
two outcome models are the same in the most typical case: when $X_{t}$ is
modeled as depending linearly on $(\overline{X}_{t-1},\overline{A}_{t-1})$.
That is, the assumption in (9) is then equivalent to

$\mbox{$\mathbb{E}$}[Y\mid\overline{X},\overline{A}]=\alpha+\beta^{\top}\overline{A}+\lambda^{\top}\overline{X}.$

Although cosmetically different, our outcome model for orthogonality
regression corresponds to a routine setup. This continues to be the case in
the nonlinear regression models that we turn to next.

###  3.4 Log-linear model

As a first extension of orthogonalization to nonlinear models, we consider
log-linear models. As before, we define the orthogonalized covariate
$\widetilde{X}$ as in (8). Additionally, we assume that it follows a location
model: $X_{t}\stackrel{{\scriptstyle
d}}{{=}}g_{t}(\overline{a}_{t-1},\overline{X}_{t-1})+\epsilon_{t}$ for some
function $g_{t}$ and independent random variable $\epsilon_{t}$. Lastly, we
assume the log of the mean of $Y$ is linear in the treatment and
orthogonalized covariates:

$\mbox{$\mathbb{E}$}[Y\mid\widetilde{X},\overline{A}]=\exp\\{\alpha+\beta^{\top}\overline{A}+\lambda^{\top}\widetilde{X}\\}.$
(11)

In this setting, regression coefficients from the log-linear regression in
(11) are causal estimates, as stated next.

######  Proposition 4 (Orthogonalized regression for log-linear models)

In the setting above, $\psi$ is log-linear:

$\psi(\overline{a})=\exp\\{\alpha^{\prime}+\theta^{\top}\overline{a}\\},$

for some $\alpha^{\prime}\in\mathbb{R}$. Moreover, $\theta$ is equal to
$\beta$ in the full log-linear model in (11).

That is, as with linear models, the g-formula simplifies such that the causal
parameters coincide with the regression coefficients.

###  3.5 Gaussian-probit model

As a second foray into nonlinear models, we consider probit models for a
binary response $Y$. As before, the orthogonalized regression coefficients
from Algorithm 1 will be valid causal estimates, although the story is
slightly more complex in this model, which will be explained soon. As before,
we define the orthogonalized covariate $\widetilde{X}$ as in (8). We now
additionally assume that that the covariates follow a Gaussian location model:
$X_{t}\stackrel{{\scriptstyle
d}}{{=}}g_{t}(\overline{a}_{t-1},\overline{X}_{t-1})+\epsilon_{t}$ for some
function $g_{t}$ and independent random variable
$\epsilon_{t}\sim\mathcal{N}(0,\sigma_{t}^{2})$. Lastly, we assume that the
binary response $Y$ follows a probit model:

$P\left(Y=1\mid\overline{A},\widetilde{X}\right)=\Phi\left(\alpha+\beta^{\top}\overline{A}+\lambda^{\top}\widetilde{X}\right),$
(12)

where $\Phi$ is the CDF of the standard normal distribution. Similar to
before, the g-formula evaluates to a probit model on $\overline{A}$, but now
with re-scaled coefficients, as stated next.

######  Proposition 5 (Orthogonalized regression for probit-Gaussian models)

In the setting above, $\psi$ follows a probit model:

$\psi(\overline{a})=\Phi(\alpha^{\prime}+\theta^{\top}\overline{a}).$

Moreover, $\theta=\beta/\sqrt{1+\sigma^{2}}$, where $\beta$ are the regression
coefficients from (11) and
$\sigma^{2}=\lambda_{1}^{2}\sigma_{1}^{2}+\dots+\lambda_{T}^{2}\sigma_{T}^{2}$
with $\lambda$ from (12).

##### Why is there re-scaling?

One may wonder why the probit regression coefficients need to be re-scaled to
recover the $g$-formula. The intuition is that with a probit regression (or a
logistic regression) the magnitude of the coefficient vector is encoding the
certainty/uncertainty. Larger magnitude corresponds to higher confidence. In
the marginal model, we have less certainty, so the magnitude of the
coefficients shrink. The direction of the coefficient vector remains
unchanged, which means that a coefficient in the marignal model is zero only
when the corresponding coefficient in the full model is zero.

##### Beyond the probit model

The calculation in this seemingly-specific model class is more general than it
first appears, in that the results will approximately hold even when the
regression in (12) is not a probit regression, and when the residuals are not
exactly Gaussian. The reason is the central limit theorem; with many time
points, we can expect $\sum_{t}\lambda_{t}\widetilde{X}_{t}$ to be
approximately Gaussian. As an example, consider logistic regression. Logistic
regression is a _threshold model_ : it can be thought of as the probability of
a logistic random variable falling below a threshold:

$P(Y=1\mid\overline{A},\overline{X})=P(Z\leq\overline{A}\beta+\overline{X}\lambda^{\prime}),$

where $Z$ is a logistic random variable. After orthogonalization and an
application of the g-formula, we obtain

$\psi(\overline{a})=P(Z-\sum_{t}\lambda_{t}\widetilde{X}_{t}\leq\overline{a}\beta)\approx\Phi(\overline{a}\beta/\sigma^{2}),$

where $\sigma^{2}$ is the variance of
$Z-\sum_{t}\lambda_{t}\widetilde{X}_{t}$. Thus, this marginalized logistic
model approaches a probit model as the variance of
$\sum_{t}\lambda_{t}\widetilde{X}_{t}$ grows large. Roughly speaking,
threshold models like logistic regression converge to probit models when
marginalizing out many time-varying covariates.

###  3.6 Cox proportional hazards model

Lastly, we consider the Cox propotional hazards model (Cox,, 1972), the most
commonly used model for modeling survival times. As before, we define the
orthogonalized covariate $\widetilde{X}$ as in (8). We additionally assume
that the the orthogonalized covariates follow a Gaussian location model. That
is, we are assuming $X_{t}\stackrel{{\scriptstyle
d}}{{=}}g_{t}(\overline{a}_{t-1},\overline{X}_{t-1})+\epsilon_{t}$ for some
function $g_{t}$ and independent random variable
$\epsilon_{t}\sim\mathcal{N}(0,\sigma_{t}^{2})$. Lastly, we assume that the
survival time $Y$ follows a proportional hazard model:

$h(t\mid\overline{A},\overline{X})=h_{0}(t)\cdot\exp\\{\beta^{\top}\overline{A}+\lambda^{\top}\widetilde{X}\\},$
(13)

where $h_{0}$ is a unknown baseline hazard rate. Recall that this implies that
the survival probabilities are

$P(Y\geq
t\mid\overline{A},\overline{X})=\exp\left\\{-H_{0}(t)\exp\\{\beta^{\top}\overline{A}+\lambda^{\top}\widetilde{X}\\}\right\\},$

where

$H_{0}(t)=\int_{0}^{t}h_{0}(t^{\prime})dt^{\prime}$

is the baseline cumulative hazard function.

We next apply the g-formula to show that the regression coefficients from
fitting the regression model above have a causal meaning. In particular, we
will show that the implied g-formula has the form of a Cox model with an
additional random effect term, and coefficients matching those from the
regression in (13), as stated next.

######  Proposition 6 (Orthogonalized Cox regression)

In the setting above, $\psi$ can be expressed as follows:

$\psi(\overline{a},t):=P(Y^{\overline{a}}\geq
t)=\int\exp\left\\{-H_{0}(t)\exp\\{\theta^{\top}\overline{A}+Z\\}\right\\}P(Z=z)dx,$
(14)

where $Z\sim\mathcal{N}(0,\sigma^{2})$. Moreover, the causal parameters
$\theta$ are identical to the coefficients $\beta$ from the full model in
(13), and $\sigma^{2}=\sum_{t}\lambda_{t}^{2}\sigma_{t}^{2}$.

The model arising from (14) is a cousin of the usual Cox model; it is a Cox
model with an additional random effects term that is marginalized out. Put
another way, notice that this is the survival curve for a population comprised
of individuals with (Cox model) hazard functions

$h(t\mid\overline{A},Z)=h_{0}(t)\exp\\{\beta^{\top}\overline{A}+Z\\},$

where $Z\sim\mathcal{N}(0,\sigma^{2})$ is a random effect varying across each
member of the population. Thus, the regression from the full Cox model in (13)
has a sensible causal interpretation: it is the causal effect in a random-
effects Cox regression model. Note however, that this is not the same estimand
as the usual marginal structural Cox model (Hernán et al.,, 2000).

###  3.7 Confidence intervals

All of the models that we are fitting are regular, finite dimensional
parametric models. The orthogonalized least squares estimators are
$M$-estimators. It follows from Theorem 5.21 of Van der Vaart, (2000) that the
estimates are asymptotically normal, converging at rate $n^{-1/2}$. The
bootstrap provides valid asymptotic confidence intervals by Theorem 23.7 of
Van der Vaart, (2000).

##  4 Numerical experiments

We next investigate the orthogonalization estimator numerically for the linear
model and the Cox regression model. In both cases, we investigate two aspects
of performance. First, we demonstrate that they are valid: they produce
estimated effects of zero (up to sampling variability) under the null. Second,
we compare the efficiency of the orthogonalization estimator to the IPW
estimator. Code reproducing these experiments is available at
https://github.com/stephenbates19/causal_orthogonalization.

###  4.1 Linear model

####  4.1.1 Validity under the null

We consider the linear regression setting of Section 3.2 with $T=5$ time
points. We simulate binary treatment assignments $A_{t}$ and covariates
$X_{t}$ whose mean depends on $A_{t-1}$. Moreover, the treatment assignment
$A_{t}$ depends on $X_{t}$, and there is a phantom $U$, so there is possible
collider bias. Lastly, the response $Y$ is heavy-tailed and asymmetric. There
is no causal effect from $A_{t}$ on $Y$. We treat the causal effect from
$A_{t-1}$ to $X_{t}$ as a control parameter, investigating the behavior across
a range of values. In this experiment, we take a sample of 1 million data
points, so the results are approximately the results of the population
regression.

The left panel of Figure 3, shows the results of a regression of $Y$ onto
$(X_{t},A_{t})_{t=1}^{5}$ (denoted _naive_) with an regression of $Y$ onto
$(\widetilde{X}_{t},A_{t})_{t=1}^{5}$ (denoted _ortho_). In particular, we
report the magnitude of the largest coefficient of $A_{1},\dots,A_{5}$ in each
setting. Because there is no causal effect, the ground-truth value is zero.
However, we observe that the naive regression produces large coefficients when
there is a large effect size from the treatment to the covariate. This is
caused by collider bias, and it was anticipated by our analytic calculations
in the Gaussian case. By contrast, the orthogonalized regression produces
coefficients that are nearly zero (they are not exactly zero because we take
the maximum of five numbers with Monte Carlo variability). Thus, the
orthogonalized regression is producing valid causal estimates in this case.

####  4.1.2 Evaluating efficiency

Continuing with this setting, we next investigate the efficiency of
orthogonalized regression compared to fitting a MSM with IPW; see van der Wal
and Geskus, (2011) for more details of the latter approach. We consider two
cases. In the first case, we have high overlap: $A_{t}$ is independent of
$X_{t-1}$, so the propensity score is always the same. In the second case, we
have poor overlap: $A_{t}$ is equal to the sign of $X_{t-1}$ with probability
95%. In this case, the IPW weights are more extreme. Lastly, we vary the
strength of the effect from $U$ to $X_{t}$ as a control parameter.

We report the results in the right panel of Figure 3. We find that as the
effect from $U$ to $X_{t}$ increases in strength, orthogonalized regression
outperforms IPW by a larger margin. This makes sense; orthogonalized
regression is able to regress out the variation from $X_{t}$ (and by
extension, part of $U$), but IPW is not able to regress this out. Thus, for
larger effects from $U$ to $X_{t}$, orthogonalized regression will have
correspondingly more precise estimates. Secondly, the performance of IPW is
worse when there is less overlap. This is a well-known problem with IPW; when
there is poor overlap the weights become extreme, leading to a small effective
sample size and less efficiency. We emphasize that the overlap in this case is
not implausible; at each step any treatment assignment can be seen with
probability at least 5%. Nonetheless, IPW has poor relative efficiency.

(a) Consistency

(b) Efficiency compared to IPW

Figure 3: Results of the linear, non-Gaussian simulation experiment. (a)
Estimates of the causal effect. The true effect is zero. (b) Relative variance
of IPW compared to orthogonalized regression in high and low overlap
scenarios.

###  4.2 Cox proportional hazards model

####  4.2.1 Validity under the null

We next demonstrate orthogonalization with the Cox model with a synthetic
example across 5 time points. The treatment assignment is binary, and is
assigned based on the covariate at the previous time point. The covariate is a
location model based on the phantom, past treatments, and independent Gaussian
noise. The phantom in this case is the baseline hazard rate for each subject,
which spans a factor of 3 in the population. (That is, the unhealthiest person
has a hazard rate three times larger than the healthiest person.) We examine
the null case where there is no treatment effect of $A_{t}$ on $Y$. We treat
the causal effect from $A_{t-1}$ to $X_{t}$ as a control parameter—larger
values correspond to more collider bias.

We compare two approaches in this setting. First, we fit a Cox model with four
covariates: $X_{t}$, $A_{t}$, the cumulative average of $X_{t}$, and the
cumulative average of $A_{t}$, calling this approach _naive_. We will report
the coefficient of the cumulative average of $A_{t}$. Since this a null case,
these should be zero. We compare this with a Cox model with covariates
$X_{t}$, $A_{t}$, the cumulative average of $\widetilde{X}_{t}$, and the
cumulative average of $A_{t}$, _ortho_. In this case, $\widetilde{X}_{t}$ is
computed by simply taking the residuals of a linear regression fit of $X_{t}$
on $A_{t-1}$.

We report the results in Figure 4, showing the average of $250$ runs with
$n=2,000$ subjects. We find that the naive regression model leads to nonzero
coefficients, even under the null. The bias gets worse as the size of the
causal effect from $A_{t-1}$ to $X_{t}$ increases, which is anticipated by our
theoretical results in the linear Gaussian model and previous simulations in
the (non-Gaussian) linear simulation. By contrast, the orthogonalized
regression yields an estimate near across all regimes.

Figure 4: Results of the Cox model experiment. We plot the regression
coefficient corresponding to the cumulative average of past treatment
assignments in the naive Cox regression model and in the orthogonalized Cox
regression model.

####  4.2.2 Evaluating power

Next, we consider the same data generating mechanism as in the previous
example, but now with a treatment effect. In particular, $A_{t}=1$ (i.e.,
receiving treatment) reduces the risk of an event during that next time step
by constant factor. In this case, our goal is to detect if the treatment
effect is present, so we will check whether the confidence intervals for the
coefficient for $A_{t}$ include zero. (Recall that the estimand for these two
methods is different, but both are zero under the null. Therefore, checking
power provides us a fair way to compare the efficiency.) We consider IPW as
implemented in van der Wal and Geskus, (2011) and orthogonalized regression,
but but not the naive method, since we know it is not a valid estimator of the
causal effect. In each trial, we simulated $n=2,000$ individuals and record
the resulting coefficient for $A_{t}$. We do this for many generated random
data sets, and report the coefficients after standardizing them to have
variance one (i.e., a z-score for the fitted coefficient).

Figure 5: The z-score for the causal effect with two valid regression methods.
The three panels correspond to the strength of the effect from $X_{t}$ to
$A_{t}$ (higher means less overlap).

We vary the strength of the treatment effect size and the strength of the
effect of $X_{t}$ on $A_{t}$ as two control parameters. The first is measured
as the fraction of reduced risk, and the second varies across three levels,
described precisely in the software accompanying this manuscript. We report
the results in Figure 5. We find that as the treatment effect increases, the
z-scores separate away from zero in the negative direction, as expected.
Turning to the sensitivity to the $X_{t}$ to $A_{t}$ effect size, we find that
orthogonalized regression is mostly unaffected by this parameter. However, as
this parameter grows, IPW becomes less effective at detecting the treatment
effect. This is expected; as the effect grows, IPW produces more extreme
weights and thus has more variability.

##  5 A realistic example

Lastly, we consider a synthetic observation HIV data set, introduced in the
ipw package to evaluate IPW-based fitting methods. The goal is to evaluate
whether a treatment regimen affects survival time. The data comprises 1200
individuals, measured for 15 time points, on average. There are two baseline
covariates, age and sex, as well as a time-varying covariate: the square root
of the CD4 cell count, which decreases as HIV progresses. In this setting,
patients initiate treatment during one of the measured time points, and then
continue the treatment until the study concludes or a death occurs.

We compare the estimate from the orthogonalized Cox regression estimator and
the IPW estimator with a Cox marginal model. For orthogonalization, we use Cox
regression with the two baseline covariates, time-varying covariate, and time-
varying treatment. For the MSM, we use a Cox regression model with the two
baseline covariates and the time-varying treatment. (Recall MSMs never include
the time-varying covariate.) Both estimates suggest a protective effect, but
the results are not statistically significant; using the bootstrap estimate of
standard error the z-score for the orthogonalized estimate is -1.37 and that
for the IPW estimate is -1.08. The basic bootstrap confidence intervals
(Davison and Hinkley,, 1997) at the 90% level are $(-1.32,.35)$ for the
orthogonalized regression estimator and $(-1.82,1.21)$ for the IPW estimator.
See Figure 6 for a visualization of the point estimates and bootstrap
replicates. A careful look at the histogram of bootstrap values indicates one
further point of concern about the IPW estimate; the bootstrap values are
biased downwards, indicating possible positive bias to the IPW point estimate.
Nonetheless, the IPW and orthogonalized regression estimates are generally
consistent with each other in this example, a good sign given the rather
different assumptions underlying the two methods.

(a) Orthogonalization

(b) IPW

Figure 6: Estimates on bootstrap subsamples of the synthetic HIV survival
data. The vertical blue line shows the point estimate from the full data.

##  6 Discussion

In this work, we have formulated a causal estimator based on regression models
with time-varying treatments and time-varying covariates. Since our estimator
is a modified regression fit, it inherits a highly-developed toolbox of
diagnostics and visualization methods. At its core, our method builds on the
connection between the causal concepts of collider bias and the classical
regression concept of collinearity. That is, we break the edges from treatment
variables to post-treatment covariates by regressing out the former. We expect
that the idea of removing collider bias by regressing out variables will be
useful in other contexts as well.

## Acknowledgements

We thank Michal Abrahamowicz, Thomas Richardson, Terry Therneau, and Lu Tian
for helpful discussions. R.T. was supported by the National Institutes of
Health (5R01 EB001988-16) and the National Science Foundation (19 DMS1208164).
S.B. was supported by the Foundations of Data Science Institute postdoctoral
fellowship.

## References

  * Bang and Robins, (2005) Bang, H. and Robins, J. M. (2005).  Doubly robust estimation in missing data and causal inference models.  Biometrics, 61(4):962–973.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal, 21(1):C1–C68.
  * Cox, (1972) Cox, D. R. (1972).  Regression models and life-tables.  Journal of the Royal Statistical Society. Series B (Methodological), 34(2):187–220.
  * Dáaz, (2019) Dáaz, I. (2019).  Statistical inference for data-adaptive doubly robust estimators with survival outcomes.  Statistics in Medicine, 38(15):2735–2748.
  * Davison and Hinkley, (1997) Davison, A. C. and Hinkley, D. V. (1997).  Boostrap methods and their applications.  Cambridge University Press, New York.
  * Drton et al., (2009) Drton, M., Eichler, M., and Richardson, T. S. (2009).  Computing maximum likelihood estimates in recursive linear models with correlated errors.  Journal of Machine Learning Research, 10:2329–2348.
  * Evans and Didelez, (2021) Evans, R. J. and Didelez, V. (2021).  Parameterizing and simulating from causal models.  arXiv preprint.  arXiv:2109.03694.
  * Frisch and Waugh, (1933) Frisch, R. and Waugh, F. V. (1933).  Partial time regressions as compared with individual trends.  Econometrica, 1(4):387–401.
  * Guo and Perković, (2021) Guo, F. R. and Perković, E. (2021).  Efficient least squares for estimating total effects under linearity and causal sufficiency.  arXiv preprint.  arXiv:2008.03481.
  * Hernán et al., (2000) Hernán, M. Á., Brumback, B., and Robins, J. M. (2000).  Marginal structural models to estimate the causal effect of zidovudine on the survival of hiv-positive men.  Epidemiology, 11(5).
  * Hernán et al., (2001) Hernán, M. A., Brumback, B., and Robins, J. M. (2001).  Marginal structural models to estimate the joint causal effect of nonrandomized treatments.  Journal of the American Statistical Association, 96(454):440–448.
  * Hernán and Robins, (2020) Hernán, M. A. and Robins, J. M. (2020).  Causal inference: what if.  Chapman & Hall/CRC, Boca Raton.
  * Keil et al., (2014) Keil, A. P., Edwards, J. K., Richardson, D. B., Naimi, A. I., and Cole, S. R. (2014).  The parametric g-formula for time-to-event data: intuition and a worked example.  Epidemiology (Cambridge, Mass.), 25(6):889–897.
  * Kennedy, (2016) Kennedy, E. H. (2016).  Semiparametric theory and empirical processes in causal inference.  In Statistical causal inferences and their applications in public health research, pages 141–167. Springer.
  * Lancaster, (2002) Lancaster, T. (2002).  Orthogonal parameters and panel data.  Review of Economic Studies, 69(3):647–666.
  * Lendle et al., (2017) Lendle, S. D., Schwab, J., Petersen, M. L., and van Der Laan, M. J. (2017).  Itmle: An R package implementing targeted minimum loss-based estimation for longitudinal data.  Journal of Statistical Software, 81(1).
  * Luedtke et al., (2017) Luedtke, A. R., Sofrygin, O., van der Laan, M. J., and Carone, M. (2017).  Sequential double robustness in right-censored longitudinal models.  arXiv preprint arXiv:1705.02459.
  * McGrath et al., (2020) McGrath, S., Lin, V., Zhang, Z., Petito, L. C., Logan, R. W., Hernán, M. A., and Young, J. G. (2020).  gformula: An r package for estimating the effects of sustained treatment strategies via the parametric g-formula.  Patterns, 1(3):100008.
  * McGrath et al., (2022) McGrath, S., Young, J. G., and Hernán, M. A. (2022).  Revisiting the g-null Paradox.  Epidemiology, 33(1):114–120.
  * Naimi et al., (2016) Naimi, A. I., Cole, S. R., and Kennedy, E. H. (2016).  An introduction to g methods.  International Journal of Epidemiology, 46(2):756–762.
  * Nandy et al., (2017) Nandy, P., Maathuis, M. H., and Richardson, T. S. (2017).  Estimating the effect of joint interventions from observational data in sparse high-dimensional settings.  Annals of Statistics, 45(2):647–674.
  * Robins, (1986) Robins, J. (1986).  A new approach to causal inference in mortality studies with a sustained exposure period—application to control of the healthy worker survivor effect.  Mathematical modelling, 7(9-12):1393–1512.
  * Robins, (2000) Robins, J. M. (2000).  Marginal structural models versus structural nested models as tools for causal inference.  In Statistical models in epidemiology, the environment, and clinical trials, pages 95–133. Springer.
  * Robins et al., (2000) Robins, J. M., Hernán, M. Á., and Brumback, B. (2000).  Marginal structural models and causal inference in epidemiology.  Epidemiology, 11(5).
  * Robins and Wasserman, (1997) Robins, J. M. and Wasserman, L. (1997).  Estimation of effects of sequential treatments by reparameterizing directed acyclic graphs.  In Proceedings of the Thirteenth Conference on Uncertainty in Artificial Intelligence, UAI’97, page 409–420, San Francisco, CA, USA. Morgan Kaufmann Publishers Inc.
  * Schomaker et al., (2019) Schomaker, M., Luque-Fernandez, M. A., Leroy, V., and Davies, M. A. (2019).  Using longitudinal targeted maximum likelihood estimation in complex settings with dynamic interventions.  Statistics in Medicine, 38(24):4888–4911.
  * Shpitser et al., (2018) Shpitser, I., Evans, R. J., and Richardson, T. S. (2018).  Acyclic linear sems obey the nested markov property.  In Uncertainty in artificial intelligence: proceedings of the Conference on Uncertainty in Artificial Intelligence, volume 2018. NIH Public Access.
  * Tian et al., (2005) Tian, L., Zucker, D., and Wei, L. J. (2005).  On the Cox model with time-varying regression coefficients.  Journal of the American Statistical Association, 100(469):172–183.
  * van der Laan and Rose, (2018) van der Laan, M. J. and Rose, S. (2018).  Targeted learning in data science: causal inference for complex longitudinal studies.  Springer.
  * Van der Vaart, (2000) Van der Vaart, A. W. (2000).  Asymptotic statistics.  Cambridge university press.
  * van der Wal and Geskus, (2011) van der Wal, W. and Geskus, R. (2011).  ipw: An r package for inverse probability weighting.  Journal of Statistical Software, 43(13):1–23.
  * Wang et al., (2017) Wang, C., Liu, H., and Gao, S. (2017).  A penalized Cox proportional hazards model with multiple time-varying exposures.  Annals of Applied Statistics, 11(1):185–201.
  * Wen et al., (2021) Wen, L., Young, J. G., Robins, J. M., and Hernán, M. A. (2021).  Parametric g-formula implementations for causal survival analyses.  Biometrics, 77(2):740–753.
  * Xiao et al., (2014) Xiao, Y., Abrahamowicz, M., Moodie, E. E., Weber, R., and Young, J. (2014).  Flexible marginal structural models for estimating the cumulative effect of a time-dependent treatment on the hazard: Reassessing the cardiovascular risks of didanosine treatment in the Swiss HIV Cohort study.  Journal of the American Statistical Association, 109(506):455–464.
  * Young and Tchetgen Tchetgen, (2014) Young, J. G. and Tchetgen Tchetgen, E. J. (2014).  Simulation from a known Cox MSM using standard parametric models for the g-formula.  Statistics in medicine, 33(6):1001–1014.
  * Zheng et al., (2016) Zheng, W., Petersen, M., and van der Laan, M. J. (2016).  Doubly Robust and Efficient Estimation of Marginal Structural Models for the Hazard Function.  International Journal of Biostatistics, 12(1):233–252.

##  Appendix A Computation

We next discuss the computational cost of the orthogonalization estimator, and
compare it to that of other methods. The orthogonalization estimator requires
$T+1$ model fits: a model for $X_{t}$ given the past for $t=1,\dots,T$ and the
final regression model with the orthogonalized covariates. Note, however, that
if $X_{t}$ is vector-valued with $d$ coordinates, then fitting the models for
$X_{t}$ given the past may be more involved than a model with a scalar outcome
variable. In the extreme case, one could fit each coordinate seperately, which
would then amount to $d\cdot T+1$ total model fits. However, regression models
that accommodate a vector-valued responses can be used to keep the total
number of model fits at $T+1$, if necessary. The IPW estimator for the MSM
also requires $T+1$ model fits.

We next compare this to the the computation required by the iterative
conditional expectations estimator (Hernán and Robins,, 2020) and the
efficient estimator of Bang and Robins, (2005). These two estimators require
order $T$ model fits, similar to the orthogonalization method. However, these
two methods only estimate $\phi(\overline{a})$ for a single pre-specified
$\overline{a}$; to estimate $\phi(\overline{a})$ for two different values of
$\overline{a}$, two different sequences of models must be fit. By contrast,
orthogonalization and IPW-fitted MSMs give an estimate for
$\phi(\overline{a})$ for all $\overline{a}$ in one pass, making them much more
efficient when the analyst is interested in understanding the effect of many
settings of $\overline{a}$.

##  Appendix B Identification conditions

The causal parameter $\psi(\overline{a}_{t})$ is identified and is given by
the $g$-formula under three conditions:

($\mathghost 1$) No interference. $Y_{t}(\overline{a}_{t})=Y_{t}$ when
$\overline{A}_{t}=\overline{a}_{t}$. This means that someone’s outcome depends
on their treatment but not other subjects’ treatments.

($\mathghost 2$) Positivity. There is an $\epsilon>0$ such that
$\pi(a_{t}|h_{t})\geq\epsilon$ for all $a_{t}$ and $h_{t}$ where $h_{t}$
denotes all variables that occur before $A_{t}$.

($\mathghost 3$) No unmeasured confounding. For every $t$, $A_{t}$ is
independent of $Y_{t}(\overline{a}_{t})$ given all past variables $H_{t}$.
This means that the data behave like a randomized experiment where the
randomization probabilities can depend on the past observed variables.

##  Appendix C Proofs

Proof. (Proposition 2) Let
$g_{t}(\overline{a}_{t-1},\overline{x}_{t-1}))=\mbox{$\mathbb{E}$}[X_{t}\mid\overline{X}_{t-1},\overline{A}_{t-1}]$.
We start from the g-formula and then manipulate the expression as follows:

$\displaystyle\psi(\overline{a})$
$\displaystyle=\int\dots\int\left[\alpha+\beta^{\top}\overline{a}+\sum_{t=1}^{T}\lambda_{t}(x_{t}-g_{t}(\overline{a}_{t-1},\overline{x}_{t-1}))\right]\
\prod_{t}p\left(X_{t}=x_{t}|\overline{A}_{t-1}=\overline{a}_{t-1},\overline{x}_{t-1}\right)dx_{1}\dots
dx_{T}$ $\displaystyle=\alpha+\beta^{\top}\overline{a}.$

$\Box$

Proof. (Proposition 4) We start from the g-formula and then manipulate the
expression as follows:

$\displaystyle\psi(\overline{a})$
$\displaystyle=\int\dots\int\exp\\{\alpha+\beta^{\top}\overline{a}+\sum_{t}\lambda_{t}(x_{t}-g_{t}(\overline{a}_{t-1},\overline{x}_{t-1}))\\}\
\prod_{t}p\left(X_{t}=x_{t}|\overline{A}_{t-1}=\overline{a}_{t-1},\overline{x}_{t-1}\right)dx_{1}\dots
dx_{T}$
$\displaystyle=\int\dots\int\exp\\{\alpha+\beta^{\top}\overline{a}+\sum_{t}\lambda_{t}x_{t}\\}\
\prod_{t}p\left(X_{t}=x_{t}+g_{t}(\overline{a}_{t-1},\overline{x}_{t-1})|\overline{A}_{t-1}=\overline{a}_{t-1},\overline{x}_{t-1}\right)dx_{1}\dots
dx_{T}$
$\displaystyle=\int\dots\int\exp\\{\alpha+\beta^{\top}\overline{a}+\sum_{t}\lambda_{t}x_{t}\\}\
\prod_{t}p\left(\epsilon_{t}=x_{t}\right)dx_{1}\dots dx_{T}$
$\displaystyle=\exp\\{\alpha+\beta^{\top}\overline{a}\\}\int\dots\int\exp\\{\sum_{t}\lambda_{t}x_{t}\\}\
\prod_{t}p\left(\epsilon_{t}=x_{t}\right)dx_{1}\dots dx_{T}$
$\displaystyle=\exp\\{\alpha^{\prime}+\beta^{\top}\overline{a}\\}$

for some new intercept $\alpha^{\prime}$. $\Box$

Proof. (Proposition 5) As usual, we will use the g-formula to derive and
expression for $\psi$, and we will relate this to the regression coefficients
from (12).

$\displaystyle\psi(\overline{a})$
$\displaystyle=\int\dots\int\Phi\left(\alpha+\beta^{\top}\overline{a}+\sum_{t}\lambda_{t}(x_{t}-g_{t}(\overline{a}_{t-1},\overline{x}_{t-1}))\right)\
\prod_{t}p\left(X_{t}=x_{t}|\overline{A}_{t-1}=\overline{a}_{t-1},\overline{x}_{t-1}\right)dx_{1}\dots
dx_{T}$
$\displaystyle=\int\dots\int\Phi\left(\alpha+\beta^{\top}\overline{a}+\lambda^{\top}\overline{x}\right)\
\prod_{t}p\left(X_{t}=x_{t}+g_{t}(\overline{a}_{t-1},\overline{x}_{t-1})|\overline{A}_{t-1}=\overline{a}_{t-1},\overline{x}_{t-1}\right)dx_{1}\dots
dx_{T}$
$\displaystyle=\int\dots\int\Phi\left(\alpha+\beta^{\top}\overline{a}+\lambda^{\top}\overline{x}\right)\
\prod_{t}p\left(\epsilon_{t}=x_{t}\right)dx_{1}\dots dx_{T}$
$\displaystyle=\int\dots\int
P(Z\leq\alpha+\beta^{\top}\overline{a}+\lambda^{\top}\overline{x})\prod_{t}p\left(\epsilon_{t}=x_{t}\right)dx_{1}\dots
dx_{T}$ $\displaystyle=P(Z-Z^{\prime}\leq\alpha+\beta^{\top}\overline{a})$
$\displaystyle=P(\sqrt{1+\sigma^{2}}Z\leq\alpha+\beta^{\top}\overline{a})$

where $Z\sim\mathcal{N}(0,1)$ and $Z^{\prime}\sim\mathcal{N}(0,\sigma^{2})$,
independently. $\Box$

Proof. (Proposition 6) We apply the g-formula to the indicator of the event
$\\{Y\geq t\\}$ for each $t$ to compute the survival times of the marginal
model:

$\displaystyle\psi(\overline{a},t)$
$\displaystyle=\int\dots\int\exp\left\\{-H_{0}(t)\exp\\{\beta^{\top}\overline{a}+\sum_{t}\lambda_{t}(x_{t}-g_{t}(\overline{a}_{t-1},\overline{x}_{t-1}))\\}\right\\}\cdot$
$\displaystyle\qquad\qquad\prod_{t}p\left(X_{t}=x_{t}|\overline{A}_{t-1}=\overline{a}_{t-1},\overline{x}_{t-1}\right)dx_{1}\dots
dx_{T}$
$\displaystyle=\int\dots\int\exp\left\\{-H_{0}(t)\exp\\{\beta^{\top}\overline{a}+\sum_{t}\lambda_{t}x_{t}\\}\right\\}\cdot$
$\displaystyle\qquad\qquad\prod_{t}p\left(X_{t}=x_{t}+g_{t}(\overline{a}_{t-1},\overline{x}_{t-1})|\overline{A}_{t-1}=\overline{a}_{t-1},\overline{x}_{t-1}\right)dx_{1}\dots
dx_{T}$
$\displaystyle=\int\dots\int\exp\left\\{-H_{0}(t)\exp\\{\beta^{\top}\overline{a}+\sum_{t}\lambda_{t}x_{t}\\}\right\\}\
\prod_{t}p\left(\epsilon_{t}=x_{t}\right)dx_{1}\dots dx_{T}$
$\displaystyle=\int\exp\left\\{-H_{0}(t)\exp\\{\beta^{\top}\overline{a}+Z\\}\right\\}\
p(Z=z)dz$

$\Box$