Deconfounding Multi-Cause Latent Confounders: A Factor-Model Approach to Climate Model Bias Correction

1 Introduction

Global Climate Models (GCMs), such as those developed by the Coupled Model Intercomparison Project phase 6 (CMIP6) [Eyring et al., 2016], are vital tools for predicting future climate changes. These models simulate the physical and chemical processes of the Earth systems-including the atmosphere, oceans, land, and ice-to provide detailed climate forecasts. Despite significant advancements in GCM, their output still exhibits systematic biases. These biases primarily come from uncertainties within the models, simplifications in parameterization processes, and inadequate representations of complex climate phenomena [Mouatadid et al., 2023]. For instance, processes such as cloud formation, precipitation, radiation, and convection are often simplified into parameterization formulas in GCMs, which may not accurately reflect real-world conditions, thus leading to biases. Consequently,

[Lafferty and Sriver, 2023] highlighted that bias corrections were particularly crucial for near-term precipitation projections, especially in cases where observational data are inconsistent with the GCM Output.

To enhance the reliability and accuracy of the GCM Outputs, numerous bias correction techniques, ranging from simple linear scaling to advanced quantile mapping, have been proposed to modify the GCM Outputs and align them more closely with actual observations [Casanueva et al., 2020; Chen et al., 2013; Dowdy, 2020; Enayati et al., 2020; Feigenwinter et al., 2018; Lafon et al., 2013; Maraun, 2016; Mehrotra et al., 2018; Miao et al., 2016; Nahar et al., 2018; Piani et al., 2010; Smitha et al., 2018; Teutschbein and Seibert, 2012; Wu et al., 2022b]. These methods typically rely on historical observation data and statistical approaches to adjust climate model outputs and correct systematic errors [Maraun and Widmann, 2018].

However, these traditional methods have been found to inflate simulated extremes, raising concerns about their use in climate change applications where extremes are significant, such as drought and flooding [Past´en-Zapata et al., 2020]. Another major limitation of most bias correct methods is their assumption that all relevant factors are known and observable, which is unrealistic in practical applications. Climate systems are complex and include many factors that may not be fully observed, such as microclimate effects, regional climate characteristics, and anthropogenic influences. These potential unobserved (confounding) factors dynamically impact time series forecast. However, since they cannot be fully observed, they are frequently overlooked by current bias correction methods. This overlook limits the potential of most existing methods in bias correction.

To deal with the challenge of extreme conditions, the study by [Nivron et al., 2024] incorporates advanced time series forecasting models into the bias correction of extreme weather events, such as heatwaves, by considering the bias correction as a time-indexed regression model with stochastic output. This provides a new perspective: adapting time series forecasting models to bias correction can improve model performance. However, although advanced models have shown significant potential in time series forecasting [Zhou et al., 2021; Wu et al., 2022a; Wu et al., 2023a; Wang et al., 2024; Gao et al., 2024], particularly in climate forecasting [Bi et al.,

2023; Wu et al., 2023b], they typically overlook the existence of unobserved confounding factors, which leads to biased results.

In the field of causal inference with time series, recent studies have started addressing the challenge of unobserved confounders in predicting the potential outcomes of treatments in time series data. Instead of assuming unconfoundedness [Pearl, 2000], these studies operate under a weaker assumption that only multi-cause confounders exist [Bica et al., 2020; Wang and Blei, 2019; Li et al., 2024; Cheng et al., 2023]. A variable is said to be a confounder if it is a common cause of both the treatment and the outcome [Hern´an and Robins, 2020]. For example, in the Directed Acyclic Graph (DAG) shown in Figure 1(a), the confounder X affects both the treatment A and outcome Y . In our bias correction problem, precipitation output is always treated as the outcome, and other related climate variables are treated as treatment variables. The most common method for predicting potential outcomes is to use all treatment variables. However, neglecting unobserved multi-cause confounders Z, such as large-scale atmospheric circulation patterns or oceanic processes, can lead to simultaneous effects on various climate variables, including precipitation. These multi-cause confounders make it possible to apply current research outcomes to bias correction in GCM.

Inspired by the work of [Bica et al., 2020], in this paper, we propose a deconfounding bias correction method. It is important to note that estimating hidden confounders in climate bias correction is much more complex. This increased complexity arises not only from the intricate nature of hidden confounders in climate science but also because these confounders need to be inferred from GCM and observation data. Our contribution can be summarized as the following:

1. We use causality to identify and understand unobserved confounders, allowing us to obtain an unbiased outcome in the presence of multi-cause confounders. By identifying latent confounders through constructing a factor model over time, our approach can capture these unobserved factors, resulting in unbiased outcomes.
2. We develop a two-phase algorithm: Deconfounding and Correction. The Deconfounding Bias Correction (BC) factor model captures the confounders from both Global Climate Model (GCM) output data and observational data, and the correction model uses these confounders as additional information for the bias correction task. Inspired by the TSD (Time Series Deconfounder) [Bica et al., 2020], our method extends to GCM and observation data.

2 ProblemFormulation

As shown in Table 1, let the random variable V represent the set of climate variables, such as humidity, temperature, pressure, precipitation and others. Among these variables, Y denote the outcome of interest, for example, precipitation. In our problem setting, because we will use both historical and current data of these variables, to avoid confusion, Xt represents the historical version of V \ Y . At represents the current version of V \ Y .

These variables are collected from two sources: Global climate models (GCMs) and observations. Moreover, we use XGt and AGt to represent historical and current variables of the GCM data respectively, while XOt and AOt for historical and current variables. The outcome from GCMs and observations are denoted as Y G t and Y O t respectively.

The data for a location, also known as the location trajectory, consist of realizations of the previously described random variables {xGt , aGt , yG t , xOt , aOt , yO t }. Let y(¯at) represent the potential outcome (precipitation), which could be either factual or counterfactual, for each possible treatment course ¯at, where ¯at = (a1, ..., at). Consequently, we have yG(¯at) and yO (¯at). The concept of potential outcomes allows us to consider what the outcome y would be under different treatment scenarios, which is essential for causal inference [Rubin, 1974].

To correct future k steps precipitation bias ∆Y = Y O - Y G which is the difference between observational precipitation yO and GCM output precipitation yG, a common way is to use all data to obtain a regression model [Nivron et al., 2024] which can be represented by Figure 1(a). However, the existence of hidden confounders (e.g. unobserved atmospheric and oceanic circulation) will result in biased result. To address the issue of hidden confounders in time series data, [Bica et al., 2020] have developed a method, assuming the presence of multi-cause hidden confounders Z, as shown in the Figure 1(b).

However, no attempt has been made to adapt their work to bias correction, including bias correction for climate models. In this paper, we extend the work for the bias correction problem as shown in the Figure 1(c) which makes the most of information we have. We assume that hidden variables Z affect both GCM and observations. Learning Z, we would estimate the potential outcome in observations for each location conditional on the location history of covariates ¯XOt = (XO1 , . . . ,XOt ) ∈ X O t , treatments ¯AOt = (AO1 , . . . ,AOt ) ∈ AOt and confounders ¯Zt = (Z1, . . . ,Zt) ∈ Zt:

3 ProposedMethod

The existence of confounders can result in the obtained association among X,A, Y not accurately representing the true relationships, potentially leading to biased results. To solve this, we (1) input the climate model data into the Deconfounding Bias Correction (BC) factor model to obtain the multicause hidden confounder which is essential for bias correction (we call this step 'Deconfounding'). Then (2) we use a hidden multicause confounder as a bias source, combine it with observational data, building a precipitation correction model to help the climate model to have a better estimate of future output (we call this step 'Correction') 1 .

3.1 Deconfounding

To address the challenge of deconfounding time series data with time-varying latent confounders, [Bica et al., 2020] proposed the Time Series Deconfounder (TSD), which extends the deconfounder methodology introduced by [Wang and Blei, 2019] to the time series domain. The fundamental principle of the TSD is that it utilizes factor model to infer substitutes for hidden confounders as treatments.

The goal of this part is to generalize the factor model proposed by [Bica et al., 2020] to a Deconfounding BC factor model for bias correction. This extension aims to extend the Time Series Deconfounder to a complex time series setting of two sources.

Deconfounding BC Factor Model

To extend time series deconfounder [Bica et al., 2020] to bias correction, We propose an enhanced multi-source factor model specially for this task, termed the Deconfounding BC factor model.

