# Automated causal inference in application to randomized controlled clinical
trials

Ji Q. Wu  Department of Pathology and Molecular Pathology, University
Hospital, University of Zurich, Zurich, Switzerland.  Corresponding authors,
<EMAIL_ADDRESS><EMAIL_ADDRESS>Nanda Horeweg  Department of
Radiation Oncology, Leiden University Medical Center, Leiden, The Netherlands.
Marco de Bruyn  Department of Obstetrics and Gynecology, University of
Groningen, University Medical Center Groningen, Groningen, The Netherlands.
Remi A. Nout  Department of Radiation Oncology, Leiden University Medical
Center, Leiden, The Netherlands.  Currently employed at department of
Radiotherapy, Erasmus MC Cancer Institute, University Medical Center
Rotterdam, Rotterdam, The Netherlands.  Ina M. Jürgenliemk-Schulz  Department
of Radiation Oncology, University Medical Center Utrecht, Utrecht, The
Netherlands.  Ludy C.H.W. Lutgens  Maastricht Radiation Oncology Clinic,
Maastricht, The Netherlands.  Jan J. Jobsen  Department of Radiotherapy,
Medisch Spectrum Twente, Enschede, The Netherlands.  Currently employed at
department of Clinical Epidemiology, Medisch Spectrum Twente, Enschede, The
Netherlands.  Elzbieta M. van der Steen-Banasik  Radiotherapiegroep, Arnhem,
The Netherlands.  Hans W. Nijman  Department of Obstetrics and Gynecology,
University of Groningen, University Medical Center Groningen, Groningen, The
Netherlands.  Vincent T.H.B.M. Smit  Department of Pathology, Leiden
University Medical Center, Leiden, The Netherlands.  Tjalling Bosse
Department of Pathology, Leiden University Medical Center, Leiden, The
Netherlands.  Carien L. Creutzberg  Department of Radiation Oncology, Leiden
University Medical Center, Leiden, The Netherlands.  Viktor H. Koelzer
Department of Pathology and Molecular Pathology, University Hospital,
University of Zurich, Zurich, Switzerland.  Corresponding authors,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Randomized controlled trials (RCTs) are considered as the gold standard for
testing causal hypotheses in the clinical domain. However, the investigation
of prognostic variables of patient outcome in a hypothesized cause-effect
route is not feasible using standard statistical methods. Here, we propose a
new automated causal inference method (AutoCI) built upon the invariant causal
prediction (ICP) framework for the causal re-interpretation of clinical trial
data. Compared to existing methods, we show that the proposed AutoCI allows to
efficiently determine the causal variables with a clear differentiation on two
real-world RCTs of endometrial cancer patients with mature outcome and
extensive clinicopathological and molecular data. This is achieved via
suppressing the causal probability of non-causal variables by a wide margin.
In ablation studies, we further demonstrate that the assignment of causal
probabilities by AutoCI remain consistent in the presence of confounders. In
conclusion, these results confirm the robustness and feasibility of AutoCI for
future applications in real-world clinical analysis.

## Introduction

Many clinical studies are driven by the research questions that are
statistical at first glance but causal by nature [1]. For instance, how safe
and efficient is a vaccine against viral infection [2]? How are
clinicopathological variables related to cancer patient survival [3]? From the
causal perspective, the common ground of these problems starts with
determining the causal variables of the outcome of interest. In clinical
medicine, randomized controlled trials (RCTs) are considered to be the gold
standard to investigate cause-effect relationships [4]. In a prototypical RCT,
a participant is randomly assigned to the experimental or control arm and the
outcome of interest is observed. In the context of causal inference, such a
randomization can be modeled with do-intervention [5].

In this article, we develop a new automated causal inference method (AutoCI)
and apply this to two large-scale, practice changing RCTs of endometrial
carcinoma patients conducted in the Netherlands from 1990-1997 (PORTEC 1 [6,
7]) and 2002-2006 (PORTEC 2 [8, 9]), with full clinicopathological datasets
and mature outcome data. Endometrial carcinoma is the most common type of
gynaecological cancer for women in developed countries [10]. The majority of
women with endometrial cancer (EC) that are diagnosed with early-stage
disease, have a favorable prognosis and are treated with surgery [11].
Approximately 15–20% of patients have an unfavorable prognosis with a high
risk of distant metastasis [11]. For those patients, different adjuvant
therapies, such as vaginal brachytherapy, external beam radiotherapy and
chemoradiation, are recommended based on their risk group [12]. The two trials
(PORTEC 1 and 2) used in this study made a key contribution to clinical
practice by investigating how these therapies impact the risk of recurrence
rates and survival [6, 8]. According to the latest ESGO/ESTRO/ESP guidelines
[12] for the management of patients with endometrial carcinoma, the risk
classification is based on a series of clinical and pathological variables
such as tumor grading (Grade), lymphovascular space invasion (LVSI),
myometrial invasion, etc. and molecular variables including but not limited to
polymerase epsilon mutant EC (POLEmut), mismatch repair deficient (MMRd) EC,
p53 abnormal EC (p53abn) and EC with no specific molecular profile (NSMP)
[12]. In a recent study [3], correlative statistical methods were used to
investigate the hazardous relevance of these variables to EC recurrence. There
is strong evidence to suggest that these variables impact EC recurrence, but a
systematic investigation from a modern causal inference perspective to support
this understanding has not been done.

Causal inference addresses the determination of cause and effect relationships
from data [13, 14, 15, 16]. Given either observational data or by inclusion of
additional interventional data, clinical studies aim to either 1) quantify the
causal effect of a treatment given the outcome [13] or 2) infer the underlying
causal structure of relationships between patient and treatment
characteristics and relevant outcomes [17, 5].

The former can be well formulated as the difference between the outcome
expectations conditioned on different treatments (Average Causal Effect (ACE))
[13]. A wide range of studies have built upon this methodology, including but
not limited to target trial specification [18], target trial emulation [19,
20] and extending inferences from randomized trials to new target populations
[21].

In comparison to the causal effect identification, we refer to 2) as (causal)
structure identification [5]. The goal of structure identification is often to
learn the entire causal structure, _i.e_., a directed acyclic graph composed
of nodes and edges that connect nodes. However, this is in general a non-
deterministic polynomial-time (NP) hard problem [22, 16]. To learn the cause-
effect relations without inferring the entire causal structure, invariant
causal prediction (ICP) [23] was proposed to determine the set of causal
variables given an outcome variable. Under the ICP framework, we utilize the
concept ‘random variable’ (informally, a function that assigns a set of
possible samples to a measurable quantity) and define ‘causal variable’ as
follows (See Tab. 3 in the supplement for summarized terminologies used in the
paper):

######  Definition 1.

Assume there exists a set of environments $u\in U$, given a collection of
random variables $\bm{X}=(X_{1},X_{2},\ldots,X_{n})$ and an outcome variable
$Y$, if there exists a $\bm{X}_{S^{*}}=(X_{S^{*}_{1}},\ldots,X_{S^{*}_{j}})$
with indices $S^{*}\subseteq\\{1,\ldots,n\\}$ such that

$Y=f^{*}(\bm{X}_{S^{*}}^{u})+\delta^{u}\;\forall u\in
U,\;\mathrm{where}\;f^{*}:\mathbb{R}^{|S^{*}|}\mapsto\mathbb{R}.$ (1)
$\delta^{u}\;\mathrm{are}\begin{cases}\text{identically distributed
(i.d.)}&\quad\text{if $\exists$ hidden confounders}\\\ \text{i.d. and
$\delta^{u}\Perp\bm{X}_{S^{*}}^{u}$}&\quad\text{else}\end{cases}$ (2)

Then $\bm{X}_{S^{*}}$ are the Plausible Causal Variables (under $U$). Here,
$\bm{X}_{S^{*}}^{u}=(X_{S^{*}_{1}}^{u},\ldots,X_{S^{*}_{j}}^{u})$ are the
corresponding random variables to
$\bm{X}_{S^{*}}=(X_{S^{*}_{1}},\ldots,X_{S^{*}_{j}})$ created under the
environment $u$. For example, the (experimental) environment $u$ can arise via
do intervention [17] ($do(X_{0}:=c)$) on $X_{0}$, then $X_{0}^{u}$ only
samples the fixed value $c$. It is worth noting that the cause-effect relation
$f^{*}$ in Eq. 1 stays invariant and independent of the environments $U$.
Next, we introduce the definition of identifiable causal variables.

######  Definition 2.

Following the specification of environments $U$, random variables
$\bm{X}=(X_{1},X_{2},\ldots,X_{n})$ and an outcome variable $Y$ in Definition
1, if there exists a
$\bm{X}_{\bar{S}}=(X_{\bar{S}_{1}},\ldots,X_{\bar{S}_{j}})$ with indices
$\bar{S}\subseteq\\{1,\ldots,n\\}$ such that

$\bar{S}:=\bigcap\\{S\subseteq\\{1,\ldots,n\\}\>|\;\bm{X}_{S}\>\text{are
plausible causal variables}\\},$ (3)

then $\bm{X}_{\bar{S}}$ are the Identifiable Causal Variables (under $U$), and
are hence referred to as such.

