\section{Introduction} \label{sec:intro}

Causal inference is a central goal of health research, and is widely understood as the task of assessing the impact of intervening on a given exposure on an outcome of interest.\cite{Hernan2020} In a perfect randomized controlled trial (i.e., with no loss to follow-up, etc.) causal effects can in principle be directly estimated by comparing the estimated expected outcome in those randomised to the exposure group with the estimated expected outcome in those randomised to be unexposed. However, estimation of causal effects in observational studies is not as straightforward and requires more sophisticated methods, in particular to adjust for potential confounding bias due to the lack of randomisation.

With continuous outcomes, causal effects are typically defined as the contrast of the mean potential outcome under exposure versus the mean under no exposure, which would be the estimand in the hypothetical target trial that we seek to emulate with observational data. However, epidemiological studies may suffer from skewed outcome data, for which the mean may no longer be a meaningful summary statistic. The mean is pulled by more extreme values in the longer tail of the distribution, and therefore is no longer a measure of centrality of the distribution. Examples of skewed outcomes are abundant in health research, particularly in areas where scale scores are used to measure aspects such as self-reported quality of life (via the PedsQL\cite{Varni1999}), self-esteem (via the RSES\cite{Rosenberg1989}), or childhood behaviour (via the SDQ\cite{Goodman1997}). Other examples include time-to-event outcomes in the absence of censoring (e.g., survival time), or duration of events (e.g., breastfeeding duration).  

When faced with this challenge in practice, the typical approach is to either "ignore or transform" - ignore the skewness in the data entirely and use the mean as a summary statistic, or transform the outcome to obtain a more symmetric distribution for which the mean is interpretable as the centre of the distribution. Both approaches have their advantages, although neither is entirely satisfactory in many practical settings. Opting to ignore the skewness in the outcome distribution may be appropriate when the mean remains the statistic of interest, and allows established confounding-adjustment methods to be applied (e.g., linear regression, g-methods, doubly robust methods). However, if the central value of the outcome distribution is of direct interest, then this approach is not suitable as the mean is not an optimal measure of centrality for skewed distributions. Transformation of the outcome to be more symmetrically distributed could be an alternative solution. However, this relies on a suitable transformation existing (which may not be feasible for highly skewed distributions), and makes interpretation of the causal effects more complex than interpretation in the original scale (e.g., log-years instead of years). 

In such a situation, defining the causal effect using a contrast of median potential outcomes could be an appealing alternative. Although the contrast of mean potential outcomes is the most commonly used estimand, the idea of defining the causal effect using the median is not a novel one; the causal effect has been generally defined as a contrast of any functional of the distributions of counterfactual outcomes under different exposure values.\cite{Hernan2004} Despite being a widely acknowledged concept, there is scarce availability of confounding-adjustment methods for estimation of this causal estimand in practice. This may be due to the mathematical properties of the median adding complexity to the development of such methods.

A handful of previous studies have acknowledged the need for such methods, and presented derivations of approaches to estimate causal effects as contrasts of distribution quantiles more generally. A study by Zhang et al. (2012) derived a number of methods - a quantile regression estimator, an inverse-probability weighted estimator, and a stratified estimator using propensity scores. Based on results from a simulation study and application to an empirical dataset, the authors suggested that the quantile regression approach was advantageous when concerned about bias or efficiency.\cite{Zhang2012} A more recent study in the field of environmental science defined a novel estimator which uses a class of balancing weights from functions of the propensity score model to weight each group to a selected target population.\cite{Sun2021} The performance of three estimation methods (propensity score regression, inverse probability weighting and overlap weighting) was explored under homogenous and heterogeneous causal effects and across a number of simulation settings. The study recommended the inverse probability weighting approach if the effect is believed to be heterogeneous, with the other two methods preferred otherwise.\cite{Sun2021} 

Despite these advances, application of these confounding-adjustment methods in health and medical research remains scarce. The current discussion and evaluation of such methods is relatively limited, which has potentially led to a lack of awareness in their existence. In addition, there is minimal understanding of how these methods perform in realistic settings across various scenarios, specifically in terms of the degree of skewness in the outcome which has only been explored minimally and not for all the methods considered here.\cite{Sun2021}

In this paper we aimed to describe and evaluate confounding-adjustment methods to estimate the causal difference in medians, with the aim to increase understanding of their utility and encourage application in practice where appropriate. The methods identified are multivariable quantile regression\cite{Sun2021}, an IPW estimator\cite{Zhang2012} and weighted quantile regression\cite{Sun2021}. We also consider adaptations of the g-computation approach. We begin this paper by defining the causal effect of interest and identifiability assumptions alongside an illustrative example from the Longitudinal Study of Australian Children\cite{Sanson2002} (LSAC). Secondly we outline the confounding-adjustment methods considered. Thirdly, we report findings from a simulation study motivated by the LSAC example, assessing the performance of each method in a realistic setting and considering varying degrees of skewness in the distribution of the outcome variable. The paper concludes by demonstration of the methods applied to the LSAC example, before summarising the key findings, strengths, limitations and practical recommendations.



\section{Defining the causal effect using medians}\label{sec:causal_effect}

Let us consider an observational study with continuous skewed outcome variable $Y$, a binary 0/1 exposure variable $A$, and a vector of $K$ confounder variables $\bm{C}$. We assume that the confounder vector includes only binary or continuous variables, noting that categorical confounders can be represented as a set of binary indicators for the categories. For simplicity of discussion, we have restricted the exposure to be binary and assume that no variable is subject to missingness. 

The example used throughout this paper involves data from a longitudinal cohort study of 4,882 children from the LSAC kindergarten cohort.\cite{Sanson2002, Christensen2017} Children aged 4-5 years were recruited in 2004 (wave 1 of LSAC; approved by the Australian Institute of Family Studies Ethics Committee), with follow-ups every two years in subsequent waves. The example investigation examined the impact of maternal mental health on a child's behaviour in early childhood. The exposure variable ($A$) was a binary indicator variable of probable serious maternal mental illness at wave 1 ($A=1$ if present, $A=0$ otherwise), with the outcome variable ($Y$) being the child's behavioural difficulties as measured by the Strengths and Difficulties Questionnaire\cite{Goodman1997} (SDQ) completed at wave 3. For SDQ scores, higher values indicate increased behavioural difficulties, with the distribution of the SDQ scores in the LSAC cohort being positively skewed (see supplementary material Figure S1). Potential confounders ($\bm{C}$) were measured at wave 1 and included demographic information about the child and mother (see Table \ref{tab:lsac_variables} for a full description of all variables used in this example). 



We define $Y^{A=a}$ to be the potential outcome when the exposure is set to level $a$, where $a=0,1$. In the LSAC example, $Y^a$ would represent the SDQ score for a child when their mother is set to have a probable serious mental illness ($a=1$) versus not ($a=0$). Here we note that the median potential outcome under $A=a$ is defined as the $50^{\mathrm{th}}$ centile from the cumulative distribution function (CDF) $F_a(y)$ of $Y^{A=a}$, assuming that the CDFs are continuous and strictly increasing for $a=0,1$.\cite{Sun2021} Here we define the shorthand notation $m_a$ to denote the median counterfactual outcome under exposure level $A=a$ for simplicity, such that $m_a = m[Y^a]=F^{-1}_{Y^{a}}(0.5)$. Therefore, the causal effect (CE) of interest, denoted by $\delta$, is defined as the difference between the median counterfactual outcomes under the two exposure levels:
\begin{eqnarray}
\delta = F^{-1}_{Y^{a=1}}(0.5) - F^{-1}_{Y^{a=0}(0.5)} = m[Y^{a=1}] - m[Y^{a=0}] = m_1 - m_0. 
\end{eqnarray}
For the LSAC example, $\delta$ represents the difference in median SDQ scores if all children were exposed to maternal mental health problems compared to if none of them were exposed. Here it useful to note that the median of the difference in individual potential outcomes is not equivalent to the difference in median potential outcomes, i.e., $m[Y^{a=1} - Y^{a=0}] \neq m[Y^{a=1}] - m[Y^{a=0}]$. This also holds for other functionals (e.g., risk ratio) as has been noted by Dawid (2000)\cite{Dawid2000}, reiterating the rationale for defining causal effects as contrasts of functionals of the population counterfactual distribution across exposure groups. That is, defining them as effects in a target trial as opposed to defining them as the average of individual causal effects, which are generally not identifiable, even in a trial.\cite{Hernan2008, Hernan2016, Moreno2021} 