For single source data, the unobserved confounder affects X,A, Y from one source, allowing us to infer the sequence of unobserved confounders zt = g(¯ht-1), where ¯ht-1 = {¯at-1, ¯xt-1, ¯zt-1} is the realization of ¯Ht-1. Specifically, factorization can be expressed as follows:

1 Code is avaliable at github: https://github.com/Wentao-Gao/ A-Factor-Model-Approach-to-Climate-Model-Bias-Correction

To extend this model for multi-source data, consider both GCM and Observations as in Figure 2. The unobserved confounder will affect X,A,Y in both sources, allowing us to infer the sequence of unobserved confounders zt = g(¯ht-1) that can be used to render both source treatments conditionally independent, where ¯ht-1 = {¯hOt-1 , ¯hGt-1 } is the realization of ¯Ht-1. Specifically, the factorization for the multisource scenario can be expressed as follows:

It allows us to infer the sequence of latent variables Zt such that, conditioned on Zt and the covariates from both sources XGt and XOt , the treatments become conditionally independent. As in one location, the same confounders Z will affect both GCM and observations. In this case, the Z we learned should render treatments in both GCM and observational data conditionally independent with covariate Xt.

The structure of the factor model depends on causality, which relies on the assumptions listed below.

Assumption 1. Consistency. If ¯ A≥t = ¯a≥t, then the potential outcomes for following the treatment ¯a≥t is the same as the factual outcome Y (¯a≥t) = Y .

Assumption 2. Positivity (Overlap)[Bica et al., 2020]: if P(¯At-1 = ¯at-1,Xt = ¯xt) ̸= 0 then P(At = at | ¯At-1 = ¯at-1,Xt = ¯xt) > 0 for all at.

Assumption 3. Sequential Single Strong Ignorability

Next, we provide a theoretical analysis for the soundness of the learned Zt by introducing the concept of sequential Kallenberg construction [Bica et al., 2020], as follows.

Definition 1. Sequential Kallenberg construction At timestep t, we say that the distribution of assigned causes (AGt1 , . . . ,AGtk ), (AOt1 , . . . ,AOtk ) admits a sequential Kallenberg construction from the random variables Zt = g(¯HOt-1 , ¯HGt-1 ) and XOt ,XGt if there exist measurable functions f O tj : ZOt ×XOt × [0, 1] →AOj and f G tj : ZGt ×XGt × [0, 1] → AGj and random variables U O tj , U G tj ∈ [0, 1], with j = 1, . . . , k such that:

where U O tj , U G tj marginally follow Uniform[0, 1] and jointly satisfy:

for all ¯aO≥t and ¯aG≥t .

We present the following theorem to guarantee that the learned Zt can serve as a substitute for the multi cause hidden confounders.

Theorem 1. The soundness of the learned Zt If at every timestep t, the two distribution of assigned causes (AGt1 , . . . ,AGtk ), (AOt1 , . . . ,AOtk ) admit a Kallenberg construction from Zt = g(¯HOt-1 , ¯HGt-1 ) and XOt ,XGt separately, then the learned Zt can serve as a substitute for the multi cause hidden confounders.

Theorem 1 establishes that the latent variable Zt can be inferred by leveraging the conditional independence structure of the treatments: within each source (GCM or Observation), the assigned treatments {At1, . . . ,Atk} are conditionally independent given Zt and their respective covariates XGt and XOt . This structure, illustrated in Figure 3, provides the necessary identifiability condition for inferring Zt through a shared factor model across both sources. The result means that, at each timestep, for both sources, the variable ¯Xt, ¯Zt,

¯At-1 contain all of the dependencies between the potential outcomes Y(¯a≥t) and the assigned causes At.

Then for both sources we have,

This equation implies that the potential outcome can be represented by E[Y | ¯a≥t, ¯At-1, ¯Xt, ¯Zt]. Therefore, our estimate of the potential outcome is unbiased given the previous treatments, current covariates, and latent variables. This unbiasedness is crucial for making valid inferences about causal effects.

Then, we will show how to implement in practical based on factorization Formula (3) and (4). The joint distribution of all collected T time steps in our enhanced Deconfounding BC Factor Model can be formulated as follows:

The same formulation holds for the GCM source by replacing O with G.

where θO 1:k and θG 1:k are the model parameters. The treatment distributions p(¯aOT ) and p(¯aGT ) are the corresponding marginal distributions.