Initially, vanilla ICP [23] was presented and verified on linear cause-effect
relations. Next, Heinze-Deml  _et al_. [24] defined an environment $u$ as a
random variable that is neither the descendant nor the parent of outcome
variable Y, and conducted multiple conditional independence tests (CIT) on
nonlinear cause-effect settings (NICP). Recently, Gamella  _et al_. [25]
suggested to investigate the stable set of variables instead, which is a
relaxation of the set of identifiable causal variables (AICP). This progress
in ICP has opened unprecedented paths to interpret complex datasets,
especially the ones collected from RCTs. Finding which variables determine
whether a treatment works or whether a patient will have a recurrence using
ICP methods has great scientific potential. In application to the clinical
domain, data interpretation by causal inference methods could improve our
understanding of disease and aid in the design of new experiments and clinical
trials. However, non-negligible efforts are required to adopt the existing
ICPs to the clinical domain. This is due to 1) the complexity and multitude of
variables that are considered relevant for treatment outcomes including
patient level characteristics (patient demographics, text data from clinical
records), information derived from images (radiology, pathology) and molecular
data (genomic sequencing) and 2) the incompatibility between the error-
tolerant implementation for the simulated dataset and safe-critical
application relevant for medical decisions. Candidate ICP methods therefore
need to be robust against noise and need to provide meaningful outputs that
can be related to clinical risk in order to inform patient stratification.

Figure 1: The overall model illustration and performance of the proposed
AutoCI. Top: Illustrative scheme of the proposed AutoCI. Bottom left: The
sampled numbers of type-safe functions vs generic functions. Here the size is
the maximum amount of nn and pred functions allowed during the program
synthesis. Bottom middle: The learning curve of the Jaccard Similarity (JS)
for top-four type-safe functions achieved in the case with pathological,
molecular and immune variables. Bottom right: The running time of determining
the causal variables for pathological variables (P), pathological and
molecular variables (PM), pathological, molecular and immune variables (PMI).
Here the proposed AutoCI utilizes the function
$\mathsf{COMP}(\mathsf{nn},\mathsf{CAT}(\mathsf{FILTER}(\mathsf{pred})))$.

## Results

### Clinical variables overview

The PORTEC 1 and 2 trials [6, 8] recruited 714 (since 1990 - 1997) and 427
(since 2000 - 2006) patients with early stage endometrial carcinoma
respectively. 305 cases from PORTEC 1 (42.7%) and 335 cases from PORTEC 2
(78.5%) with complete clinicopathological datasets were aligned and used in
the experiments. Clinicopathologic characteristics of these subgroups were
similar to the original trial populations (less than or equal to 17.3%
absolute difference in frequency of any variable). Importantly, there was no
statistically significant difference in the variable of interest (mean and
5-year RFS) for causal variable identification between the excluded and
included patient datasets, supporting that our analysis is representative of
the overall study population (Supplementary Fig. 4, Tab. 4, 5). Considering
PORTEC 1 and 2 as the two experimental environments, we aimed to determine the
causal pathological, molecular and immune-related variables of EC recurrence
status.

Pathological variables (P). Pathological criteria tumor grading (Grade) [26,
27], lymphovascular space invasion (LVSI) [28] and myometrial invasion [26,
27] are examined in the study, all of which are important indicators for an
elevated risk of EC recurrence (See also the guideline [12]). All variables
were re-evaluated on formalin-fixed paraffin-embedded (FFPE) tumor material by
specialized gynecopathologists to guarantee variable consistency for the two
environments (trials).

Molecular variables (M). The molecular classification of endometrial cancer
distinguishes four subtypes with validated prognostic impact: i) ultra-mutated
EC with DNA-polymerase epsilon exonuclease domain mutations (POLEmut) with an
excellent prognosis; ii) hypermutated EC with mismatch-repair deficiency
(MMRd) with an intermediate prognosis; iii) copy-number-high EC with frequent
TP53 mutations (p53abn) with an unfavorable prognosis; and iv) copy-number-low
EC without a specific molecular profile (NSMP) with an intermediate prognosis
[29, 27]. Pathogenic POLE mutations were detected by next generation
sequencing of POLE hotspot exons [30]. MMR and p53 status were determined by
immunohistochemistry [31]. Cases with more than one classifying feature were
classified according to the dominant molecular feature on the basis of
pathogenicity [32]. Overexpression of L1CAM by tumor cells was assessed by
immunohistochemistry using a cut-off of $>=10\%$ for positivity, and is
associated with an increased risk of metastasis and death [33, 34].

Immune variable (I). (Intraepithelial) CD8+ T cell infiltration is an
independent favorable prognostic indicator in early stage endometrial cancer
[3]. To quantity CD8+ T cell infiltration in tissue samples of the PORTEC 1
and 2 trials, we compute CD8+ cell density derived from tissue microarrays
(TMA) by immunohistochemistry and image analysis [35, 3]. Specifically, tissue
microarrays capture cancer tissue samples from each patient in a highly
standardized manner, allowing for the highly accurate evaluation of tumor and
microenvironment-related factors in cancer samples for investigation with
clinical outcomes [36].

Based on existing domain knowledge and biological understanding [29, 37, 31,
38, 39, 3], we consider the pathological ($\mathsf{P}$), molecular
($\mathsf{M}$) and immune ($\mathsf{I}$) variables to be the proxies of causal
variables ($S_{\mathsf{prox}}$) (See Tab. 5 in the supplement for more
characteristics details).

Sanity check variables. To investigate the robustness of the causal inference
models, we intentionally include a randomized number as Patient ID and the
vital tissue area of each TMA core (tissue area = sum of randomly sampled
tumor and stroma areas from each case) in our subsequent analysis as non-
causal variables. Based on prior domain knowledge, the variables "Patient ID"
and "Tissue area" should not be causally related to cancer outcomes. The
inclusion of these two non-causal variables thus serves as an important
benchmark for comparison of the methods presented in this study. Importantly,
as the two variables are expected not to impact clinical outcome and we
attempt to verify that they do not impact the outcome, this strategy can also
be considered as a typical use case of negative control [40] in RCT studies.

### PORTEC experiments overview

As presented in Fig. 1 (top), our proposed AutoCI is composed of two key
components: 1) A program synthesis language that searches the type-safe
function candidates automatically; 2) A novel causal differentiable learning
scheme that determines the causal variables. In the following, we first
present the experimental results with emphasis on the individual components 1)
and 2). We then demonstrate the overall performance on the main and ablation
studies for the complete AutoCI method, confirming that the integration of
program synthesis with a causal differentiable learning scheme is a critical
step towards automated causal inference for clinical applications.

Automated type-safe function search. Type-safe functions are favorable
candidates to strengthen the security and reliability of critical software
applications [41]. However, the manual verification of type-safe properties
can be cumbersome, especially when examining thousands of function candidates.
Fig. 1 (bottom left) shows that the AutoCI can efficiently filter out a subset
of promising type-safe differentiable functions (See also similar results in
Valkov  _et al_. [42]), and this can be done in a short period of time,
$42.26s,118.56s,119.44s$ for $\text{size}=\\{3,4,5\\}$. By automatically
excluding a large amount of generic functions that are not type-safe, it can
greatly improve the development efficiency and algorithm safety compared to
manual function design. After obtaining the type-safe candidates, we execute
the causal differentiable learning scheme (Fig. 1 (top)) on the candidates and
determine the set of predicted causal variables $S_{\mathsf{pred}}$. Here we
utilize the Jaccard Similarity (JS) [25] as the key metric to measure the
prediction accuracy,

$\mathsf{JS}(S_{\mathsf{pred}},S_{\mathsf{prox}})=\frac{|S_{\mathsf{pred}}\cap
S_{\mathsf{prox}}|}{|S_{\mathsf{pred}}\cup S_{\mathsf{prox}}|}.$ (4)

Fig. 1 (bottom middle) demonstrates the JS accuracy with the growing epochs
for the top-four type-safe candidates (See also Supplementary Tab. 6). We
conclude that the function
$f=\mathsf{COMP}(\mathsf{nn},\mathsf{CAT}(\mathsf{FILTER}(\mathsf{pred})))$
achieves the optimal JS score, _e.g_., $91.9\pm 0.06\%$, and thus is used as
the default function for further analysis. As pointed out in Valkov  _et al_.
[42], HOUDINI allows us to transfer high-level modules across learning tasks.
More specifically, the type-safe candidates discovered via the search
algorithm are agnostic of disease-specific features for hazard analysis.
Independent of the hazard analysis conducted for cancer studies, these type-
safe candidates can therefore be re-used and fine-tuned to perform causal
variable identification given data on the survival outcome for each patient.
Depending on the JS score achieved by the candidates we determine the optimal
learned type-safe model. As a result, the proposed AutoCI approach can pave
the way towards an efficient causal analysis for many of the real-world cancer
studies.

Table 1: The comparison of causal variable determination for the PORTEC
dataset among ICP, NICP and the proposed method. Left column: The causal
variable determination of P, PM and PMI. Here, we report Yes (Causal) or No
(Non-causal) for ICPs and causal probability for the proposed AutoCI. Right
column: The hazard analysis including the results of Hazard ratio (HR), 95%
Confidence interval (CI) and p-value for P, PM and PMI, where the p-value is
computed from the chi-square test. The Proposed (Warm-up) refers to the Step
1-2 of the pseudo code in Fig. 1.

