A short user’s guide for QUALYPSO

Univ. Grenoble Alpes, INRAE, IGE - Institut des Géosciences de l’Environnement
guillaume.evin@inrae.fr

Univ. Grenoble Alpes, INRAE, IGE - Institut des Géosciences de l’Environnement
benoit.hingray@univ-grenoble-alpes.fr

Abstract

This R package aims at characterizing the overall variability in an ensemble of climate projections (or an ensemble of projections of indicators impacted by climate change) given that different sources of uncertainty are present. It provides estimates of climate change responses of all simulation chains and of all uncertainty variables. It additionally propagates uncertainty due to missing information in the estimates.

1 Introduction

The package provides an easy-to-use statistical tool that performs analysis of variance (ANOVA) of ensembles of climate projections. In particular, different regression models are proposed to estimate the climate response of a given chain and filter out fluctuations resulting from internal variability. The corresponding “Quasi-Ergodic Analysis of Climate Projections Using Data Augmentation” approach (QUALYPSO) was first designed to address this challenge and the treatment of incomplete ensemble using data augmentation techniques (Evin et al. 2019). It has been applied to different ensembles of projections: seasonal precipitation and temperature in Evin et al. (2019), and Evin et al. (2021); hydrological variables in Evin et al. (2026), Lemaitre-Basset et al. (2021), and Aitken et al. (2023); photovoltaic energy in Bichet et al. (2019).

Most of the studies about climate change impacts rely on future projections of some indicators describing the evolution of the different components of the earth system (e.g. atmosphere, oceans, cryosphere). Regularly, Coupled Model Intercomparison Projects (CMIPs) provide large ensembles of projections simulated from different global climate models (GCMs) and regional climate models (RCMs) driven by a small selection of emission scenarios (RCP, SSP) in order to improve knowledge on climate change and its impacts. Impact models (e.g. hydrological models), driven by climate projections, provide valuable information about climate change impacts in many sectors (agriculture, water resources, energy production, infrastructures, transport, economy, etc.).

Multi-model projections of climate change impacts can be used to assess mean future trends, but characterising the overall variability of the ensemble and identify the main sources of uncertainty is arguably as important. Using an ANOVA implemented in QUALYPSO, Evin et al. (2026) provide the contribution of each factor to the total variance of the ensemble. In this study, up to five sources of uncertainty are considered:

• uncertainty on future emissions (i.e. scenario uncertainty), assessed using the simulations based on different emission scenarios,
• global and regional climate model uncertainty, arising from the diversity of representation of the climate system by the different GCMs and RCMs,
• hydrological model uncertainty, as a result of different representations of the hydrological system in four hydrological models,
• climate internal variability (IV), which represent the inter-annual to decadal fluctations around the long-term response of climate change as a result of the chaotic nature of the climate system (Lehner and Deser 2023).

Depending on the application, additional sources uncertainties can be considered. For example, Lemaitre-Basset et al. (2021) considers different calibration schemes to evaluate the potential contribution of the hydrological model calibration to the total uncertainty.

As a result of IV, climate projections include large fluctuations. The magnitude of IV depends on 1. the temporal (e.g. interannual or decadal variability) and/or spatial aggregation and 2. the indicators under consideration: IV is typically small for temperature (or indicators derived from temperature) and large for hydroclimate indicators (e.g. precipitation, flood), see Hawkins and Sutton (2009), Evin et al. (2026). Due to IV, the estimation of the climate response, i.e. long-term evolutions of climate or hydroclimatic variable, can be challenging. QUALYPSO embeds different statistical models to estimate these climate responses. In addition, QUALYPSO also addresses:

• incomplete ensembles: Some combinations of models and scenarios can be missing. This is usually the case when climate projections are simulated by RCMs.
• unbalanced ensembles: Some models of a given model category (GCM, RCM) can be often over- or under-represented in large ensembles, see, e.g., the large EUROCORDEX MME assessed in Evin et al. (2021).

2 Concepts

To characterize how a given indicator is projected to change, QUALYPSO focuses on the climate change response (CCR) projected for the considered indicator. The CCR projected by a given simulation chain defines, for any future time, the absolute or relative difference of the long-term climate response for this future time compared to the reference time. QUALYPSO aims to estimate, for any future time, the CCR of each modeling chain, the mean CCR of the MME, the dispersion between the CCRs of individual chains, the sources of uncertainty explaining this dispersion, and the additional potential dispersion in future realizations of climate indicators that can result from IV.

The QUALYPSO package contains a set of average winter temperature (DJF) climate projections for a region covering most of continental Europe (the SREX Central-Eastern Europe “CEU” region). The dataset Y_DJFTas contains projections for the period 1971-2099 (129 years) for 20 combinations of GCMs (CMIP5) and RCMs models, obtained with the RCP8.5 scenario, and is illustrated in Fig. 2.1. The main steps applied in QUALYPSO will be illustrated using this dataset and the name of the climate models (type of GCM and RCM) used to simulate the temperatures and provided by the data.frame scen_DJFTas (see Table 2.1), so-called ``factors” in ANOVA analyses.

Table 2.1: Two factors characterise the climate projections for the ensemble of average winter temperatures: the GCM and the RCM.

GCM	RCM
CNRM-CM5	RACMO22E
CNRM-CM5	CCLM4-8-17
CNRM-CM5	HIRHAM5
CNRM-CM5	RCA4
CNRM-CM5	REMO
EC-EARTH	RACMO22E

# plot of the 20 projections
vecGCM = unique(scen_DJFTas$GCM)
X_time_vec=1971:2099
plot(-1,-1,xlim=range(X_time_vec),ylim=range(Y_DJFTas),xlab="Years",ylab="DJF Temperature (°C)")
for(i in 1:nrow(Y_DJFTas)){
  lines(X_time_vec,Y_DJFTas[i,],col=match(scen_DJFTas$GCM[i], vecGCM))
}
legend("bottomright",legend = vecGCM, col=1:length(vecGCM),lty=1)