The causal difference in medians is identifiable from observational data under the usual assumptions\cite{Greenland1986} and has been shown elsewhere.\cite{Sun2021} Firstly we require the consistency assumption, which states that the exposure $A$ corresponds to a well-defined intervention that in turn corresponds to the versions of the exposure in the data.\cite{Hernan2020} Under these conditions, the potential outcome $Y^a$ is well-defined and would be equal to $Y$ if an individual received exposure level $A=a$ (assumption 1). Secondly, we require the conditional exchangeability assumption given the selected set of confounders, which states that the potential outcome $Y^a$ is independent of the received exposure $A$ given $\bm{C}$, i.e., $Y^a \perp\!\!\!\perp A|\bm{C}$ (assumption 2). Thirdly, we require the positivity assumption, which states that every individual in the population has a positive probability of being exposed or unexposed, that is $Pr[A=a|\bm{C}=\bm{c}] > 0$ for all $\bm{c}$ with positive probability of occurring (assumption 3). Whether these assumptions hold in practice is a matter of debate, however, for the remainder of this paper we assume these conditions to hold.

\begin{center}
	\begin{table}[t!]%
		\caption{Overview of variables from the Longitudinal Study of Australian Children (LSAC) example\cite{Sanson2002, Christensen2017} examining the impact of maternal mental health on a child's behavioural difficulties in early childhood. The exposure and confounder variables are recorded at wave 1 (2004), with the outcome variable recorded at wave 3 (2008). \label{tab:lsac_variables}}
		\centering
		\begingroup \fontsize{9}{14}\selectfont
		\begin{tabular} {lll
			\toprule
			\textbf{Role} & \textbf{Variable} & \textbf{Values and additional details} \\ 
			\midrule
			Outcome (wave 3) & Behavioural difficulties score & Range 0-40; Strengths \& Difficulties Questionnaire \\
			Exposure (wave 1) & Probable serious maternal mental illness & Yes/No; Yes: K10 score $<$ 4\cite{Kessler2002, Andrews2001} \\
			Confounders (wave 1) & Sex of child & Male/Female \\
			& Whether the child has siblings & Yes/No \\
			& Child's physical functioning score & Range 0-100; Pediatric Quality of Life Inventory \\
			& Behavioural difficulties score (baseline) & Range 0-40; Strengths \& Difficulties Questionnaire \\
			& Maternal age & Recorded in years \\
			& Maternal smoking status & Yes/No \\
			& Maternal risky alcohol consumption & Yes/No; Yes: $>$ 2 standard alcoholic drinks per day \\
			& Maternal completion of high school & Yes/No \\
			& Family financial hardship score & Range 0-6 \\
			& Consistent parenting score & Range 1-5 \\
			\bottomrule
		\end{tabular}
		\endgroup
	\end{table}
\end{center}


\section{Confounding-adjustment methods}\label{sec:methods}
Under the aforementioned assumptions, the causal difference in medians $\delta$ can in principle be estimated from observational data using methods that adjust for the confounder variables $\bm{C}$.  Within this section we will introduce the proposed confounding-adjustment methods investigated in our study.


\subsection{Multivariable quantile regression} \label{sec:quant_reg}
In the presence of confounding, multivariable regression is a common approach to estimate the average causal effect (i.e., the causal difference in means in the case of continuous outcomes) under the assumption of a constant causal effect across confounder strata. This approach adjusts for confounding through conditioning on the confounders. To estimate the causal difference in medians $\delta$, a natural adaptation of this approach is to use quantile regression (QR); a common method for modelling the quantiles of the distribution of a random variable conditional on a set of covariates.\cite{Koenker1978} Within this approach, a QR model is fitted to model the outcome variable $Y$ conditional on both the exposure $A$ and $K$ confounder variables $\bm{C}$. Under the fitted QR model, the $\tau^{th}$ quantile is modelled as
\begin{eqnarray}
Q_\tau(Y) = \beta_0(\tau) + \beta_1(\tau)A + \bm{\beta_3}(\tau)\bm{C}, 
\end{eqnarray}
where $\bm{\beta_3}$ is a vector of length $K$. By setting $\tau=0.5$, the coefficient of the exposure variable $\beta_1$ encodes the difference in the conditional median outcome between exposure groups for every substrata, $m[Y|A=1, C=c] - m[Y|A=0, C=c]$. Under assumptions (1-3), the assumption of a constant causal effect across confounder strata and that the regression is correctly specified, the exposure coefficient estimator $\hat{\beta}_1$ has been shown to consistently estimate the CE $\delta$.\cite{Sun2021}

Quantile regression is a widely-applied method in practice and has implementation readily available in statistical software (e.g., the \textit{quantreg} package in R \cite{quantreg2019}; the \textit{qreg} function in STATA) making it an accessible approach. It is important to note however, that the assumption of a constant causal effect across confounder strata may be too simplistic. Therefore less-restrictive confounding-adjustment methods may be required in practice, in particular g-methods. 


\subsection{IPW estimator}
One family of g-methods are those using the framework of inverse probability weighting (IPW).\cite{Horvitz1952, Hernan2000} The motivation behind the IPW framework is to create a pseudo-population in which the distribution of covariates $\bm{C}$ are balanced between exposure groups, such that the association between $A$ and $Y$ in the pseudo-population provides an unbiased estimate for the causal effect of $A$ on $Y$. To create the pseudo-population, observations are re-weighted in a way that is inversely proportional to the probability of their observed exposure group conditional on the confounding variables, with these probabilities estimated from a model for the propensity score (PS), $\pi(\bm{c}) = Pr[A=1|\bm{C}=\bm{c}]$ (e.g., a logistic regression model for $A$ conditional on $\bm{C}$).  

A study by Zhang et al. (2012) derived an IPW estimator for the causal effect defined as the difference in a given quantile of the counterfactual outcome distribution, in particular the median which is our focus here.\cite{Zhang2012} The method can be used to estimate $\delta$ as follows. Recall that the median of the potential outcome under exposure level $a$ is defined as the $50^{th}$ centile from the cumulative counterfactual distribution function $F_a(y)=P(Y^a \leq y)$. Given assumptions (1-3) hold, then for a given exposure level $a$ this is equal to\cite{Zhang2012}
\begin{eqnarray}
F_a(y) = \mathrm{E} \left[ \frac{I(A=a)I(Y \leq y)}{Pr(A=a|\bm{C})} \right]. 
\end{eqnarray}
Under each exposure level $a\in\{0,1\}$, the expectations can be estimated using sample averages\cite{Horvitz1952, Robins1994} and can be regarded as weighted empirical distributions of $Y$.\cite{Zhang2012} The IPW estimator of the median outcome value $m_a$ is therefore defined as the solution to the equation
\begin{eqnarray}
\sum_{i=1}^n W_{a,i} I(Y_i \leq m_a) = 0.5,
\label{eq:ipw_estimator}
\end{eqnarray}
where $W_{a,i}$ denotes the weight for observation $i=1,\ldots,n$ under exposure level $a$ that could be estimated as $1/\hat{\pi}(\bm{c})$ for those who were exposed, and $1/(1-\hat{\pi}(\bm{c}))$ for those who were not exposed, where $\hat{\pi}(\bm{c})$ is an estimate of the propensity score (i.e., probability of being exposed conditional on the observed set of confounders). Following the suggestion of Zhang et al. (2012), normalised weights are preferred to improve finite-sample performance\cite{Zhang2012}, and are calculated by dividing each IP weight by the sum of IP weights in the associated exposure group, such that
\begin{eqnarray}
W_{a,i} = \left. \frac{I(A_i=a)}{Pr(A_i=a|\bm{c}_i)} \ \right/ \sum_{i=1}^{n}\frac{I(A_i=a)}{Pr(A_i=a|\bm{c}_i)}
\end{eqnarray}
To obtain the estimates, $\hat{m}_1$ and $\hat{m}_0$, of the median values, equation \ref{eq:ipw_estimator} is solved for each exposure level $a$, with the weights $W_{a,i}$ substituted by their estimates $\hat{W}_{a,i}$ obtained as above using estimates of the propensity score. Under assumptions (1-3) and assuming that the PS model is correctly specified, the difference between these two values consistently estimates the causal difference in medians $\delta$. 



\subsection{Weighted quantile regression}
An alternative implementation of the IPW approach uses the IP weights to fit an appropriate weighted regression model, that is, weight the score equations of the regression as opposed to the observed outcomes. This can be applied using weighted quantile regression (WQR), which estimates the causal difference in medians by modelling via regression $Y$ on $A$. Using the fitted WQR model, the $\tau^{th}$ quantile is obtained as
\begin{eqnarray}
Q_\tau(Y) = \beta_0^*(\tau) + \beta_1^*(\tau)A, 
\end{eqnarray}
where the coefficient $\beta_1^*(\tau)$ is estimated by solving a weighted sum of check functions.\cite{Sun2021} The weighting is implemented using the same IP weights as outlined for the previous IPW estimator. By setting $\tau=0.5$, the coefficient of the exposure variable encodes the difference in medians between each exposure level in the pseudo-population, $m[Y|A=1] - m[Y|A=0]$. Under assumptions (1-3) and assuming the propensity score model is correctly specified, the estimate $\hat{\beta_1^*}$ has been shown to be a consistent estimator for the causal difference in medians $\delta$.\cite{Sun2021}

Unlike the previous IPW estimator, which needs to be hand-coded at present, implementation of WQR is readily available, with many of the aforementioned software packages enabling specification of weights within the QR implementation (e.g., the \textit{weights} argument within the \textit{rq} function in R\cite{R2019}). 



\subsection{Adaptation of g-computation}

G-computation is another popular method to adjust for confounding when estimating average causal effects.\cite{Robins1986} The method arises from the g-formula which identifies the mean potential outcome under an exposure value as the weighted average of the means across confounder strata given that exposure value, with weights corresponding to the total-sample distribution of the strata. To our knowledge, no similarly intuitive identification formula has been derived for the median potential outcome under an exposure value, and therefore no methods of this type have been examined for the causal difference in medians. 

Here we explored an adaptation of the g-computation method. This was motivated purely from an exploratory stance, considering how practitioners might adapt the g-computation approach in an attempt to estimate the causal difference in medians. For this approach, we consider the definition of the causal effect as essentially aiming to contrast statistics (e.g., medians) that summarise the outcome distribution under each exposure value within a target trial.\cite{Hernan2008, Hernan2016, Moreno2021} In an observational study, where there is lack of randomisation, we cannot assume that summary statistics across observed exposure groups can be directly compared given the imbalance in distribution of outcome predictors. Therefore, by replicating our sample twice and obtaining predicted potential outcomes for each individual under each exposure level we recreate two synthetic exposure groups with the same confounder distribution and thus in which the summary statistics can be compared. Given the outcome is skewed, we assume each conditional distribution (given $\bm{C}$) is skewed also, and therefore obtain those predictions as predicted medians rather than means. The final step requires summarising those outcomes in each synthetic exposure group. 

The method is implemented as follows. Firstly quantile regression is used to model the outcome $Y$ conditional on the exposure $A$ and confounders $\bm{C}$. Using the fitted QR model, predictions of the median outcome are obtained for every observation twice - firstly assuming that each observation was exposed to estimate $m[Y_i|A_i=1, \bm{C}_i=\bm{c}_i]$, and then assuming every observation was unexposed to estimate $m[Y_i|A_i=0, \bm{C}_i=\bm{c}_i]$. Finally, to estimate the median counterfactual outcome for each exposure group separately, and thereby estimate the causal difference in medians, predictions needed to be appropriately aggregated across each respective exposure group. The natural approach would be to take the median value, such that $m_a$ is approximated using the \textit{median} of predictions for under each exposure value separately. Under this approach, $\delta$ was estimated as
\begin{eqnarray}
\hat{m}\left[\hat{m}[Y|A=1, \bm{C}=\bm{c}]\right] - \hat{m}\left[\hat{m}[Y|A=0, \bm{C}=\bm{c}]\right]. 
\end{eqnarray}
However, if the distribution of these predictions is in fact symmetric, then using the median or mean will produce similar estimates. Therefore, we also considered estimating $m_a, a \in (0,1)$ as the \textit{mean} of predictions for each exposure group separately.  The causal difference in medians $\delta$ was therefore estimated as
\begin{eqnarray}
\hat{\mathrm{E}}\left[\hat{m}[Y|A=1, \bm{C}=\bm{c}]\right] - \hat{\mathrm{E}}\left[\hat{m}[Y|A=0, \bm{C}=\bm{c}]\right].
\end{eqnarray}
We considered both variations in our simulation studies to determine their performance (denoted as \textit{g-comp (med)} and \textit{g-comp (mean)} in subsequent sections).

Unlike the previous quantile regression method in Section \ref{sec:quant_reg}, the fitted QR model within this method can incorporate exposure-confounder interactions as the regression coefficients are not interpreted directly, but rather the model is used for prediction. Therefore, like other g-computation methods, this method does not rely on the assumption of a constant causal effect across strata. Given assumptions (1-3) hold and that the fitted QR model is correctly specified, we conjectured that the causal difference in medians $\delta$ could be approximated using this adaptation of the g-computation approach. We empirically evaluated this in the simulation study below.  



\section{Simulation Study}\label{sec:simulation_study}


A simulation study was conducted motivated by the LSAC example, to investigate the performance of each confounding-adjustment method in a realistic setting and under varying degrees of skewness in the outcome variable.  We explored four different skewness scenarios (i.e., different skewed distributions for $Y$), with the simulated datasets generated as detailed below. A full outline of each variable and their generating distribution is provided in the supplementary material (Table S1). 

We generated 1000 samples consisting of 1000 observations for each of the four skewness scenarios. For each scenario, dataset and record, five confounder variables $C_k$ for $k=1,\dots,5$ (three binary and two continuous) were generated sequentially in that order based on variables in the LSAC data set. Observations for the binary confounders $C_1, C_3$ and $C_4$ (based on sex, maternal education and financial hardship) were generated from a binomial distribution with the success probability defined by a logistic regression model using all previous confounders as predictors. For the continuous confounders $C_2$ and $C_5$ (maternal age and log-transformed baseline SDQ) a similar approach was used; values were generated from normal distributions with the means defined by a linear regression model including all previous confounders as predictors. 

Secondly, a binary exposure $A$ was generated for each record conditional on the values of the five confounder variables. The value for the exposure was generated from a binomial distribution with success probability defined by a logistic regression model including all confounders as predictors. A skewed continuous outcome $Y$ was then generated for each record. Values for log($Y$) were generated from a normal distribution with the mean defined by a linear regression model including the exposure and all confounders as predictors. The linear regression model was specified to include all main effects and two exposure-confounder interaction terms. Different skewed distributions in the outcome variable were established by setting the standard deviation in the generating normal distribution for log($Y$) to $\sigma=0.75, 1, 1.25, 1.5$ for each of the four increasing skewness scenarios, respectively. Values for log($Y$) were then exponentiated to obtain the outcome value $Y$, with the corresponding distribution of $Y$ being positively skewed (see Figure S2 in the supplementary material for a visualisation of these distributions). Here we note that $\sigma$ characterises other properties of the outcome distribution beyond the measure of skewness. Therefore, findings across skewness scenarios are not expected to behave in a monotonic pattern, for example, increasing bias with increasing skewness. 

Aside from varying skewness, a further consideration of this study was the strength of confounding bias in the simulated datasets and whether this impacted the methods' performance. Therefore the data generation process was conducted under two different settings - weak confounding bias (approximately 10\% relative bias in the unadjusted estimate relative to the true value) and strong confounding bias (approximately 20\% relative bias). Modifications to coefficients in the outcome-generating model were used to achieve this.

The true causal difference in medians $\delta$ under this data generation set-up was computed by empirical methods\cite{Sun2021}, with further details about how we implemented this procedure and the true values provided in the supplementary material (Section S2.3 and Table S2). Coefficients in the outcome-generating model were modified to ensure the true value was large enough to be estimable with the given sample size with adequate power (approximately 80\%) in an unadjusted analyses.

The confounding-adjustment methods outlined in Section \ref{sec:methods} (alongside an unadjusted approach contrasting the medians across exposure groups in the observed data) were applied to each simulated dataset to estimate the causal difference in medians $\delta$. For the multivariable quantile regression approach, the outcome model included $A$ and $\bm{C}$ as predictors with main effect terms only. For both the IPW estimator and the weighted quantile regression approaches, the propensity score model regressed $A$ on the confounders, including main effect terms only. For both g-comp (med) and g-comp( mean), the outcome model was specified as $Y$ conditional on $A$ and $\bm{C}$, including all two-way interaction terms between $A$ and each confounder in $\bm{C}$. 

Based on previous recommendations\cite{Zhang2012} and for ease of implementation and comparability, standard errors and confidence intervals (CIs) were estimated using bootstrap procedures (using the percentile method for CIs) for all methods. For each confounding strength and skewness scenario, metrics assessing the performance of each method (e.g., bias, coverage, empirical and model-based standard errors for $\hat{\delta}$) were calculated and averaged over the 1000 samples, with key focus on bias in the point estimate. Monte Carlo standard errors were estimated for each metric. All analysis was conducted in R 4.0.2\cite{R2019}, using the \textit{quantreg}\cite{quantreg2019} and \textit{boot} packages within self-developed code (available at \url{https://github.com/daisyshep/CI-medians.git}), with the IPW estimator implemented using R code supplied in the method's source paper.\cite{Zhang2012}    




\subsection{Results from simulation study}\label{sec:results}

Within the simulation study, estimates for $\delta$ differed depending on the confounding-adjustment method used for analysis, with the variation in estimates increasing with increasing skewness for both confounding strengths (Figure \ref{fig:sim_results}). Estimates obtained using quantile regression (QR) were biased across all skewness scenarios (relative bias range: 6.6\% to 8.1\% weak confounding, 9.0\% to 14.4\% strong confounding; Table \ref{tab:sim_results}); an expected result given the method's strict assumption of a constant causal effect across stratum. Estimates were heavily biased for adaptions of the g-computation method (median and mean variations) despite relaxing this assumption, with the average bias exceeding that of the unadjusted estimate in the majority of settings. Notably both g-comp (med) and g-comp (mean) had similar performance to one another across all skewness scenarios and both confounding settings. In contrast, both the IPW estimator and weighted quantile regression performed well and had similar performance to one another, with minimal bias in estimates for $\delta$ (relative bias < 5\% in the majority of skewness scenarios). These observations were consistent across both the weak and strong confounding settings. 

There was little bias in estimating the model standard error (SE) with the bootstrap when using QR, the IPW estimator or weighted QR, although adaptations of g-computation (both variations) consistently overestimated the model SE by a larger margin (relative error in model SE > 10\% in the majority of skewness scenarios). The coverage probability was close to nominal level for all confounding-adjustment methods (range: 93.70\% to 97.70\%) across all skewness scenarios and both confounding settings. Monte Carlo SEs of performance estimates were all acceptable (supplementary material Table S3), and therefore we feel comfortable drawing conclusions about the methods based on the 1000 samples used per skewness scenario. 

\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.65]{sim_results_plot}
	\caption{Calculated bias in causal difference estimates obtained under each method, skewness scenario (strength of skewness increases with scenario number, with scenario 1 and 4 corresponding to the weakest and highest skewness, respectively) and confounding setting in the simulated datasets (1000 datasets per skewness scenario) analysed under each method. \label{fig:sim_results}}
\end{figure}


\begin{table}[htbp]
	\caption{Performance of confounding-adjustment methods across confounding and skewness scenarios in the simulation study with maximum Monte Carlo standard errors (SE) provided in the table footnote (see supplementary material Table S3 for all Monte Carlo SEs). \label{tab:sim_results}}
	\centering
	\begingroup \fontsize{9}{12}\selectfont
	\begin{tabular}{p{1.9cm} p{1.4cm} r >{\raggedleft\arraybackslash}p{1.2cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.25cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm}}
		\toprule
		\textbf{Confounding} & \textbf{Skewness scenario} & \textbf{Method} & \textbf{Bias} & \textbf{Relative bias} (\%) & \textbf{Empirical SE} & \textbf{Model SE} & \textbf{Error SE} (\%) & \textbf{Coverage} (\%) \\ 
		\toprule
		
		Weak & 1 & Unadjusted & 0.090 & 10.07 & 0.369 & 0.381 & 3.21 & 94.60 \\ 	
		& & QR & 0.059 & 6.63 & 0.384 & 0.375 & -2.24 & 94.80 \\ 
		& & IPW estimator & -0.004 & -0.48 & 0.442 & 0.454 & 2.78 & 94.90 \\ 		
		& & Weighted QR & -0.023 & -2.53 & 0.440 & 0.453 & 2.98 & 95.20 \\ 
		& & G-comp (med) & 0.083 & 9.24 & 0.428 & 0.467 & 9.09 & 97.10 \\ 				
		& & G-comp (mean) & 0.072 & 8.10 & 0.419 & 0.454 & 8.55 & 96.90 \\ \hline
		
		Weak & 2 & Unadjusted & 0.123 & 10.08 & 0.526 & 0.528 & 0.32 & 94.90 \\ 
		& & QR & 0.099 & 8.13 & 0.540 & 0.517 & -4.35 & 94.20 \\ 
		& & IPW estimator & 0.050 & 4.09 & 0.625 & 0.640 & 2.43 & 93.70 \\ 
		& & Weighted QR & 0.023 & 1.91 & 0.621 & 0.636 & 2.54 & 94.20 \\ 		
		& & G-comp (med) & 0.164 & 13.43 & 0.604 & 0.665 & 10.10 & 96.40 \\ 				
		& & G-comp (mean) & 0.151 & 12.35 & 0.594 & 0.650 & 9.31 & 95.40 \\ \hline

		Weak & 3 & Unadjusted & 0.160 & 10.00 & 0.684 & 0.727 & 6.27 & 95.60 \\ 
		& & QR & 0.114 & 7.14 & 0.693 & 0.695 & 0.31 & 94.70 \\ 
		& & IPW estimator & 0.019 & 1.17 & 0.806 & 0.870 & 7.94 & 95.00 \\ 		
		& & Weighted QR & -0.020 & -1.23 & 0.803 & 0.864 & 7.62 & 95.20 \\ 		
		& & G-comp (med) & 0.246 & 15.39 & 0.779 & 0.925 & 18.74 & 97.70 \\ 		
		& & G-comp (mean) & 0.239 & 14.95 & 0.763 & 0.899 & 17.85 & 97.30 \\ \hline
		
		Weak & 4 & Unadjusted & 0.193 & 10.09 & 0.833 & 0.894 & 7.24 & 96.80 \\ 
		& & QR & 0.133 & 6.95 & 0.864 & 0.842 & -2.62 & 95.50 \\ 
		& & IPW estimator & 0.123 & 6.45 & 1.032 & 1.098 & 6.40 & 95.30 \\ 		
		& & Weighted QR & 0.081 & 4.22 & 1.017 & 1.091 & 7.24 & 95.50 \\ 
		& & G-comp (med) & 0.324 & 16.94 & 1.014 & 1.173 & 15.66 & 95.70 \\ 				
		& & G-comp (mean) & 0.319 & 16.68 & 1.015 & 1.153 & 13.54 & 95.00 \\ \midrule
		
		Strong & 1 & Unadjusted & 0.171 & 20.17 & 0.391 & 0.398 & 1.665 & 93.20 \\ 
		& & QR & 0.098 & 11.58 & 0.398 & 0.390 & -2.10 & 94.30 \\ 
		& & IPW estimator & 0.003 & 0.31 & 0.441 & 0.466 & 5.75 & 95.40 \\ 		 
		& & Weighted QR & -0.017 & -1.95 & 0.440 & 0.464 & 5.34 & 95.30 \\ 
		& & G-comp (med) & 0.129 & 15.16 & 0.433 & 0.484 & 11.76 & 96.80 \\
		& & G-comp (mean) & 0.110 & 12.95 & 0.418 & 0.468 & 11.99 & 97.10 \\ \hline

		Strong & 2 & Unadjusted & 0.240 & 20.05 & 0.544 & 0.568 & 4.40 & 94.00 \\ 
		& & QR & 0.117 & 9.81 & 0.556 & 0.556 & -0.03 & 95.00 \\ 
		& & IPW estimator & 0.030 & 2.55 & 0.621 & 0.673 & 8.34 & 96.00 \\ 
		& & Weighted QR & 0.004 & 0.33 & 0.617 & 0.670 & 8.61 & 95.40 \\ 
		& & G-comp (med) & 0.186 & 15.58 & 0.620 & 0.710 & 14.53 & 97.20 \\ 
		& & G-comp (mean) & 0.172 & 14.36 & 0.604 & 0.692 & 14.61 & 97.00 \\ \hline

		Strong & 3 & Unadjusted & 0.307 & 20.11 & 0.729 & 0.734 & 0.72 & 93.50 \\ 
		& & QR & 0.219 & 14.39 & 0.726 & 0.709 & -2.32 & 93.70 \\ 
		& & IPW estimator & 0.117 & 7.68 & 0.847 & 0.897 & 5.90 & 95.40 \\ 		
		& & Weighted QR & 0.082 & 5.35 & 0.843 & 0.891 & 5.78 & 95.50 \\ 		
		& & G-comp (med) & 0.339 & 22.20 & 0.851 & 0.953 & 11.98 & 95.10 \\ 
		& & G-comp (mean) & 0.318 & 20.83 & 0.821 & 0.930 & 13.32 & 95.50 \\ \hline

		Strong & 4 & Unadjusted & 0.420 & 20.02 & 1.004 & 1.052 & 4.76 & 95.10 \\ 
		& & QR & 0.189 & 8.98 & 1.002 & 0.993 & -0.88 & 94.50 \\ 
		& & IPW estimator & 0.038 & 1.81 & 1.189 & 1.252 & 5.37 & 95.50 \\ 		
		& & Weighted QR & -0.009 & -0.43 & 1.180 & 1.241 & 5.19 & 95.50 \\ 		
		& & G-comp (med) & 0.384 & 18.31 & 1.151 & 1.361 & 18.28 & 96.80 \\ 
		& & G-comp (mean) & 0.363 & 17.28 & 1.126 & 1.322 & 17.35 & 95.90  \\ 
		\bottomrule
	\end{tabular}
	\begin{tablenotes
		\item Maximum Monte Carlo SE (performance measure): 0.038 (bias), 0.010 (empirical SE), 
		0.029 (model SE), 6.956\% (relative error in model SE), 0.796\% (coverage). 
	\end{tablenotes}
	\endgroup
\end{table}




\section{Application to LSAC}\label{sec:lsac}

Methods were applied to the LSAC data set using complete cases only (n=3,245). For the multivariable quantile regression approach, the outcome model included $A$ and $\bm{C}$ (consisting of 9 confounders; Table \ref{tab:lsac_variables}) as predictors with main effect terms only. The propensity score model (used for the IPW estimator and weighted quantile regression) regressed $A$ on $\bm{C}$ including main effect terms only. For both g-comp (med) and g-comp (mean), the outcome model was specified as $Y$ conditional on $A$ and $\bm{C}$, including two-way interaction terms between $A$ and three confounders (child's sex, maternal completion of high school, consistent parenting).  

Applying the confounding-adjustment methods to the LSAC example yielded a lower estimate of the causal difference in medians for all methods comparatively to the unadjusted approach (Figure \ref{fig:lsac_results}). Estimated effects were more consistent across the methods than observed within the simulation study, although the quantile regression approach and both g-comp (med) and g-comp (mean) yielded slighty lower estimates of $\delta$ than the IPW estimator and weighted quantile regression approaches. Notably both the IPW estimator and the weighted quantile regression approach produced the same point estimates and bootstrap confidence intervals as one another when applied to the LSAC example.  


\begin{figure}[t!]
	\centering
	\includegraphics[scale=0.8]{LSAC_results_CCA}
	\caption{Estimated causal difference in medians $\hat{\delta}$ for the Longitudinal Study Of Australian Children (LSAC) example\cite{Sanson2002, Christensen2017} obtained under each confounding-adjustment method. Point estimates and their corresponding 95\% confidence intervals are presented alongside the figure.\label{fig:lsac_results}}
\end{figure}


\section{Discussion}\label{sec:conclusion}

In the presence of skewed outcome data, the common approach to "ignore or transform" may not be optimal, and defining the causal effect using a contrast of median potential outcomes may be more appropriate. Despite being a widely acknowledged concept, there is scarce availability of confounding-adjustment methods to estimate the causal difference in medians in practice. A handful of previous studies have proposed potential approaches to estimate causal effects as contrasts of distribution quantiles more generally\cite{Sun2021, Zhang2012}, although the application of these methods in health and medical studies remains scarce. 

In this paper we described and evaluated confounding-adjustment methods to estimate the causal difference in medians. We identified methods from previous literature (multivariable quantile regression\cite{Sun2021}, an IPW estimator\cite{Zhang2012} and weighted quantile regression\cite{Sun2021}) alongside adaptations of the g-computation approach. Methods were selected and described with a key focus on their accessibility and ease of implementation in a practical setting, to encourage their use in practice where applicable. 

Results from the simulation study indicated varied performance of the confounding-adjustment methods when estimating the causal difference in medians. Results showed support for the IPW estimator and weighted quantile regression methods, with both producing estimates with minimal bias, aligning with results seen elsewhere.\cite{Sun2021} Meanwhile, as anticipated, the multivariable quantile regression was too simplistic for the nature of the simulated datasets (e.g., assumption of a constant causal effect across stratum), and therefore resulted in biased estimates for the causal difference in medians. The adaptations of the g-computation approach also showed poor performance, even though they relaxed the strict assumption of a constant causal effect across stratum. When applied to the LSAC case study, estimated effects were more consistent across methods than observed within the simulation study, with the IPW estimator and weighted quantile regression producing slightly higher estimates than the other confounding-adjustment methods. 

These findings need to be interpreted in light of the fact that our data generation approach may have favoured the IPW-based approaches (the IPW estimator and weighted quantile regression) over other approaches. The data generation mechanism used within the simulation study ensured the propensity score model was correctly specified (a key assumption of the method), which may have contributed to these methods' good performance. This aligned with the previous findings of Zhang et al. (2012), which indicated the necessity of the propensity model to be correctly specified.\cite{Zhang2012} In contrast, the data-generation process used within our simulation study did not guarantee correct specification of the outcome model used in the multivariable quantile regression or adaptations of the g-computation approach, which is an assumption of each method.

In order to address this limitation, we explored a second data generation process (see supplementary material Section S3 for full details and results). Comparatively to the approach presented in this paper, this alternative data-generation approach did not guarantee correct specification of the propensity score model. Results under this approach indicated the g-comp (med) method performed well, while the IPW estimator and weighted quantile regression methods did not; an opposite observation to the findings under our main simulation design. We consider the findings for the adaptations of the g-computation approach to be less conclusive, due to no guarantee that the outcome model was correctly specified in either simulation design. We also note that the bootstrap consistently overestimated the model standard error for both g-comp (med) and g-comp (mean), which may suggest that the bootstrap procedure is not appropriate for these approaches.

A strength of this work was the design of our simulation study which was motivated by the LSAC example, allowing us to investigate the performance of these methods in a range of realistic settings. These settings considered a range of varying skewness distributions in the outcome variable, alongside two different strengths of confounding bias present in the datasets, covering more complex and realistic scenarios than those explored in previous papers.\cite{Zhang2012} 

A potential limitation of our study was the restriction to singly robust methods only, which rely on correct specification of the respective model. However, previous studies have presented a handful of promising doubly robust methods.\cite{Zhang2012,Xie2020,Xu2018,Diaz2017} Zhang et al. (2012) derived an adaptation of the augmented inverse probability weighted estimator and a further adaptation incorporating stratification based on the propensity score, with results suggesting both estimators are preferable over the singly robust alternative if precision is of concern.\cite{Zhang2012} A later study by Xu et al. (2018) compared the performance of the doubly robust estimators of Zhang et al.\cite{Zhang2012} and a proposed targeted maximum likelihood estimator (TMLE) for quantiles by D{\'i}az (2017)\cite{Diaz2017} with a Bayesian non-parametric (BNP) approach for causal inference on quantiles.\cite{Xu2018} Results indicated the BNP performed consistently well, whilst the frequentist approaches suffered with under-coverage and large bias under some simulation scenarios, and the TMLE approach had high computation times and bias under some simulation settings.\cite{Xu2018} Xie et al.\cite{Xie2020} proposed an additional multiply robust estimator, using multiple outcome and propensity score models, and found the estimator to be less biased under model misspecification compared to the doubly robust estimators of Zhang et al.\cite{Zhang2012}. Results from these studies indicated the proposed doubly robust methods performed well, with the natural advantage being the method's robustness to model misspecification. However, the doubly robust methods have not yet been evaluated in a range of complex and realistic scenarios as is done here for singly robust methods, with their implementation in software not as readily available. Both of these factors would be useful to pursue in future work.  

Within this study, we restricted our simulation design to explore a binary exposure only. Therefore, potential future directions could explore the performance of methods in the presence of a continuous, or multi-level exposure variables. In addition, the outcome distributions explored within this study were all heavily right-skewed, so results may represent an extreme representation of the methods' performance under heavily skewed outcome data. Furthermore, we also reiterate the need for the causal estimand of interest to be driven by the research question at hand. Even in the presence of skewed outcome data, it may be appropriate or preferred to estimate the average causal effect defined using a contrast of mean potential outcomes. Therefore we caution the reader to not be purely data-driven, but to consider the usefulness of the target estimand when drawing conclusions about their study. 

In conclusion, when estimating the causal difference in medians the IPW estimator or weighted quantile regression present promising approaches, provided a richly specified model is used such that correct specification of the propensity score model is likely. Implementation of both methods is readily available and therefore accessible to implement (e.g., \textit{rq} function in R\cite{R2019}, open-source code of the IPW estimator\cite{Zhang2012}). These methods provide appealing alternatives to the common "ignore or transform" approach or the stringent constant-effect assumption of multivariable quantile regression, enhancing our capability to obtain meaningful causal effect estimates with skewed outcome data. 



\subsection*{Author contributions}

Development of the project and study design was conducted by DS and MMB. DS conducted the simulation study and data analysis, with regular input from MMB. The manuscript was developed by DS, with critical input provided by MMB. 

\subsection*{Financial disclosure}
MMB was supported by an Australian Research Council Discovery Early Career Researcher Award (project number DE190101326) funded by the Australian Government and a NHMRC Investigator Grant (ID 2009572). 

\subsection*{Data availability statement}

Data used for the illustrative example used in this paper are available from the DSS. Restrictions apply to the availability of these data, which were used under license for this study. Data are available at \href{https://dataverse.ada.edu.au/dataset.xhtml?persistentId=doi:10.26193/F2YRL5} with the permission of the DSS.

\subsection*{Conflict of interest}

The authors declare no potential conflict of interests.





\section{The LSAC case study}

\begin{figure}[htbp]
	\centering
	\includegraphics[scale=0.45]{LSAC_outcome_distribution}
	\caption{Distribution of the Strengths \& Difficulties questionnaire (SDQ) scores at wave 3 for the Longitudinal Study Of Australian Children (LSAC) case study. \cite{Sanson2002, Christensen2017}\label{fig:LSAC_outcome}}
\end{figure}

\bigskip


\section{Additional simulation study details}


\subsection{Data generation models}
\begin{table}[htbp]
	\caption{Detail on the models used to generate variables for the simulation study, with the structure based on the LSAC dataset unless otherwise specified in Section 4 of the main text.  \label{tab:datagen_tab}}
	\begingroup\fontsize{9}{14}\selectfont
	\begin{tabular}{p{1.5cm} l l p{7cm}}
		\hline
		\textbf{Variable} & \textbf{LSAC variable} & \textbf{Generating distribution} & \textbf{Additional details} \\
		\hline
		$C_1$ & Sex & $C_1 \sim$ Binomial(1, 0.51) & \\
		$C_2$ & Maternal age & $C_2 \sim$ Normal(35.17, 5.47) & \\
		$C_3$ & Maternal education & $C_3 \sim$ Binomial(1, $Pr(C_3=1)$) & logit$[Pr(C_3=1)] = -1.41 + 0.78C_1 + 0.04C_2$ \\
		$C_4$ & Financial hardship$^a$ & $C_4 \sim$ Binomial(1, $Pr(C_4=1)$) & logit$[Pr(C_4=1)] = -1.55 + 0.47C_1 + 0.03C_2 + 0.80C_3$ \\
		$C_5$ & Baseline SDQ (logged) & $C_5 \sim$ Normal($\mu_{C_5}$, 0.63) & $\mu_{C_5} = 1.91 + 0.03C_1 + 0.01C_2 + 0.05C_3 + 0.12C_4$ \\
		$A$ & Maternal mental health & $A \sim$ Binomial(1, $Pr(A=1)$) & logit$[Pr(A=1)] = -2.39 + 0.04C_1 - 0.05C_2 - 0.09C_3 + 0.51C_4 + 1.07C_5$\\
		log($Y$) & SDQ score & log($Y$) $\sim$ Normal($\mu_Y, \sigma$) & $\mu_Y = 1.40 + 0.49A + 0.03C_1 - 0.01C_2 + 0.01C_3 + 0.03C_4 + 0.26C_5 + 0.12AC_1 - 0.01AC_2$ \\
		& & & $\sigma=0.75, 1, 1.25, 1.5$\\
		\hline
	
	\end{tabular} \\
	$^a$Dichotomised: No (score = 0 on original scale), Yes (score $>$ 0 on original scale)
	\endgroup
	
\end{table}


\newpage



\subsection{Distribution of the outcome variable}

\begin{figure}[htbp]
	\centering
	\includegraphics[scale=0.7]{outcome_dist}
	\caption{Distribution of the outcome variable ($Y$) under each skewness scenario and confounding bias strength used within the simulation study. \label{fig:skew_true}}
\end{figure}



\subsection{True difference in medians} \label{app:true}

Within the simulation study, the true difference in medians $\delta$, with respect to which assess bias, was computed by empirical methods outlined in the supplementary material of Sun et al. (2021).\cite{Sun2021} Specific details of our implementation are as follows. Initially we generated a large dataset (1,000,000 observations) per skewness scenario and strength of confounding bias (weak or strong). For each large dataset, a quantile regression model was fitted under different quantiles $\tau$; we used 200 different values for $\tau$ equally distributed across the range [0.05, 0.95]. For our vector of $m_a$ support values, we covered the range [3,8] in increments of 0.005. Different values of $\delta$ were obtained for each skewness scenario and under the two different strengths of confounding bias, as outlined below in Table \ref{app:true_value}.\\

\begin{table}[htbp!]
	\centering
	\caption{True values for the causal difference in medians $\delta$ obtained under each skewness scenario and confounding setting. \label{app:true_value}}
	\begin{tabular}{ccccc}
		\toprule
		\textbf{Confounding bias} & \multicolumn{4}{c}{\textbf{Skewness scenario}} \\ 
		\textbf{strength} & \textbf{1} & \textbf{2} & \textbf{3} & \textbf{4} \\
		\midrule
		Weak & 0.895 & 1.220 & 1.600 & 1.910 \\
		Strong & 0.850 & 1.195 & 1.525 & 2.100 \\
		\bottomrule
	\end{tabular}
\end{table}


\newpage
\subsection{Monte Carlo standard errors}

\begin{table}[htbp!]
	\centering
	\caption{Monte Carlo standard errors of performance estimates calculated over the 1000 simulated datasets per skewness scenario and confounding bias strength under each of the confounding-adjustment methods. \label{tab:mc_results}}	
	\begingroup\fontsize{9}{11.3}\selectfont
	\begin{tabular}{p{2cm} p{1.5cm} r >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm}}
		\toprule
		\textbf{Confounding} & \textbf{Skewness scenario} & \textbf{Method} & \textbf{Bias} & \textbf{Empirical SE} & \textbf{Model SE} & \textbf{Relative error SE (\%)} & \textbf{Coverage (\%)} \\ 
		\midrule
		Weak & 1 & Unadjusted & 0.012 & 0.008 & 0.002 & 2.387 & 0.715 \\ 
		& & QR & 0.012 & 0.009 & 0.002 & 2.234 & 0.702 \\ 
		& & IPW estimator & 0.014 & 0.010 & 0.003 & 2.414 & 0.696 \\ 
		& & Weighted QR & 0.014 & 0.010 & 0.003 & 2.420 & 0.676 \\ 
		& & G-comp (med) & 0.014 & 0.010 & 0.003 & 2.510 & 0.531 \\ 
		& & G-comp (mean) & 0.013 & 0.009 & 0.003 & 2.502 & 0.548 \\ \hline
		
		Weak & 2 & Unadjusted & 0.017 & 0.008 & 0.004 & 2.607 & 0.715 \\ 
		& & QR & 0.017 & 0.009 & 0.003 & 2.364 & 0.702 \\ 
		& & IPW estimator & 0.020 & 0.010 & 0.007 & 2.768 & 0.696 \\ 
		& & Weighted QR & 0.020 & 0.010 & 0.007 & 2.750 & 0.676 \\ 
		& & G-comp (med) & 0.019 & 0.010 & 0.005 & 2.721 & 0.531 \\ 
		& & G-comp (mean) & 0.019 & 0.009 & 0.005 & 2.711 & 0.548 \\ \hline
		
		Weak & 3 & Unadjusted & 0.022 & 0.008 & 0.009 & 3.352 & 0.715 \\ 
		& & QR & 0.022 & 0.009 & 0.007 & 2.783 & 0.702 \\ 
		& & IPW estimator & 0.025 & 0.010 & 0.013 & 3.735 & 0.696 \\ 
		& & Weighted QR & 0.025 & 0.010 & 0.013 & 3.744 & 0.676 \\ 
		& & G-comp (med) & 0.025 & 0.010 & 0.011 & 3.605 & 0.531 \\ 
		& & G-comp (mean) & 0.024 & 0.009 & 0.011 & 3.567 & 0.548 \\ \hline
		
		Weak & 4 & Unadjusted & 0.026 & 0.008 & 0.014 & 4.412 & 0.715 \\ 
		& & QR & 0.027 & 0.009 & 0.010 & 3.405 & 0.702 \\ 
		& & IPW estimator & 0.033 & 0.010 & 0.022 & 5.563 & 0.696 \\ 
		& & Weighted QR & 0.032 & 0.010 & 0.022 & 5.590 & 0.676 \\ 
		& & G-comp (med) & 0.032 & 0.010 & 0.019 & 5.139 & 0.531 \\ 
		& & G-comp (mean) & 0.032 & 0.009 & 0.019 & 5.210 & 0.548 \\ \midrule
		
		Strong & 1 & Unadjusted & 0.012 & 0.009 & 0.002 & 2.357 & 0.796 \\ 
		& & QR & 0.013 & 0.009 & 0.002 & 2.242 & 0.733 \\ 
		& & IPW estimator & 0.014 & 0.010 & 0.003 & 2.491 & 0.662 \\ 
		& & Weighted QR & 0.014 & 0.010 & 0.003 & 2.480 & 0.669 \\ 
		& & G-comp (med) & 0.014 & 0.010 & 0.002 & 2.559 & 0.557 \\ 
		& & G-comp (mean) & 0.013 & 0.009 & 0.002 & 2.561 & 0.531 \\ \hline
		
		Strong & 2 & Unadjusted & 0.017 & 0.009 & 0.005 & 2.582 & 0.796 \\ 
		& & QR & 0.018 & 0.009 & 0.004 & 2.397 & 0.733 \\ 
		& & IPW estimator & 0.020 & 0.010 & 0.007 & 2.879 & 0.662 \\ 
		& & Weighted QR & 0.020 & 0.010 & 0.007 & 2.869 & 0.669 \\ 
		& & G-comp (med) & 0.020 & 0.010 & 0.006 & 2.852 & 0.557 \\ 
		& & G-comp (mean) & 0.019 & 0.009 & 0.006 & 2.877 & 0.531 \\ \hline
		
		Strong & 3 & Unadjusted & 0.023 & 0.009 & 0.008 & 3.106 & 0.796 \\ 
		& & QR & 0.023 & 0.009 & 0.007 & 2.734 & 0.733 \\ 
		& & IPW estimator & 0.027 & 0.010 & 0.017 & 4.559 & 0.662 \\ 
		& & Weighted QR & 0.027 & 0.010 & 0.017 & 4.532 & 0.669 \\ 
		& & G-comp (med) & 0.027 & 0.010 & 0.012 & 3.780 & 0.557 \\ 
		& & G-comp (mean) & 0.026 & 0.009 & 0.012 & 3.826 & 0.531 \\ \hline
		
		Strong & 4 & Unadjusted & 0.032 & 0.009 & 0.020 & 5.633 & 0.796 \\ 
		& & QR & 0.032 & 0.009 & 0.014 & 4.172 & 0.733 \\ 
		& & IPW estimator & 0.038 & 0.010 & 0.029 & 6.956 & 0.662 \\ 
		& & Weighted QR & 0.037 & 0.010 & 0.028 & 6.831 & 0.669 \\ 
		& & G-comp (med) & 0.036 & 0.010 & 0.026 & 6.438 & 0.557 \\ 
		& & G-comp (mean) & 0.036 & 0.009 & 0.024 & 6.367 & 0.531 \\
		\bottomrule
	\end{tabular}
	\endgroup
\end{table}


\newpage

\section{Alternative data generation approach} \label{app:datagen2}

Alongside the data generation approach presented in the manuscript, we explored an additional method to generate datasets within the simulation study. This second approach was motivated by the lack of guarantee of the outcome model being correctly specified under the current data generation approach, and attempted to overcome this. The majority of the simulation study design was the same as that outlined in the manuscript (Section 2), except for how the outcome variable $Y$ was generated. The process of generating $Y$ and the corresponding simulation study results are presented below.   


\subsection{Generating the outcome variable}
A key characteristic when generating the simulated datasets was ensuring the true difference in medians $\delta$ with respect to which assess bias was known. Within this data generation approach, the DAG in Figure \ref{fig:dag} depicts the structure of the data generation process which enabled $\delta$ to be directly set, instead of needing to use empirical methods to derive it as for the other generating approach (i.e., Section \ref{app:true}). Initially we generated 1000 samples consisting of 1200 observations for each skewness scenario. The five confounder variables $\bm{C}$ and binary exposure variable $A$ were generated as outlined in Section 2 and using the same data generation models detailed in Table \ref{tab:datagen_tab}. 

\begin{figure}[htbp!]
	\begin{center}
		\includegraphics[scale=0.6]{DAG_datagen}
		\caption{Directed acyclic graph (DAG) depicting the data generation process used within the simulation study, where $\bm{C}$ denotes the confounder variables, $A$ the binary exposure, $Y$ the skewed outcome variable and $M$ is an indicator of missing information in $Y$. \label{fig:dag}}
	\end{center}
\end{figure}

A skewed continuous outcome $Y$ was generated for each of the 1200 observations within a sample. Values for $Y$ were generated from a log-normal distribution conditional on the exposure level only. Different skewed distributions in the outcome variable were established by varying the parameters of the log-normal distribution as outlined in Table \ref{tab:lognormal} (see Figure \ref{fig:skew_true_alt} for a visualisation of these distributions). At this stage there was no confounding in the data and the estimated difference in medians between the two groups would be unbiased for the causal difference in medians. 

In order to introduce confounding, i.e. open backdoor paths between $A$ and $Y$, missing values were introduced into the outcome variable and then records with missing outcome were deleted (missingness indicator denoted as $M$ in the DAG, Figure \ref{fig:dag}). This in effect means we are conditioning on a collider, opening the path between $\bm{C}$ and $Y$, and thus a backdoor path between $A$ and $Y$  that could be closed by  adjusting for $\bm{C}$. To generate those missing values, each observation was assigned a probability that the value $Y$ was missing, with the probability defined by a logistic regression model including all confounders and the outcome, such that

\begin{eqnarray}
\mathrm{logit}(Pr[M=1]) = \beta_0 + \beta_1C_1 + \beta_2C_2 + \beta_3C_3 + \beta_4C_4 + \beta_5C_5 + \beta_6Y_{std}.
\end{eqnarray}
A standardised version of the outcome $Y_{std}$ was used instead of $Y$ as a covariate in this model to improve comparability between scenarios with differing degrees of skewness in the outcome variable. Coefficients $\beta_1, \dots, \beta_5$ were estimated as outlined in Table \ref{tab:datagen_tab}. The value of $\beta_0$ was fixed such that the probability of $Y$ missing was approximately 0.2. In actuality, the 1000 observations with the smallest missingness probabilities were retained in the dataset, with the remaining observations discarded. The value of $\beta_6$ was selected to obtain a certain strength of confounding bias (weak or strong; defined in Section 4 of the manuscript) in the simulated datasets. 

\begin{table}[t!]
	\centering
	\caption{Parameters for the generating distribution of the outcome variable $Y$. \label{tab:lognormal}}
	\begin{tabular}{lcccc}
		\toprule
		\multirow{2}{*}{\textbf{Skewness scenario}} & \multicolumn{2}{c}{\textbf{Exposure:} $A=0$} & \multicolumn{2}{c}{\textbf{Exposure: }$A=1$}                \\
		& $\mu_1$ & $\sigma_1$ & $\mu_2$ & $\sigma_2$ \\
		\midrule
		1 & log(7) & 0.25 & log(8) & 0.25 \\
		2 & log(7) & 0.5 & log(8.5) & 0.5 \\
		3 & log(7) & 0.75 & log(9.5) & 0.75 \\
		4 & log(7) & 1 & log(10)) & 1 \\
		\bottomrule       
	\end{tabular}
\end{table}

\begin{figure}[htbp!]
	\centering
	\includegraphics[scale=0.6]{skew_dist_old_approach2}
	\caption{Distribution of the outcome variable ($Y$) under each skewness scenario used within the alternative data generation approach outlined in Section \ref{app:datagen2}. \label{fig:skew_true_alt}}
\end{figure}




\subsection{Results}
Metrics were used to assess the performance of each confounding-adjustment method, with a key focus on the bias in the point estimate. Results under the alternative data geenration approach are presented in Table \ref{tab:datagen2_results}. The g-comp (med) approach consistently produced the least biased estimates across all skewness scenarios and both confounding settings. Meanwhile the multivariable quantile regression, IPW estimator and weighted quantile regression approaches exhibited large biases. 

Under this data generation approach, there was no guarantee that the propensity score model was correctly specified given the additional selection into exposure groups through the deletion of records with missing outcome; this may contribute to the poor performance of these methods under this data generation process as correct model specification for the IPW-based methods is key as acknowledged elsewhere.\cite{Zhang2012} A couple of variations in the calculation of weights were explored for both the IPW estimator and weighted quantile regression approaches - trimming of propensity scores and incorporating the probability of missingness in calculations of the weights. However, neither of these approaches offered any improvement to the results presented here. In addition, under this data generation approach there still remained no guarantee that the outcome model was correctly specified; a limitation also of the other data generation approach in the main text. 

\clearpage 

\newpage
\begin{table}[htbp]
	\centering
	\caption{Performance of confounding-adjustment methods across confounding and skewness scenarios in the simulation study using the alternative data generation approach. \label{tab:datagen2_results}}	
	\begingroup\fontsize{9}{11.3}\selectfont
	\begin{tabular}{p{2cm} p{1.5cm} r >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm} >{\raggedleft\arraybackslash}p{1.5cm}}
		\toprule
		\textbf{Confounding} & \textbf{Skewness scenario} & \textbf{Method} & \textbf{Bias} & \textbf{Relative bias} (\%) & \textbf{Empirical SE} & \textbf{Model SE }& \textbf{Error SE} (\%) \\ 
		\midrule
		Weak & 1 & Unadjusted & -0.100 & -10.01 & 0.189 & 0.200 & 6.232 \\ 
		& & QR & -0.159 & -15.93 & 0.175 & 0.184 & 5.518  \\ 
		& & IPW estimator & -0.103 & -10.29 & 0.215 & 0.235 & 9.175  \\
		& & Weighted QR & -0.113 & -11.25 & 0.214 & 0.235 & 9.385 \\ 		
		& & G-comp (med) & 0.009 & 0.90 & 0.238 & 0.261 & 10.011 \\ 		 		
		& & G-comp (mean) & -0.138 & -13.80 & 0.205 & 0.228 & 11.216 \\ \hline
		
		Weak & 2 & Unadjusted & -0.150 & -10.00 & 0.399 & 0.399 & 0.037 \\ 
		& & QR & -0.189 & -12.60 & 0.387 & 0.386 & -0.267 \\ 
		& & IPW estimator & -0.148 & -9.86 & 0.461 & 0.477 & 3.427 \\ 		
		& & Weighted QR & -0.171 & -11.41 & 0.459 & 0.476 & 3.577 \\ 		
		& & G-comp (med) & 0.022 & 1.49 & 0.509 & 0.531 & 4.302 \\ 
		& & G-comp (mean) & -0.147 & -9.81 & 0.459 & 0.486 & 5.908 \\ \hline
		
		Weak & 3 & Unadjusted & -0.250 & -10.01 & 0.624 & 0.664 & 6.394  \\ 
		& & QR & -0.282 & -11.28 & 0.624 & 0.633 & 1.480 \\ 
		& & IPW estimator & -0.253 & -10.10 & 0.737 & 0.798 & 8.353 \\ 		
		& & Weighted QR & -0.291 & -11.62 & 0.734 & 0.798 & 8.701 \\ 
		& & G-comp (med) & -0.027 & -1.08 & 0.819 & 0.870 & 6.228 \\ 				
		& & G-comp (mean) & -0.205 & -8.20 & 0.739 & 0.818 & 10.621 \\ \hline
		
		Weak & 4 & Unadjusted & -0.300 & -10.00 & 0.889 & 0.904 & 1.686 \\ 
		& & QR & -0.329 & -10.96 & 0.913 & 0.873 & -4.338 \\ 
		& & IPW estimator & -0.279 & -9.32 & 1.066 & 1.102 & 3.322 \\ 
		& & Weighted QR & -0.331 & -11.05 & 1.063 & 1.096 & 3.050 \\ 
		& & G-comp (med) & 0.009 & 0.31 & 1.133 & 1.211 & 6.896 \\ 						
		& & G-comp (mean) & -0.154 & -5.13 & 1.062 & 1.159 & 9.127 \\ \midrule
		
		Strong & 1 & Unadjusted & -0.200 & -20.01 & 0.186 & 0.191 & 2.823 \\ 
		& & QR & -0.242 & -24.22 & 0.175 & 0.173 & -0.931 \\ 
		& & IPW estimator & -0.170 & -17.03 & 0.207 & 0.222 & 7.190  \\
		& & Weighted QR & -0.180 & -17.98 & 0.207 & 0.222 & 6.998 \\ 		
		& & G-comp (med) & -0.042 & -4.17 & 0.238 & 0.256 & 7.497 \\ 		 		
		& & G-comp (mean) & -0.178 & -17.83 & 0.203 & 0.222 & 9.197 \\ \hline
		
		Strong & 2 & Unadjusted & -0.300 & -20.00 & 0.349 & 0.382 & 9.395  \\ 
		& & QR & -0.356 & -23.72 & 0.342 & 0.353 & 3.285 \\ 
		& & IPW estimator & -0.265 & -17.64 & 0.412 & 0.448 & 8.702 \\ 		
		& & Weighted QR & -0.285 & -19.03 & 0.411 & 0.447 & 8.806 \\ 		
		& & G-comp (med) & -0.030 & -2.01 & 0.491 & 0.518 & 5.613 \\ 
		& & G-comp (mean) & -0.217 & -14.46 & 0.418 & 0.458 & 9.398 \\ \hline
		
		Strong & 3 & Unadjusted & -0.501 & -20.00 & 0.627 & 0.617 & -1.599  \\ 
		& & QR & -0.559 & -22.34 & 0.610 & 0.573 & -5.983 \\ 
		& & IPW estimator & -0.469 & -18.74 & 0.725 & 0.731 & 0.821 \\ 		
		& & Weighted QR & -0.501 & -20.03 & 0.722 & 0.729 & 0.960 \\ 
		& & G-comp (med) & -0.149 & -5.97 & 0.807 & 0.844 & 4.592 \\ 				
		& & G-comp (mean) & -0.392 & -15.68 & 0.696 & 0.764 & 9.800 \\ \hline
		
		Strong & 4 & Unadjusted & -0.600 & -20.00 & 0.813 & 0.842 & 3.545 \\ 
		& & QR & -0.646 & -21.55 & 0.803 & 0.796 & -0.831 \\ 
		& & IPW estimator & -0.585 & -19.50 & 0.935 & 1.013 & 8.361 \\ 
		& & Weighted QR & -0.630 & -21.03 & 0.933 & 1.008 & 8.101 \\ 
		& & G-comp (med) & -0.148 & -4.93 & 1.093 & 1.167 & 6.738 \\ 						
		& & G-comp (mean) & -0.374 & -12.47 & 0.983 & 1.086 & 10.432 \\ 
		\bottomrule
	\end{tabular}
	\endgroup
\end{table}


\clearpage
\newpage