Determining causal variables with clear differentiation. We compare the AutoCI
with the state-of-the-art ICP methods, ICP and NICP. Although competitive
results are achieved by AICP on the toy experiments, AICP requires the re-
generation of additional interventional data in each learning step. Such a
learning scheme is incompatible to the real-world RCTs setting, hence AICP is
not applicable (N/A) for the PORTEC experiments. As ICP and NICP explicitly
accept or reject the variable of interest, we report YES or NO in Tab. 1. For
the proposed AutoCI, we report the mean of the causal probabilities for each
variable. Overall, the proposed AutoCI outperforms the ICP and NICP in terms
of differentiating causal and non-causal variables ($\geq 50\%$ vs $<25\%$ for
PMI), demonstrating its advantages over the SOTA methods with a clear margin
(See Tab. 1 (left column)). When examining the individual variables of
interest, we can see that ICP fails to determine the proxy variables to be the
causal ones, while for NICP all the variables including the sanity check ones
are considered to be causal. These results clearly contrast to the
methodological comparison of the ICP methods on the toy data respecting the
normal distributions. If we decompose the proposed causal learning scheme, we
witness the suppression of the causal probability on the sanity variables over
the warm-up stage (Step 1-2 of the pseudo code in Fig. 1), while the causal
variables do not show deterioration of performance. This clear differentiation
of causal and non-causal variables can aid the definition of meaningful cut-
offs by AutoCI on a given cohort guided by clinical expertise.

Hazard analysis of the individual variables. In parallel to the causal
variable determination, the corresponding hazard ratio (HR) analysis on
endometrial cancer (EC) recurrence is also performed. In the scenario when
unknown spurious (non-causal) variables are included in the hazard analysis,
the causal cut-off can help reducing the noise introduced by non-causal
variables. For instance, Tab. 1 (right column) reports the decreased hazard of
Tissue area (0.90 (0.88-0.92)), Patient ID (0.92 (0.91-0.94)) achieved in the
warm-up stage for PMI case. Without causal analysis, one may falsely conclude
that larger tissue area leads to a slightly lower risk of cancer recurrence.
Besides, the hazard ratio achieved within the warm-up and complete learning
stage remains stable and consistent to the standard clinical interpretation,
_i.e_., the values assigned to the variable indeed correspond to either poor
or favorable outcome correctly. Therefore, this learning scheme can help
delivering reliable outputs that are understandable for clinical experts.

Ablation study with hidden confounders. To elaborate the robustness of AutoCI,
we conducted ablation studies with the influence of confounding. This is
achieved by step-wise inclusion of pathological variables (P), pathological
and molecular variables (PM). As shown in Tab. 1 (top and middle rows), the
proposed AutoCI presents consistent advantages over the existing ICP methods
in terms of learning meaningful causal probabilities for both non-causal and
causal variables. For instance, the Tissue area, Patient ID are determined to
be 35.96%, 20.61% for P and 24.94%, 17.48% for PM, while the causal
probability of all the proxy variables remain close to or above 50%. In the
hazard analysis, the results of the confounding studies P and PM are also
consistent with the results reported for the main study (PMI). For instance,
L1CAM and p53abn are assigned with increased hazards for PM and PMI,
indicating poor prognostic outcome. These results are in agreement with
clinical understanding [29, 33, 27]. The numerical improvements from P to PMI
in the accuracy of probabilistic predictions (See Fig. 2) provide
complementary evidence confirming the robustness of AutoCI. When comparing to
the warm-up stage of AutoCI, we further observe a reduction in variance in the
evaluation of the complete causal-aware AutoCI. Lastly, the running time of
the P, PM and PMI studies grows only mildly for the proposed method, in
contrast to the dramatic increase in time complexity for ICP and NICP, where
the bottleneck of the ICPs lies in the exponential increase $o(2^{|S|})$ in
search space with the growing number of variables $S$ (Eq. 3).

Figure 2: The evaluation metrics of hazard analysis conducted on the PORTEC
dataset. The box plots of Concordance Index (Left), Brier Score (Middle), and
Binomial Log-likelihood (Right) that are derived from n=640 patients included
in the PORTEC dataset, where the box bounds the interquartile range (IQR)
divided by the median, and whiskers extend to $\pm 1.5\times$IQR beyond the
box.

## Discussion

In this study, we proposed a novel automated causal inference algorithm
(AutoCI). Taking two large RCTs as the experimental environments, we re-
interpret the clinical variables of interest from a causal perspective.
Compared to the existing ICP methods, the proposed AutoCI demonstrates
consistent advantages in determining causal variables with the presence of
hidden confounders. Complementary to the standard hazard analysis, AutoCI
provides an automatic tool for medical data analysis to investigate causal
association of patient variables from a new perspective, offering informative
and critical evidence to support clinical interpretation. Specifically, the
accurate determination and exclusion of the spurious (non-causal) variables is
a key step to enable more precise patient stratification in the future.

Design choices of AutoCI. Dissimilar to generic algorithms, clinical
algorithms must deal with unique challenges in terms of ensuring the safety
[43] and robustness. Error-prone algorithms can potentially lead to critical
errors in medical care. Driven by the need to develop safe-critical
applications, AutoCI is carried out with a type-safe program synthesis method
[42]. By further incorporating a newly proposed causal-aware module into this
framework, we are able to synthesize a subset of differentiable type-safe
candidates well suited for causal-aware learning. Compared to the laborious
and error-prone manual function design, this implementation improves the
efficiency and safety of AutoCI. Moreover, to achieve the robustness on the
real clinical tasks, we introduce a novel causal differentiable learning
scheme that utilizes the Fréchet inception distance (FID) [44]. As a whole,
the proposed AutoCI is the seamless integration of both components.

Comparison to existing ICPs. Application of the prior ICP methods has
confirmed the feasibility of causal variable learning on toy experiments (See
Tab. 2). This is substantiated by the outstanding results in the absence of
confounders. However, the error-tolerant implementations of prior ICPs on the
synthesized experiments are not well-tailored for real clinical applications,
especially in the presence of hidden confounders. Compared to the ICP, AICP
and NICP, the AutoCI presents robust results on both toy and PORTEC
experiments. With the inclusion of confounders, AutoCI demonstrated a robust
differentiation between causal and non-causal variables for PORTEC, and
achieves superior quantitative scores on both the finite sample and ABCD
settings.

Clinical interpretations. Importantly, the hazard analysis and ranking of
clinicopathological and molecular variables using AutoCI (See Tab. 1) is
generally consistent with the common biological and clinical interpretation.
Taking pathological variables as an example, studies [26, 27, 28] indeed show
that grade, deep myometrial invasion and LVSI are significant independent
predictors of early EC recurrence. The causal probabilities provided by AutoCI
thereby give additional information on the relevance of each variable for the
determination of outcome, and the likelihood of each variable is consistent
with domain expertise. LVSI is considered to be a significant predictor
independent of molecular subgroup and is ranked with the highest causal
probability, while grade and myometrial invasion are indeed weaker but
independent prognostic indicators.

Further, AutoCI correctly identifies the prognostic associations of the
molecular variables of endometrial cancer, assigns the appropriate hazards for
outcome, and ranks the molecular subgroups in the order of causal probability
that would be expected by expert’s domain knowledge. Specifically, the
molecular factors with the highest adverse risk, p53 abnormality (3.01 (PMI),
3.22 (PM)) and L1CAM over-expression (2.16 (PMI), 2.33 (PM)) [33] are
recognised as such, while the POLEmut is consistently associated with a
reduced risk of disease relapse as confirmed in previous studies [30]. Adding
an immune variable further refines the model, as expected by domain expertise
[3], and highlights a causal relationship between cytotoxic T-cell responses
and EC recurrence in early-stage endometrial cancer. In summary, AutoCI
correctly quantifies and ranks causal pathological, molecular and immune
variables for patient outcomes in the clinical trial setting.

Going beyond academic toy models, the proposed AutoCI extends the researchers
current statistical toolbox with a new causal-driven method that can assign
the causal likelihood to prognostic and predictive variables. As such, this
method will enable identification of clinically relevant variables among the
ever-increasing number of biomarkers in cancer research that show statistical
correlation with clinical outcome. Hence, while the direct real-world
application of this method is primarily scientific, subsequent clinical
validation and development may enable better selection of (bio)markers to
stratify patients for cancer treatments and prediction of prognosis.

Limitation. Despite the clear cut-off between non-causal and proxy variables
provided by AutoCI, some of the proxy variables present borderline hazard
ratios, for instance MMRd. Due to small effect size, clinical studies [45, 46]
usually associate Mismatch Repair Deficiency (MMRd) with intermediate patient
prognosis, not significantly different from the prognosis of NSMP EC. However,
from the biological perspective, MMRd is highly relevant for a well-defined
cascade of molecular changes in cancer cells with favorable prognostic impact
as proven by a large number of well-designed experimental and translational
studies [31, 47]. The loss of DNA mismatch repair capability in cancer cells
leads to a strong increase in tumor mutational burden caused by mismatch,
frameshift and insertion/deletion mutations. Due to the structure of
eucaryotic DNA, frameshift mutations frequently lead to the translation of
truncated peptides that are highly immunogenic, and contribute to the
induction of an effective anti-tumoral immune response [47]. In consistency
with biological understanding, AutoCI identifies MMRd as causally related to
outcome in the present study, though the lack of a statistically significant
association with EC recurrence requires further investigation.

Conclusion. In this study, we investigate the causal variable determination
among multiple RCTs and present its advantages over the compared methods on
both toy and PORTEC experiments. For clinical application, further validation
of this methodology in independent clinical trial datasets will be needed to
ensure generalisation.

Table 2: The comparison of causal variable determination for the toy datasets
between ICPs and the proposed method. Top: The illustration of the Structural
Causal Model (SCM), from which the data are sampled for the cases of finite
sample and ABCD setting. Middle: The results of the compared methods for the
Finite sample setting. Bottom: The results of the compared methods for the
ABCD setting.