Figure 2.1: Ensemble of average winter temperature (DJF) climate projections for the SREX Central-Eastern Europe CEU region, obtained with 20 combinations of four GCMs and five RCMs. The different colors correspond to the different GCMs.

2.1 Main assumptions

Following Hawkins and Sutton (2009) and Hingray and Saïd (2014), QUALYPSO is based on two main hypotheses expressed as follows:

• We consider that the projection of a simulation chain \(i\) combines two components: the climate response and internal variability. For a future time \(t\), the climate response corresponds to the long-term trend of the simulated projection. We hypothesize that it has a temporal variation that is inherently gradual and smooth. The high- to mid-frequency fluctuations around the climate responses are assumed to result solely from IV.
• The CCR of a simulation chain \(i\) can be expressed as a linear sum of the ensemble mean and the main effects of the different components of the chain (e.g. main effects of the scenario \(s\), GCM \(g\), …).

The first assumption implies that it is reasonable to consider a trend model to estimate the climate response of a chain and, in turn, fluctuations due to IV (deviations from the climate response). Hingray et al. (2019) have shown that this assumption allows for providing more precise estimates for all uncertainty components than those obtained with standard time-slice approaches, i.e. using an ANOVA applied to different future time windows, see, e.g. Yip et al. (2011) and many following works.

The additive decomposition model of the second assumption can be derived for any future time using a fixed-effects ANOVA model. In the case of a complete MME (no missing combinations), the main effect of a given scenario or model is easily interpreted as the mean difference between 1) the CCR of all the projections using this scenario or model and 2) the mean CCR of the ensemble.

2.2 Estimation of the climate response

In this vignette, we adopt the notations used in Evin et al. (2026). In QUALYPSO, the climate response \(CR_i(t)\) of each modelling chain \(i\) is estimated by applying a trend model to the corresponding projection \(Y_i(t)\). Several regression models are available (smoothing splines, polynomial regressions, local linear regressions) and corresponding parameters can be modified to provided more or less flexibility. This degree of smoothing is arbitrary, but a moderate flexibility limits the influence of IV. Scherrer et al. (2024) compare several approaches and advise to apply local linear regressions (loess) with a large span (42 years).

QUALYPSO provides different options for the trend model. Smoothing splines can be fitted using smooth.spline, where the argument df indicates the desired equivalent number of degrees of freedom. poly fits a polynomial model and the parameter degree is an integer indicating the degree of the polynomial. lowess fits locally-weighted polynomial regressions and the parameter f in (0,1) is the proportion of points which influence the smooth at each value.

Figure 2.2 illustrates these different options, for one projection of temperatures. For smoothing splines, small df values almost fit a linear curve (see plain red curve, ‘SS: df=3’), whereas larger df values provide more flexibility but might overfit decadal oscillations, which should attributed to internal variability. Similarly, a fourth-degree polynomial as in Hawkins and Sutton (2009) (see dotted blue curve, ‘poly: degree=4’) or a small span parameter for lowess (green plain curve, ‘lowess: f=0.1’) lead to overfitted climate trends.

# illustration for the first simulation
iSimu=1
X_time_vec=1971:2099
dfCR = data.frame(x=X_time_vec,y=Y_DJFTas[iSimu,])
plot(-1,-1,xlim=range(dfCR$x),ylim=range(dfCR$y),xlab="Years",ylab="DJF Temperature (°C)")
lines(dfCR$x,dfCR$y,col="black")
lableg = c()
colleg = c()
ltyleg = c()

# fit smoothing splines, df can take values between 1 and the length of the simulations (here, 129)
vec_df_SS = c(3,6,15)
for(idf in 1:length(vec_df_SS)){
  ySS = stats::smooth.spline(x = dfCR$x,y = dfCR$y, df=vec_df_SS[idf])$y
  lines(X_time_vec,ySS,col="red",lty=idf,lwd=2)
}
lableg = c(lableg,paste0("SS: df=",vec_df_SS))
colleg = c(colleg,rep("red",length(vec_df_SS)))
ltyleg = c(ltyleg,1:length(vec_df_SS))

# fit polynomial curves of order 1,2, and 3
vec_poly = c(1,2,4)
for(iorder in 1:length(vec_poly)){
  yPOLY = stats::lm(y ~ stats::poly(x, vec_poly[iorder]), data = dfCR)
  lines(X_time_vec,yPOLY$fitted.values,col="blue",lty=iorder,lwd=2)
}
lableg = c(lableg,paste0("poly: degree=",vec_poly))
colleg = c(colleg,rep("blue",length(vec_poly)))
ltyleg = c(ltyleg,1:length(vec_poly))

# fit LOWESS with different degrees of smoothness
vec_lowess = c(0.2,0.4,0.6)
for(ilowess in 1:length(vec_lowess)){
  yLOWESS <- stats::lowess(x=dfCR$x,y=dfCR$y, f=vec_lowess[ilowess], iter=0, delta=0)
  lines(yLOWESS$x,yLOWESS$y,col="green",lty=ilowess,lwd=2)
}
lableg = c(lableg,paste0("lowess: f=",vec_lowess))
colleg = c(colleg,rep("green",length(vec_lowess)))
ltyleg = c(ltyleg,1:length(vec_lowess))

legend("bottomright",legend = lableg,lty=ltyleg,col=colleg)

Figure 2.2: Climate responses obtained with different trend models (smoothing splines, ‘SS’; polynomial models, ‘poly’, and lowess) and different smoothing parameters.

2.3 Estimation of the climate change response and of internal variability

