Study design
This is a retrospective cohort study conducted at Tygerberg Hospital (TBH) during the first two waves of the COVID-19 pandemic between 27 March 2020 and 10 February 2021. The TBH is a 1380-bed hospital that serves as the main teaching hospital for Stellenbosch University, Faculty of Medicine and Health Sciences. TBH was designated as a centre for COVID-19 management with additional critical care services. It provides tertiary services to around 3.5 million people.
Study population and sample size
The main study had 490 participants and our study included data from 168 adult patients admitted with severe COVID-19 pneumonia. These participants were on either hydrocortisone or dexamethasone. The diagnosis was confirmed with a positive SARS-CoV-2 polymerase chain reaction (PCR). Details regarding admission criteria to ICU are documented in the Western Cape Government’s provincial guidelines[16].
Data collection
Clinical data were extracted from ICU clinical notes and entered a REDCap® (Research Electronic Data Capture, Stellenbosch, South Africa) database, a secure web application. Laboratory data were imported from the National Health Laboratory Service (NHLS) Laboratory Information System (TrakCare® Lab Enterprise) onto the REDCap database. Data quality assurance was done by the research assistants and verified by the supervisor of the research team before analysis. Detailed information about the clinical parameters is defined in the previously published articles [17, 18].
Targeted maximum likelihood estimation.
TMLE incorporates a two-stage approach. The first step involves estimating an outcome regression function then finally, a targeting step which updates the initial estimate to return an unbiased and efficient estimator of the target parameter [19].
Consider a binary treatment \(A\), outcome \(Y\) (either binary or continuous) and measured covariates \(W\) from an observational dataset; \({O}_{obs}=\left(W,A,Y\right)\).
Causal inference from observational studies using TMLE is dependent on three assumptions. Consistency ( \({Y}_{i}^{a}={Y}_{i }\forall i, with A=a\)), Conditional Exchangeability (\({Y}^{a}⫫A|W\) ) and Positivity (\(Pr \left[A=a|W=w\right]>0 \forall w, with Pr [W=w]\ne 0\)) [2].
Consider ATE as a parameter of interest. It is defined as \({\psi }_{0}=\widehat{{E}_{W}}\left(\widehat{E}\left(Y|A=1,W\right)-\widehat{E}\left(Y|A=0,W\right)\right)\), and interpreted as the difference in mean outcome between treatment groups averaging over the measured covariates[20, 21].
Under TMLE, the initial step estimates the outcome regression function \({\text{E}}_{0}\left[\text{Y}|\text{A},\text{W}\right]\) represented as \({\overline{Q}}_{0}\left(A,W\right)\), whose estimate is represented as \({\overline{Q}}_{n}\left(A,W\right)\)[20].
There are various methods that can be implemented to obtain \({\overline{Q}}_{n}\left(A,W\right)\), the outcome regression function.
Finally, an updated estimate \({\overline{Q}}_{n}^{\text{*}}\left(A,W\right)\) is returned by the targeting step using the initial estimate \({\overline{Q}}_{n}\left(A,W\right)\) and an estimate of the propensity score \({g}_{n}\left(A|W\right)=P\left(A=1|W\right)\) [20].
During the updating step, there too is estimation of the clever covariate \({H}_{n}\left(A,W\right)\) and \({\text{ϵ}}_{\text{n}}\), an estimate of the coefficient \(ϵ\) of the clever covariate\({H}_{n}\left(A,W\right)\) per formula, \(f\left({\overline{Q}}_{n}^{\text{*}}\left(A,W\right)\right)=f\left({\overline{Q}}_{n}\left(A,W\right)\right)+ϵ.{H}_{n}\left(A,W\right)\), with \(f\left(.\right)\) being a suitable link function[20].
\({H}_{n}\left(A,W\right)\) takes the form of \({H}_{n}\left(A,W\right)=\frac{A}{{g}_{n}\left(A\mid W\right)}-\frac{A}{1-{g}_{n}\left(A\mid W\right)}\) [20].
The TMLE estimator is given by the equation[21],
$${\psi }_{n}^{\mathcal{T}MLE}=\frac{1}{n}\sum _{i=1}^{n}{\overline{Q}}_{n}^{\text{*}}\left(1,{W}_{i}\right)-{\overline{Q}}_{n}^{\text{*}}\left(0,{W}_{i}\right)$$
Comprehensive mathematical formulation and detailed step-by-step guide on implementing TMLE can be found in the following studies: Schuler et al, 2017, VanderLaan et al, 2022, Berkeley et al, 2009 and Luque-Fernandez et al, 2018. [3, 20–22].
Super learner (SL) and TMLE
The initial step involved creating a node list where variable roles were defined into W (covariates), A (treatment) and Y (outcome) [20].
The next step involved defining learners which took the form of a list of sl3 learners. Instead of selecting a learner for each likelihood factor to be estimated i.e. Q and g, as illustrated by VanderLaan et al, 2022[20], a stack of suitable base learners taking into consideration the data type of the outcome (continuous) and treatment (binary) was defined based off criteria from Phillips et al,2022[23]. Two meta-learners (Nonnegative linear least squares, Nonlinear Optimization via Augmented Lagrange) were included in the learner library instead of assembling ensemble learners. This was due to lack of computational resources to enable the building of a learner library where these meta-learners would be accommodated as described in literature.
Table 1 shows the defined learners and their respective modelling specifications. Discrete super learners with specified loss functions, i.e., mean squared error for the continuous outcome and binomial log likelihood loss function for the binary treatment, were defined as meta learners to select the best performing algorithm in fitting the outcome regression function \({\overline{Q}}_{n}\left(A,W\right)\) and propensity score function \({g}_{n}\left(A|W\right)\), respectively.
Table 1
Prediction algorithms in their respective code formats R Software.
Algorithm as per R code specification. | Model |
lrnr_mean <- make_learner (Lrnr_mean) | Mean model |
lrnr_rf <- make_learner (Lrnr_ranger) | Faster random forests |
lrnr_glm<-make_learner (Lrnr_glm) | Generalized linear models |
lrn_lasso <- make_learner (Lrnr_glmnet) #(alpha = 1) | Lasso-penalized regressions |
lrn_nnls <- make_learner (Lrnr_nnls) | Nonnegative linear least squares. |
lrn_solnp<-make_learner (Lrnr_solnp, loss_function = loss_loglik_binomial, learner_function = metalearner_logistic_binomial) | Nonlinear Optimization via Augmented Lagrange |
lrn_gam<-make_learner (Lrnr_gam) | Generalised additive models |
lrn_xgb<-make_learner (Lrnr_xgboost) | Extreme gradient boosting |
BART AND TMLE
It involved defining a list of confounders, treatment variable, outcome variable, method (TMLE) to fit the outcome regression function \({\overline{Q}}_{n}\left(A,W\right)\), method (BART) to fit the treatment assignment mechanism, a common support rule (none) to exclude any observations and the estimand ATE.
The method to fit \({\overline{Q}}_{n}\left(A,W\right)\); TMLE fit the outcome regression function using BART and further adjusted with TMLE. For the case of the treatment assignment mechanism, it was fit using BART[24]. Elaborate explanations on the background processes of the package can be found at the webpage1.
Parametric regression and TMLE
Initial step involved transforming the continuous outcome to be bounded between 0 and 1.
Linear regression models were implemented to obtain the outcome regression function and potential outcomes. Thereafter prediction of propensity scores was done by fitting a logistic regression model to the binary treatment.
Clever covariates \({H}_{1}\left(A,W\right)\) and \({H}_{0}\left(A,W\right)\) which corresponded to treatment assignment 1 (dexamethasone) and O (hydrocortisone) were estimated and the fluctuation parameter ϵ estimated.
These were used to update the initial predictions of potential outcomes and a mean difference calculated and rescaled to give the average treatment effect.
The efficient influence curve (EIF) was calculated, and 95% confidence intervals constructed.
The manual implementation of TMLE followed the illustrations from Luque-Fernandez et al, 2018[22, 25]. Figure 1 showed Flow diagram illustrating the implementation of TMLE variations.
Statistical analysis
Descriptive statistics, such as frequency, percentages, and median with interquartile range (IQR), were done to summarize the patient characteristics. The chi-square test, Fisher’s exact test, and Wilcoxon rank sum test compared the patient characteristics between the dexamethasone and hydrocortisone groups. The statistical significance level for these tests was 5%. Super learner, BART, and parametric regression methods were used to obtain the initial outcome regression function estimate \({\overline{Q}}_{n}\left(A,W\right)\). Then, the respective estimates were used to implement TMLE.
SL and TMLE were conducted cohesively under the tmle3 and sl3 libraries in R software. The node list, tmle3_specification object, and defined learner list were passed to the tmle3 function that returned the average treatment effect estimate.
BART and TMLE were done using the R package bartCause. The defined arguments were passed to the bartC function, which returned the average treatment effect estimate. Missing data were imputed using multiple imputation with chained equations through the MICE package. Additional file 1.docx in Table 1 showed the characteristics of original data and imputed data. All analyses were done in R version 4.2.0 (R Core Team, Austria) with R Studio V.1.4 (R Studio Team, Austria) statistical software.