## Methods

As embedded in the AutoCI abbreviation, automated and causal are the two
building blocks of the proposed method. Concretely, the automation component
is implemented with a type-safe program synthesis language HOUDINI [42].
Further, we introduce a novel differentiable causal learning scheme that is
built upon ICP.

### Type-safe program synthesis

HOUDINI [42] is a typed language with a rich set of pythonic higher-order
functions such as MAP, FOLD, COMP (Pythonic:
$\mathsf{map()},\mathsf{reduce()},\mathsf{lambda\;x:f(g(x))}$), etc (See Fig.
1 (top left)). Relying on the built-in method for program search, it allows us
to efficiently search promising type-safe differentiable program candidates.
Compared to other program synthesis languages [48, 49, 50, 51], HOUDINI rules
out the error-prone functions that undermine the software safety and presents
itself as an ideal candidate for our task (See Supplementary Tab. 7 for
further comparison).

Despite of rich built-in functions provided by the HOUDINI, it lacks an
explicit flow control mechanism. Driven by the need of integrating the causal-
aware learning, we introduce the predicate module ($\bm{\mathsf{PRED}}$)
containing the function (cau) (See also Fig. 1 (top left)),

$\mathsf{cau}(\bm{x};\bm{\theta})=\mathsf{mask}\odot\mathsf{sigmoid}(\bm{\theta})\odot\bm{x},$
(5)

where $\bm{\theta}$ are the learnable weights (normalized by
$\mathsf{sigmoid}$), $\odot$ is the element-wise multiplication,
$\mathsf{mask}$ is the vector containing $0$ or $1$ manipulated in Step 5 of
Fig. 1 (top), $\mathsf{mask}\odot\mathsf{sigmoid}(\bm{\theta})$ presents the
causal probability for each variable. Together with the newly introduced
higher-order function $\bm{\mathsf{FILTER}}$, we are able to synthesize type-
safe causal-aware programs.

### Causal differentiable learning

In Definition 1, the ICP does not make assumptions on the function $f^{*}$
(Eq. 1). In real-world applications, it is reasonable to specify the search
space of $f^{*}$. If we assume that $f^{*}$ is differentiable, then we have a
trivial extension:

$\displaystyle f:$ $\displaystyle\quad\mathbb{R}^{n}\mapsto\mathbb{R}$ (6)
$\displaystyle\quad\underbrace{X_{i_{0}},\ldots,X_{i_{|S^{*}|}}}_{\bm{X}_{S^{*}}},\underbrace{X_{i_{|S^{*}|+1}},\ldots,X_{i_{n}}}_{\bm{X}_{S^{*c}}}\rightarrow
f^{*}(\bm{X}_{S^{*}}),\quad\forall u\in U,$

where $f$ remains differentiable _w.r.t_. all the variables $\bm{X}$. More
importantly, the gradient norms _w.r.t_. the non (plausible) causal variables
$\bm{X}_{S^{*c}}$ should vanish, _e.g_., $\|\nabla_{S^{*c}}f\|=0$. Motivated
by the extension, we first assume $f^{*}$ to be differentiable in Definition
1, then we have the claim:

######  Claim 1.

Following the specification of environments $U$, random variables
$\bm{X}=(X_{1},X_{2},\ldots,X_{n})$ and an outcome variable $Y$ in Definition
1, if $\bm{X}_{\hat{S}}=(X_{\hat{S}_{1}},\ldots,X_{\hat{S}_{j}})$ with indices
$\hat{S}\subseteq\\{1,\ldots,n\\}$ are the identifiable causal variables, then
there exists a differentiable function
$f(\bm{X}):\mathbb{R}^{n}\mapsto\mathbb{R}$ satisfying (1, 2) such that $f$
has the maximum amount $|\hat{S}^{c}|$ of variables with
$\|\nabla_{\hat{S}^{c}}f\|=0$.

From this perspective, we can reduce the ICPs to learning an invariant
differentiable function $f$, where $f$ has the (maximum amount of) vanishing
gradient norms on the non-causal variables. Such reduction enables us to
smoothly integrate the ICP into the modern differentiable learning framework.

Algorithm design. To impose the vanishing gradient norms, we seek for the
$\mathsf{mask}$ vector in Eq. 5 as the solution. Initially, we assign
$\mathsf{mask}=1$ for all the variables. Assume we mask a causal variable
$X_{i}$ by flipping $\mathsf{mask_{i}}=0$, then it should greatly disturb the
learning errors across multiple environments. If this is the case, we restore
$\mathsf{mask_{i}}=1$ and take $X_{i}$ as the causal variable. Otherwise we
reject the variable $X_{i}$ and $\mathsf{mask_{i}}$ remains to be $0$, which
imposes the zero gradient with respect to $X_{i}$. To quantify the disturbance
when variables of interest are missing, existing ICPs use several statistical
tests [52, 53, 54]. These tests suffer from capturing the nuance of
distributions related to higher-dimensional clinical data. Motivated by the
recent success in complex vision data [55], we utilize the Fréchet inception
distance (FID) [44] that is derived from the Wasserstein distance [56]. More
specifically, the square root of FID between Gaussian distributions is exactly
$W_{2}$ Wasserstein distance [56], _i.e_., satisfying three axioms: identity,
symmetry, triangular inequality, whereas the statistical tests used in [23,
24, 25] are generally not mathematical metrics. As a result, the maximum FID
(mFID) of different environments $U$ is proposed to measure the distribution
difference,

$\mathsf{mFID}=\max_{u\in U}\mathsf{FID}\,(\bm{\mu}_{u},\bm{\mu}_{u^{c}}),$
(7)

where $\bm{\mu}_{u},\bm{\mu}_{u^{c}}$ are the distributions w.r.t the data
sampled from the environment(s) $\\{u\\}$ and $\\{u\\}^{c}$. In the
supplement, Tab. 8 presents side by side comparisons between $\mathsf{mFID}$,
F-test + t-test [23, 25] and Levene-test + Wilcoxon-test [24], all of which
are applied for training the same type-safe function
$\mathsf{COMP}(\mathsf{nn},\mathsf{CAT}(\mathsf{FILTER}(\mathsf{pred})))$
under the proposed causal differentiable learning scheme. As displayed in the
supplementary Tab. 8, we conclude that the proposed $\mathsf{mFID}$
outperforms the compared statistical tests with a clear margin. For the pseudo
code of proposed algorithm please check the top plot of Fig. 1 (See the box of
‘causal differential learning’ for more details).

Figure 3: The learning performance and time of the optimal type-safe
$f=\mathsf{COMP}(\mathsf{nn},\mathsf{CAT}(\mathsf{FILTER}(\mathsf{pred})))$ on
toy datasets. Left: The learning curve of the Jaccard Similarity(JS) for top-
four type-safe functions (Finite sample setting). Middle: The learning curve
of the Jaccard Similarity(JS) for top-four type-safe functions (ABCD setting).
Right: The running time of the compared methods for the Finite sample and ABCD
settings. Here the proposed AutoCI utilizes the function
$\mathsf{COMP}(\mathsf{nn},\mathsf{CAT}(\mathsf{FILTER}(\mathsf{pred})))$.

### Proof of concept

For the sake of concept validation, we first conduct experiments on toy
datasets. We compare the proposed AutoCI to the SOTA methods ICP, NICP and
AICP. Specifically, we follow the two experimental protocols presented in AICP
[25]: Finite sample setting and ABCD setting [57]. The former presents the
ideal scenario where the same amount of data (1000) are sampled from both
observational and experimental (interventional) environments, the latter
simulates a more realistic case where limited experimental data (10) are
collected in conjunction with a large amount of observational data (1000). The
data of both settings are generated from randomly chosen linear structural
causal models (SCM). In our experiments, 400 SCM models are tested to
guarantee the reliability of our results. For the compared ICP methods, we
applied the optimal strategies discussed in the paper and parameters are fine-
tuned to the experiments. Specifically, careful parallelization and code
optimization is also performed for ICP methods. We use 16 cores of CPU
Intel(R) Core(TM) i7-7820X CPU @ 3.60GHz to train the ICP methods in parallel
and the GPU NVIDIA TITAN V (12GB) to train the AutoCI. For the proposed
AutoCI, we use the standard Adam optimization [58] with the learning rate
$0.02$ throughout the experiments. For the warm-up stage (Step 1-2 of Fig. 1)
of causal differentiable learning we adopt eight epochs. The batch size is set
to be $64$ for all the experiments. We calibrate the $\lambda=5,1$ for toy and
PORTEC on a small subset of unused data. For the toy experiment, we apply the
MSE loss to supervise the learning process and report the result obtained by
training the AutoCI one time, where the non-causal variable is determined to
be the one with $0$ causal probability (Eq. 5). For the PORTEC experiment, we
utilize the partial likelihood to learn the hazard coefficient [59]. To fully
utilized the PORTEC patient data and incorporate into the differentiable cox
model [59], the molecular subtype variables POLEmut, MMRd and p53abn are
assigned with 1 if present else (including NSMP) 0. To guarantee the
representativeness of the PORTEC results, we independently train the AutoCI 64
times and average the causal probability for each variable. Complementary to
JS score, we also report the Family-wise error rate
$\mathsf{FWER}=P(S_{\mathsf{pred}}\nsubseteq S_{\mathsf{prox}})$ (type-I
error).