The climate change response (CCR) of each simulation \(i\), denoted by \(CCR_i(t)\), is obtained by computing the differences in the climate response \(CR_i(t)\) between future times and the reference time \(t=c\). Absolute changes \(CCR_i(t)=CR_i(t) - CR_i(c)\) or relative changes \(CCR(t)= CR_i(t) / CR_i(c) - 1\) can be considered. Corresponding deviations \(D_i(t)\) between the raw projection \(Y_i(t)\) and the climate response \(CR_i(t)\) are obtained as \(D_i(t)=Y_i(t) - CR_i(t)\) and \(D_i(t)= [ Y_i(t) - CR_i(t) ] / CR_i(c)\), respectively. By construction, the ensemble of CCRs has a null variability for the reference time \(c\) and this variability generally increases for further time horizons.

Internal variability \(IV(t)\) is then estimated by the time variance of these deviations \(Var(D_i(t))\) by fitting the following gamlss model: \(D_i(t) \sim N(0, \sigma_t)\), where \(\log(\sigma_t)\) is a linear function of time. The log-link is the default option in gamlss and avoids obtaining negative values for \(\sigma_t\).

2.4 ANOVA model: Estimation of the ensemble mean and the main effect of each scenario/model

To partition uncertainty sources, QUALYPSO assumes an additive model for the climate change response \(CCR_i(t)\), i.e. a sum of the main effects corresponding to the different factors of dispersion (e.g. GCM, RCM) associated to this chain. In our example, the simulation \(i=1\) has been obtained with the GCM CNRM-CM5 and the RCM RACMO22E: \[\begin{equation} CCR_1(t) = M(t) + eG(CNRM, t) + eR(RACMO, t) + \epsilon_1(t), \tag{2.1} \end{equation}\] where \(M(t)\) is the ensemble mean for a future period \(t\), and \(eG(CNRM, t)\) and \(eR(RACMO, t)\) are, respectively, the main effect of the GCM CNRM-CM5 and the main effect of the RCM RACMO22E for the time \(t\), and where \(\epsilon_1(t)\) is the residual term, i.e. the part of the climate change response that is not explained by the ANOVA model (2.1), i.e. the sum of the individual main effects. The residuals of the ANOVA model arise from the imperfect estimation of the climate responses and the limitations of the additive assumption. For example, as a result of the additive assumption, the residuals can potentially include interactions between scenarios and models and/or between models of different categories, see, e.g. Hingray and Saïd (2014); Evin et al. (2021).

For each future time horizon \(t\), the estimation of the grand mean and the main effects is carried out simultaneously using the function lm of the library stats (as in Evin et al. 2026) or a Bayesian model (Evin et al. 2019). When lm is applied, the least squares of the residuals \(\epsilon_i(t), i=1,\dots,N\), where \(N\) is the number of simulations, are minimized. The Bayesian approach gives similar results but provides a more comprehensive description of the uncertainty associated to the parameter estimates. However, it is computationally expensive and lm is the default (preferred) option. By construction, for both approaches, the sum of the main effects of a given factor is zero. Both approaches also provide unbiased estimates even if the ensemble is unbalanced (or incomplete). The uncertainty corresponding to the estimation of the main effects (standard deviations and confidence intervals) is directly related to the number of simulations available for one particular factor.

For a future time horizon \(t\), the uncertainties associated with each factor (GCM, RCM, ..) correspond to the variance of the corresponding main effects. For example, GCM uncertainty is estimated by the variance of the main effects \(eG(g,t)\), \(g=1..N_g\), estimated for the \(N_g\) GCMs. Finally, the residual variability \(RV(t)\) is the variance of the residuals \(\epsilon_i(t)\) of the ANOVA model.

The total uncertainty variance \(CCRV(t)\) of the ensemble is obtained as the sum of the variance of each factor, and residual variability. In our example: \[\begin{equation} CCRV(t) = Var(eG(g,t)) + Var(eR(r,t)) + RV(t). \end{equation}\] These different estimates can be used to estimate the fraction of the total uncertainty resulting from each source of uncertainty: \[\begin{equation} FV_x(t) = Var_x(t)/CCRV(t), \end{equation}\] where \(FV_x(t)\) is the fraction of total uncertainty variance \(CCRV(t)\) explained by the source of uncertainty associated with \(x\) (in the example \(x\) = GCM, RCM, or RV).

3 The package

The QUALYPSO package has one main function QUALYPSO which applies the different steps described above and has two mandatory arguments: Y and scen:

• Y is an ensemble of projections, i.e. a matrix nS x nY where nS projections are provided for nY years or time steps,
• scen is a data.frame that provides nFactor factors for each projection (e.g. GCM, RCM).

3.1 Optional arguments

Additional features are provided if needed:

• X: Projections are often assessed for different future time horizons (e.g. 2050, 2070) or given warming levels (e.g. +2\(^\circ\)C, +3\(^\circ\)C). In QUALYPSO, these horizons correspond to the predictors during the estimation of the climate responses and the argument X can be used to provide specific predictors corresponding to the projections, e.g. time or warming level. It can be a vector if the predictor is the same for all projections (e.g. X=2001:2100) or a matrix of the same size as Y if these predictors are different for the different projections (see case study where X are global temperatures). By default, a vector 1:nY is created.

• Xfut: The argument Xfut can be used to define the nF values of the predictor (time horizon or warming level) over which the ANOVA will be applied and for which QUALYPSO outputs will be produced. It must be a vector of values within the range of values of X to avoid extrapolation of the fitted climate responses. By default, it corresponds to X if X is a vector, 1:nY if X is NULL or a vector of 10 values equally spaced between the minimum and maximum values of X if X is a matrix. The first value defines the reference time or warming level c used to obtain \(CR_i(c)\) when climate change responses are computed.

