%% BioMed_Central_Tex_Template_v1.06
%%                                      %
%  bmc_article.tex            ver: 1.06 %
%                                       %

%%IMPORTANT: do not delete the first line of this template
%%It must be present to enable the BMC Submission system to
%%recognise this template!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                     %%
%%  LaTeX template for BioMed Central  %%
%%     journal article submissions     %%
%%                                     %%
%%          <8 June 2012>              %%
%%                                     %%
%%                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% For instructions on how to fill out this Tex template           %%
%% document please refer to Readme.html and the instructions for   %%
%% authors page on the biomed central website                      %%
%% http://www.biomedcentral.com/info/authors/                      %%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%% BioMed Central currently use the MikTex distribution of         %%
%% TeX for Windows) of TeX and LaTeX.  This is available from      %%
%% http://www.miktex.org                                           %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% additional documentclass options:
%  [doublespacing]
%  [linenumbers]   - put the line numbers on margins

%%% loading packages, author definitions

%\documentclass[twocolumn]{bmcart}% uncomment this for twocolumn layout and comment line below
\documentclass{bmcart}
\usepackage[utf8]{inputenc}
\usepackage{tabularx}
\usepackage{apacite}
%\usepackage[numbers]{natbib}
\usepackage{natbib} % generate Lee et al. (2010) format
\usepackage{url}
\usepackage{amsmath}
\usepackage{latexsym}
\usepackage{amssymb}
\usepackage{graphicx}
\makeatletter
\newcommand*{\rom}[1]{\expandafter\@slowromancap\romannumeral #1@}
\makeatother
\usepackage{xcolor}
\usepackage{todonotes}
\usepackage{soul}
\usepackage{microtype}
\newcommand{\tabincell}[2]{\begin{tabular}{@{}#1@{}}#2\end{tabular}}  
\newcommand\independent{\protect\mathpalette{\protect\independenT}{\perp}}
\def\independenT#1#2{\mathrel{\setbox0\hbox{$#1#2$}%
\copy0\kern-\wd0\mkern4mu\box0}} 
\usepackage{threeparttable}
\usepackage{pifont}% http://ctan.org/pkg/pifont
\newcommand{\cmark}{\ding{51}}%
\usepackage{float}
\usepackage{fixltx2e}
\usepackage{makecell}
\usepackage{longtable}
\usepackage{adjustbox}
% Hide table cells:
\usepackage{array}
\newcolumntype{H}{>{\setbox0=\hbox\bgroup}c<{\egroup}@{}}
\usepackage{appendix}
\usepackage{cleveref}
\usepackage{bbm}
\usepackage{multirow}

%%% Load packages
%\usepackage{amsthm,amsmath}
%\RequirePackage{natbib}
%\RequirePackage[authoryear]{natbib}% uncomment this for author-year bibliography
%\RequirePackage{hyperref}
%\usepackage[utf8]{inputenc} %unicode support
%\usepackage[applemac]{inputenc} %applemac support if unicode package fails
%\usepackage[latin1]{inputenc} %UNIX support if unicode package fails


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                             %%
%%  If you wish to display your graphics for   %%
%%  your own use using includegraphic or       %%
%%  includegraphics, then comment out the      %%
%%  following two lines of code.               %%
%%  NB: These line *must* be included when     %%
%%  submitting to BMC.                         %%
%%  All figure files must be submitted as      %%
%%  separate graphics through the BMC          %%
%%  submission process, not included in the    %%
%%  submitted article.                         %%
%%                                             %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%\def\includegraphic{}
%\def\includegraphics{}


%%% Put your definitions there:
\startlocaldefs
\endlocaldefs


%%% Begin ...
\begin{document}

%%% Start of article front matter
\begin{frontmatter}

\begin{fmbox}
\dochead{Research}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the title of your article here     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Estimating the effect of adjuvant chemo-therapy for colon-cancer using registry data: a method comparison and validation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors here                   %%
%%                                          %%
%% Specify information, if available,       %%
%% in the form:                             %%
%%   <key>={<id1>,<id2>}                    %%
%%   <key>=                                 %%
%% Comment or delete the keys which are     %%
%% not used. Repeat \author command as much %%
%% as required.                             %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\author[
   addressref={aff1},                   % id's of addresses, e.g. {aff1,aff2}
   corref={aff1},                       % id of corresponding address, if any
   %noteref={n1},                       % id's of article notes, if any
   email={L.SHEN@uvt.nl}   % email address
]{\inits{SLJ}\fnm{Lingjie} \snm{Shen}}
\author[
   addressref={aff2},
   email={Visser.Erik@kpmg.nl}
]{\inits{EV}\fnm{Erik} \snm{Visser}}
\author[
   addressref={aff4},
   email={hans.dewilt@radboudumc.nl}
]{\inits{JHW}\fnm{Hans} \snm{de Wilt}}
\author[
   addressref={aff5},
   email={Henk.Verheul@radboudumc.nl}
]{\inits{HV}\fnm{Henk} \snm{Verheul}}
\author[
   addressref={aff3},
   email={F.vanErning@iknl.nl}
]{\inits{F.N.}\fnm{Felice} \snm{van Erning}}
\author[
   addressref={aff3},
   email={G.Geleijnse@iknl.nl}
]{\inits{GG}\fnm{Gijs} \snm{Geleijnse}}
\author[
   addressref={aff2},
   email={M.C.Kaptein@uvt.nl}
]{\inits{MCK}\fnm{Maurits} \snm{Kaptein}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter the authors' addresses here        %%
%%                                          %%
%% Repeat \address commands as much as      %%
%% required.                                %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\address[id=aff1]{%                           % unique id
  \orgname{Department of Methodology and Statistics, Tilburg University}, % university, etc
  \street{Warandelaan 2},                     %
  \postcode{5000 LE}                                % post or zip code
  \city{Tilburg},                              % city
  \cny{The Netherlands}                                    % country
}
\address[id=aff2]{%
  \orgname{Jheronimus Academy of Data Science},
  \street{ Sint Janssingel 92},
  \postcode{5211 DA}
  \city{'s-Hertogenbosch},
  \cny{The Netherlands}
}
\address[id=aff3]{%
  \orgname{IKNL (Integraal Kankercentrum Nederland)},
  \street{Godebaldkwartier 419},
  \postcode{3511 DT}
  \city{Utrecht},
  \cny{The Netherlands}
}
\address[id=aff4]{%
  \orgname{Department of Surgery, Radboud University Medical Center},
  \street{Geert Grooteplein Zuid 10},
  \postcode{6525 GA}
  \city{Nijmegen},
  \cny{The Netherlands}
}
\address[id=aff5]{%
  \orgname{Department of Medical Oncology, Radboud University Medical Center},
  \street{Geert Grooteplein Zuid 10},
  \postcode{6525 GA}
  \city{Nijmegen},
  \cny{The Netherlands}
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Enter short notes here                   %%
%%                                          %%
%% Short notes will be after addresses      %%
%% on first page.                           %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{artnotes}
%\note{Sample of title note}     % note to the article
%\note[id=n1]{Equal contributor} % note, connected to author
\end{artnotes}

\end{fmbox}% comment this for two column layout

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Abstract begins here                 %%
%%                                          %%
%% Please refer to the Instructions for     %%
%% authors on http://www.biomedcentral.com  %%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstractbox}

\begin{abstract} % abstract
\parttitle{Background} 
Although randomized controlled trials (RCT) are the gold standard to estimate treatment effects, they are often criticized in terms of generalizability. Observational data might alleviate this problem by being readily available in large quantities. However, observational data are potentially confounded. In this methodological study we use a large-scale RCT as the gold standard to examine the performance of various statistical methods to control for potential confounding in observational data.
\parttitle{Methods}
In this paper we compare three types of methods that allow researchers to correct for such potential confounding: direct methods, inverse probability weighting (IPW) methods and doubly robust (DR) methods. We uniquely compare estimates obtained from the population-wide Netherlands Cancer Registry (NCR) on colon cancer ($n=52621$) with estimates obtained from a large-scale RCT. As the RCT differs from the observational data both in its sampling mechanism and in its treatment assignment mechanism, we first resample the NCR data to reflect the distribution in RCT data. Next, we correct for potential confounding using three alternative types of methods and consequentially evaluate their estimates to those obtained in the RCT. 
\parttitle{Results}
We find that while all estimators qualitatively approximate to findings in the RCT, methods that can flexibly model the response (i.e., direct estimation and DR estimation) performed consistently superior to the inverse propensity score method. Subgroup analysis indicates that relatively simple models allow us to properly estimate the treatment effect. However, these simple models do not properly identify heterogeneous treatment effects in stage \rom{2} colon cancer. Careful sensitivity analysis using more flexible models demonstrates both the uncertainty and the potential heterogeneous treatment effect in stage \rom{2} cancer and provides robust estimation of treatment effect in stage \rom{3} cancer.   
\parttitle{Conclusions}
Our results suggest that both the direct method and the DR method, when executed with care, can be used to reliably estimate treatment effects based on registry data. This methodological validation opens the door to more extensive use of registry data for the estimation of (individual) treatment effects. Additionally, our identification of potentially meaningful subgroups of stage \rom{2} colon cancer patients who, based on our analysis seem to benefit from chemotherapy, should be further explored. 
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The keywords begin here                  %%
%%                                          %%
%% Put each keyword in separate \kwd{}.     %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{keyword}
\kwd{average treatment effect}
\kwd{methods comparison}
\kwd{confounding}
\kwd{guideline}
\end{keyword}

% MSC classifications codes, if any
%\begin{keyword}[class=AMS]
%\kwd[Primary ]{}
%\kwd{}
%\kwd[; secondary ]{}
%\end{keyword}

\end{abstractbox}
%
%\end{fmbox}% uncomment this for twcolumn layout

\end{frontmatter}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% The Main Body begins here                %%
%%                                          %%
%% Please refer to the instructions for     %%
%% authors on:                              %%
%% http://www.biomedcentral.com/info/authors%%
%% and include the section headings         %%
%% accordingly for your article type.       %%
%%                                          %%
%% See the Results and Discussion section   %%
%% for details on how to create sub-sections%%
%%                                          %%
%% use \cite{...} to cite references        %%
%%  \cite{koon} and                         %%
%%  \cite{oreg,khar,zvai,xjon,schn,pond}    %%
%%  \nocite{smith,marg,hunn,advi,koha,mouse}%%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%% start of article main body
% <put your article body there>

%%%%%%%%%%%%%%%%
%% Background %%
%%

\section{Background} 
\label{sec: Background}
Colorectal cancer is the third most commonly diagnosed cancer and the second leading cause of cancer death worldwide in 2018 \cite{bray2018global}. The majority of patients are curatively treated with surgical removal of the tumor, surrounding tissues, and lymph nodes leading to excellent long term survival \cite{brouwer2019gap}. For stage \rom{3} colon cancer, Randomized Control Trials (RCTs) \cite{labianca1995efficacy, andre2009improved, yothers2011oxaliplatin} have shown an improved overall survival when surgery is followed by adjuvant chemotherapy compared to surgery only. However, results from RCTs were not supportive of prescribing adjuvant chemotherapy to stage \rom{2} patients \cite{gill2004pooled} due to lack of convincing evidence. For example, the meta-analysis by Erlichman and Charles (1999) \cite{erlichman1999efficacy} showed in a pooled analysis a hazard ratio (HR) of 0.83 (90\% CI 0.72, 1.07) for disease-free survival and a HR of 0.86 (90 \% CI 0.68, 1.07) for overall survival for fluoropyrimidine monotherapy compared to observation. The same finding was also reported \cite{jongeneel2019estimating}: while both RCT and real world data suggest a clinically relevant benefit of adjuvant chemotherapy for stage \rom{2} colon cancer, the estimates do not reach statistical significance. A borderline significantly improved overall survival was found in stage \rom{2} cancer in QUASAR trial (Relative Risk 0.82 95\% CI 0.7-0.95) and consequently the authors suggest adjuvant chemotherapy for all stage \rom{2} patients. Without critical pathological and clinical characteristics to identify the heterogeneous treatment effect in subgroups (e.g., emergency resection, tumor size, grade, number of positive lymph nodes), conclusions from the QUASAR trial may however not be applicable to all stage \rom{2} patients. 

These mixed and inconsistent results from current trials may be partly due to the heterogeneity between populations resulting from selection rules of eligible patients. The strict selection rules in RCTs make baseline populations in trials differ from daily clinical practice, which can threaten the external validity of studies. Registry data, without strict patient selection, may truly reveal the clinical relevance and may provide additional evidence regarding the effect of the adjuvant chemotherapy on a population outside the context of trials. Consequently, researchers have recently tuned their interest to such registry data to estimate treatment effects because it covers a large, heterogeneous and country-wide population \cite{hugen2013prognosis,hugen2015colorectal,verhoeff2016adjuvant,jongeneel2019estimating}. 

Estimating the effect of adjuvant chemotherapy using registry data, however, is challenging. First, variations in the type, the dose, and the length of the chemotherapy commonly occur across individuals, which make comparison of treatment effect across studies difficult. For instance, although capecitabine and oxaliplatin (CAPOX) and fluorouracil, leucovorin, and oxaliplatin (FOLFOX) are regarded as current standard treatments regarding their improved survival compared to fluorouracil and leucovorin (FU/LV) \cite{andre2009improved}, the detailed information of the chemotherapy regimen for each individual is not always properly recorded. This hinders the estimation of the long-term effect of CAPOX and FOLFOX using existing registry data. To fully utilize the registry data regardless of the limitations, researchers  uses a binary variable to indicate the chemotherapy (FU/LV and FOLFOX) and control, and consequentially compares treatment effect using registry data with the estimation from trials (FU/LV vs. control) \cite{jongeneel2019estimating}; an approach we also adopt in the current study. 

Second, unlike RCT data, in which the observed and unobserved confounding variables are balanced by random treatment assignment, the observational registry data runs the risk of being confounded. More precisely, in observational data, the treatment assignment mechanism is not exactly known, and thus confounding cannot be ruled out \cite{austin2011introduction}. Understanding the performance of statistical methods to correct for such potential confounding is an extremely important methodological question: if we are able to properly control for confounding in registry data we open up a treasure throve of new data to be used to estimate (individual) treatment effects.

Many statistical methods to control for potential confounding have been suggested over the past few decades, most often under the assumptions of positivity and unconfoundedness which we describe in detail in section \ref{sec: controlling for confounding}. These methods can roughly be classified into one of three classes \cite{atan2018deep}:
\begin{enumerate}
    \item \textit{Direct methods:} Direct methods aim to directly model the outcome as a function of the treatment and potential covariates. Some methods consider the treatment itself as a covariate and fit the outcome directly, while others fit the separate outcome for the treatment and control groups \cite{dorie2019automated}. Using multivariate regression to directly model the associations between the outcome and the treatment while adjusting the covariates is the most commonly applied implementation of the direct method \cite{prentice1976use}. Other, more modern, implementations use more flexible models like Bayesian additive regression trees (BART) \cite{chipman2010bart, pratola2016, mohammadi2019continuous} or Causal Forests \cite{wager2018estimation}. 
    \item \textit{Inverse probability weighting (IPW) methods:} IPW methods control for potential confounding by estimating, and controlling for, the conditional probability that an individual unit receives the treatment (i.e., the propensity score, PS) \cite{austin2011introduction}. IPW methods estimate the treatment effect by weighting each sample with the inverse PS \cite{austin2015moving}, thereby generating the pseudo-population whose covariates are independent of the treatment assignment (as would be the case in an RCT).  
    \item \textit{Doubly robust (DR) methods:} DR methods combine a model for the treatment assignment mechanism and a model for the outcome, which can practically yield accurate estimation when one of the two models is good (but not necessarily consistent) \cite{dudik2011doubly,bang2005doubly}.
\end{enumerate}

Note that for any of the methods above researchers can make a variety of modeling choices. Historically, parametric models have been commonly used. For instance, in the majority of published studies about the treatment effect estimation in the field of colorectal oncology \cite{algra2012effects,nishihara2013long,jongeneel2019estimating} the direct method, using standard logistic regression, is employed. Although these parametric models are efficient and well-understood, they impose strong parametric assumptions about the relationship between potential confounders and outcomes. Recently, more flexible semi-parametric and non-parametric models have demonstrated good performance \cite{dorie2019automated}. In particular, BART seems well suited to uncover the true treatment effect when the treatment assignment is confounded \cite{dorie2019automated,hahn2017bayesian,wendling2018comparing,hill2011bayesian}.

When comparing direct estimation with IPW and DR methods in simulation studies, an improved performance of direct estimation over IPW methods is found \cite{dorie2019automated}, while in other study DR methods outperform both direct estimation and IPW methods \cite{lunceford2004stratification}. Comparisons between these different classes of methods are however hard to carry out in non-simulated scenarios: often when observational data (registry data) is available there is no ground truth of the estimate of interest available (i.e., there is no RCT data available). The unavailability of evaluations of different methods to control for confounding on real-world data hinders the adoption of these methods in general. 

The case of adjuvant chemotherapy for colon cancer is different, however: here we uniquely have access to a large RCT \cite{quasar2007adjuvant} (QUASAR trial) as an experimental benchmark and we have access to the observational data collected by the NCR. This provides us with the opportunity to evaluate the performance of different (classes of) methods to control for confounding on real-world registry data. To our best knowledge, the QUASAR trial is the best possible gold standard to compare the effect of chemotherapy verses control as it has large samples, long follow up years, and good completion of six months of chemotherapy.

To conclude, observational registry data and RCT data potentially differ in two respects which requires us to be aware of when estimating and comparing the treatment effect \cite{Imbens2018ACB}:
\begin{enumerate}
    \item \textit{the sampling mechanisms}, i.e., the mechanism by which subjects are included into studies, can differ substantially. While the absence of strict inclusion criteria is an argument in favor of the use of observational data, the sampling process will ultimately affect the distributions of the covariates and subsequently will affect the estimated average treatment effect. The sampling mechanism may also lead to differences in effect estimation between the RCT and the observational study. 
    \item \textit{the treatment assignment mechanism} differs: while in the RCT the assignment mechanism is known – and often uniformly random – for observational studies, the assignment mechanism is unknown and confounding is potentially present.
\end{enumerate}

In this study, of which we provide an overview in Figure \ref{fig:flowchart}, we attempt to control for both sources of divergence and aim to explore the effects of three classes of correction methods of controlling for the second source of divergence. This study thus helps us understand how these different methods perform on real-world registry data, potentially opening up their use to inform the current treatment guidelines. 

In the next section, we first describe the two data sources used in this study. In section \ref{sec: controlling for confounding} we describe the different methods for controlling for confounding that we examine in this paper: direct methods, IPW methods, and DR methods. In section \ref{sec: study population sampling} we explain our approach to ensure that, next to controlling for confounding, the sampling mechanism playing in both data sources are comparable: effectively, we estimate the distribution of the main covariates in the QUASAR trial as $\mathcal{D}_{QUASAR'}$, and subsequently subsample the NCR data accordingly ($X_{sub} \sim \mathcal{D}_{QUASAR'}$). Note that we explicitly order these two steps such that the models used to control for confounding are estimated using all the available data ($X_{NCR} \sim \mathcal{D}_{NCR}$). In section \ref{sec: results} we present our results: we compare several estimates obtained using both direct, IPW, and DR methods with the estimates of the average treatment effect (ATE) originating from the RCT, and explore the performance of IPW and DR with different trimming values. The results of the subgroup analysis and sensitivity analysis are also presented. Finally, in section \ref{sec: discuss} and \ref{sec:conclusion} we discuss our findings and provide recommendations for the analysis of observational data. 

\section{Overview of data sources} 
\label{sec: data sources} 
In this study we effectively used two data sources: on one hand we used the estimates resulting from the QUASAR trial \cite{quasar2007adjuvant} as the gold standard. Note that we did not have access to the patient-level data of the QUASAR trial and relied on the estimates and descriptive statistics provided in the original paper. On the other hand, we used registry data from the NCR ($N=52621$) as our observational dataset of choice. The impact of the discrepancy of the treatments in two data sources on the estimation will be further discussed and interpreted in section \ref{sec: discuss}.

\subsection{The gold standard: the QUASAR trial}
\label{subsec: QUASAR trial}
The QUASAR trial was an RCT designed to test the effects of adjuvant chemotherapy for colorectal cancer patients. Between May 25, 1994, and Dec 24, 2003, $N=3239$ patients from 150 centres in 19 countries were entered into the trial ($2963$ [91.5\%] with stage \rom{2} disease, $260$ [8\%] with stage \rom{3}, 16 [0.5\%] with stage 1, $2291$ [71\%] with colon cancer, $948$ [29\%] with rectum or colorectal cancer, median age 63 [IQR 56-68] years), of whom 82\% of patients were entered from UK.

In the trial, eligible patients had undergone the resection and had no evidence of distant metastases. Furthermore, they had no definite contraindications to chemotherapy. $1622$ patients were randomly assigned to receive chemotherapy consisting of FU/LV. The primary outcome was overall survival. We refer readers to \cite{quasar2007adjuvant} for more details regarding the exact protocol of the QUASAR trial.

\subsection{The NCR}
\label{subsec: NCR}
The NCR is a registry containing all incidences of cancer in the Netherlands managed by the Netherlands Comprehensive Cancer Organisation (Integraal Kankercentrum Nederland, IKNL). We analyzed a subset of patients who were diagnosed with colon cancer between 2006 and 2015 and who underwent surgical resection \footnote{We only focus on the colon cancer patients instead of colorectal cancer patients in NCR for two reasons: first in NCR, the treatment regimen for rectum cancer patients are much more complicated than colon cancer; second, reported by \cite{quasar2007adjuvant}, the treatment effect of adjuvant chemotherapy didn’t differ significantly by tumor site, namely, colon or rectum.}. Follow-up was completed until January 31, 2019 . The total study population consisted of pathological stage \rom{2} and \rom{3} ($N=52621$) patients . We excluded $100$ patients who had the missing value of treatment and incidence. The primary outcome is overall survival.

\subsubsection{Covariates description and data processing}
\label{subsec: covariates in NCR }
Table \ref{tab:covariates} presents the distribution of the covariates in the current colon-cancer guideline and some important prognostic covariates as well as the basic demographic covariates, in which we exclude $4586$ patients with clinical metastasis. For height and weight, a large number of values were not recorded (82.3 \% and 81.3 \%) which are imputed with median. The BMI is then computed with height and weight. For categorical covariates, the missing value is replaced with category as ’not recorded’. To explore the impact of covariates included on the estimation, we will consider four groups of covariates for analysis: 
\begin{enumerate}
    \item \textit{3 covariates:} age, sex, and stage;
    \item \textit{5 covariates:} 2 covariates (age, sex) in addition to 3 covariates in the guideline (stage, MSI, high risk \footnote{risk factors implied in the 2014 Dutch guideline for colon cancer \cite{oncoguide}: perforation at diagnosis = yes, pT = 4, lymphatic invasion = yes, angio invasion = EMVI or IMVI, grade = poor to undifferentiated or unknown, the number of lymph nodes assessed less than 10. As long as one of risk factors occurs in stage \rom{2} patients, patients are identified as ‘high risk = high’. For stage \rom{3}, 'high risk = not applicable'.});
    \item \textit{13 covariates:} 2 covariates (age, sex) in addition to 11 covariates in the guideline (stage, cM, pT, pN, high risk, MSI, number of lymph nodes assessed, colon perforation, lymphomatic invasie, agio invasie, grade);
    \item \textit{All covariates in Table \ref{tab:covariates}}
\end{enumerate}

\section{Controlling for confounding: Direct, IPW, and DR methods} 
\label{sec: controlling for confounding}
In this section, we will introduce three groups of methods for controlling for confounding in detail. Note that this controls for confounding; we cover controlling for sampling mechanism in section \ref{sec: study population sampling}.

To formalize our methods, we first formalize the problem of estimating treatment effects. We use the potential outcome framework by Rubin \cite[see][chapter 1]{rubin2005causal}. We have a population with distribution $\mathcal{D}$ and size $n$, let $X$ denote a vector of covariates with length $d$ sampled according to $\mathcal{D}$, $Z$ denote the binary treatment variable ($1$=adjuvant chemotherapy, $0$=no adjuvant chemotherapy) and $Y$ denote a binary response variable ($1=$=death, $0=$=alive). Capital roman letters denote random variable, while lower case letters denote realized values. We use $Y_{i}(1)$ to denote the potential outcome for subject $i$ under the treatment $Z=1$ and $Y_{i}(0)$ denote the potential outcome under the treatment $Z=0$. The observed outcome under the realized treatment for the subject $i$ can be denoted as
\begin{equation}\label{equ: observed outcome}
    Y_{i}=Z_{i}Y_{i}(1)+(1-Z_{i})Y_{i}(0).
\end{equation}
The population average treatment effect (ATE) with covariates $X \sim \mathcal{D}$ can be formulated as:
\begin{equation}
    \text{ATE} = \mathbb{E}[Y_{i}(1)-Y_{i}(0)] =  \mathbb{E}[Y_{i}(1)]-\mathbb{E}[Y_{i}(0)] 
\end{equation}
However, by definition we cannot learn the true treatment effect because only one of the $Y_{i}(1)$ and $Y_{i}(0)$ is observed, which is the fundamental problem of causal inference. To estimate the average treatment effect, assumptions of unconfoundedness
\begin{equation} \label{assumption: unconfoundedness}
    \{Y(0),Y(1)\} \independent Z | X
\end{equation}
and positivity (or overlap) 
\begin{equation}
    0 < P (Z = 1|X ) < 1 .
\end{equation} 
are made \cite{rosenbaum1983central}. In RCTs, the known assignment mechanism causes the treatment assignment $Z$ independent of both observed and unobserved covariates $X$, consequentially leading to the treatment $Z$ independent of potential outcomes under the assumption of unconfoundedness. In this case, the average treatment effect (ATE) can be formulated as 
\begin{equation}
    \begin{split}
        \text{ATE} & = \mathbb{E}[Y_{i}(1)-Y_{i}(0)] \\
                   & = \mathbb{E}[Y_{i}(1)]-\mathbb{E}[Y_{i}(0)]  \\
                   & = \mathbb{E}[Y_{i}(1)|Z_{i}]-\mathbb{E}[Y_{i}(0)|Z_{i}] \\
                   & = \mathbb{E}[Y_{i}|Z_{i}=1]-\mathbb{E}[Y_{i}|Z_{i}=0] \\
    \end{split}
\end{equation}
In observational study where the treatment assignment $Z$ is correlated with the covariates $X$, the ATE can be formulated as \citep{rubin1974estimating}:
\begin{equation}\label{equation: ATE}
    \begin{split}
        \text{ATE} & = \mathbb{E}[Y_{i}(1)-Y_{i}(0)] \\
                   & = \mathbb{E}[Y_{i}(1)]-\mathbb{E}[Y_{i}(0)]  \\
                   & = \mathbb{E}_{X \sim \mathcal{D}}\{\mathbb{E}[Y_{i}(1)|X_{i}]\} - \mathbb{E}_{X \sim \mathcal{D}}\{\mathbb{E}[Y_{i}(1)|X_{i}]\} \  \text{(law of iterated expectation)} \\
                   & = \mathbb{E}_{X \sim \mathcal{D}}\{\mathbb{E}[Y_{i}(1)|Z_{i},X_{i}]\} - \mathbb{E}_{X \sim \mathcal{D}}\{\mathbb{E}[Y_{i}(0)|Z_{i},X_{i}]\} \  \text{(unconfoundedness)} \\
                   & = \mathbb{E}_{X \sim \mathcal{D}}[\mathbb{E}(Y_{i}|Z_{i}=1,X_{i})] - \mathbb{E }_{X \sim \mathcal{D}}[\mathbb{E}(Y_{i}|Z_{i}=0,X_{i})] \  \text{(equation 1)}\\ 
    \end{split}
\end{equation}

Note that the population treatment effect is defined as the average treatment effect on the target population $\mathcal{D}$. To make the estimator of the treatment effect comparable across the QUASAR ($\mathcal{D}_{QUASAR}$) and NCR ($\mathcal{D}_{NCR}$), we need to resample the subjects $X_{sub}$ in NCR according to $\mathcal{D}_{QUASAR}$ ($X_{sub} \sim \mathcal{D}_{QUASAR} $)such that the population in NCR is matched to QUASAR. Although the true $\mathcal{D}_{QUASAR}$ is unknown, provided by the descriptive statistics of the covariates in QUASAR, we can estimate the $\mathcal{D_{QUASAR'}}$ by assuming the covariates are independent, which will be introduced in section \ref{sec: study population sampling}. 

Now that we have formalized the problem of estimating the treatment effect, we will introduce different methods for modeling the $\mathbb{E}[Y_{i}(1)]$ and $\mathbb{E}[Y_{i}(0)]$ in observational study, namely, direct methods, IPW methods, and DR methods. In each case we have multiple implementations: we discuss implementations relying on relatively simple regression models and implementations relying on the flexible, non-parametric BART model. All the models are fitted based on the population (N=$52521$). Then the estimators of the treatment effect of the target population in NCR can be obtained by averaging the treatment effect of subjects sampled according to $\mathcal{D_{QUASAR'}}$.

\subsection{Direct methods}
\label{subsec: direct estimation}
The direct method, also known as regression adjustment, uses the regression method to directly fit the outcome model $\mathbb{E}[Y|Z,X]$. If the unconfoundedness assumption holds, $\mathbb{E}[Y_{i}(1)]$ and $\mathbb{E}[Y_{i}(0)]$ can be obtained by marginalizing the conditional outcome of $Y$ given covariates $X$ across the population under treatment $Z=1$ and $Z=0$, respectively, and can be formalized as follows: 
\begin{equation}
    \begin{split}
        \mathbb{E}[Y(z)]  = \mathbb{E}_{X \sim \mathcal{D}}[\mathbb{E}(Y_{i}|Z_{i}=z,X_{i})], \text{for} \ z = 0,1  \\ 
    \end{split}
\end{equation}
which have been proved in equation \ref{equation: ATE}. Therefore, by fitting the model for potential outcomes conditional on the covariates to adjust for the direct and indirect associations between covariates and outcomes, we can obtain the potential outcomes for each subject and subsequentially the ATE of the population. In the following section, we will introduce two methods to fit the model for potential outcomes $\mathbb{E}(Y|Z,X)$ based on observed outcomes, namely, a parametric model of logistic regression (Direct/LReg) and a non-parametric model of BART (Direct/BART). 

\subsubsection{Logistic regression (LReg)} 
Logistic regression (LReg) is a commonly used parametric method to control for confounding. When the covariates are strongly predictive of outcomes, we can efficiently identify the relation between treatment and outcome by making use of the information about relations between covariates and outcomes. By specifying the main effect and interaction effect between treatment and each independent covariate, $\mathbb{E}[Y(1)]$ and $\mathbb{E}[Y(0)]$ can be expressed as:
\begin{equation}
    \begin{split}
        \mathbb{E}[Y(z)] &= \mathbb{E}_{X \sim \mathcal{D}}[\mathbb{E}(Y_{i}|Z_{i}=z,X_{i})] \\
                        &= n^{-1}\sum_{i=1}^{n}\frac{1}{1+exp(\hat{\beta}_{0}+\boldsymbol{\hat{\beta}_{1}}^{\prime}X_{i} + \hat{\alpha}z + \boldsymbol{\hat{\beta}_{2}}^{\prime}X_{i}z)^{-1}} \ \ \ \text{for} \ \  z=1,0
    \end{split}
\end{equation}
where $\hat{\beta}_{0}$, $\boldsymbol{\hat{\beta}_{1}}$, $\boldsymbol{\hat{\beta}_{1}}$ and $\hat{\alpha}$ are estimated by maximum likelihood. We perform lasso regularization to prevent overfitting \citep{zou2005regularization}.

\subsubsection{Bayesian additive regression tree (BART)} 
\label{subsub: direct/BART}
BART is a sum-of-trees model where each tree is constrained by a regularization prior to a weak learner, and fitting and inference are obtained by Bayesian Markov Chain Monte Carlo (MCMC) algorithm that generates samples from a posterior \cite{chipman2010bart}. \cite{hill2011bayesian} and \cite{green2012modeling} first promoted the use of BART to obtain the heterogeneous and average treatment effect. \cite{wendling2018comparing} also compared the BART with other methods using simulated observational data in healthcare and results show that BART has a low bias. Unlike many other non-parametric models, BART is flexible and it can handle non-linear main effects and multi-way interactions without much input from researchers \cite{tan2019bayesian}. We refer the readers to \cite{chipman2010bart} for details. For binary outcomes, $\mathbb{E}(Y|Z,X)$ can be estimated using a probit model and then $\mathbb{E}[Y(1)]$ and $\mathbb{E}[Y(0)]$ can be obtained by averaging the outcome conditional on $X$ for $Z=1$ and $Z=0$ across the population. 

In our study, we use the default setting in \cite{chipman2010bart} because of BART's robust performance with respect to various hyperparameter settings: the number of trees = $50$ (we found that the results show no gain with more number of trees incorporated), the number of posterior draws = $1000$, and the prior distribution of $\sigma$ is inverse gamma distribution $IG(\alpha,\beta)$ where shape parameter(splitting probability) of $\alpha$ is $0.95$ and rate parameter (depth penalty) $\beta$ is $2$, respectively. 

\subsection{Inverse probability weighting estimation (IPW) methods}
\label{subsec: IPW estimation}
IPW is one of the PS methods for correcting the mismatch of data distribution produced by the treatment assignment. We define PS as follows: 
\begin{equation}
    \pi(X) = \mathbb{E}(Z=1|X)
\end{equation}
We refer readers to \cite{rosenbaum1983central,lunceford2004stratification} for details about the PS. The IPW method is considered to have some methodological advantages because it imitates the same treament assignment mechanism as the RCT. Weighting by the PS effectively generates the pseudopopulation in which the treatment assignment is independent of observed covariates, generating the similar distribution of covariates between two groups. Hence, in this pseudopopulation, we can obtain $\mathbb{E}[Y(1)]$ and $\mathbb{E}[Y(0)]$ simply by averaging the outcome in treatment and control groups after weighting by the inverse of PS. Then the $\mathbb{E}[Y(1)]$ and $\mathbb{E}[Y(0)]$ with respect to the distribution of $X \sim \mathcal{D}$ in the target population with $n$ subjects can be obtained via 
\begin{equation}
        \mathbb{E}[Y(1)] = \left[\sum_{i=1}^{n}\frac{\mathbbm{1}(Z_{i}=1)}{\pi(X_{i})}\right]^{-1}\sum_{i=1}^{n}Y_{i}\frac{\mathbbm{1}(Z_{i}=1)}{\pi(X_{i})} 
\end{equation}
and
\begin{equation}
        \mathbb{E}[Y(0)] = \left[\sum_{i=1}^{n}\frac{\mathbbm{1}(Z_{i}=0)}{1-\pi(X_{i})}\right]^{-1}\sum_{i=1}^{n}Y_{i}\frac{\mathbbm{1}(Z_{i}=0)}{1-\pi(X_{i})}.
\end{equation}
Then the ATE can be estimated by directly comparing the difference between $\mathbb{E}[Y(1)]$ and $\mathbb{E}[Y(0)]$. In this study, we compute the PS using two regression models, namely, LReg with lasso regularization (IPW/LReg), and BART (IPW/BART) with default setting. We trim the weights at $0.1$ and $0.9$, which is implemented in both IPW and DR methods, and explore the effect of different trimmed values on the ATE estimation with these methods, which we will further illustrate in section \ref{sec: results}.

\subsection{Doubly robust methods (DR)}
\label{subsec: DR estimation}
The DR estimator combines a model for the potential outcomes as a function of covariates $\mathbb{E}(Y|Z, X)$ and a model for assignment mechanism as a function of covariates $\mathbb{E}(Z|X)$ (i,e, $\pi(X)$). Particularly, we fit a separate outcome model $m_{z}= \mathbb{E}(Y_{i}|Z_{i}=z,X_{i})$ for the treatment group $Z=1$ and control group $Z=0$. $m_{z}(X_{i})$ is the predicted values for the subject $i$ given $X_{i}$ and the treatment $Z_{i}=z$. Then $\mathbb{E}[Y_{i}(1)]$ and $\mathbb{E}[Y_{i}(0)]$ with respect to the distribution of $X \sim \mathcal{D}$ in the target population can be estimated via \cite{lunceford2004stratification}: 
\begin{equation}
    \begin{split}
        \mathbb{E}[Y(1)] & = n^{-1}\sum_{i=1}^{n}\left\{\frac{\mathbbm{1}(Z_{i}=1)}{\pi(X_{i})}\left[Y_{i}-m_{1}(X_{i})\right] + m_{1}(X_{i}) \right\} \\
        \mathbb{E}[Y(0)] & = n^{-1}\sum_{i=1}^{n}\left\{\frac{\mathbbm{1}(Z_{i}=0)}{1-\pi(X_{i})}\left[Y_{i}-m_{0}(X_{i})\right] + m_{0}(X_{i}) \right\} 
    \end{split}
\end{equation}

The name DR is derived from the fact that the estimation remains unbiased when either a model for treatment assignment mechanism or a model for the potential outcomes is correctly specified \citep{bang2005doubly}. Combining a direct estimator and an IPW estimator under the form of a DR estimator leads to lower bias than the direct estimator alone, and lower variance than the IPW estimator alone \citep{bibaut2019more}.See \cite[][p148-149]{tsiatis2007semiparametric} for more details of the proof of the DR property. In this study, four DR estimators are constructed based on two models for potential outcomes $\mathbb{E}[Y|Z,X]$ fitted by LReg and BART and two models for assignment mechanism $\mathbb{E}[Z|X]$ fitted by LReg and BART. Table \ref{table:methods} shows the details of these four combinations.

\subsection{Additional methods}
\label{subsec: additional methods}
Next to the three groups of methods described above, we also include two other alternatives:
\begin{enumerate}
    \item To explore the effect of sampling mechanism and treatment assignment on estimation, we roughly model potential outcomes $\mathbb{E}(Y_{i}|Z_{i})$ on the whole population before the subsampling ($(X_{NCR} \sim \mathcal{D}_{NCR}$) and the study population after subsampling ($X_{sub} \sim \mathcal{D}_{QUASAR'}$)respectively. 
%because the Crude LReg is the model for the association observed when comparing subpopulations (defined by the levels of treatment) of the study population \citep{robins2000marginal}. 
    \item Furthermore, in addition to covariates, we also include the PS in the direct method $\mathbb{E}[Y_{i}|Z_{i},X_{i},\pi(X_{i})]$ (Direct/ps-BART). Both the model for potential outcomes and PS (i.e., $\pi(X_{i})=\mathbb{E}(Z_{i}|X_{i})$) are fitted by BART with the default setting. Direct/ps-BART incorporates the estimate of PS in the specification of the outcome model, implicitly inducing a covariate-dependent prior to the regression function, which even outperforms the BART model (\cite{hahn2017bayesian, santos2018tree}.
\end{enumerate}

\subsection{Overview of methods}
Table \ref{table:methods} provides an overview of the methods introduced in section \ref{subsec: direct estimation}, \ref{subsec: IPW estimation}, \ref{subsec: DR estimation} and \ref{subsec: additional methods}. Parametric LReg and non-parametric BART with the default setting introduced in section \ref{subsub: direct/BART} are used to fit both models for the potential outcomes $\mathbb{E}[Y_{i}|Z_{i},X_{i}]$ and the treatment assignment $\mathbb{E}(Z_{i}|X_{i})$.

When fitting the model using LReg, we need to specify the relationship between each variable: the main effect of covariates and the interaction effect between the treatment and each covariate are specified for the model of $\mathbb{E}(Y_{i}|X_{i}, Z_{i})$ in the direct method; the main effect of all covariates and two-way interaction effects of all covariates are specified for the model of $\mathbb{E}(Y_{i}|X_{i}, Z_{i}=z)$ in DR method and the model of $\mathbb{E}(Z_{i}|X_{i})$ in IPW method. To prevent overfitting, lasso regularization is implemented on all the models fitted by LReg, the loss function is to minimize the mean squared error, and the ten-fold cross-validation is implemented to choose the regularization parameter. %When fitting all the models with BART, due to BART’s robustness to change in the exact specification \cite{chipman2010bart} and its advantage of handling the non-linear association, the cross-validation and the manual pre-specifications of associations between each variable are left out.

Note that all the models are initially fitted based on the whole population ($X_{NCR} \sim \mathcal{D}_{NCR}$) and the estimator of ATE for comparison with the QUASAR are then obtained from the subsamples ($X_{sub} \sim \mathcal{D}_{QUASAR'}$) from the whole population. We quantify the sampling uncertainty by stratified sampling the study population for 5000 times, which we will introduce in section \ref{sec: study population sampling}.

\section{Study population sampling}
\label{sec: study population sampling}
In this section, we will select the study population from the whole population ($N=52621$) to ensure that the subsamples of the NCR ($X_{sub} \sim \mathcal{D}_{QUASAR'}$) is similar to the samples in the QUASAR trial ($X_{QUASAR} \sim \mathcal{D}_{QUASAR}$). We first select the eligible patients ($N=45539$) by complying with the inclusion criteria of QUASAR. Then we subsample the target population according to the probablity estimated from the joint distribution of $\mathcal{D}_{QUASAR'}$ over important covariates by assuming independence.

\subsection{Eligible participates}
To select the eligible participates for the study population, we first sought to comply with the same inclusion criteria as the QUASAR trial \cite{quasar2007adjuvant}. From a total of $N=52621$ patients with colon cancer, $100$ patients missed the information of treatment and incidence, $45442$ patients were without distant metastases (cM==1), diagnosed stage \rom{2} or stage \rom{3} and aged 23 to 86. All of these $45442$ remaining patients are included for a two-year follow-up analysis and $n=34832$ for a five-year follow-up analysis in which patients diagnosed after 2013 are excluded. Figure \ref{fig:pateints selection} provides an overview of this first selection step \footnote{It is worth noting that in our study the models for potential outcomes and assignment mechanism were fitted based on a very large sample ($N=52521$), utilizing the information from the general population which is beyond the scope of most RCTs. Thus, it is reasonable to, in the future, estimate the treatment effect of patients who are not included in the RCT using these methods.}. 

\subsection{Sampling}
\label{subsec: sampling mechanism}
Next to following the inclusion criteria in the QUASAR trial, we also estimate the distribution $\mathcal{D}_{QUASAR'}$ over covariates by correcting for the distribution employed in the NCR. Effectively, we ensured that our subsample $X_{sub} \sim \mathcal{D}_{QUASAR'}$ resulting from the NCR data corresponds – in terms of the distribution of several important covariates – as closely as possible to the QUASAR trial data. By correcting for the sampling mechanism such that the distribution of covariates are the same in both data sources (i.e., $\mathcal{D}_{QUASAR'}$ and $\mathcal{D}_{QUASAR}$), the estimator obtained by averaging the outcomes across subsamples of NCR ($X_{sub} \sim \mathcal{D}_{QUASAR'}$) can be reliably compared to the QUASAR in further analysis. To do so, we divide the eligible population into $16$ subgroups by stage, sex, and age, and sampled from each subgroup according to the joint probability of their cancer stage, their sex, and their age as observed in the QUASAR sample (see Appendix) \footnote{No other covariate distributions are reported in the paper of QUASAR trial. Note that only the marginals are reported, and we estimate the joint distribution by assuming independence supported by \cite{hayne2001current}.}. 

Then we implemented this subsampling such that the resulting sample was as large as possible while maintaining the desired joint distribution of the covariates. The stratified subsampling leads to the selection of $7704$ and $6378$ patients from the NCR for two-year follow-up and five-year follow-up analysis; Table \ref{tab:subsampling matching} demonstrates that our subsampling procedure effectively ensured the desired joint distribution of the NCR subsample comparable to the original QUASAR sample. Note that, to incorporate the variability induced by our subsampling procedure in our remaining procedure, we create $m=5000$ stratified (sub)samples.

\section{Results}
\label{sec: results} 
In this section, we compare the performance of the different methods to control for confounding. We examine the ATE in two and five years of follow-up. The main results are presented in Figure \ref{fig:ATE}, where the ATEs in each method are estimated with four groups of covariates. Numerical details are provided in Table \ref{tab:Results2} and \ref{tab:Results5}. In addition, to assess the magnitude of violation of positivity assumption in IPW and DR methods, we visualize the PS distribution between treatment and control groups and explore the variation of IPW and DR to different trim values – an approach to improve the performance in face of the positivity violation – which are shown in Figure \ref{fig:trim and ps}. We also show the results of subgroup analysis to identify the heterogeous treatment effect and sensitivity analysis to assess the assumption of unconfoundedness. %In the following section, we will interpret the results of direct, IPW, and DR methods in details and give a short discussion of comparison.

\subsection{Results of the estimation of treatment effect}
\label{subsec: ate}
To explore the impact of different covariates sets used on the estimation, we compare the results of three methods using four sets of covariates, which are shown in Figure \ref{fig:ATE}. We find that the methods using the least covariates are the closest to the ground truth of RCT while the results using more covariates substantially overestimate the ATE. 

We then explore the results of the different methods after the two-year follow- up. In the RCT the estimated ATE -0.021 (95\% CI -0.039, -0.001). The result is presented in Figure \ref{fig:ATE} using the black dotted line (point estimate) and the grey band (CI). Using the NCR data without adjusting for any confounders — i.e., the Crude method — we obtain an ATE of -0.069 (95\% CI -0.076, -0.062) for $X_{NCR} \sim \mathcal{D}_{NCR}$ and a median ATE of -0.019 (range from -0.041, 0.003) for $X_{sub} \sim \mathcal{D}_{QUASAR'}$. Thus, if neither control for sampling mechanism nor treatment assignment, the Crude method overestimates the treatment effect. This overestimation could be alleviated by controlling for the sampling mechanism, although the estimation still deviates from the RCT and shows large variance produced by the sampling process. 

Regarding the impact of the number of covariates controlled for on the ATE estimation, it seems that using sex, age, and stage can properly obtain the ATE as close as possible to the RCT. When high risk and more covariates added (shown in the second, the third and the fourth row in Figure \ref{fig:ATE}), the estimate deviates far from the results from the RCT. 

Regarding the direct methods, the direct/BART and direct/ps-BART are very close to the ground truth while the direct/LReg seems to overestimate the ATE. Notably, the variation over different subsamples obtained to correct for the sampling mechanism in direct methods is very small: the direct method seems very robust to this sampling variation. Compared to direct methods, the IPW methods, however, seem to have a huge variation. The DR methods seem to – as expected – combine the best of both worlds: the estimates are close to the ground truth, and they exhibit small sampling variance compared to IPW methods. For each of the DR methods, irrespective of using parametric or non-parametric models, the majority of the point estimates are within the confidence bound of the original ICR. The DR using 13 and all covariates all show the same results, which probably means that unconfoundedness is hold in both cases. 

In the five-year follow-up analysis, the ATE in RCT is -0.025 (95\% CI -0.059, 0.008). In NCR, ATEs using crude method for $X_{NCR} \sim \mathcal{D}_{NCR}$ and $X_{sub} \sim \mathcal{D}_{QUASAR'}$  are -0.073 (95\% CI -0.083, -0.063) and a median of -0.031 (range from -0.065, 0.005), respectively. These results are qualitatively similar to those obtained after two years: the Crude method overestimates the treatment effect, all the DR methods are very close to the ground truth. The majority of the point estimates obtained from the DR method – whether using a parametric or a non-parametric approach – fall within the confidence bound obtained in the RCT. 

To conclude, our results show that a) in general all the estimates are consistent in terms of direction with the results of RCT; b) strength of the association between treatment and outcome varies across statistical methods and covariates included, which is also indicated in \cite{cnossen2018adjusting}. Clearly, without cautious controlling for the confounding and sampling mechanism, results using the crude method will not be reliable in the aspect of both magnitude and directionality of the treatment effect, as reported in Crude methods on population $X_{NCR} \sim \mathcal{D}_{NCR}$ and $X_{sub} \sim \mathcal{D}_{QUASAR'}$. However, when correction methods are used carefully, we do obtain estimates that are close to the ground truth; c) regarding the performance of the four groups of methods, the DR estimations, particularly modeled by the non-parametric BART, seem more robust in obtaining the treatment effect than IPW. Within the direct methods the choice of a model, the non-parametric model, namely, ps-BART, is quite flexible and performs properly. 

\subsection{The results of the estimation with different trimmed values}
\label{subsec: results of trimmed values}
Regarding the effect of covariates on the estimation of ATE, we find that the estimate using the least number of covariates is closest to the ground truth from the RCT while estimates from models using more covariates deviate from the ground truth. The right panel of Figure \ref{fig:trim and ps} shows the distribution of PS between treatment and control, indicating that adjusting for more covariates leading to more extreme PS. Although IPW estimator is theoretically unbiased, it is sensitive to extreme weights, and hence weights are commonly truncated or bounded to shrink the variance of the estimator. %However, there is no common practice of which trim value to choose, therefore, we need to carefully explore the performance of IPW and DR corresponding to various trim values.

In our study, weights are trimmed at a fixed level, namely, at 0.1 and 0.99, 0.5 and 0.95, 0.1 and 0.9, and 0.2 and 0.8 \cite{austin2015moving} and the results are shown in the left panel of Figure \ref{fig:trim and ps}. We found that across the results, trimming the IPW weights decreases the variance of the estimator but makes the point estimate substantially vary across trimmed values. Besides, we also find that IPW methods benefit from the trimming, which is quite clear when looking at complex models using 5, 13 and all covariates. This impact on the complex model can also be validated by looking at the estimates with only three covariates: in this case, the point estimation doesn’t vary across different trim values probably due to the comparable specification of PS model with the ground truth. Overall, in our study, trimming the extreme weight can help the IPW method approach to the ground truth.

\subsection{Subgroup analysis}
\label{subsec: subgroup analysis}
From the previous results, we find that high risk and other pathological characteristics are sensitive to the estimation. To understand the reason that leads to this overestimation, we explore the estimation in the subgroups of stage \rom{2} and stage \rom{3} classified by the high-risk level using Direct/ps-BART. The results are shown in Figure \ref{fig:subgroup analysis}. We found that by adjusting the high risk the method can identify the heterogeneous treatment effect among the subgroups in stage \rom{2} cancer. This finding is supported by a) the treatment effect on stage \rom{3} patients are the same in the model with and without adjusting for high risk covariate ; b) the treatment effect for stage \rom{2} colon cancer who are not identified as high risk are the same and is approximately equal to zero in all the models, regardless of  high risk being adjusted or not; c) the heterogeneity of treatment effect in stage \rom{2} cancer is supported by many studies indicating that part of stage \rom{2} high-risk colon cancer. for example, pT4, can achieve the survival benefit from the adjuvant chemotherapy \cite{babcock2018high,verhoeff2016adjuvant}. 

Based on these findings, although the models using three covariates can estimate the treatment effect comparable to the RCT, we can’t rule out the possibility that the true treatment effect of the NCR data is very close to the estimate from the models using 5, 13 and all covariates, because no pathological information was provided in the RCT to identify the heterogeneous treatment effect of the high risk stage \rom{2} patients. Overall, the model using three covariates without adjusting for high risk can estimate the treatment effect as close as possible to the ground truth from the RCT while at the cost of inability to identify the heterogeneity among the population, in particular, the subgroup of stage \rom{2} cancer. 

\subsection{Sensitivity analysis}
The subgroup analysis implicitly indicates that the estimation for stage \rom{3} patients are robust to various covariates included while stage \rom{2} patients are quite sensitive to the type and the number of covariates. Hence, we can use sensitivity analysis in subgroups to test for the impact of an unmeasured confounder and assess the extent to which we allow the unconfoundedness assumption to be violated in subgroups.

In our study, a two-parameter sensitivity analysis based on flexible non-parametric methods BART is implemented \cite{dorie2016flexible} and the results are presented in Figure \ref{fig:senstivity analysis}. It implies that the estimate of stage \rom{2} colon cancer, particularly those are unknown risk, is sensitive to the unmeasured confounder and the sign of the estimate can even change in the case where, for instance, some high risk factors are not controlled for. This implies that more information about this subgroup should be collected to identify the high risk stage, otherwise the estimation might not truly reveal the treatment effect. However, for stage \rom{3} patients, the estimation of the treatment effect is quite robust to the unmeasured confounders and seemingly deterministic.  

\section{Discussion}
\label{sec: discuss}
In this study, we exploited the unique opportunity of having access to large-scale colon cancer registry data to compare three different methods of controlling for potential confounding in observational data, namely, direct method, IPW method, and DR method. Our study, by comparing the resulting estimates to those obtained in a large scale RCT, validated that these methods can be used to obtain the estimates of treatment effects close to the ground truth. This opens up the door to use the large-scale, high-quality registry data regarding colon cancer to further strengthen the (Dutch) guidelines by adding, as a third source of information next to RCT data and expert consensus, the available registry data. It has to be noted, however, that adding this third source requires care: methods using naive estimators — such as the crude method evaluated in this study — can be highly biased. We can only use the registry data when both the sampling process and the assignment process are properly controlled for.

Next to demonstrating the utility of observational data in this specific case, we also find evidence of the advantages of modeling the outcome (as done in the direct method and the DR method) over solely modeling the treatment assignment process (IPW method). This latter result might be because it is hard to properly estimate this process. Besides, the supreme performance of non-parametric BART over simple parametric LReg is quite promising, which implies its use for further analysis. In the subgroup analysis, we find that when adjusting for the high-risk in models, the heterogeneous treatment effect could be identified, although the estimation might deviate from the result of the RCT. This dilemma of whether or not to include more covariates as to obtain both the accurate and heterogeneous estimation could be solved by access to a variety of data with complete information in each subgroup, which could probably shrink the variation of estimation from the model using all the covariates as shown in Figure \ref{fig:subgroup analysis}. The unconfoundedness assumption is also explored using sensitivity analysis and the result is qualitatively reasonable. 

In general, the estimates we obtained from the NCR observational data — as can be seen in Figure \ref{fig:ATE} — are often higher than those in the ground truth Quasar trial. While we have discussed several sources of this discrepancy due to our explored correction methods, there are possible two alternative explanations. On one hand, although we have resampled the registry data to as much as possible match the population in two data sources, there is still uncorrectable heterogeneity between populations at baseline. For example, more stage \rom{2} patients are identified as high risk in registry data than RCT, which consequentially increases the overall estimation. On the other hand, it is now known that CAPOX and FOLFOX are more effective than FU/LV \cite{andre2009improved}. The fact that in the NCR data more patients received CAPOX/FOLFOX as opposed to FU/LV, while in the QUASAR trial all patients received FU/LV, might truly increase the average treatment effect of chemotherapy overall in the NCR.

Regarding the subsampling process, we made two assumptions a) the assumption of independence among gender, age and stage when estimating the distribution of population in QUASAR based on the population-level statistics provided in the end of the enrollment; b) the assumption of no change of the estimated distribution over time when subsampling the target population accordingly from the NCR data. These two assumptions, fortunately, can approximately be checked by the paper published during that time. Some studies reported that males and females have a similar distribution of incidence by age from 1971 to 1999 \cite{quinn2001cancer,mistry2011cancer}. Therefore, the variation of the distribution of age and gender in colon cancer over time could be ignorable and hence doesn't violate the assumption. However, all these two assumptions are still untestable without patient-level data with covariate information. %Finally, the diagnose system, surgery and other unknown variables in two data sources might change over time, which we assume doesn’t affect drawing the conclusion. 

\section{Conclusion}
\label{sec:conclusion}
To conclude, our study shows that it is possible to obtain treatment effects from observational data when analyzed carefully. Without cautious controlling for confounders, the results from observational data may mislead people to make incorrect decisions. Our study also shows that the colon cancer registry data, when handled with care, can potentially be used to inform future guidelines: at the very least the ATE obtained from this data source – with the appropriate controls – provides estimates that are similar to those obtained in an RCT.


%\nocite{oreg,schn,pond,smith,marg,hunn,advi,koha,mouse}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                          %%
%% Backmatter begins here                   %%
%%                                          %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{backmatter}

\section*{List of abbreviations}
RCT: randomized controlled trial

IPW: inverse probability weighting

DR: doubly robust

NCR: Netherlands cancer registry 

HR: hazard ratio

QUASAR: quick and simple and reliable

CAPOX: capecitabine and oxaliplatin

FOLFOX: fluorouracil, leucovorin, and oxaliplatin

FU/LV: fluorouracil and leucovorin

BART: bayesian additive regression trees

IKNL: Integraal Kankercentrum Nederland

BMI: body mass index

ATE: average treatment effect

LReg: logistic regression



\section*{Declarations}

\subsection*{Ethics approval and consent to participate}
Permission to access the Netherlands Cancer registry (NCR) and the databases of Statistics Netherlands (CBS) was granted by both IKNL (Netherlands Comprehensive Cancer Care Organisation) and the CBS (Statistics Netherlands). According to the Dutch Medical Research Involving Human Subjects Act, this study did not require medical-ethical approval, as was confirmed in writing by the medical ethical committee of the AMC, 10 July 2013. We followed the ethical principles for medical research involving human subjects as laid down in the Declaration of Helsinki and adopted by the World Medical Association (WMA Declaration of Helsinki, 2000).

\subsection*{Consent for publication}
Not applicable.

\subsection*{Availability of data and materials}
The data that support the findings of this study are available from IKNL and CBS but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available.

\subsection*{Competing interests}
The authors declare that they have no competing interests.

\subsection*{Funding}
The project was funded by the grant from the China Scolarship Council Program, grant number 201806860047. The funding body had no role in the design of the study and collection, analysis, and interpretation of data, nor in writing the manuscript.

\subsection*{Authors' contributions}
SLJ and EV carried out the statistical analyses. SLJ and MCK drafted the manuscript. SLJ participated in the interpretation of the statistical analyses. All authors critically read the drafts of this paper, and read and approved its final version.

\subsection*{Acknowledgements}
The authors thank the registration team of the Netherlands Comprehensive Cancer Organisation (IKNL) for the collection of data for the Netherlands Cancer Registry, the Prospective National Colon Rectum Cancer (PLCRC) working group for the collection of Patient Reported Outcome data as well as IKNL staff for scientific advice.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  The Bibliography                       %%
%%                                                         %%
%%  Bmc_mathpys.bst  will be used to                       %%
%%  create a .BBL file for submission.                     %%
%%  After submission of the .TEX file,                     %%
%%  you will be prompted to submit your .BBL file.         %%
%%                                                         %%
%%                                                         %%
%%  Note that the displayed Bibliography will not          %%
%%  necessarily be rendered by Latex exactly as specified  %%
%%  in the online Instructions for Authors.                %%
%%                                                         %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% if your bibliography is in bibtex format, use those commands:
\bibliographystyle{bmc-mathphys} % Style BST file (bmc-mathphys, vancouver, spbasic).
\bibliography{References}      % Bibliography file (usually '*.bib' )
% for author-year bibliography (bmc-mathphys or spbasic)
% a) write to bib file (bmc-mathphys only)
% @settings{label, options="nameyear"}
% b) uncomment next line
%\nocite{label}

% or include bibliography directly:
% \begin{thebibliography}
% \bibitem{b1}
% \end{thebibliography}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Figures                       %%
%%                               %%
%% NB: this is for captions and  %%
%% Titles. All graphics must be  %%
%% submitted separately and NOT  %%
%% included in the Tex document  %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Do not use \listoffigures as most will included as separate files

\section*{Figures}

\begin{figure}[h!]
   %\centering
   %\includegraphics[width=0.9\columnwidth]{flowchart.pdf}
   \caption{the flowchart of this study}
   \label{fig:flowchart}
\end{figure}

\begin{figure}[h!]
    %\centering
    %\includegraphics[width=0.9\columnwidth]{patient_selection.pdf}
    \caption{The first step of selecting eligible patients from the NCR; patients without follow up, and patients who did not meet the inclusion criteria for the QUASAR trial are removed.}
    \label{fig:pateints selection}
  \end{figure}

\begin{figure}[h!]
    %\centering
    %\includegraphics[width=0.9\columnwidth]{ATE.pdf}
    \caption{ATE in two (left panel) and five (right panel) years of follow-up. It is worth noting that the treatment effect on colon cancer in the QUASAR trial should be smaller than the dotted line in the plot because in the QUASAR 29\% rectum cancer patients were included for whom the treatment effect is higher than colon cancer based on the results reported in the paper.}
    \label{fig:ATE}
\end{figure}

\begin{figure}[h!]
    %\centering
    %\includegraphics[width=0.9\columnwidth]{trim_ps.pdf}
    \caption{left panel: ATE of IPW and DR methods with different trimming values in two years follow-up. Right panel: the distribution of propensity score modeled by BART across the treatment and control groups.}
    \label{fig:trim and ps}
\end{figure}

\begin{figure}[h!]
    %\centering
    %\includegraphics[width=0.9\columnwidth]{subgroup.pdf}
    \caption{Subgroup analysis of stage \rom{2} (risk=normal, risk=high) and stage \rom{3} in two and five years follow up.}
    \label{fig:subgroup analysis}
\end{figure}

\begin{figure}[h!]
    %\centering
    %\includegraphics[width=0.9\columnwidth]{SA.pdf}
    \caption{The sensitivity analysis implemented using treatSens package in stage \rom{2} and stage \rom{3} cancer, controlling for gender and age. NS denotes ‘Not significant’.}
    \label{fig:senstivity analysis}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Tables                        %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Use of \listoftables is discouraged.
%%

%\section*{Tables}
%\begin{table}[h!]
%\caption{Sample table title. This is where the description of the table should go.}
%      \begin{tabular}{cccc}
%        \hline
%           & B1  &B2   & B3\\ \hline
%       A1 & 0.1 & 0.2 & 0.3\\
%        A2 & ... & ..  & .\\
%       A3 & ..  & .   & .\\ \hline
%      \end{tabular}
%\end{table}

\section*{Tables}
\begin{table}[h!]
    \caption{Baseline characteristics of stage \rom{2} and \rom{3} patients who underwent surgery between 2006 And 2015
    } \label{tab:covariates}
    \begin{tabular}{l l l l} 
    \hline
    \textbf{Characteristics}& \makecell[l]{\textbf{Whole population}\\\textbf{(n=47935)}}& \makecell[l]{\textbf{Chemotherapy}\\ \textbf{(n=15712)}} & \makecell[l]{\textbf{No chemotherapy}\\\textbf{(n=32223)}} \\
    \hline
    \textbf{Age}             & 70.6     ( 11.2          )& 64.3      ( 9.8           )& 73.7     ( 10.4          )\\
    \textbf{Sex}\\
    male                         & 24383    ( 50.9          )& 8336      ( 53.1          )& 16047    ( 49.8          )\\
    female                       & 23552    ( 49.1          )& 7376      ( 46.9          )& 16176    ( 50.2          )\\
    \textbf{BMI}                         & 26.4     ( 6.3           )& 26.5      ( 7.2           )& 26.4     ( 5.9           )\\
    \textbf{Number of co-morbidity}\\
    1                & 9159     ( 19.1          )& 3646      ( 23.2          )& 5513     ( 17.1          )\\
    $\geq$2                 & 6101     ( 12.7          )& 1528      ( 9.7           )& 4573     ( 14.2          )\\
    not recorded           & 32675    ( 68.2          )& 10538     ( 67.1          )& 22137    ( 68.7          )\\
    \textbf{Perforation at diagnosis}\\
    no          & 6095     ( 12.7          )& 2248      ( 14.3          )& 3847     ( 11.9          )\\
    yes          & 330      ( 0.7           )& 125       ( 0.8           )& 205      ( 0.6           )\\
    unknown    & 508      ( 1.1           )& 181       ( 1.2           )& 327      ( 1.0             )\\
    not recorded    & 41002    ( 85.5          )& 13158     ( 83.7          )& 27844    ( 86.4          )\\
    \textbf{Location of the tumor}\\
    proximal                     & 27472    ( 57.3          )& 8017      ( 51.0            )& 19455    ( 60.4          )\\
    distal                       & 19479    ( 40.6          )& 7401      ( 47.1          )& 12078    ( 37.5          )\\
    other/not otherwise specified(nos)               & 984      ( 2.1           )& 294       ( 1.9           )& 690      ( 2.1           )\\
    \textbf{pT:pathological size of tumor}\\
    T0                          & 3        ( 0.0             )& 3         ( 0.0             )& 0        ( 0.0             )\\
    T1                          & 515      ( 1.1           )& 330       ( 2.1           )& 185      ( 0.6           )\\
    T2                          & 1890     ( 3.9           )& 1204      ( 7.7           )& 686      ( 2.1           )\\
    T3                          & 36970    ( 77.1          )& 10423     ( 66.3          )& 26547    ( 82.4          )\\
    T4                          & 8557     ( 17.9          )& 3752      ( 23.9          )& 4805     ( 14.9          )\\
    \textbf{pN:pathological spread to lymph nodes}\\
    N0                          & 25622    ( 53.5          )& 2139      ( 13.6          )& 23483    ( 72.9          )\\
    N1                          & 14687    ( 30.6          )& 8647      ( 55.0            )& 6040     ( 18.7          )\\
    N2                          & 7327     ( 15.3          )& 4884      ( 31.1          )& 2443     ( 7.6           )\\
    NX                          & 299      ( 0.6           )& 42        ( 0.3           )& 257      ( 0.8           )\\
    \textbf{Number of lymph nodes assessed}              & 17     ( 10.5  )& 18      ( 10.9 ) & 17     ( 10.3 )\\
    \textbf{cM:clinical presence of metastasis}\\
    M0                          & 45790    ( 95.5          )& 15157     ( 96.5          )& 30633    ( 95.1          )\\
    MX                          & 2145     ( 4.5           )& 555       ( 3.5           )& 1590     ( 4.9           )\\
    \textbf{Stage}\\
    \rom{2}                      & 25801    ( 53.8          )& 2114      ( 13.5          )& 23687    ( 73.5          )\\
    \rom{3}                      & 22134    ( 46.2          )& 13598     ( 86.5          )& 8536     ( 26.5          )\\
    \textbf{Morphology}\\
    mucinous                     & 6991     ( 14.6          )& 2108      ( 13.4          )& 4883     ( 15.2          )\\
    adenocarcinoom               & 39703    ( 82.8          )& 13189     ( 83.9          )& 26514    ( 82.3          )\\
    signet ring cell             & 611      ( 1.3           )& 238       ( 1.5           )& 373      ( 1.2           )\\
    neuroendocrine tumors        & 211      ( 0.4           )& 33        ( 0.2           )& 178      ( 0.6           )\\
    other/nos                    & 419      ( 0.9           )& 144       ( 0.9           )& 275      ( 0.9           )\\
    \textbf{Grade}\\
    good to moderate          & 35561    ( 74.2          )& 11314     ( 72.0            )& 24247    ( 75.2          )\\
    poor to undifferentiated  & 8078     ( 16.9          )& 3035      ( 19.3          )& 5043     ( 15.7          )\\
    unknown                   & 4296     ( 9.0             )& 1363      ( 8.7           )& 2933     ( 9.1           )\\
    \textbf{Microsatellite instable (MSI) present}\\
    no: microsatellite stable                        & 1205     ( 2.5           )& 703       ( 4.5           )& 502      ( 1.6           )\\
    yes: microsatellite unstable                     & 328      ( 0.7           )& 117       ( 0.7           )& 211      ( 0.7           )\\
    unknown                      & 4158     ( 8.7           )& 1355      ( 8.6           )& 2803     ( 8.7           )\\
    not recorded                 & 42244    ( 88.1          )& 13537     ( 86.2          )& 28707    ( 89.1          )\\
    \textbf{Lymphatic invasion}\\
    not found      & 4213     ( 8.8           )& 1351      ( 8.6           )& 2862     ( 8.9           )\\
    yes     & 1079     ( 2.3           )& 607       ( 3.9           )& 472      ( 1.5           )\\
    suspect & 4        ( 0.0             )& 2         ( 0.0             )& 2        ( 0.0             )\\
    unknown & 273      ( 0.6           )& 102       ( 0.6           )& 171      ( 0.5           )\\
    not recorded & 42366    ( 88.4          )& 13650     ( 86.9          )& 28716    ( 89.1          )\\
    \textbf{Angio invasion}\\
    not found            & 2848     ( 5.9           )& 961       ( 6.1           )& 1887     ( 5.9           )\\
    extramural venous invasion (EMVI)           & 753      ( 1.6           )& 408       ( 2.6           )& 345      ( 1.1           )\\
    intramural venous invasion (IMVI)           & 220      ( 0.5           )& 118       ( 0.8           )& 102      ( 0.3           )\\
    no EMVI, intramural unknown        & 2505     ( 5.2           )& 843       ( 5.4           )& 1662     ( 5.2           )\\
    not applicable             & 6        ( 0.0             )& 3         ( 0.0             )& 3        ( 0             )\\
    unknown        & 667      ( 1.4           )& 245       ( 1.6           )& 422      ( 1.3           )\\
    not recorded        & 40936    ( 85.4          )& 13134     ( 83.6          )& 27802    ( 86.3          )\\
    \textbf{High risk*} \\
    high                & 11532    ( 24.1          )& 1638      ( 10.4          )& 9894     ( 30.7          )\\
    normal               & 643      ( 1.3           )& 12        ( 0.1           )& 631      ( 2.0             )\\
    unknown             & 13626    ( 28.4          )& 464       ( 3.0             )& 13162    ( 40.8          )\\
    not applicable      & 22134    ( 46.2          )& 13598     ( 86.5          )& 8536     ( 26.5          )\\
    \textbf{Maximum length of stay}   & 9.4      ( 20.8          )& 8.3       ( 31.3          )& 9.9      ( 13.0            )\\ 
    \hline 
    \end{tabular}
    \\ Data is mean (sd) or numbers (percentage) 
\end{table}


\begin{table}[h!]
    \centering
    \caption{Methods description}
    \begin{tabular}{ l c c c c c}
    \hline
    \textbf{Methods}  & \multicolumn{5}{c}{\textbf{Models}} \\
    \cline{2-6}
    & $\mathbb{E}(Y_{i}|Z_{i})$ & \multicolumn{2}{c}{$\mathbb{E}(Y_{i}|X_{i},Z_{i})$} & \multicolumn{2}{c}{$\mathbb{E}(Z_{i}|X_{i})$} \\
    \cline{2-6}
    &LReg&LReg&BART&LReg&BART\\   
    \cline{2-5}
    \hline
    Crude          & \cmark & -         & -         &-      &-             \\
    Direct/LReg    &  -     & \cmark    & -         &-      &-             \\
    Direct/BART    &  -     & -         & \cmark    &-      & -            \\
    IPW/LReg       &  -     & -         & -         &\cmark & -            \\
    IPW/BART       &  -     & -         & -         &-      & \cmark       \\
    DR/LReg-LReg   &  -     & \cmark    & -         &\cmark & -            \\
    DR/LReg-BART   &  -     & \cmark    & -         &-      & \cmark       \\
    DR/BART-LReg   &  -     & -         & \cmark    &\cmark & -            \\
    DR/BART-BART   &  -     & -         & \cmark    &-      & \cmark       \\
    Direct/ps-BART &  -     & -         & \cmark    &-      & \cmark       \\
    \hline
    \end{tabular}
    \label{table:methods}
\end{table}


\begin{table}[h!]
    \centering
    \caption{Baseline characteristics of the study population of in the NCR after subsampling and QUASAR}
    \label{tab:subsampling matching}
    \begin{tabular}{l l l l} 
    \hline
    \textbf{Covariates}&\makecell[l]{\textbf{the two-year follow-up}\\ \textbf{NCR after subsampling}}&\makecell[l]{\textbf{the five-year follow-up}\\ \textbf{NCR after subsampling}} &\textbf{QUASAR} \\
    & (n=7704)&(n=6378)&(n=3239) \\[0.5ex]
    \hline 
    \textbf{Stage}\\
    \uppercase\expandafter{\romannumeral2} & 7011(91\%)&5803(91\%) & 2963(91\%)\\ 
    \uppercase\expandafter{\romannumeral3} & 693(9\%) &575(9\%)& 260(9\%)\\
    \textbf{Sex}\\
    Male & 4699(61\%)&3890(61\%)&1979(61\%) \\
    Female & 3005(39\%)&2488(39\%) &1260(39\%) \\
    \textbf{Age}\\
    \textless 50 &847(11\%)&702(11\%)&370(11\%) \\ 
    50-59 & 1996(26\%))&1651(26\%)&855(26\%)\\
    60-69 & 3236(42\%)&2679(42\%)&1351(42\%)\\
    \textgreater 70 & 1625(21\%)&1632(21\%) & 663(21\%)\\ 
    \hline
    \end{tabular}
\end{table}


\begin{table}[h!]
    \begin{threeparttable}
    \centering
    \caption{Stratified sampling profile}
    \label{tab: Stratified sampling profile}
    \begin{tabular}{l l l l l l} 
    \hline
    \textbf{index}&\textbf{stage; sex; age} &\makecell{\textbf{\# people}\\ \textbf{in \textit{Stratum}}\\ (N= 45442)}  &\textbf{P(\textit{Stratum})} & \makecell{\textbf{\# people}\\ \textbf{in \textit{STRATUM}}\\(N=7704)} &\textbf{P(\textit{STRATUM})}\\[0.5ex]
    \hline
    1&\rom{2}; male; \textless 50   &  470            &0.010 &    470 &           0.061 \\[1pt]
    2&\rom{2}; male; 50-59          &  1259           &0.028 &   1110 &           0.144 \\[1pt]
    3&\rom{2}; male; 60-69          &  3740           &0.082 &   1795 &           0.233 \\[1pt]
    4&\rom{2}; male; 70+            &  7010           &0.154 &    901 &           0.117 \\[1pt]
    5&\rom{2}; female; \textless 50 &   476           &0.010 &    300 &           0.039 \\[1pt]
    6&\rom{2}; female; 50-59        &  1163           &0.026 &    709 &           0.092 \\[1pt]
    7&\rom{2}; female; 60-69        &  3036           &0.067 &   1148 &           0.149 \\[1pt]
    8&\rom{2}; female; 70+          &  7148           &0.157 &    578 &           0.075 \\[1pt]
    9&\rom{3}; male; \textless 50   &   557           &0.012 &     46 &           0.006 \\[1pt]
    10&\rom{3}; male; 50-59         &  1385           &0.030 &    108 &           0.014 \\[1pt]
    11&\rom{3}; male; 60-69         &  3614           &0.080 &    177 &           0.023 \\[1pt]
    12&\rom{3}; male; 70+           &  5407           &0.119 &     92 &           0.012 \\[1pt]
    13&\rom{3}; female; \textless 50&   574           &0.013 &     31 &           0.004 \\[1pt]
    14&\rom{3}; female; 50-59       & 1323            &0.029 &     69 &           0.009 \\[1pt]
    15&\rom{3}; female; 60-69       & 2961            &0.065 &    116 &           0.015 \\[1pt]
    16&\rom{3}; female; 70+         & 5319            &0.117 &     54 &           0.007 \\[1pt]     
    \hline 
    \end{tabular}
    P(\textit{Stratum}): the probability of stratum in the population \textit{P} with sample size $45485$ \\
    P(\textit{STRATUM}): the probability of stratum in the study population \textit{S} with sample size N=7704, which is consistent with that of QUASAR \\
    For the first \textit{STRATUM}, $P(STRATUM)= P(stage\rom{2} \cap  male \cap age\textless50)= P(stage \rom{2}) \times P(male) \times P(age\textless50)=
    0.91 \times 0.61 \times 0.11=0.061 $ (see Table 3 for more details of distributions of these covariates in the QUASAR).
\end{threeparttable}
\end{table}


\begin{table}[h!]
\centering
\caption{Results for 2-year follow-up}
\label{tab:Results2}
\begin{tabular}{l l c c}
\hline
covariates &Methods  & ATE(range) &overlap(\%) \\
\hline
\multirow{10}{*}{3 cov}&Crude($X_{sub} \sim \mathcal{D}_{QUASAR'}$)&-0.019 (-0.041,  0.003)& 99.48\\
&     Direct/LReg &-0.029  (-0.03, -0.029)&   100\\
&     Direct/BART  &-0.02  (-0.02, -0.019)&   100\\
&        IPW/LReg &-0.029 (-0.057,  0.006)& 89.16\\
&        IPW/BART &-0.034  (-0.06,      0)& 76.72\\
&    DR/LReg-LReg &-0.028 (-0.053,  0.008)&  92.9\\
&    DR/LReg-BART &-0.028 (-0.051,  0.004)& 94.44\\
&    DR/BART-LReg &-0.021 (-0.046,  0.016)& 98.72\\
&    DR/BART-BART &-0.022 (-0.045,   0.01)& 99.14\\
& Direct/ps-BART &-0.022 (-0.022, -0.022)&   100\\
\hline
\multirow{10}{*}{5 cov}&Crude($X_{sub} \sim \mathcal{D}_{QUASAR'}$)&-0.019 (-0.041,  0.006)& 99.36\\
&    Direct/LReg &-0.047 (-0.047, -0.046)&     0\\
&    Direct/BART &-0.034 (-0.035, -0.034)&   100\\
&       IPW/LReg  &-0.04 (-0.068, -0.013)& 47.54\\
&       IPW/BART &-0.046 (-0.074, -0.022)& 19.02\\
&   DR/LReg-LReg &-0.044 (-0.068, -0.022)& 21.92\\
&   DR/LReg-BART &-0.043 (-0.066, -0.023)& 25.96\\
&   DR/BART-LReg  &-0.04 (-0.064, -0.018)& 44.22\\
&   DR/BART-BART &-0.042 (-0.064, -0.022)& 34.22\\
& Direct/ps-BART &-0.041 (-0.041,  -0.04)&     0\\
\hline
\multirow{10}{*}{13 cov}&Crude($X_{sub} \sim \mathcal{D}_{QUASAR'}$)&-0.019 (-0.041,  0.006)& 99.36\\
&    Direct/LReg &-0.055 (-0.057, -0.054)&     0\\
&    Direct/BART &-0.056 (-0.057, -0.055)&     0\\
&       IPW/LReg &-0.046 (-0.073,  -0.02)& 20.18\\
&       IPW/BART &-0.053  (-0.08, -0.027)&  4.24\\
&   DR/LReg-LReg  &-0.05 (-0.074, -0.029)&  4.44\\
&   DR/LReg-BART  &-0.05 (-0.073, -0.031)&  3.02\\
&   DR/BART-LReg  &-0.05 (-0.072,  -0.03)&  4.46\\
&   DR/BART-BART &-0.052 (-0.074, -0.033)&  1.14\\
& Direct/ps-BART &-0.051 (-0.052,  -0.05)&     0\\
\hline
\multirow{10}{*}{all cov}&Crude($X_{sub} \sim \mathcal{D}_{QUASAR'}$)&-0.019 (-0.041,  0.006)& 99.36\\
&    Direct/LReg &-0.053 (-0.055, -0.052)&     0\\
&    Direct/BART &-0.045 (-0.046, -0.044)&     0\\
&       IPW/LReg &-0.043 (-0.072, -0.017)& 31.16\\
&       IPW/BART &-0.049 (-0.075, -0.021)& 11.54\\
&   DR/LReg-LReg &-0.046 (-0.069, -0.026)& 12.32\\
&   DR/LReg-BART &-0.046 (-0.068, -0.025)& 12.46\\
&   DR/BART-LReg &-0.045 (-0.065, -0.026)& 17.58\\
&   DR/BART-BART &-0.048 (-0.068, -0.028)&   7.3\\
& Direct/ps-BART &-0.053 (-0.055, -0.052)&     0\\
\hline       
\end{tabular}
\\RCT:-0.021 (95\% CI -0.039, -0.002), Crude ($X_{NCR} \sim \mathcal{D}_{NCR}$): -0.069 (95\% CI -0.076, -0.062)
\end{table}


\begin{table}[h!]
\centering
\caption{Results for 5-year follow-up}
\label{tab:Results5}
\begin{tabular}{l l c c}
\hline
covariates &Methods  & ATE(range) &overlap(\%) \\
\hline
\multirow{10}{*}{3 cov}&Crude($X_{sub} \sim \mathcal{D}_{QUASAR'}$)&-0.031 (-0.065, 0.005)& 99.76\\
  &     Direct/LReg &-0.035 (-0.036, -0.034)&   100\\
  &     Direct/BART &-0.025 (-0.026, -0.025)&   100\\
  &        IPW/LReg  &-0.05 (-0.095, -0.011)&  78.8\\
  &        IPW/BART  &-0.06 (-0.105,  -0.02)&  46.7\\
  &    DR/LReg-LReg &-0.041 (-0.085, -0.002)&  94.8\\
  &    DR/LReg-BART &-0.038 (-0.079, -0.003)& 98.18\\
  &    DR/BART-LReg &-0.034 (-0.078,  0.005)& 98.66\\
  &    DR/BART-BART &-0.032 (-0.072,  0.003)& 99.54\\
  &  Direct/ps-BART &-0.021 (-0.022, -0.021)&   100\\
  \hline
  \multirow{10}{*}{5 cov} &Crude($X_{sub} \sim \mathcal{D}_{QUASAR'}$)&-0.031  (-0.07,  0.008)& 99.76\\
  &    Direct/LReg &-0.065 (-0.066, -0.065)&     0\\
  &    Direct/BART &-0.053 (-0.053, -0.052)&   100\\
  &       IPW/LReg &-0.064 (-0.103, -0.026)& 33.14\\
  &       IPW/BART &-0.075 (-0.114, -0.039)&  7.94\\
  &   DR/LReg-LReg &-0.065 (-0.097,  -0.03)& 27.14\\
  &   DR/LReg-BART &-0.061 (-0.093, -0.029)& 39.66\\
  &   DR/BART-LReg  &-0.06 (-0.092, -0.026)& 44.62\\
  &   DR/BART-BART &-0.058 (-0.089, -0.025)&  53.8\\
  & Direct/ps-BART &-0.049 (-0.049, -0.049)&   100\\
  \hline
  \multirow{10}{*}{13 cov} &Crude($X_{sub} \sim \mathcal{D}_{QUASAR'}$)&-0.031  (-0.07,  0.008)& 99.76\\
  &    Direct/LReg &-0.081 (-0.082, -0.078)&     0\\
  &    Direct/BART &-0.093 (-0.095, -0.092)&     0\\
  &       IPW/LReg  &-0.08 (-0.121, -0.043)&  3.92\\
  &       IPW/BART  &-0.09 (-0.129, -0.054)&  0.24\\
  &    DR/LReg-LReg &-0.078 (-0.114, -0.049)&  1.42\\
  &    DR/LReg-BART &-0.077 (-0.111,  -0.05)&  1.54\\
  &    DR/BART-LReg &-0.077 (-0.112, -0.048)&   2.2\\
  &    DR/BART-BART &-0.077 (-0.111, -0.049)&  1.52\\
  & Direct/ps-BART  &-0.08 (-0.081, -0.078)&     0\\
  \hline
  \multirow{10}{*}{all cov} &Crude($X_{sub} \sim \mathcal{D}_{QUASAR'}$)&-0.031  (-0.07,  0.008)& 99.76\\
  &    Direct/LReg &-0.078  (-0.08, -0.076)&     0\\
  &    Direct/BART  &-0.07 (-0.071, -0.069)&     0\\
  &       IPW/LReg &-0.078 (-0.117, -0.041)&  5.66\\
  &       IPW/BART &-0.086 (-0.124, -0.052)&  0.76\\
  &   DR/LReg-LReg &-0.075 (-0.109, -0.046)&   3.4\\
  &   DR/LReg-BART &-0.073 (-0.105, -0.045)&  5.22\\
  &   DR/BART-LReg &-0.073 (-0.107, -0.044)&  4.82\\
  &   DR/BART-BART &-0.073 (-0.106, -0.044)&  4.34\\
  & Direct/ps-BART &-0.076 (-0.078, -0.075)&     0\\  
\hline       
\end{tabular}
\\RCT:-0.025 (95\% CI -0.059, 0.008),Crude ($X_{NCR} \sim \mathcal{D}_{NCR}$): -0.073(95\% CI -0.083, -0.063)
\end{table}








%\subsection*{The summary of chemotherapies which patients in the NCR received}
%  \label{app: 1}
%  \begin{table}[H]
%      \caption{The summary of chemotherapies which patients in the NCR received}
%      \begin{tabular}{ l l l l}
%          \hline
%          \textbf{ID} &\textbf{chemotherapy regimen}                  & \textbf{number}  &\textbf{description} \\
%          \hline                  
%          1 &no chemotherapy                                             &   34073          &\\
%          2 &        systemic chemotherapy                               &   9993           &\\    
%          3 &platinum-containing chemotherapy                            &   5057           &\\
%          4 &                        CAPOX                               &   2072           &capecitabine and oxalipatin\\
%          5 &                Capecitabine                                &   638            \\
%          6 &                        FOLFOX                              &   244            &5-FU,folinic acid, and exaliplatin\\
%          7 &\makecell[l]{chemotherapy with levamisole \\or leucovorin / leervorin} & 241            &\\
%          8 &                taxane-containing chemotherapy              &   55             &\\
%          9 &                                        5-FU                &   40             &\\
%          10 &                            chemotherapy nno                &   33             &\\
%          11 &                                 oxaliplatin                &   18             &\\
%          12 &    chemotherapy with angiogenesis inhibitor                &   16             &\\
%          13 &                                  irinotecan                &   14             &\\
%          14 &                                    FOLFIRI                 &   5              &\makecell[l]{5-FU,folinic acid,\\ and irinotecan}\\
%          15 &                    etoposide + carboplatin                 &   4              &\\
%          16 &    anthracycline-containing chemotherapy                   &   3              &\\
%          17 &                     etoposide + cisplatine                 &   3              &\\
%          18 &                                FOLFOXIRI                   &   2              &\makecell[l]{5-FU, folinic acid, oxaliplatin,\\ and irinotecan}\\
%          19 &\makecell[l]{Tegafur + gimeracil + \\oteracil (S1,Teysuno)}         &   2              &\\
%          20 &\makecell[l]{anthracycline 357 277 275n \\taxane-containing chemotherapy}  &   2          &\\
%          21 &                                    CAPIRI                  &   1              &capecitabine and irinotecan\\
%          22 &                                Carboplatin                 &   1              &\\
%          23 &            Paclitaxel + capecitabine/xeloda                &   1              &\\
%          24 &                         Uracil and tegafur                 &   1              &\\
%          25 &            chemotherapy with COX2 inhibitor                &   1              &\\
%          26 &                                epirubicin                  &   1              &\\
%      \hline
%      \end{tabular}
%      \label{tab: chemotherapy regimen}
% \end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                               %%
%% Additional Files              %%
%%                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\clearpage
\section*{Appendix}
\label{appendix}

\subsection*{Details regarding the subsampling procedure}
\label{app:subsampling}
Here we describe our stratified (sub) sampling approach to estimate the sampling distribution of $\mathcal{D}_{QUASAR'}$ by correcting for the sampling mechanisms in the NCR as introduced in section \ref{sec: study population sampling} in the context of two-year follow-up analysis. Our aim is by correcting for the sampling mechanism, the subsamples in NCR ($X_{sub} \sim \mathcal{D}_{QUASAR'}$) is similar with QUASAR ($X_{QUASAR} \sim \mathcal{D}_{QUASAR}$) in terms of distributions of all the known covariates. 
  
To do so, first, due to the fact that only three covariates, namely, stage, sex, and age and only their marginals were reported in QUASAR, we can compute the joint probability of each stratum according to the marginal probabilities of these three covariates by assuming independence and the joint probability is given by $P(stage \cap sex \cap age) = P(stage) \times P(sex) \times P(age) $ which is presented in the column $P(STRATUM)$ in Table \ref{tab: Stratified sampling profile}.  
  
Second, subsampling the population in NCR according to the joint probability of each stratum such that the probability of each stratum in the study population after subsampling is the same as QUASAR. Here we subsample the target population for two-year follow-up analysis. Let the eligible population in NCR be \textit{P} with sample size $45442$, the study population after sampling be \textit{S} with sample size $N$, the first stratum in \textit{P} be \textit{Stratum1} with sample size $n1$, the first stratum in the \textit{S} be \textit{STRATUM1} with sample size $N1$. To maximize the $N$ of \textit{S}, presume all the units in \textit{Stratum1} (n1=470) are included in \textit{S}, then the sample size of \textit{S} should be $ N=\frac{n1}{P(STRATUM1)} = \frac{470}{0.061} \approx 7704$. Based on this study population ($N=7704$), the sample size of \textit{STRATUM2} should be $N2=N \times P(STRATUM2)=1259$; likewise, the sample size of other \textit{STRATUMs} in the \textit{S} can be obtained in the same way. Under this scenario, all the units in each \textit{STRATUM} in \textit{S} can be sampled from \textit{Stratum} of the population \textit{P} without replacement. Next, we continue another round of stratified sampling in the same way. Different from the first round, this time we presume all the units in \textit{Stratum2} ($n2=1259$) are included in the \textit{S} and the sample size of the \textit{S} should be $ N=\frac{n2}{P(STRATUM2)}=\frac{1259}{0.144}=8743$. However, if we use this scenario, then the sample size of the \textit{STRATUM1} will be $N1 = N \times P(STRATUM1) = 8743 \times 0.061 = 533$, which is much larger than the sample size of \textit{Stratum1} $(n1=470)$. Hence we can't sample the units from \textit{Stratum1} without replacement and we won't use this sampling scenario. Next, we replicate this sampling method for the third time, and include all the samples of \textit{Stratum3} into \textit{S}.  Then the sample size of \textit{S} is equal to $ N=\frac{n3}{P(STRATUM3)} = \frac{3740}{0.233}=16051$. On the basis of the sample size $N$ of \textit{S}, we check if we can stratified sample from each $Stratum$ of $P$ without replacement. We replicate this sampling method in total for 16 times in which each \textit{Stratum} is chosen in turn to generate the sample size of \textit{S}. Finally, only one method is adopted: all the units in \textit{Stratum1} are included into \textit{S} to generate the sample size of \textit{S} ($N=7704$), and subsample from \textit{P} according to $P(STRATUM)$ without replacement on the basis of the sample size of \textit{S}, as details shown in Table \ref{tab: Stratified sampling profile}. After subsampling, the distributions of the stage, sex and age of the study population \textit{S} in the NCR are consistent with the QUASAR, as shown in Table \ref{tab:subsampling matching} in section \ref{sec: study population sampling}. In further analysis, to quantify the uncertainty of sampling procedure, we will subsample the \textit{P} using this stratified sampling method for 5000 times.  

\end{backmatter}
\end{document}