As shown in Tab. 2, our AutoCI achieved competitive JS and FWER scores
compared to the ICP methods. Clearly, the proposed method is more resistant to
the influence of hidden confounders and all the results reach $>90\%$ JS
accuracy. This is achieved by the optimal type-safe function
$\mathsf{COMP}(\mathsf{nn},\mathsf{CAT}(\mathsf{FILTER}(\mathsf{pred})))$
(Fig. 3 (left and middle), Supplementary Tab. 9, 10). Such advantages also
confirm the effectiveness of the proposed causal learning scheme with the
utilization of mFID metric. Similar to the PORTEC experiments, because of the
exhaustive subset research required in Eq. 3, the time complexity of ICP and
NICP raises dramatically from ABCD to Finite sample settings (Fig. 3 (right)).

## Data availability

The code used to generate the data of the toy experiments is available via the
link https://github.com/juangamella/aicp. TransPORTEC complies with the FAIR
guiding principles for scientific data management and stewardship. TransPORTEC
is oriented towards a knowledge and data sharing policy. In other words,
expertise and structure of different groups will be made available for all the
proposals and the projects presented by the member groups. In the knowledge-
sharing policy of the consortium, the possibility and the power offered by the
collected databank of samples and data can be used by every participating
member. All member researchers/groups are welcome to introduce new proposals.
In principle, all the anonymized pathological and molecular data already
included in the TransPORTEC biobank will be available for every member of the
consortium upon request. In the future, after completion of the primary joint
projects, the biobank data will be opened towards use by researchers worldwide
upon request.

## Code availability

Our code is implemented with PyTorch and publicly accessible via
https://github.com/CTPLab/AutoCI, which is released under the MIT licence.

## Acknowledgments

The authors would like to convey the gratitude to all clinicians and
technicians participated in the PORTEC 1 and 2 trials (registration number
ISRCTN16228756), and all scientists, pathologists and patients involved in the
data processing and analysis. The PORTEC 1 and PORTEC 2 trials were supported
by the grants from the Dutch Cancer Society (CKTO 90–01 and CKTO 2001–04,
respectively). The molecular profiling was supported by the grants from the
Dutch Cancer Society (KWF UL2012-5447 and KWF/YIG 10418, respectively). Prof.
Koelzer reports grants from Promedica Foundation (F-87701-41-01) during the
conduct of the study. Dr. Horeweg reports grants from the Dutch Cancer Society
(KWF-2021-13400, KWF-2021-13404) during the conduct of the study.

## Author Contributions

J.W. and V.H.K. conceived the research idea. J.W. implemented the algorithm
and ran the experiments. J.W., V.H.K., and N.H. analysed the data and results.
J.W. and V.H.K. wrote the manuscript. N.H. and V.H.K. reviewed the manuscript
and supervised this study. R.A.N., I.M.J.-S., J.J.J., L.C.H.W.L.,
E.M.v.d.S.-B., T.B., C.L.C., M.d.B., H.W.N., T.B., and V.T.H.B.M.S. provided
the data (PORTEC1 and PORTEC2 trials). C.L.C. and V.T.H.B.M.S. conceptualized
and designed the PORTEC 1 and PORTEC 2 study. R.A.N., T.B., C.L.C., and
V.T.H.B.M.S. supervised the PORTEC 1 and PORTEC 2 study.

## Competing interests

The authors declare no competing interests.

## Ethical approval

The PORTEC study protocols were approved by the Dutch Cancer Society and by
the medical ethics committees at participating centers. Both studies were
conducted in accordance with the principles of the Declaration of Helsinki.
All patients provided informed consent for study participation. The PORTEC 1
trial was registered at the Daniel Den Hoed Cancer Center (DDHCC) Trial
Office. The PORTEC 2 trial was registered in ClinicalTrials.gov under the
identifier NCT00376844.