• listOption: named list where each item provides a different option. In particular:

  • typeClimateResponse is the type of climate response fitted to the projections. It can be either smooth.spline where the parameter df in (1,nY] is the equivalent number of degrees of freedom (value by default), poly, a polynomial model fitted with lm where the parameter degree is an integer indicating the degree of the polynomial, lowess where the parameter f in (0,1) is the proportion of points which influence the smooth at each value, or tweedie, a Generalized Linear Model with the Tweedie distribution (Jørgensen 1987) where the parameter var.power is the index of power variance function (positive value). The tweedie model is the only model that handles zeros in the climate projections, typically with var.power in (1,2), see case study with projections of snow maxima hereafter.
  • parClimateResponse: parameter of the model applied for the extraction of the climate response (i.e. df, degree, f, or var.power). Default values are df=4, degree=1, f=0.4 or var.power=1.5 for the smoothing spline, polynomial, lowess, or tweedie models, respectively.
  • typeChangeVariable: type of climate change response \(CCR(t)\): abs (absolute, value by default) or rel (relative). It can be interesting to assess relative changes when relative increases/decreases are more informative than absolute changes (e.g. for seasonal precipitation).

3.2 Outputs

QUALYPSO returns a list containing different outputs given as time series corresponding to Xfut (except CHANGEVAR_t and D_t):

• CR_tfut: Fitted climate responses [nS x nF],
• CCR_tfut: Climate change responses [nS x nF],
• D_t: Deviation from the climate change responses [nS x nY],
• CHANGEVAR_t: Change variable from the projections [nS x nY],
• M_tfut: List of estimates for the grand mean (mean, standard deviation, and credible/confidence intervals).
• E_tfut: List of estimates for the main effects \(eX(x,t)\) (mean, standard deviation, and credible/confidence intervals).
• VARE_tfut: For each factor, provide the corresponding uncertainty \(Var_x(t)\), i.e. the variance of the corresponding main effects (i.e. variability between the different RCMs, GCMs,..).
• CONTRIB_VARE_tfut: Contribution of each individual effect to its component (as a percentage), e.g. what is the contribution of GCM1 to the GCM uncertainty.
• RV_tfut: List of estimates for the variance of the residual errors (mean, standard deviation, and credible/confidence intervals).
• IV_tfut: Internal variability \(IV(t)\).
• CCRV_tfut: total variability \(CCRV(t)\).
• FV_CCRV_tfut: decomposition of the total variability for each component.
• RESERR: \(\epsilon_i(t)\), the differences between the climate change responses and the additive ANOVA formula (sum of the grand mean and the main effects).

In addition to CCRV_tfut and FV_CCRV_tfut, we also provide the corresponding estimates including the interval variability, i.e. CCRV_IV_tfut and FV_CCRV_IV_tfut, where CCRV_IV(t) is the sum of the total uncertainty CCRV(t) and the internal variability IV(t).

3.3 Plotting facilities

Different plotting functions are provided to diagnose the outputs of the function QUALYPSO:

• plotQUALYPSO_CCR plots the ensemble of climate change responses \(CCR_i(t)\),
• plotQUALYPSO_CR plots one simulation and the corresponding fitted climate response \(CR_i(t)\),
• plotQUALYPSO_effect plots the estimated main effects for one factor \(Ex(t)\) and the corresponding estimation (sampling) uncertainties,
• plotQUALYPSO_M plots the estimates of the grand mean \(M(t)\) and its uncertainty,
• plotQUALYPSO_IV plots the deviations \(D_i(t)\) for all projections and an interval based on the estimates of the internal variability,
• plotQUALYPSO_FV_CCRV_IV plots the fractional variance of \(CCRV(t) + IV(t)\).

4 Case studies

4.1 Example 1: Synthetic scenarios

This first example aims at illustrating the application of QUALYPSO when the true values of the main effects are known. It relies on a synthetic ensemble of 5 projections for 3 fictive GCMs and 2 fictive RCMs. The projections consider different combinations of these GCMs/RCMs for which the main effects are known. One combination is missing (GCM3/RCM2). Simple linear evolutions of time are considered in this toy model, for 21 future pseudo-years. Gaussian noise with a constant standard deviations are generated around these main effects and the grand mean is zero.

# create nS=3 fictive climate scenarios with 2 GCMs and 2 RCMs, for a period of nY=20 years
n=20
t=0:n/n

# GCM effects (sums to 0 for each t)
effGCM1 = t*2
effGCM2 = t*1
effGCM3 = t*-3

# RCM effects (sums to 0 for each t)
effRCM1 = t*1
effRCM2 = t*-1

# These climate scenarios are a sum of effects and a random gaussian noise
scenGCM1RCM1 = effGCM1 + effRCM1 + rnorm(n=n+1,sd=0.5)
scenGCM1RCM2 = effGCM1 + effRCM2 + rnorm(n=n+1,sd=0.5)
scenGCM2RCM1 = effGCM2 + effRCM1 + rnorm(n=n+1,sd=0.5)
scenGCM2RCM2 = effGCM2 + effRCM2 + rnorm(n=n+1,sd=0.5)
scenGCM3RCM1 = effGCM3 + effRCM1 + rnorm(n=n+1,sd=0.5)
Y_synth = rbind(scenGCM1RCM1,scenGCM1RCM2,scenGCM2RCM1,scenGCM2RCM2,scenGCM3RCM1)

# Here, scen_Synth indicates that the first scenario is obtained with the combination of the
# GCM "GCM1" and RCM "RCM1", the second scenario is obtained with the combination of
# the GCM "GCM1" and RCM "RCM2", etc.
scen_Synth = data.frame(GCM=c('GCM1','GCM1','GCM2','GCM2','GCM3'),
RCM=c('RCM1','RCM2','RCM1','RCM2','RCM1'))