In practical applications, as shown Figure 3, we need to build a model to obtain Z with predicting treatments as the constraint. To manage time-varying treatments AG,AO , we follow a pragmatic approach by utilizing a recurrent neural network (RNN) with multitask output for implementation.

The recurrent model infers the latent variables Zt based on location history:

where L consists of randomly initialized trainable parameters. Note that in each time step t, the assignments AOt = [AOt1 , . . . ,AOtk ] and AGt = [AGt1 , . . . ,AGtk ] are conditionally independent given Zt, XtO and Zt, XtG separately:

for all j and t, where θG j and θO j are the parameters of the fully connected layers. Using these as multitask output, treatments is conditionally independent given Xt,Zt for both sources, this can be ensured by inferring each treatment AG or AO as different tasks with only XGt ,Zt or XOt ,Zt, if we condition on XGt ,Zt or XOt ,Zt, they are conditionally independent.

Training

The factor model is trained using gradient descent methods on the observational and GCM dataset. The architecture leverages the dependencies between multiple treatments and the location history to infer the latent variables.

3.2 Correction

After obtaining hidden multi-cause confounders Z from both GCM and observational data in one area. In our problem formulation, we aim to correct the bias between Y G and Y O while considering unobserved confounders. A key question is how to utilize this information for bias correction. Advanced time series forecasting models such as iTransformer [Liu et al., 2024] have demonstrated significant potential in capturing temporal information, and integrating these models into our bias correction approach promises substantial improvements. Compared to models like TimeMixer and PatchTST, iTransformer balances complexity and performance.

We propose to use a probability model for bias correction. As shown in Figure 4, we build a prediction model based on the observational data and the multi-cause confounder Z to predict future precipitation.

The probability model employed in our framework is designed to generate predictions of ∆Y , the correction term for the GCM output. Formally, the model learns to approximate the conditional distribution: Pθ(∆Yt+k | AG, AO , Z), where AG and AO are the current covariates from GCM and observations, and Z is the inferred latent confounder.

To achieve high performance in forecasting tasks involving complex temporal dependencies, we choose to use iTransformer [Liu et al., 2024] as the prediction model. iTransformer is a lightweight and effective transformer variant specifically designed for long-sequence time series forecasting, offering a strong balance between expressiveness and efficiency.

After we obtain the ∆Yt+k, we can derive our corrected GCM output with:

4 Experiments And Results

The objectives of the experiments are as follows: 1) Use synthetic data sets to evaluate the correctness of the latent confounder captured by our deconfounding method. 2) Build a correction model to correct precipitation predictions made by the climate model. Implementation code is attached.

4.1 Experiments On Synthetic Data

To assess the effectiveness of our deconfounding method, we conduct controlled experiments on synthetic data, since the true influence of hidden confounders cannot be measured directly in real-world settings [Wang and Blei, 2019]. We generate data from a two-source autoregressive model based on the causal summary graph in Figure 1(c). At each timestep t, we sample time-varying covariates Xsource1 t and Xsource2 t , together with a multi-cause latent confounder Zt. We simulate N = 500 locations over T = 3650 timesteps, with k = 3 covariates and treatments and autoregressive order p = 5. The resulting dataset is split into 80% training, 10% validation, and 10% testing.

Simulated Dataset

Our evaluation focuses on two properties. First, the inferred latent variable Z, together with the observed covariates X, should accurately predict the treatment vector A = [A1, . . . ,Ak]. Concretely, we model p(A | Z,X) = ∏kj=1 p(Aj | Z,X).

Second, the learned Z should match the ground-truth confounders used in simulation. A low mean squared error (MSE) between inferred and true Z indicates successful recovery of the hidden structure. The details for data generation can be found in our GitHub repository.

Figure ?? summarizes our predictive checks. The Deconfounding BC method yields an MSE of 0.06157 when predicting the treatment vector, demonstrating strong alignment between model predictions and the simulated treatments. Comparing the inferred Z against the true latent confounders gives an MSE of 0.00187, confirming that our factor-model approach accurately captures the hidden variables. Together, these results validate that the treatments become conditionally independent given the inferred confounders and covariates, and that the learned latent representation faithfully recovers the ground-truth confounding signal.