## References

  * [1] Pearl, J.  Causal inference in the health sciences: a conceptual introduction.  _Health services and outcomes research methodology_ 2, 189–220 (2001).
  * [2] Voysey, M. _et al._ Safety and efficacy of the chadox1 ncov-19 vaccine (azd1222) against sars-cov-2: an interim analysis of four randomised controlled trials in brazil, south africa, and the uk.  _The Lancet_ 397, 99–111 (2021).
  * [3] Horeweg, N. _et al._ Prognostic integrated image-based immune and molecular profiling in early-stage endometrial cancer.  _Cancer Immunology Research_ 8, 1508–1519 (2020).
  * [4] Hariton, E. & Locascio, J. J.  Randomised controlled trials—the gold standard for effectiveness research.  _BJOG: an international journal of obstetrics and gynaecology_ 125, 1716 (2018).
  * [5] Peters, J., Janzing, D. & Schölkopf, B.  _Elements of causal inference_ (The MIT Press, 2017).
  * [6] Creutzberg, C. L. _et al._ Surgery and postoperative radiotherapy versus surgery alone for patients with stage-1 endometrial carcinoma: multicentre randomised trial.  _The Lancet_ 355, 1404–1411 (2000).
  * [7] Creutzberg, C. L. _et al._ Fifteen-year radiotherapy outcomes of the randomized portec-1 trial for endometrial carcinoma.  _International Journal of Radiation Oncology* Biology* Physics_ 81, e631–e638 (2011).
  * [8] Nout, R. A. _et al._ Vaginal brachytherapy versus pelvic external beam radiotherapy for patients with endometrial cancer of high-intermediate risk (portec-2): an open-label, non-inferiority, randomised trial.  _The Lancet_ 375, 816–823 (2010).
  * [9] Wortman, B. _et al._ Ten-year results of the portec-2 trial for high-intermediate risk endometrial carcinoma: improving patient selection for adjuvant therapy.  _British journal of cancer_ 119, 1067–1074 (2018).
  * [10] Sung, H. _et al._ Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries.  _CA: a cancer journal for clinicians_ 71, 209–249 (2021).
  * [11] van den Heerik, A. S. V., Horeweg, N., de Boer, S. M., Bosse, T. & Creutzberg, C. L.  Adjuvant therapy for endometrial cancer in the era of molecular classification: radiotherapy, chemoradiation and novel targets for therapy.  _International Journal of Gynecologic Cancer_ 31 (2021).
  * [12] Concin, N. _et al._ Esgo/estro/esp guidelines for the management of patients with endometrial carcinoma.  _International Journal of Gynecologic Cancer_ 31 (2021).
  * [13] Hernán, M. A. & Robins, J. M.  Causal inference: what if (2020).
  * [14] Zenil, H., Kiani, N. A., Zea, A. A. & Tegnér, J.  Causal deconvolution by algorithmic generative models.  _Nature Machine Intelligence_ 1, 58–66 (2019).
  * [15] Prosperi, M. _et al._ Causal inference and counterfactual prediction in machine learning for actionable healthcare.  _Nature Machine Intelligence_ 2, 369–375 (2020).
  * [16] Luo, Y., Peng, J. & Ma, J.  When causal inference meets deep learning.  _Nature Machine Intelligence_ 2, 426–427 (2020).
  * [17] Pearl, J. _et al._ Causal inference in statistics: An overview.  _Statistics surveys_ 3, 96–146 (2009).
  * [18] Hernán, M. A., Sauer, B. C., Hernández-Díaz, S., Platt, R. & Shrier, I.  Specifying a target trial prevents immortal time bias and other self-inflicted injuries in observational analyses.  _Journal of clinical epidemiology_ 79, 70–75 (2016).
  * [19] Dickerman, B. A., García-Albéniz, X., Logan, R. W., Denaxas, S. & Hernán, M. A.  Avoidable flaws in observational analyses: an application to statins and cancer.  _Nature medicine_ 25, 1601–1606 (2019).
  * [20] Caniglia, E. C. _et al._ Emulating a target trial of antiretroviral therapy regimens started before conception and risk of adverse birth outcomes.  _AIDS (London, England)_ 32, 113 (2018).
  * [21] Dahabreh, I. J., Robertson, S. E., Steingrimsson, J. A., Stuart, E. A. & Hernan, M. A.  Extending inferences from a randomized trial to a new target population.  _Statistics in medicine_ 39, 1999–2014 (2020).
  * [22] Zhu, S., Ng, I. & Chen, Z.  Causal discovery with reinforcement learning.  _International Conference on Learning Representations_ (2019).
  * [23] Peters, J., Bühlmann, P. & Meinshausen, N.  Causal inference by using invariant prediction: identification and confidence intervals.  _Journal of the Royal Statistical Society. Series B (Statistical Methodology)_ 947–1012 (2016).
  * [24] Heinze-Deml, C., Peters, J. & Meinshausen, N.  Invariant causal prediction for nonlinear models.  _Journal of Causal Inference_ 6 (2018).
  * [25] Gamella, J. L. & Heinze-Deml, C.  Active invariant causal prediction: Experiment selection through stability.  _arXiv preprint arXiv:2006.05690_ (2020).
  * [26] Scholten, A. N. _et al._ Postoperative radiotherapy for stage 1 endometrial carcinoma: long-term outcome of the randomized portec trial with central pathology review.  _International Journal of Radiation Oncology* Biology* Physics_ 63, 834–838 (2005).
  * [27] Stelloo, E. _et al._ Improved risk assessment by integrating molecular and clinicopathological factors in early-stage endometrial cancer—combined analysis of the portec cohorts.  _Clinical cancer research_ 22, 4215–4224 (2016).
  * [28] Bosse, T. _et al._ Substantial lymph-vascular space invasion (lvsi) is a significant risk factor for recurrence in endometrial cancer–a pooled analysis of portec 1 and 2 trials.  _European journal of cancer_ 51, 1742–1750 (2015).
  * [29] Kandoth, C. _et al._ Mutational landscape and significance across 12 major cancer types.  _Nature_ 502, 333–339 (2013).
  * [30] Church, D. N. _et al._ Prognostic significance of pole proofreading mutations in endometrial cancer.  _JNCI: Journal of the National Cancer Institute_ 107 (2015).
  * [31] Stelloo, E. _et al._ Refining prognosis and identifying targetable pathways for high-risk endometrial cancer; a transportec initiative.  _Modern pathology_ 28, 836–844 (2015).
  * [32] Vermij, L., Smit, V., Nout, R. & Bosse, T.  Incorporation of molecular characteristics into endometrial cancer management.  _Histopathology_ 76, 52–63 (2020).
  * [33] Bosse, T. _et al._ L1 cell adhesion molecule is a strong predictor for distant recurrence and overall survival in early stage endometrial cancer: pooled portec trial results.  _European journal of cancer_ 50, 2602–2610 (2014).
  * [34] Van Gool, I. C. _et al._ Prognostic significance of l1cam expression and its association with mutant p53 expression in high-risk endometrial cancer.  _Modern pathology_ 29, 174–181 (2016).
  * [35] Koelzer, V. H., Sirinukunwattana, K., Rittscher, J. & Mertz, K. D.  Precision immunoprofiling by image analysis and artificial intelligence.  _Virchows Archiv_ 474, 511–522 (2019).
  * [36] Zlobec, I., Koelzer, V. H., Dawson, H., Perren, A. & Lugli, A.  Next-generation tissue microarray (ngtma) increases the quality of biomarker studies: an example using cd3, cd8, and cd45ro in the tumor microenvironment of six different solid tumor types.  _Journal of translational medicine_ 11, 1–7 (2013).
  * [37] Creutzberg, C. L. _et al._ Nomograms for prediction of outcome with or without adjuvant radiation therapy for patients with endometrial cancer: a pooled analysis of portec-1 and portec-2 trials.  _International Journal of Radiation Oncology* Biology* Physics_ 91, 530–539 (2015).
  * [38] Karnezis, A. N. _et al._ Evaluation of endometrial carcinoma prognostic immunohistochemistry markers in the context of molecular classification.  _The Journal of Pathology: Clinical Research_ 3, 279–293 (2017).
  * [39] Talhouk, A. _et al._ Molecular subtype not immune response drives outcomes in endometrial carcinoma.  _Clinical Cancer Research_ 25, 2537–2548 (2019).
  * [40] Lipsitch, M., Tchetgen, E. T. & Cohen, T.  Negative controls: a tool for detecting confounding and bias in observational studies.  _Epidemiology (Cambridge, Mass.)_ 21, 383 (2010).
  * [41] Yang, J. & Hawblitzel, C.  Safe to the last instruction: automated verification of a type-safe operating system.  _Proceedings of the 31st ACM SIGPLAN Conference on Programming Language Design and Implementation_ 99–110 (2010).
  * [42] Valkov, L., Chaudhari, D., Srivastava, A., Sutton, C. & Chaudhuri, S.  Houdini: Lifelong learning as program synthesis.  _Advances in Neural Information Processing Systems_ 8687–8698 (2018).
  * [43] Allen, B.  The role of the fda in ensuring the safety and efficacy of artificial intelligence software and devices.  _Journal of the American College of Radiology_ 16, 208–210 (2019).
  * [44] Heusel, M., Ramsauer, H., Unterthiner, T., Nessler, B. & Hochreiter, S.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  _Advances in neural information processing systems_ 6626–6637 (2017).
  * [45] Smyth, E. C. _et al._ Mismatch repair deficiency, microsatellite instability, and survival: an exploratory analysis of the medical research council adjuvant gastric infusional chemotherapy (magic) trial.  _JAMA oncology_ 3, 1197–1203 (2017).
  * [46] León-Castillo, A. _et al._ Molecular classification of the portec-3 trial for high-risk endometrial cancer: impact on prognosis and benefit from adjuvant therapy.  _Journal of clinical oncology: official journal of the American Society of Clinical Oncology_ 38, 3388–3397 (2020).
  * [47] Kloor, M. & von Knebel Doeberitz, M.  The immune biology of microsatellite unstable cancer.  _Hereditary Colorectal Cancer_ 367–384 (2018).
  * [48] Gaunt, A. L., Brockschmidt, M., Kushman, N. & Tarlow, D.  Differentiable programs with neural libraries.  _International Conference on Machine Learning_ 1213–1222 (2017).
  * [49] Mao, J., Gan, C., Kohli, P., Tenenbaum, J. B. & Wu, J.  The neuro-symbolic concept learner: Interpreting scenes, words, and sentences from natural supervision.  _International Conference on Learning Representations_ (2018).
  * [50] Vedantam, R. _et al._ Probabilistic neural symbolic models for interpretable visual question answering.  _International Conference on Machine Learning_ 6428–6437 (2019).
  * [51] Ellis, K. _et al._ Dreamcoder: Growing generalizable, interpretable knowledge with wake-sleep bayesian program learning.  _arXiv preprint arXiv:2006.08381_ (2020).
  * [52] Pfanzagl, J. & Sheynin, O.  Studies in the history of probability and statistics xliv a forerunner of the t-distribution.  _Biometrika_ 891–898 (1996).
  * [53] Levene, H.  Robust tests for equality of variances.  _Contributions to probability and statistics. Essays in honor of Harold Hotelling_ 279–292 (1961).
  * [54] Wilcoxon, F.  Individual comparisons by ranking methods.  _Breakthroughs in statistics_ 196–202 (1992).
  * [55] Lucic, M., Kurach, K., Michalski, M., Gelly, S. & Bousquet, O.  Are GANs created equal? a large-scale study.  _arXiv preprint arXiv:1711.10337_ (2017).
  * [56] Villani, C.  _Optimal transport: old and new_ , vol. 338 (Springer Science & Business Media, 2008).
  * [57] Agrawal, R., Squires, C., Yang, K., Shanmugam, K. & Uhler, C.  Abcd-strategy: Budgeted experimental design for targeted causal structure discovery.  _The 22nd International Conference on Artificial Intelligence and Statistics_ 3400–3409 (2019).
  * [58] Kingma, D. & Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * [59] Kvamme, H., Borgan, Ø. & Scheel, I.  Time-to-event prediction with neural networks and cox regression.  _arXiv preprint arXiv:1907.00825_ (2019).

Supplementary material

Abbreviation | Definition
---|---
RCT | Randomised controlled trials
PORTEC | Post operative radiation therapy in endometrial carcinoma
EC | Endometrial carcinoma (cancer)
ESGO | European Society of Gynaecological Oncology
ESTRO | European Society for Radiotherapy and Oncology
ESMO | European Society of Medical Oncology
Grade | Tumor grading
LVSI | Lymphovascular space invasion
POLEmut | Polymerase epsilon mutant EC
MMRd | Mismatch repair deficient EC
p53abn | p53 abnormal EC
NSMP | EC with no specific molecular profile
L1CAM | L1 cell adhesion molecule
P | Pathological variables
PM | Pathological and molecular variables
PMI | Pathological, molecular and immune variables
HR | Hazard ratio
CI | Confidence interval
i.d. | identically distributed
NP | Non-deterministic polynomial-time
SCM | Structural causal model
ICP | Invariant causal prediction
PRED | Predicate module
FID | Fréchet inception distance
JS | Jaccard similarity
FWER | Family-wise error rate
ABCD | Active budgeted causal design strategy
Table 3: The abbreviation table of clinical, statistical and causal definitions. Figure 4: The consort diagram presenting the process of patient selection. Abbreviations: QC - quality control, IHC - immunohistochemistry, EEC- endometrioid endometrial carcinoma, NEEC- non-endometrioid endometrial carcinoma. Table 4: The characteristics comparison of excluded and included patients. * After a posteriori central review 2 cases were classified as stage II, 2 as stage IIIA and 1 as stage IIIB. ** The p-values of RFS and OS are computed from log-rank test. Table 5: The characteristics of study participants for PORTEC 1 and 2. * After a posteriori central review 2 cases were classified as stage II, 2 as stage IIIA and 1 as stage IIIB. Table 6: The JS score and its standard deviation of compared type-safe candidates for PORTEC study (PMI). | Task | Typed | Functional | Code availability
---|---|---|---|---
NTPT[48] | Misc. | ✗ | ✗ | ✗
NS-CL[49] | VQA | ✗ | ✓ | ✓
Prob-NMN[50] | VQA | ✗ | ✓ | ✓
DreamCoder[51] | Misc. | ✓ | ✓ | ✗
HOUDINI[42] | Misc. | ✓ | ✓ | ✓
Table 7: The comparison between existing program synthesis languages. Table 8:
The comparison of statistical measurements for toy experiments. Top: The
results of the compared statistical measurements for the Finite sample
setting. Bottom: The results of the compared statistical measurements for the
ABCD setting. Here, all the measurements are applied for training the same
type-safe function
$\mathsf{COMP}(\mathsf{nn},\mathsf{CAT}(\mathsf{FILTER}(\mathsf{pred})))$
under the proposed causal differentiable learning scheme. F-test + t-test is
used in ICP and AICP. Levene-test + Wilcoxon-test is used in NICP. Table 9:
The JS score and its standard deviation of compared type-safe candidates for
toy experiments (Finite sample setting). Table 10: The JS score and its
standard deviation of compared type-safe candidates for toy experiments (ABCD
setting).


# PnP for Two-Dimensional Pose Estimation