We can run QUALYPSO for these projections available for 20 pseudo-years between 2000 to 2100. By construction, the long-term evolution of each synthetic climate projection is linear. Here, a linear trend model is adequate to estimate the climate response of each projection (option typeClimateResponse="poly" with degree=1 by default).

# run QUALYPSO
X = seq(from=2000,to=2100,by=5)
Q_synth = QUALYPSO(Y=Y_synth, scen=scen_Synth, X=X,listOption = list(typeClimateResponse="poly"))

We can check that climate responses are well fitted and compare them to the “true” climate response used to generate the projections (red line):

plotQUALYPSO_CR(Q_synth,i=1)
lines(X, effGCM1 + effRCM1,lty=2)
legend("topleft",inset=c(0,0.1), legend = c("True response"), lty = 2, bty = "n")

Figure 4.1: Raw projection (plain curve), fitted climate response (thick line) and true climate response (dotted line)

For GCM and RCM effects, we can superimpose the estimated effects and the true known main effects (dotted lines). Note that they are adequately estimated even if the projection for the combination GCM3/RCM2 is missing.

plotQUALYPSO_effect(Q_synth,nameEff="GCM",xlab="Years")
lines(X, effGCM1,lty=2,lwd=2,col="black")
lines(X, effGCM2,lty=2,lwd=2,col="red")
lines(X, effGCM3,lty=2,lwd=2,col="green")
legend("topleft",inset=c(0,0.1), legend = c("Truth"), lty = 2, bty = "n")

Figure 4.2: Median estimates of the main effects (plain lines) and corresponding 90% credible intervals, and true main effects (dotted lines)

4.2 Example 2: Uncertainty Analysis for different time horizons

The second example aims at illustrating QUALYPSO with real projections using the set of average winter temperature (DJF) climate projections described above. The dataset Y_DJFTas contains projections for the period 1971-2099 (129 years) for 20 combinations of GCMs (CMIP5) and RCMs models, obtained with the RCP8.5 scenario. We apply QUALYPSO with default options: the climate response of each projection is estimated using smoothing splines with df=4 and the climate change response of each projection \(i\) is defined by absolute differences \(CR_i(t)-CR_i(c)\).

# call QUALYPSO
Q_time = QUALYPSO(Y=Y_DJFTas,scen=scen_DJFTas,X=1971:2099)

# diagnostic of the climate responses
for(iS in 1:2) plotQUALYPSO_CR(QUALYPSOOUT = Q_time, iS = iS)

Figure 4.3: Fitted climate response for two projections

Figure 4.4: Fitted climate response for two projections

In this example, we also illustrate how the assumption about the evolution of internal variability can be checked. We remind that the logarithm of the standard deviation of \(D_i(t)\) is a linear function of the predictor (e.g. time or warming level). The function plotQUALYPSO_IV provides a diagnostic plot where the deviations \(D_i(t)\) are plotted as a function of the predictor (here the year), along with an interval [-1.645\(\sqrt{IV(t)}\) , +1.645\(\sqrt{IV(t)}\)] that represent the model for \(D_i(t)\) and should contain about 90% of the points according to the gaussian assumption. Figure 4.5 shows the output of this plot for this example. We see that internal variability decreases as a function of time. In this case, the interval contains most of the points and assumption about the evolution of internal variability seems reasonable.

# diagnostic of the internal variability
plotQUALYPSO_IV(QUALYPSOOUT = Q_time)

Figure 4.5: Diagnostic plot for internal variability. Points represent all \(D_i(t)\) estimates as a function of \(X\). The red curves correspond to [-1.645\(\sqrt{IV(t)}\) , +1.645\(\sqrt{IV(t)}\)] and should contain about 90% of the points.

4.3 Example 3: Uncertainty Analysis for different Global Warming Levels

The third example performs the ANOVA decomposition with the same ensemble of winter temperatures but using warming level as the predictor. X_DJFTas_WL contains smoothed warming levels at the planetary scale simulated by four different CMIP5 GCMs for the same period 1971-2099. All GCMs exceed +4 \(^\circ\)C at the end of the century (see Figure 4.7). HadGEM2-ES is warmer than the three other GCMs.

vecGCM = unique(scen_DJFTas$GCM)
X_time_vec = 1971:2099
plot(-1,-1,xlim=range(X_time_vec),ylim=range(X_DJFTas_WL),xlab="Years",ylab="Global warming levels (°C)")
for(i in 1:nrow(X_DJFTas_WL)){
  lines(X_time_vec,X_DJFTas_WL[i,],col=match(scen_DJFTas$GCM[i], vecGCM))
}
legend("bottomright",legend = vecGCM, col=1:length(vecGCM),lty=1)

Figure 4.6: Global warming levels from the four GCMs.

We can apply QUALYPSO as a function of these global warming levels simply by replacing the argument X. Here, we ask to produce the outputs for different warming levels between +1 \(^\circ\)C and +3 \(^\circ\)C.

# call QUALYPSO
Q_wl = QUALYPSO(Y=Y_DJFTas,scen=scen_DJFTas,X=X_DJFTas_WL,Xfut=seq(from=1,to=3,by=0.1))

# variance decomposition
plotQUALYPSO_effect(QUALYPSOOUT = Q_wl,nameEff="GCM",xlab="lobal warming levels (°C)")

Figure 4.7: Fraction of total variance for absolute changes in mean winter temperature as a function of global warming levels.

4.4 Example 4: Projections with many zero values and Tweedie functions

In many studies about climate change impacts, indicators can be bounded by physical thresholds. For example, many indicators are greater than or equal to zero (e.g. precipitation, low flows, snow intensities). When many zero values are present in the projections, extracting the climate responses using standard regression methods (linear models, smoothing splines, lowess) is not adapted because there is a probability mass in zero and residuals are not Gaussian. In addition, in this case, non-zero values usually follow skewed distributions, e.g. gamma distributions.

