\section*{Abstract}
Inverse probability weighting is increasingly used in causal inference, but the g-computation constitutes a promising alternative. We aimed to compare the performances of these methods with time-to-event outcomes. Given the limitations of the interpretability of the hazard ratio for causal inference, our target estimand of interest is the difference in restricted mean survival times. We report the findings of an extensive simulation study showing that both inverse probability weighting and g-computation are unbiased under a correct model specification, but g-computation is generally more efficient. We also analyse two real-world datasets to illustrate these results. Finally, we update the {\tt{R}} package {\tt{RISCA}} to facilitate the implementation of g-computation.\\

\section*{Keywords} Causal inference, Parametric g-formula, Propensity score, Restricted mean survival time, Simulation study.



\section{Introduction}
Real-world evidence is scientific evidence obtained from data collected outside the context of randomised clinical trials \cite{Sherman_RWE_2016}. The absence of randomisation complicates the estimation of the marginal causal effect (hereafter referred to merely as causal effect) of exposure (including treatment or intervention) due to a potential risk of confounding \cite{Hernan_2004}. Rosenbaum and Rubin \cite{rosenbaum_central_1983} introduced the propensity score (PS) as a tool for causal inference in the presence of measured confounders. In a binary exposure setting, it has been shown that the estimated PS is a balancing score, meaning that conditional on the estimated PS, the distribution of covariates is similar for exposed and unexposed patients. Following this property, the PS can be used in four ways to provide estimates of the causal exposure effect: matching, stratification, adjustment, and inverse probability weighting (IPW) \cite{robins_marginal_2000}. Stratification leads to residual confounding and adjustment relies on strong modelling assumptions \cite{lunceford_stratification_2004, Vansteelandt_Daniel_2014}. Although matching on PS has long been the most popular \cite{ali_reporting_2015}, IPW appears to be less biased and more precise in several studies \cite{le_borgne_comparisons_2015, Hajage__2016, austin_performance_2013}. Moreover, King and Nielsen \cite{King_PSM} argued for halting the use of PS matching for many reasons, including covariate imbalance, inefficiency, model dependence, and bias. 

Causal effects can also be estimated using the g-computation (GC), a maximum likelihood substitution estimator of the g-formula \cite{Robins_1986, Snowden_2011, Keil_2014}. While IPW is based on exposure modelling, the GC relies on the prediction of the potential outcomes for each subject under each exposure status. 

Several studies (see \cite{Chatton_2020} and references therein) compared the IPW and GC in different contexts. They reported a lower variance for the GC than the IPW. Nevertheless, to the best of our knowledge, no study has focused on time-to-event outcomes. The presence of right-censoring and its magnitude is of prime importance since a small number of observed events due to censoring may impact the estimation of the outcome model involved in the GC. By contrast, the IPW may perform well as long as the number of exposed patients is sufficient to estimate the PS and the total sample size is sufficiently large to limit variability in the estimated weights.
In addition to the lack of simulation-based studies comparing the GC and the IPW in a time-to-event context, the GC-based estimates, such as hazard ratio or restricted mean survival, are not straightforward to obtain. Continuous-time data result in methodological difficulties, mainly due to the time-varying distribution of individual baseline characteristics.

In the present paper, we aimed to detail the statistical framework for using the GC in time-to-event analyses. We restricted our developments to time-invariant confounders because the joint modelling of time-to-event and longitudinal data requires for controlling the time-varying confounding by GC or IPW \cite{Hernan_2001}, which is beyond the scope of this paper. We also compared the performances of the GC and IPW. The rest of this paper is structured as follows. In section 2, we detail the methods. Section 3 presents the design and findings of a simulation study. In section 4, we propose a practical comparison with two real-world applications related to treatment evaluations in multiple sclerosis and kidney transplantation. Finally, we discuss the results and provide practical recommendations to help analysts to choose the appropriate analysis method.




\section{Methods}


\subsection{Notations}