Joshua Wang

###### Abstract

We propose a PnP algorithm for a camera constrained to two-dimensional
movement (applicable, for instance, to many wheeled robotics platforms).
Leveraging this assumption allows performance improvements over 3D PnP
algorithms due to the reduction in search space dimensionality. It also
reduces the incidence of ambiguous pose estimates (as, in most cases, the
spurious solutions fall outside the plane of movement). Our algorithm finds an
approximate solution using geometric criteria and refines its prediction
iteratively. We compare this algorithm to existing 3D PnP algorithms in terms
of accuracy, performance, and robustness to noise.

##  I Introduction

The perspective-$n$-point problem, or PnP, asks us to find the pose of a
perspective camera given $n$ 3D points $\textbf{p}_{i}$ and their
corresponding projections onto the image plane, $\textbf{q}_{i}$. Thus, wish
to find a rigid transformation matrix R from the camera frame to the world
frame such that:

$z_{i}\begin{bmatrix}\textbf{q}_{i}\\\ 1\end{bmatrix}=\textbf{I}_{3\times
4}\textbf{R}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\ 1\end{bmatrix}$

Here, the $z_{i}$ are the depth factors and we use known camera intrinsics to
remap the $\textbf{q}_{i}$ so that they model the normalized perspective
camera. This is exactly solvable in the case $n=3$, giving rise to the class
of P3P algorithms as in [1]. P4P and P5P algorithms have also been developed
in this vein [2]. However, using all available correspondences is desirable to
reduce the impact of error, in which case the problem is overconstrained and
we must settle for an approximate solution. One approach is to minimize the
reprojection error, the sum of squared distances between each $\textbf{q}_{i}$
and the projection of its $\textbf{p}_{i}$ onto the predicted image plane.
This is given by:

$\text{Error}(\textbf{R})=\sum_{i}\left|\begin{bmatrix}\textbf{q}_{i}\\\
1\end{bmatrix}-\frac{1}{z_{i}}\textbf{I}_{3\times
4}\textbf{R}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\ 1\end{bmatrix}\right|^{2}$

Note that this does not necessarily measure the distance in pixels due to our
remapping of the $\textbf{q}_{i}$. Minimization of the reprojection error is a
simple, geometric criterion commonly used in iterative PnP solutions; there
are many alternative errors, such as the object-space error of Lu et al. [3].
Non-iterative solutions also exist, such as the EPnP algorithm presented by
Lepetit et al., which utilizes four weighted control points to reduce time
complexity to $O(n)$ [4], or the OPnP algorithm of Zheng et al. which
analytically minimizes an algebraic error by solving a system of polynomial
equations [5]. At the time of publication, state-of-the-art PnP solutions were
mostly non-iterative polynomial system solvers, but iterative solutions such
as Lu et al.’s LHM still remain competitive, and many non-iterative solvers,
including [4, 5], have an iterative phase at the end to refine the predicted
pose (for this purpose, a simple minimization algorithm like Gauss-Newton
suffices).

PnP algorithms have excellent accuracy for large numbers of correspondences
spread widely in 3D space, but many behave suboptimally when the points are
coplanar. Some algorithms like EPnP include a special case for planar
configurations, and assuming coplanarity allows the use of homography-based
solvers such as the IPPE solver of Collins and Bartoli [6]. However, planar
configurations present a challege even for specialized solvers because they
can produce ambiguity; the images seen from multiple camera positions can be
very similar (or even indistinguishable given the presence of noise). The
problem worsens at low resolutions and long distances, or when the plane
containing the points is almost parallel to the image plane. We illustrate the
possible ambiguity with a square marker, the corners of which were used for
pose estimation:

Figure 1: Ambiguity in a square’s estimated pose

This problem is well-known [7], and many methods have been proposed to
mitigate it. Iterative algorithms can be initialized with the last known pose
to improve stability (although this means that they can become persistently
stuck in the wrong pose). Some algorithms, such as IPPE, attempt to compute
all possible solutions, either selecting the best fit or returning the
solutions to the user. The outputs of existing algorithms can be filtered as
in [8]. However, as robust as these solutions are, they solve only a facet of
the deeper issue of imprecision for planar configurations.

This paper deals with the case of a camera constrained to two-dimensional
movement (for instance, mounted to a wheeled robot). In this situation, using
PnP algorithms to estimate 3D pose loses valuable information, as the returned
3D poses must be projected onto the plane of movement. An algorithm which
takes this information into account in the first place would likely exhibit
improved accuracy because it can account for error while remaining in the
plane of motion instead of ’escaping’ into disallowed poses; by reducing the
dimensionality of the search space from 6 to 3, performance improvements could
also be realized. We derive a simple example of such an algorithm, called
2DPnP, which finds an approximate solution and then iteratively minimizes the
reprojection error.

##  II 2DPnP

Assume that the camera origin is on the world $xy$-plane and that the camera
is constrained to move on this plane. Then,
$\textbf{R}=\textbf{P}\begin{bmatrix}\textbf{C}&\textbf{0}\\\
\textbf{0}^{\intercal}&1\end{bmatrix}$ where P is an unknown transformation
preserving the $xy$-plane and C is a known rotation matrix. Let
$\textbf{C}=\textbf{R}_{\alpha}\textbf{R}_{\beta}\textbf{R}_{\gamma}$ be C’s
decomposition into the classical Euler angles $\alpha,\beta,\gamma$. Then,
note that $\textbf{R}_{\alpha}$ can be “absorbed” into P, so we can assume
that $\alpha=0$. We have:

$\frac{1}{z_{i}}\textbf{I}_{3\times
4}\textbf{R}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\
1\end{bmatrix}=\frac{1}{z_{i}}\textbf{I}_{3\times
4}\begin{bmatrix}\textbf{R}_{\gamma}^{-1}\textbf{R}_{\beta}^{-1}&\textbf{0}\\\
\textbf{0}^{\intercal}&1\end{bmatrix}\textbf{P}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\
1\end{bmatrix}$ $=\frac{1}{z_{i}}\textbf{R}_{\gamma}^{-1}\textbf{I}_{3\times
4}\begin{bmatrix}\textbf{R}_{\beta}^{-1}&\textbf{0}\\\
\textbf{0}^{\intercal}&1\end{bmatrix}\textbf{P}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\
1\end{bmatrix}$

Note that the rotation $\textbf{R}_{\gamma}$ does not change $z_{i}$ because
it rotates about the camera’s axis, so:

$\left|\textbf{R}_{\gamma}\begin{bmatrix}\textbf{q}_{i}\\\
1\end{bmatrix}-\frac{1}{z_{i}}\textbf{I}_{3\times
4}\begin{bmatrix}\textbf{R}_{\beta}^{-1}&\textbf{0}\\\
\textbf{0}^{\intercal}&1\end{bmatrix}\textbf{P}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\
1\end{bmatrix}\right|$ $=\left|\begin{bmatrix}\textbf{q}_{i}\\\
1\end{bmatrix}-\frac{1}{z_{i}}\textbf{I}_{3\times
4}\textbf{R}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\ 1\end{bmatrix}\right|$

Thus, letting $\begin{bmatrix}\textbf{q}_{i}^{\prime}\\\
1\end{bmatrix}=\textbf{R}_{\gamma}\begin{bmatrix}\textbf{q}_{i}\\\
1\end{bmatrix}$, we can minimize the reprojection error using the
$\textbf{q}_{i}^{\prime}$ while assuming that $\gamma=0$. Assume that the
camera has 2D pose $(x,y,\theta)$. Then, letting $\Delta x_{i}=x-p_{i,x}$ and
similarly for $\Delta y_{i}$, we have:

$\textbf{P}=\begin{bmatrix}\cos(\theta)&-\sin(\theta)&0&x\\\
\sin(\theta)&\cos(\theta)&0&y\\\ 0&0&1&0\\\
0&0&0&1\end{bmatrix}\Rightarrow\begin{bmatrix}\textbf{R}_{\beta}^{-1}&\textbf{0}\\\
\textbf{0}^{\intercal}&1\end{bmatrix}\textbf{P}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\
1\end{bmatrix}$ $=\begin{bmatrix}-(\Delta x_{i}\cos(\theta)+\Delta
y_{i}\sin(\theta))\cos(\beta)-p_{i,z}\sin(\beta)\\\ \Delta
x_{i}\sin(\theta)-\Delta y_{i}\cos(\theta)\\\ -(\Delta
x_{i}\cos(\theta)+\Delta y_{i}\sin(\theta))\sin(\beta)+p_{i,z}\cos(\beta)\\\
1\end{bmatrix}$
$\Rightarrow\text{Error}(\textbf{R})=\sum_{i}\left|\begin{bmatrix}\textbf{q}_{i}^{\prime}\\\
1\end{bmatrix}-\frac{1}{z_{i}}\textbf{I}_{3\times
4}\begin{bmatrix}\textbf{R}_{\beta}^{-1}&\textbf{0}\\\
\textbf{0}^{\intercal}&1\end{bmatrix}\textbf{P}^{-1}\begin{bmatrix}\textbf{p}_{i}\\\
1\end{bmatrix}\right|^{2}$
$=\sum_{i}\left|\begin{bmatrix}q_{i,x}^{\prime}-\frac{(\Delta
x_{i}\cos(\theta)+\Delta
y_{i}\sin(\theta))\cos(\beta)+p_{i,z}\sin(\beta)}{(\Delta
x_{i}\cos(\theta)+\Delta y_{i}\sin(\theta))\sin(\beta)-p_{i,z}\cos(\beta)}\\\
q_{i,y}^{\prime}-\frac{\Delta y_{i}\cos(\theta)-\Delta
x_{i}\sin(\theta)}{(\Delta x_{i}\cos(\theta)+\Delta
y_{i}\sin(\theta))\sin(\beta)-p_{i,z}\cos(\beta)}\end{bmatrix}\right|^{2}$