The last example illustrates this type of application with a dataset containing annual maxima of snow water equivalent (SWE) projections for an area located at a low mean elevation (0 m) and with mild and rainy winters. Due to the rare snow events, these projections contain many zero values and distributions of non-zero values are skewed (Evin et al. 2025). This ensemble is comprised of 18 projections obtained with nine different combinations of GCMs and RCMs and two emission scenarios (RCP4.5 and RCP8.5). Obviously, at larger warming levels, snow events become even rarer and the chance of having zero values increases.

plot(-1,-1,xlim=range(X_SWE_WL),ylim=range(Y_SWE),
xlab="Global warming levels (°C)",ylab = "SWE maxima (mm w.e.)")
for(i in 1:nrow(Y_SWE)){
  lines(X_SWE_WL[i,],Y_SWE[i,],col=i)
}

For this kind of data, extracting the climate responses using standard regression methods (linear models, smoothing splines, lowess) are not adapted because residuals are not Gaussian. They provide poor fits when there is a probability mass in zero. In our example, they provide negative SWE values when the warming levels exceed 4\(^\circ\)C.

An alternative model is provided by the tweedie function of the package statmod, which provides access to a wide range of generalized linear model (GLM) response distributions, including response distributions with mass at zero and continuous distributions for positive values.

# illustration for the first simulation
iSimu=1

# df can take values between 1 and the length of the simulations (here, 129)
vec_df = c(3,6,15)

# take one projection
dfpred = data.frame(x=X_SWE_WL[4,],y=Y_SWE[4,])

plot(-1,-1,xlim=range(dfpred$x),ylim=range(dfpred$y),
xlab="Global warming levels (°C)",ylab = "SWE maxima (mm w.e.)")

# observations
points(dfpred$x,dfpred$y,pch=20,cex=1.3)

# smoothing splines
smooth.spline.out<- stats::smooth.spline(x=dfpred$x,y=dfpred$y,df=4)
lines(dfpred$x,smooth.spline.out$y,col="red",lwd=2)

# polynomial function
lm_out = stats::lm(y ~ stats::poly(x, 1), data = dfpred)
lines(dfpred$x,lm_out$fitted.values,col="green",lwd=2)

# lowess
lowess_out <- lowess(x=dfpred$x,y=dfpred$y, f=0.6, iter=0, delta=0)
lines(dfpred$x,lowess_out$y,col="blue",lwd=2)

# tweedie
glm_out <- stats::glm(y ~ x, 
        family = statmod::tweedie(var.power = 1.5, link.power = 0), 
        data = dfpred)
lines(dfpred$x,glm_out$fitted.values,col="orange",lwd=2)
        
legend("topright",legend = c("smooth.spline","poly","lowess","tweedie"),
lwd=2,col=c("red","green","blue","orange"),title="climate response model")

For this ensemble of projections, climate responses show a clear decrease at larger warming levels, up to -3 mm w.e. at +3\(^\circ\)C. These changes are small in comparison of the large inter-annual variability, SWE maxima exceeding 10 mm w.e. every few years. In this ensemble, we can also notice that residual variability is larger than the uncertainty due to RCMs, GCMs, and RCP emission scenarios, indicating that the uncertainty of climate responses is not well characterized by the additive sum of RCM, GCM, and RCP main effects.

# list of options
listOption = list(typeChangeVariable='abs',typeClimateResponse="tweedie")

# call QUALYPSO
Q_SWE = QUALYPSO(Y=Y_SWE,scen=scen_SWE,X=X_SWE_WL,Xfut=seq(from=0.5,to=3,by=0.1),
listOption=listOption)

# plot climate responses
plotQUALYPSO_CCR(Q_SWE)

5 Treatment of unbalanced ensembles

As indicated in the introduction, the main motivation for the development of this package was to provide an easy access to statistical approaches that are adapted to unbalanced ensembles of climate projections. As discussed in Evin et al. (2021), this issue is sometimes ignored in some uncertainty assessments, leading to biased estimates when direct averages are taken to estimate the main effects. Sometimes, subselection of available projections are done to obtain a balanced ensemble, leading to a dramatic waste of information. Evin et al. (2019) proposes a Bayesian approach that applies augmentation techniques to estimate the main effects of incomplete ensembles. The main downside of this approach is the computational burden, since it relies on sampling the posterior distributions with Markov Chain Monte Carlo (MCMC) techniques. Large samples are needed to represent adequately the posterior distributions. Comparatively, direct estimates using ordinary least squares, as implemented in the function lm, provide similar results almost instantaneously. In the following example, we provide a comparison of both approaches where we downsampled the illustrative ensemble of winter temperature projections, such that the GCMs HadGEM2-ES and MPI-ESM-LR appear only once in the ensemble. The two estimation approaches provide very similar results in terms of mean (plain curves) and credible/confidence intervals (see the larger intervals for the GCMs HadGEM2-ES and MPI-ESM-LR). In terms of running time, the application of QUALYPSO using the lm approach takes about 0.1 s, whereas the Bayesian approach, with 5,000 MCMC samples (2,000 burn-in samples and 3,000 kept samples) takes about 18 s. In addition, while the lm approach leads to a smooth temporal evolution for all estimates (e.g. climate responses, main effects, uncertainty components), the Bayesian approach does not, due to the limited sampling of the posterior distributions.

# remove scenario for one GCM
Y_sel = Y_DJFTas[-c(12,13,14,15,17,18,19,20),]
scen_sel = scen_DJFTas[-c(12,13,14,15,17,18,19,20),]

# lm
listOption = list(typeChangeVariable='abs',typeClimateReponse = "smooth.spline", parClimateResponse = 3)
system.time({Q_lm = QUALYPSO(Y=Y_sel,scen=scen_sel,X=1971:2099,Xfut=1990:2099,listOption=listOption)})