Let $(T_i,\delta_i,A_i,L_i)$ be the random variables associated with subject $i$ ($i=1,...,n)$. $n$ is the sample size, $T_i$ is the participating time, $\delta_i$ is the censoring indicator (0 if right-censoring and 1 otherwise), $A_i$ is a binary time-invariant exposure initiated at time $T=0$ (1 for exposed subjects and 0 otherwise), and $L_i=(L_{1i},...,L_{pi})$ is a vector of the $p$ measured time-invariant confounders. Let $S_a(t)$ be the survival function of group $A=a$ at time $t$, and let $\lambda_a(t)$ be the corresponding instantaneous hazard function. Suppose $D_a$ is the number of different observed times of event in group $A=a$. At time $t_j$ ($j=1,... , D_a$), the number of events $d_{ja}$ and the number of at-risk subjects $Y_{ja}$ in group $A=a$ can be defined as $d_{ja} = \sum_{i:t_i = t_j} \delta_i \mathds{1}(A_i=a)$ and $Y_{ja} = \sum_{i:t_i \geq t_j} \mathds{1}(A_i=a)$.


\subsection{Estimands}
\label{estimand}


The hazard ratio ($HR$) has become the estimand of choice in confounder-adjusted studies with time-to-event outcomes. However, it has also been contested \cite{hernan_hazards_2010, Aalen_2015}, mainly because the baseline distribution of confounders varies over time. For instance, consider a population in which there is no confounder at baseline, \textit{i.e.}, the distribution of $L$ is identical regardless of exposure level $A$. Suppose additionally that the quantitative covariate $L_1$ and the exposure $A$ are independently associated with a higher risk of death. In this situation, as illustrated in Figure \ref{L1}, the difference between the $L_1$ mean values of the two groups increases over time. Then, even when the subject-specific proportional hazard (PH) assumption is true, \textit{i.e.}, the $HR$ between an exposed and an unexposed patient with the same characteristics $L$ is constant over time, the population-average $HR$ varies over time. Instead of $HR$, one can estimate the average over time of the different time-specific $HRs$: $AHR = \int [ \lambda_1(t) / \lambda_0(t) ] f(t)$ \cite{Schemper_2009}.


Nevertheless, Aalen \textit{et al.} \cite{Aalen_2015} concluded that it is difficult to draw causal conclusions from such a relative estimand. Hern\'an\cite{hernan_hazards_2010} advocated the use of the adjusted survival curves and related differences. For instance, the restricted mean survival time (RMST) allows us to summarise a survival curve for a specific time-window and to compare two curves by looking at the difference in RMST \cite{Royston_Parmar_2013}. The RMST difference up to time $\tau$ is formally defined as :
\begin{equation}
\Delta (\tau) = \int_{0}^{\tau} [S_1(t) - S_0(t)]dt
\label{eq_rmst}
\end{equation}

This value corresponds to the difference in terms of mean event-free time between two groups of exposed and unexposed individuals followed up to time $\tau$. A further advantage of the RMST difference is its usefulness for public health decision making \cite{Poole_2010}. Note that other alternatives that might avoid this problem exist, such as the attributable fraction or the number needed to treat \cite{Sjolander_2018}.

Hereafter, we considered $AHR$ and $\Delta (\tau)$.


\subsection{Weighting on the inverse of propensity score}


Formally, the PS is $p_{i} = P(A_i=1|L_i)$, \textit{i.e.}, the probability that subject $i$ is exposed according to her/his characteristics $L_i$. In practice, analysts often use logistic regression. Let $\omega_i$ be the weight of subject $i$. Xu \textit{et al.} \cite{xu_use_2010} defined $\omega_i = A_i P(A_i=1) / p_i + (1-A_i) P(A_i=0) / (1-p_i)$ to obtain a pseudo-population in which the distribution of covariates is balanced between exposure groups, enabling estimation of the causal effect in the entire population \cite{Hernan_2004}. The use of stabilised weights has been shown to produce a suitable estimate of the variance even when there are subjects with extremely large weights \cite{robins_marginal_2000, xu_use_2010}. The weighted numbers of events and at-risk subjects at time $t_j$ in group $A=a$ are $d_{ja}^\omega = \sum_{i:t_i = t_j} \omega_i \delta_i \mathds{1}(A_i=a)$ and $Y_{ja}^\omega = \sum_{i:t_i \geq t_j} \omega_i \mathds{1}(A_i=a)$, respectively. Cole and Hern\'an \cite{cole_adjusted_2004} proposed a weighted Kaplan-Meier estimator defined as:

\begin{equation}
\hat S_a(t)= \prod_{t_j \leq t}  \Big  [ 1- d_{ja}^\omega / Y_{ja}^\omega \Big ]
\label{eq1}
\end{equation}