4.2 Case Study: South Australia

In a bias correction (BC) process, data from climate models are compared to actual observations (or their proxies, like reanalysis products) to adjust for biases. Since global climate models (GCMs) usually have a lower resolution compared to observations or reanalysis reference data, BC often involves downscaling the resolution of GCMs. We adopt this methodology as well. For the climate model data, we selected the 15 initial condition runs from the Institut Pierre-Simon Laplace (IPSL) climate model as part of the sixth Coupled Model Intercomparison Project (CMIP6) historical experiment. The data, which is available on a monthly basis, includes climate variables such as tmax (maximum temperature 2 meters above the surface) and prate (precipitation rate), with our bias correction efforts centered on prate.

The IPSL 2 model is run at a 250km nominal resolution and is not re-gridded. We selected the closest geographical point to South Australia for the case study, covering the period from 1948 to 2014.

2 IPSL data portal: https://aims2.llnl.gov/search/cmip6/

For observational reference data, NCEP-NCAR Reanalysis 1 3 , provided by the National Oceanic and Atmospheric Administration (NOAA), is utilized. The dataset encompasses the same variables as the IPSL model (e.g., tmax & prate). NCEP-NCAR reanalysis data is available at a monthly frequency, covering the period from 1948 to 2014, with a global resolution of 2.5 degrees in both the latitudinal and longitudinal directions, and has not been re-gridded. The nearest geographical point to South Australia was selected.

This research focuses on bias correction of the precipitation rate to improve the accuracy of precipitation predictions.

Data Preprocessing

We extracted and converted climate data from 1948 to 2014 for the South Australia region. The data, originally in NetCDF format, were transformed into CSV files for further processing. Subsequently, we split the data into a training period (1948-1992) and a testing period (1993-2014). The IPSL dataset contains 30 initial conditions; to generalize our model, we selected 15 initial conditions for the experiment. In this paper, we present results from three of these initial conditions(r5i1p1f1, r6i1p1f1 and r7i1p1f1) along with their average results.

Experiment Settings

To evaluate our proposed method, we compare our method with several bias correction baseline methods, including linear scaling [Teutschbein and Seibert, 2012], variance scaling [Teutschbein and Seibert, 2012], quantile mapping [Cannon et al., 2015], quantile delta mapping [Tong et al., 2021] and Temporal BC [Nivron et al., 2024] methods. The metric we are using to compare the performance is the Mean Squared Error (MSE) and the Mean Absolute Error (MAE).

4.3 Results

Table 2 illustrates the three-step-ahead predictions using our deconfounding bias correction method on the dataset with different GCM initial settings. After applying the deconfounding step and augmenting the dataset with substitutes for the hidden confounders, we evaluated the performance of various bias correction methods using MSE and MAE. The results in Table 2 support the effectiveness of the Deconfounding BC method, which consistently achieves the lowest MSE and MAE across all experimental conditions. This indicates its superior performance in reducing prediction errors and aligning closely with the observed data. Figure 6 further visualizes the daily precipitation trends from IPSL, observations, and various correction methods, clearly highlighting the superior alignment achieved by Deconfounding BC.

3 NCEP-NCAR data portal: https://psl.noaa.gov/data/gridded/data. ncep.reanalysis.html

The QQ plots (Figure 7) reveal that among the evaluated bias correction methods, Deconfounding BC shows the closest alignment with the red dashed line (y=x) across the entire range of observed precipitation values. This indicates that Deconfounding BC is the best in accurately capturing the observed precipitation levels, outperforming other methods in reducing biases and providing reliable predictions.

To demonstrate that the latent confounder learned by our Deconfounding BC factor model contains essential information, we did ablation study to compare our correction model's results with and without the latent confounder, as shown in Table 3. The performance improved with the inclusion of the hidden confounder.

5 Conclusion

In this paper, we propose the deconfounding bias correction method for multi-cause confounders. By integrating climate bias correction techniques with causality-based time series deconfounding, our approach provides a novel perspective for future studies, emphasizing the importance of not assuming all variables are observed.

Acknowledgments

This work was supported by the ARC Discovery Project DP230101122 and the University of South Australia Research Training Program (RTP) Scholarship. We gratefully acknowledge the continued support from the CSIRO Environment Research Unit and the Data61 Business Unit.