##    user  system elapsed 
##    0.24    0.01    0.23

# Bayesian
listOption = list(typeChangeVariable='abs',parClimateResponse=3, ANOVAmethod="Bayesian",nBurn=2000,nKeep=5000)
system.time({Q_Bayes = QUALYPSO(Y=Y_sel,scen=scen_sel,X=1971:2099,Xfut=1990:2099,listOption=listOption)})

##    user  system elapsed 
##   44.72    1.28   46.90

# main GCM effects
lQ = list(lm=Q_lm,Bayesian=Q_Bayes)
for(method in names(lQ)) plotQUALYPSO_effect(lQ[[method]],nameEff="GCM",xlab="Years",main=method)

6 Estimating climate responses and internal variability before calling QUALYPSO

QUALYPSO embeds different regression models (smoothing splines, linear models, polynomial models) to estimate the climate response, and one glm model (tweedie distributions). In some cases, these options might be inadequate or too limited. As an alternative, it is possible to estimate the climate responses and internal variability using other approaches and provide the different corresponding objects (the climate change responses, the climate responses, the deviations from the climate response, and internal variability) in listOptions. The different steps can be summarized as follows:

Outside QUALYPSO

1. estimate \(CR(t)\) and \(CCR(t)\) for each projection (e.g. using robust regression),
2. extract fluctuations \(D(t)\) from IV for each projection and model the temporal evolution of \(IV(t)\),

With QUALYPSO

3. estimate the main effects and all uncertainty components,
4. use plotting functions and statistics.

Here is an example based on the function rlm of the package MASS that applies robust linear mmodels. This approach is less sensitive to outliers and can be interesting if some rare values are abnormally large in the projections. It also shows the piece of code used to estimate internal variability using gamlss.

library(MASS)
library(gamlss)

# define dimensions
Xs = 1971:2099
Xfut = 1990:2099
Xref = Xfut[1]
nY = length(Xs)
nFut = length(Xfut)
nS = nrow(Y_DJFTas)

# initialise matrices
CCR_tfut = CR_tfut = matrix(nrow=nS,ncol=nFut)
D_t = CHANGEVAR_t = matrix(nrow=nS,ncol=nY)

for(iS in 1:nS){
  Ys = Y_DJFTas[iS,]

  # data.frame with Xs and Ys
  dfpred = data.frame(x=Xs,y=Ys)

  # apply robust regression
  rlm_out = rlm(y~x,data = dfpred)

  # fitted responses at the points of the fit (for Dt), at Xfut and at Xref
  CR_t_s = predict(rlm_out, newdata = data.frame(x = Xs))
  CR_tfut_s = predict(rlm_out, newdata = data.frame(x = Xfut))
  CR_tref = predict(rlm_out, newdata = data.frame(x = Xref))     

  # store climate response for this simulation chain
  CR_tfut[iS,] = CR_tfut_s

  # fill matrices
  CCR_tfut[iS,] = CR_tfut_s-CR_tref
  D_t[iS,] = Ys-CR_t_s
  CHANGEVAR_t[iS,] = Ys-CR_tref
}

# Heteroscedasticity: assume gaussian deviations from the climate responses
# with a mean of zero
Xmat = matrix(rep(Xs,nS),byrow=T,nrow=nS,ncol=nY)
df_fit = data.frame(x = as.vector(Xmat), y = as.vector(D_t))
gamlss_out = gamlss(y ~ x, sigma.formula = ~x, family = gamlss.dist::NO(), data = df_fit,mu.fix = TRUE, 
mu.start = 0,control = gamlss.control(trace = FALSE))

# predict IV(t) for Xfut
sigma_pred = predict(object = gamlss_out, what = "sigma", type = "response", newdata = data.frame(x=Xfut),
data = df_fit)
IV_tfut = sigma_pred^2

# return objects
ClimateResponse = list(CR_tfut=CR_tfut,CCR_tfut=CCR_tfut,D_t=D_t,CHANGEVAR_t=CHANGEVAR_t,IV_tfut=IV_tfut)

# call QUALYPSO with this list of custom climate responses
Q_customCR = QUALYPSO(Y=Y_DJFTas,scen=scen_DJFTas,X=Xs,Xfut=Xfut,
listOption=list(ClimateResponse=ClimateResponse))

7 Limitations and caveats

As discussed in section concepts, QUALYPSO relies on several assumptions that must be verified to ensure that outputs provide a fair assessment of the main uncertainties and their contribution to the total uncertainty. In particular:

• Estimation of the climate response: The estimation of the climate responses is an important step that must not be overlooked. The fitted climate responses must reflect the expected reponses to climate change. The function plotQUALYPSO_CR provides a visualisation tool that can help to guide the choice of the smoothing parameter or the family of climate response models. As has been illustrated in the fourth example with SWE maxima, default options can lead to inadequate fitted climate responses.

• internal variability: QUALYPSO assumes that \(D_i(t)\) (deviations from the climate reponses) follow a normal distribution where the mean is zero and the logarithm of the standard deviation is a linear function of \(t\). This simple assumption is generally reasonable but must be verified using the function plotQUALYPSO_IV. In the example with SWE maxima, it can be seen that the normal distribution is not verified, as \(D_i(t)\) reach high positive values but not large negative values.

plotQUALYPSO_IV(Q_SWE)