To estimate the corresponding $AHR$, they suggested the use of a weighted univariate Cox PH model, in which exposure is the single explanatory variable. We use equation (\ref{eq_rmst}) to estimate the corresponding $\Delta(\tau)$.


\subsection{G-computation} 


Akin to the IPW, the GC involves two steps. The first consists of estimating the outcome model, called the Q-model \cite{Snowden_2011}. When suitable, it can consist of a PH regression: $h(t|A_i,L_i) = h_0(t)\exp(\gamma A_i + \beta L_i)$ where $h(t|A_i,L_i)$ is the conditional hazard function of subject $i$ at time $t$, $h_0(t)$ the baseline hazard function, and $\gamma$ and $\beta$ are regression coefficients. Estimates of the cumulative baseline hazard $\hat{H}_0(t)$ and the regression coefficients ($\hat \gamma, \hat \beta$) can be obtained by the joint likelihood approach proposed by Breslow \cite{Breslow_1972}. The second step consists of predicting the counterfactual mean survival function if all subjects would have been exposed ($do(A=1)$) or unexposed ($do(A=0)$):

\begin{equation}
\hat{S}_a(t) = n^{-1} \sum_{i=1}^n \exp \Big [- \hat{H}_0(t) \times \exp(\hat\gamma \times do(A_i=a) + \hat\beta L_i)  \Big ]
\label{eq2}
\end{equation}

Then, $\widehat{AHR}$ can be computed as the mean of the individual counterfactual hazard ratios at the observed event times \cite{Schemper_2009}:

\begin{equation}
\widehat{AHR} = \Big [ \sum_{i=1}^n \delta_i  \Big ]^{-1} \sum_{i=1}^n \delta_i \Big [\hat{\lambda}_1(t_i)/ \hat{\lambda}_0(t_i)  \Big ],
\label{eq3}
\end{equation}
where $\hat{\lambda}_a(t) = - \partial \log\hat{S}_a(t)/\partial t$, which is obtained from equation (\ref{eq2}) by numerical differentiation. We use equation (\ref{eq_rmst}) to estimate the corresponding $\Delta(\tau)$.


\subsection{Software} 

We performed all the analyses using {\tt{R}} version 3.6.0. \cite{R} Source code to reproduce the results is available as Supporting Information on the journal’s web page. To facilitate their use in practice, we have implemented the previous methods in the {\tt{R}} package entitled {\tt{RISCA}} (versions $\geq$ 0.8.1), which is available at {\tt{cran.r-project.org}}. We obtained the variances by bootstrap, as recently recommended by Austin \cite{Austin_boot}. 


\subsection{Identifiability conditions}


As for standard regression models, the IPW and the GC require assumptions of non-informative censoring, no measurement error, no model misspecification, and no interference \cite{Hudgens_Halloran_2008}. Three additional assumptions, called \emph{identifiability conditions}, are necessary for causal inference \cite{Naimi_2016}. (i) The values of exposure under comparisons correspond to well-defined interventions that, in turn, correspond to the versions of exposure in the data. (ii) The conditional probability of receiving every value of exposure depends only on the measured covariates. (iii) The conditional probability of receiving every value of exposure is greater than zero. These assumptions are known as \emph{consistency}, \emph{conditional exchangeability} and \emph{positivity}, respectively. 


\section{Simulation study}


\subsection{Data generation}


We generated data in three steps. (a) We simulated three covariates ($L_1$ to $L_3$) from a Bernoulli distribution with parameter equal to 0.5 and three covariates ($L_4$ to $L_6$) from a standard normal distribution (Figure \ref{DAG}). (b) We generated the exposure $A$ according to a Bernoulli distribution with probability obtained by the logistic model with the following linear predictor: $-0.5 + \log(2) \cdot L_2 + \log(1.5) L_3 +\log(1.5) L_5 +\log(2) L_6$. We set the intercept to obtain the prevalence of exposed individuals at 50\%. (c) We generated the times-to-event from a Weibull PH model. We set the scale and shape parameters to 40.0 and 2.0, respectively. Based on a random variable $U_i$ drawn from a standard uniform distribution, we then computed the time-to-event from a Weibull PH model as $40.0 * [(1 - \log(1-U_i) \exp(-\gamma A_i - \log(1.8) L_1 - \log(1.3) L_2 - \log(1.8) L_4 - \log(1.3) L_5 )) - 1]^{-2.0}$, where $\gamma = \log(1.0)$ under the null hypothesis or $\log(1.3)$ under the alternative hypothesis. We subsequently censored the times-to-event using a uniform distribution on [0,70] or [0,15], leading to approximately 40\% and 90\% censored observations, respectively. 