Because minimizing $\text{Error}(\textbf{R})$ is a nonlinear least-squares
problem, we can employ iterative algorithms such as Levenberg-Marquardt (used
by our implementation) or Gauss-Newton. These algorithms require the Jacobian
matrix J, which is computed using the above expression for
$\text{Error}(\textbf{R})$.

To initialize the iterative phase of the algorithm, we choose three of the
$\textbf{p}_{i}$, say $\textbf{p}_{1},\textbf{p}_{2},\textbf{p}_{3}$. We find
poses $(x,y,\theta)$ (with corresponding $2\times 2$ rotation matrix
$\textbf{R}_{\theta}$) such that, for $i=1,2,3$, some real number $f_{i}$
makes the following true:

$\begin{bmatrix}\textbf{R}_{\theta}\\\
\textbf{0}^{\intercal}\end{bmatrix}\textbf{R}_{\beta}\begin{bmatrix}\textbf{q}_{i}^{\prime}\\\
1\end{bmatrix}=f_{i}\begin{bmatrix}x-p_{i,x}\\\ y-p_{i,y}\end{bmatrix}$

This states that the projections of $\begin{bmatrix}\textbf{q}_{i}^{\prime}\\\
1\end{bmatrix}$ and $\textbf{p}_{i}^{\prime}$ onto the world $xy$-plane should
be collinear with the camera center, as illustrated in the figure:

$\textbf{q}_{3}^{\prime}$$\textbf{p}_{3}$$\textbf{q}_{2}^{\prime}$$\textbf{p}_{2}$$\textbf{q}_{1}^{\prime}$$\textbf{p}_{1}$$(x,y)$
Figure 2: Initialization scheme

Letting $\begin{bmatrix}r_{i}\cos(\theta_{i})\\\
r_{i}\sin(\theta_{i})\end{bmatrix}=\textbf{I}_{2\times
3}\textbf{R}_{\beta}\begin{bmatrix}\textbf{q}_{i}\\\ 1\end{bmatrix}$, the
collinearity condition is equivalent to:

$(x-p_{i,x})\sin(\theta+\theta_{i})-(y-p_{i,y})\cos(\theta+\theta_{i})=0$

We solve for $x,y$ in terms of $\theta$ using two of the equations, which are
linear in $x,y$, and substitute into the third. This gives two solutions for
$\theta$, satisfying (with indices taken modulo 3):

$\tan(\theta)=\frac{\displaystyle\sum_{i=1}^{3}(p_{i,y}\cos(\theta_{i})-p_{i,x}\sin(\theta_{i}))\sin(\theta_{i+2}-\theta_{i+1})}{\displaystyle\sum_{i=1}^{3}(p_{i,x}\cos(\theta_{i})+p_{i,y}\sin(\theta_{i}))\sin(\theta_{i+2}-\theta_{i+1})}$

We take the solution which results in the lowest reprojection error (using
those three points). In terms of choosing points, our implementation takes the
$\textbf{q}_{i}^{\prime}$ with minimal and maximal $x$-coordinates (in camera
space), and then chooses the third index $j$ such that the distance from
$(p_{j,x},p_{j,y})$ to the line formed by the projections of the two
$\textbf{p}_{i}$ is maximal. This choice prioritizes the robustness of the
estimate for $\theta$.

We chose this initialization strategy because of its $O(n)$ time complexity
and simple closed-form solution. Another plausible strategy is to find the
pose minimizing the following closely related error (over all $n$ points):

$\sum_{i}((x-p_{i,x})\sin(\theta+\theta_{i})-(y-p_{i,y})\cos(\theta+\theta_{i}))^{2}$

This is the sum of the squared distances from the projection of each
$\textbf{p}_{i}$ to the line formed by the camera center and the projection of
the corresponding $\textbf{q}_{i}^{\prime}$. One may easily find the optimal
$x,y$ in terms of $\theta$, but the equation for $\theta$ appears to be
intractable. Other initialization strategies were also of interest, such as
the method of Lu et al., which uses the weak-perspective approximation to
linearize an optimization problem [3].

(a) (a) Average translational error vs. number of points

(b) (b) Average rotational error vs. number of points

(c) (c) Average run time vs. number of points

##  III Testing

Using MATLAB, we tested our algorithm against existing PnP algorithms
calculating 3D pose for general and planar configurations of points. For the
general case, points were randomly generated in the cuboid
$[-2,2]\times[-2,2]\times[4,8]$ in camera space, projected onto the image
plane assuming a focal length of 800 pixels, and each image point was
perturbed with Gaussian noise. This follows the methodology of [9], and the
implementations that we tested came from their MATLAB toolbox. The returned 3D
poses of each algorithm were then projected onto the $xy$-plane. Each data
point represents the average of 250 tests.

We tested the translational and rotational errors of each algorithm for
numbers of points from 10 to 200. The ground truth pose was $(0,0,0)$ for each
test, and the camera rotation offset C was randomized uniformly within the
range of pitch angles $\beta\in[10^{\circ},170^{\circ}]$ to avoid inherent
singularities in the conversion of 3D rotation to 2D heading. Gaussian noise
with a standard deviation of $1$ pixel was applied to each image point. The
same conditions were used to time each algorithm for numbers of points from 50
to 1000.

We also tested the impact of noise on each algorithm. For these tests, we used
50 points with randomized camera rotation, and each image point was perturbed
by Gaussian noise with a standard deviation from 1 to 10 pixels.

The tested algorithms were 2DPnP (ours), LHM [3], EPnP with Gauss-Newton
refinement [4], OPnP [5], and MLPnP without covariance information [9]. This
is by no means a comprehensive list; because none of the 3D algorithms can
accomodate for constraints to 2D motion, they would not be expected to perform
optimally, and our intent was only to show how our method compares to a
variety of approaches to the 3D problem. All tests were done on a 2 GHz Intel
i5 processor. The MATLAB code and data files are available at
https://github.com/25wangj/2DPnPToolbox.

(a) (d) Average translational error vs. amount of noise

(b) (e) Average rotational error vs. amount of noise

##  IV Conclusion

In this paper, we presented a PnP algorithm for a camera constrained to two-
dimensional motion. We found an initial pose using a geometric criterion
amounting to the collinearity of the planar projections of the object and
image points with the camera center, and used Levenberg-Marquardt optimization
to minimize the reprojection error. In testing against currently used PnP
algorithms which find three-dimensional pose, our algorithm demonstrated
significantly better performance due to the reduced dimensionality of the
search space, increased translational and comparable rotational accuracy, and
increased robustness to noise. This algorithm shows promise for applications
involving 2D motion such as warehouse and educational robotics, especially
where processing power or resolution is limited.

## References

  * [1] Laurent Kneip, Davide Scaramuzza and Roland Siegwart  “A novel parametrization of the perspective-three-point problem for a direct computation of absolute camera position and orientation”  In _CVPR 2011_ , 2011, pp. 2969–2976  DOI: 10.1109/CVPR.2011.5995464
  * [2] Martin Bujnak, Zuzana Kukelova and Tomas Pajdla  “A general solution to the P4P problem for camera with unknown focal length”  In _2008 IEEE Conference on Computer Vision and Pattern Recognition_ , 2008, pp. 1–8  DOI: 10.1109/CVPR.2008.4587793
  * [3] C.-P. Lu, G.D. Hager and E. Mjolsness  “Fast and globally convergent pose estimation from video images”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 22.6, 2000, pp. 610–622  DOI: 10.1109/34.862199
  * [4] Vincent Lepetit, Francesc Moreno-Noguer and Pascal Fua  “EPnP: An accurate O(n) solution to the PnP problem”  In _International Journal of Computer Vision_ 81, 2009  DOI: 10.1007/s11263-008-0152-6
  * [5] Yinqiang Zheng et al.  “Revisiting the PnP Problem: A Fast, General and Optimal Solution”  In _2013 IEEE International Conference on Computer Vision_ , 2013, pp. 2344–2351  DOI: 10.1109/ICCV.2013.291
  * [6] Toby Collins and Adrien Bartoli  “Infinitesimal Plane-Based Pose Estimation”  In _International Journal of Computer Vision_ 109, 2014, pp. 252–286  URL: https://api.semanticscholar.org/CorpusID:8674118
  * [7] Gerald Schweighofer and Axel Pinz  “Robust Pose Estimation from a Planar Target”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 28, 2006, pp. 2024–2030  URL: https://api.semanticscholar.org/CorpusID:2630612
  * [8] Mohammad Amin Mehralian and Mohsen Soryani  “EKFPnP: extended Kalman filter for camera pose estimation in a sequence of images”  In _IET Image Processing_ 14.15, 2020, pp. 3774–3780  DOI: https://doi.org/10.1049/iet-ipr.2020.0606
  * [9] S. Urban, J. Leitloff and S. Hinz  “MLPNP – A REAL-TIME MAXIMUM LIKELIHOOD SOLUTION TO THE PERSPECTIVE-N-POINT PROBLEM”  In _ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences_ III–3 Copernicus GmbH, 2016, pp. 131–138  DOI: 10.5194/isprs-annals-iii-3-131-2016