• ANOVA model: Similarly, QUALYPSO assumes a simple linear model for the decomposition of the main sources of uncertainty, where the climate change reponse is assumed to be a sum of the main fixed effects and a residual error. The residual errors are assumed to be identically distributed (same variance) and for small set of projections, to follow a Gaussian distribution. More complex ANOVA decomposition could involve interactions between the different factors and/or random effects (Northrop and Chandler 2014). The motivation for this simple additive model is to provide a framework that can be applied for small ensembles of projections, which is generally the case of most applications. However, in some cases, and in particular with large ensembles Evin et al. (2026), the residual variability can be large, suggesting that potential interactions between scenarios and models and/or between models of different categories are included in this component. In other cases, a large residual variability indicates that climate responses are poorly fitted.

8 Conclusion

The package QUALYPSO is designed to tackle the main challenges of uncertainty partitioning in multi-model multi-scenario ensembles, and leverages large ensemble of projections used in climate change impact assessements. QUALYPSO enables the characterization of projected climate change responses (CCRs) across various climate and climate-related indicators. Specifically, it provides insights into the mean and the dispersion of the CCRs obtained from different modeling chains. Moreover, it facilitates the estimation of all sources of uncertainty – e.g. scenario, GCM, RCM, and impact model uncertainties – as well as potential additional dispersion in future changes arising from internal variability. Lastly, it enables the identification of the main effects associated with each model within each model category, providing a clearer understanding of how individual models compare and where the primary sources of uncertainty arise.

References

Aitken, Gordon, Lindsay Beevers, Simon Parry, and Katie Facer-Childs. 2023. “Partitioning Model Uncertainty in Multi-Model Ensemble River Flow Projections.” Climatic Change 176 (11): 153. https://doi.org/10.1007/s10584-023-03621-1.

Bichet, Adeline, Benoit Hingray, Guillaume Evin, Arona Diedhiou, Cheikh Mouhamed Fadel Kebe, and Sandrine Anquetin. 2019. “Potential Impact of Climate Change on Solar Resource in Africa for Photovoltaic Energy: Analyses from CORDEX-AFRICA Climate Experiments.” Environmental Research Letters 14 (12): 124039. https://doi.org/10.1088/1748-9326/ab500a.

Evin, G., E. Le Roux, E. Kamir, and S. Morin. 2025. “Estimating Changes in Extreme Snow Load in Europe as a Function of Global Warming Levels.” Cold Regions Science and Technology 231 (March): 104424. https://doi.org/10.1016/j.coldregions.2025.104424.

Evin, Guillaume, Benoit Hingray, Juliette Blanchet, Nicolas Eckert, Samuel Morin, and Deborah Verfaillie. 2019. “Partitioning Uncertainty Components of an Incomplete Ensemble of Climate Projections Using Data Augmentation.” Journal of Climate 32 (8): 2423–40. https://doi.org/10.1175/JCLI-D-18-0606.1.

Evin, Guillaume, Benoit Hingray, Guillaume Thirel, et al. 2026. “Uncertainty Sources in a Large Ensemble of Hydrological Projections: Regional Climate Models and Internal Variability Matter.” Hydrology and Earth System Sciences 30 (4): 1023–51. https://doi.org/10.5194/hess-30-1023-2026.

Evin, Guillaume, Samuel Somot, and Benoit Hingray. 2021. “Balanced Estimate and Uncertainty Assessment of European Climate Change Using the Large EURO-CORDEX Regional Climate Model Ensemble.” Earth System Dynamics 12 (4): 1543–69. https://doi.org/10.5194/esd-12-1543-2021.

Hawkins, E., and R. Sutton. 2009. “The Potential to Narrow Uncertainty in Regional Climate Predictions.” Bulletin of the American Meteorological Society 90 (8): 1095–107. https://doi.org/10.1175/2009BAMS2607.1.

Hingray, Benoit, Juliette Blanchet, Guillaume Evin, and Jean-Philippe Vidal. 2019. “Uncertainty Component Estimates in Transient Climate Projections.” Climate Dynamics 53 (5): 2501–16. https://doi.org/10.1007/s00382-019-04635-1.

Hingray, B., and M. Saïd. 2014. “Partitioning Internal Variability and Model Uncertainty Components in a Multimember Multimodel Ensemble of Climate Projections.” Journal of Climate 27 (17): 6779–98. https://doi.org/10.1175/JCLI-D-13-00629.1.

Jørgensen, Bent. 1987. “Exponential Dispersion Models.” Journal of the Royal Statistical Society: Series B (Methodological) 49 (2): 127–45. https://doi.org/10.1111/j.2517-6161.1987.tb01685.x.

Lehner, Flavio, and Clara Deser. 2023. “Origin, Importance, and Predictive Limits of Internal Climate Variability.” Environmental Research: Climate 2 (2): 023001. https://doi.org/10.1088/2752-5295/accf30.

Lemaitre-Basset, Thibault, Lila Collet, Guillaume Thirel, Juraj Parajka, Guillaume Evin, and Benoît Hingray. 2021. “Climate Change Impact and Uncertainty Analysis on Hydrological Extremes in a French Mediterranean Catchment.” Hydrological Sciences Journal 66 (5): 888–903. https://doi.org/10.1080/02626667.2021.1895437.

Northrop, Paul J., and Richard E. Chandler. 2014. “Quantifying Sources of Uncertainty in Projections of Future Climate.” Journal of Climate 27 (23): 8793–808. https://doi.org/10.1175/JCLI-D-14-00265.1.

Scherrer, Simon C., Cees de Valk, Michael Begert, Stefanie Gubler, Sven Kotlarski, and Mischa Croci-Maspoli. 2024. “Estimating Trends and the Current Climate Mean in a Changing Climate.” Climate Services 33 (January): 100428. https://doi.org/10.1016/j.cliser.2023.100428.

Yip, S., C. A. T. Ferro, D. B. Stephenson, and E. Hawkins. 2011. “A Simple, Coherent Framework for Partitioning Uncertainty in Climate Predictions.” Journal of Climate 24 (17): 4634–43. https://doi.org/10.1175/2011JCLI4085.1.