We explored three sample sizes: $n=$ 100, 500, and 2000. When the censoring rate was 90\%, we did not investigate the smallest sample size due to the reduced number of events. For each scenario, we randomly generated 10,000 datasets. To compute the difference in $\Delta (\tau)$, we defined $\tau$ in each dataset as the time at which at least 10\% of the individuals in each group (exposed or unexposed) were still at risk. 

We furthermore considered two sets of covariates: $(L_1, L_2, L_4, L_5)$ the risk factors of the outcome, or $(L_2, L_5)$ the true confounders.


\subsection{Performance criteria}


We computed the theoretical values of the $AHR$ and $\Delta(\tau)$ by averaging the estimations obtained, respectively, from univariate Cox PH models ($A$ as the only explanatory covariate) and by equation (\ref{eq_rmst}) where the survival functions were estimated by the Kaplan-Meier estimator, fitted from datasets simulated as above, except $A$ was simulated independently of $L$ \cite{gayat_propensity_2012}. We reported the following criteria: (i) the percentage of datasets without convergence ; (ii) the mean absolute bias $MAB=E(\hat \theta)-\theta$, where $\theta$ is the estimand of interest; (iii) the mean square error $MSE = E[(\hat\theta - \theta)^2]$; (iv) the variance estimation bias $VEB = 100 \times (E[\widehat{SD}(\hat\theta)]/SD(\hat\theta) - 1)$, where $\widehat{SD}(\bullet)$ is the asymptotic standard deviation and $SD(\bullet)$ is the empirical standard deviation; (v) the empirical coverage rate of the nominal 95\% confidence interval (CI), defined as the percentage of 95\% CIs including $\theta$; (vi) the type I error, defined as the percentage of times the null hypothesis is rejected when the null hypothesis is true; and (vii) the statistical power, defined as the percentage of times the null hypothesis is rejected when the alternative hypothesis is true. We computed the Monte Carlo standard errors for each performance measure \cite{Morris_simulation}.


\subsection{Results}


The Monte Carlo errors were weak, and we did not encounter any convergence problems. Figures \ref{resAHR} and \ref{resRMST} present the results under the alternative hypotheses for $AHR$ and $\Delta$, respectively. The results under the null hypothesis were comparable and can be found in the supplementary material available online.


The MAB values associated with IPW and GC were similar to zero in all scenarios, regardless of the considered covariate set (all the risk factors or only the true confounders). Nevertheless, the MAB under the alternative hypothesis was slightly lower with the GC when considering all the risk factors. For instance, in terms of $\Delta$ when $n=500$ with a censoring rate of 40\%, the MAB was 0.033 for GC versus 0.061 for IPW. When considering only true confounders, these values were 0.147 and 0.066, respectively.


The GC, when considering all outcome causes, produced the best results in terms of MSE, especially for small sample sizes. For instance, when $n=100$ with a censoring rate of 40\%, the MSE related to the $AHR$ was 0.056 for GC versus 0.068 for IPW. When considering only true confounders, these values were 0.076 and 0.079, respectively. 


Regarding the VEB, the results were close for the GC and IPW. One can nevertheless note lower values for the IPW when $n=100$ to estimate the $AHR$. For instance, for a censoring rate of 40\%, the VEB values were between 4.1\% and 4.7\% for GC versus 2.4\% and 3.7\% for IPW.


All scenarios broadly respected the nominal coverage value of 95\%. The power was the highest for the GC, especially when considering all the risk factors, with a gain between 10\% and 13\%.

With a censoring rate of 90\%, the two methods produced similar results, irrespective of the considered scenario.


\section{Applications}


We used data from two studies performed for multiple sclerosis and for kidney transplantation \cite{Laplaud_2019,Masset_2019}. We conducted these studies following the French law relative to clinical noninterventional research. Written informed consent was obtained. Moreover, the French commission for data protection approved the collection (CNIL decisions DR-2014-327 and 914184). To guide variable selection, we asked experts which covariates were causes of the exposure or the outcome prognosis to define the causal structure \cite{VanderWeele_Shpitser_2011}. We checked the positivity assumption and the considered covariates balance (see supplementary materials available online). The log-linearity hypothesis of continuous covariates was confirmed in the univariate analysis if the Bayesian information criterion \cite{BIC} was not reduced using natural spline transformation compared to the inclusion of the covariate in its natural scale. In case of violation, we used a natural spline transformation. We also assessed the PH assumption via the Grambsch-Therneau test at a significance level of 5\% \cite{Grambsch_Therneau_1994}. For simplicity, we performed complete case analyses.


\subsection{Dimethylfumarate versus Teriflunomide to prevent relapse in multiple sclerosis}


With the increasing number of available drugs for preventing relapses in multiple sclerosis and the lack of head-to-head randomised clinical trials, Laplaud \textit{et al.} \cite{Laplaud_2019} aimed to compare Teriflunomide (TRF) and Dimethylfumarate (DMF) using data from the multicentric cohort OFSEP. We reanalysed the primary outcome, defined as the time-to-first relapse. We presented the cohort characteristics of 1770 included patients in Table 1: 1057 patients were in the DMF group (59.7\%) versus 713 in the TRF group (40.3\%). Approximately 39\% of patients (40\% in the DMF group versus 38\% in the TRF group) had at least one relapse during follow-up. 


We presented the confounder-adjusted results in the left panel of Figure \ref{forCVC}. The GC and IPW were equally robust to the considered set of covariates. For the difference in RMST, the width of the 95\% CI was lower for the GC. For instance, when we considered all the risk factors, the CI of GC had a width of 44.2 days versus 50.9 days for IPW.


The conclusion of no significant difference between TRF and DMF was unaffected by the method or the choice of the estimand. 


\subsection{Basiliximab versus Thymoglobulin to prevent post-transplant complications}


Amongst non-immunised kidney transplant recipients, one can expect similar rejection risk between Thymoglobulin (ATG) and Basiliximab (BSX), two possible immunosuppressive drugs proposed as induction therapy. However, ATG may be associated with higher serious events, especially in the elderly. We aimed to reassess the difference in cardiovascular complications in ATG versus BSX patients \cite{Masset_2019}. Table 2 describes the 383 included patients from the multicentric DIVAT cohort: 204 patients were in the BSX group (53.3\%) versus 179 in the ATG group (46.7\%). Approximately 30\% of patients (29\% in the BSX group and 31\% in the ATG group) had a least one cardiovascular complication during follow-up. The median follow-up time was 1.8 years (min: 0.0; max: 8.2).


In the right panel of Figure \ref{forCVC}, we present the confounder-adjusted RMST differences for a cohort followed up to 3 years. First, the GC-based results obtained were slightly less sensitive to the considered set of covariates than those obtained by the IPW. Indeed, GC resulted in a $\hat{\Delta}(3)$ of 2.6 (IC 95\% from -1.2 to 7.0) and 4.0 (IC 95\% from 0.2 to 8.5) months with all the risk factors and only the true confounders, respectively. In contrast, the IPW resulted in a $\hat{\Delta}(3)$ of 1.5 (IC 95\% from -2.5 to 5.8) and 3.4 (IC 95\% from -0.5 to 7.9) months, respectively. Likewise, the log of $\widehat{AHR}$ obtained by GC was between -0.272 (IC 95\% from -0.703 to 0.123) and -0.440 (IC 95\% from -0.898 to -0.021) versus -0.188 (IC 95\% from -0.639 to 0.226) and -0.384 (IC 95\% from -0.830 to 0.002) for IPW. This variability in the estimations, according to the set of covariates, illustrates the importance of this choice. Second, we obtained the smallest 95\% CIs using the GC with all risk factors, in line with the corresponding lower variance previously reported in the simulations. Third, the conclusion differed depending on the method: the IPW-based results did not point to a significant difference in effect between BSX and ATG, in contrast to the GC with the true confounders.


\section{Discussion}


We aimed to explain and compare the performances of the GC and IPW for estimating causal effects in time-to-event analyses with time-invariant confounders. We focused on the average HR and the RMST difference. The results of the simulations showed that the two methods performed similarly in terms of MAB, VEB, coverage rate, and type I error rate. Nevertheless, the GC outperformed the IPW in terms of statistical power, even when the censoring rate was high. Furthermore, the simulations showed that the Q-model in the GC approach should preferentially include all the risk factors to ensure a smaller mean bias. Therefore, the main advantage of using the GC is the gain in statistical power. 

In the two applications, the 95\% CIs were also narrower when using the GC with all the risk factors. Moreover, while the first application in multiple sclerosis did not highlight relevant differences between the GC and IPW, the second one in kidney transplantation illustrated the importance of the set of covariates to consider. On the basis of the GC with all risk factors, we concluded there was a significant difference between Basiliximab and Thymoglobulin. In contrast, when using the GC with only the true confounders, we reported a non-significant difference. Two arguments can explain why the consideration of all risk factors in the Q-model was associated with the best performance. First, this approach reduces the residual variance and increases the precision of the predictions of potential outcomes. Second, even if a risk factor is balanced between the exposure groups at baseline, it could become unbalanced over time (as illustrated in Figure~\ref{L1}). 


Nevertheless, the higher power of the GC is counterbalanced by the inability of GC to evaluate balance of characteristics between exposure groups over time and the need for bootstrapping to estimate the variance, analytic estimators that are available for the IPW \cite{Conner_Trinquart_2019, Hajage_2018}. In practice, we must emphasise that bootstrapping the entire estimation procedure has the advantage of valid post-selection inference \cite{Efron_2014}. For instance, data-driven methods for variables selection, such as the super learner \cite{vdl_2007}, have recently been developed and may represent a promising perspective when such clinical knowledge is unavailable \cite{Blakely_2019}. 


The IPW and GC are not the only available methods to estimate the causal effect. For instance, Conner \textit{et al.} \cite{Conner_Trinquart_2019} compared the performance of IPW with that of other regression-based methods. Overall, the statistical performance was similar. However, the advantage of the IPW and GC compared to other methods is the visualisation of the confounder-adjusted results in terms of the survival curve or an indicator such as RMST. The use of doubly robust estimators (DREs) could be an extension of our work \cite{Benkeser_2018, Cai_vdl_2019}. DREs combine both the PS and GC approaches for a consistent estimation of the exposure effect if a least one model has been correctly specified, circumventing the aforementioned potential model misspecification \cite{Lendle_2013}. Unfortunately, DREs can increase bias when the two models are misspecified \cite{Kang_Schafer_2007}. Moreover, several studies have reported a lower variance for the GC than DREs \cite{Chatton_2020, Kang_Schafer_2007, Lendle_2013}. The use of the GC also represents a partial solution to prevent the selection of instrumental variables \cite{Myers_2011} since it is independent of the exposure modelling.


Our study has several limitations. First, the results of the simulations and applications are not a theoretical argument for generalising to all situations. Second, we studied only logistic and Cox PH regression: other models could be applied. Keil and Edwards \cite{Keil_Edwards_2018} proposed a review of possible Q-models with a time-to-event outcome. Third, we considered only a reduced number of covariates, which could explain the abovementioned equivalence between the GC and the IPW with the extreme censoring rate. Last, we did not consider competing events or time-varying confounders that require specific estimation methods \cite{Young_2020, Daniel_2013}.


To conclude, by means simulation and two applications on real
datasets, this study tended to show the higher power of the
GC compared to IPW to estimate the causal effect with time-to-event outcomes. All the risk factors should be considered in the Q-model. Our work is a continuation of the emerging literature that questions the near-exclusivity of propensity score-based methods in causal inference.

\section*{Acknowledgements}
The authors would like to thank the members of DIVAT and OFSEP Groups for their involvement in the study, the physicians who helped recruit patients and all patients who participated in this study.  We also thank the clinical research associates who participated in the data collection. We also thank David Laplaud and Magali Giral for their clinical expertise. The analysis and interpretation of these data are the responsibility of the authors. This work was partially supported by a public grant overseen by the French National Research Agency (ANR) to create the Common Laboratory RISCA (Research in Informatic and Statistic for Cohort Analyses, {\tt{www.labcom-risca.com}}, reference: ANR-16-LCV1-0003-01) involving the development of Plug-Stat software. 


\section*{Competing interest}
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Arthur Chatton obtained a grant from IDBC for this work.
Other authors received no financial support for the research, authorship, and/or publication of this article.


\section*{Supplemental material and code}
Supplemental material for this article is available online.


\bibliographystyle{vancouver}
