Parametric Quantile Copula Regression for Bivariate Median Based Unit Rayleigh

1. Introduction

In the realm of multivariate regression analysis, have you ever found yourself navigating the complex landscape where multiple response variables are documented for each individual or experimental unit? Imagine the interplay of various economic indicators, such as the delicate dance between income and expenditure, or the intricate relationships in biomedical research where systolic and diastolic blood pressure readings are closely intertwined. Similarly, consider environmental metrics like temperature and humidity, whose fluctuations may influence one another in subtle yet profound ways.

In these scenarios, it's striking to observe how the response variables often display a profound statistical dependence, even when we account for the effects of explanatory variables. Ignoring this intricate web of relationships can lead to inefficiencies in estimation, biases in inference, and potentially misleading conclusions regarding the true influence of predictors. Each overlooked connection risks distorting our understanding of the underlying phenomena we seek to unravel.

Traditional regression methods often involve fitting individual univariate models for each response variable. However, these methods typically assume that the responses are conditionally independent given the predictor variables, a premise that is seldom met in real-world scenarios. In practice, correlations between outcomes can stem from various sources, including shared latent factors, common measurement processes, or intrinsic structural relationships among the responses themselves. When such dependencies are present, relying on separate univariate models fails to effectively represent the complex joint behavior of the responses. Consequently, the marginal predictions produced may not align with realistic joint distributions, leading to potentially misleading conclusions.

To overcome this limitation, copula-based regression presents a flexible and theoretically robust framework for modeling the joint distribution of multiple responses, while simultaneously maintaining distinct regression models for each marginal distribution. Sklar's theorem provides a foundational framework for understanding multivariate distributions by allowing them to be expressed as a combination of their marginal distributions and a copula function. This copula function serves to describe the dependence structure among the different variables. By utilizing this decomposition, researchers can develop tailored regression models for each marginal outcome, potentially employing various families or link functions. Furthermore, the copula effectively captures the strength and nature of the relationships between these outcomes, whether they are symmetric, asymmetric, or exhibit tail dependence, thereby enhancing the analysis of complex data sets.

When two response variables $Y_1$ and $Y_2$ depend on a common predictor $X$ , the parametric copula regression model expresses their conditional joint distribution as

$$F_{Y_1,Y_2|X}\left(y_1,y_2|x\right)=C_{\theta \left(x\right)}\left(F_1\left(y_1|x\right), F_2\left(y_2|x\right)\right)$$

where $F_1$ and $F_2$ are the conditional marginal distributions and $C_{\theta \left(x\right)}$ is a parametric copula with a dependency parameter $\theta$ that may itself vary with the predictor $x$.

This formulation empowers both the marginal behaviors and the strength of association between responses to evolve dynamically in response to predictors. Such adaptability is crucial across a multitude of applications. For example, in the realm of finance, the relationship between asset returns can intensify dramatically during times of market turbulence. In environmental studies, the interplay between rainfall and humidity can vary significantly based on factors like altitude or the changing seasons. In biomedical research, the connection between two biomarkers may fluctuate with age or the prescribed treatment dosage. A copula regression model adeptly captures these intricate, covariate-dependent patterns of dependence in a manner that traditional correlation-based or linear models simply cannot achieve.

The parametric copula framework opens a door to a sophisticated realm of understanding the intricate nature of dependence between variables. By dynamically estimating theta as a function of the predictor, researchers can rigorously test important hypotheses. For instance, one can determine whether the degree of dependence remains stable across all levels of predictor or if it exhibits a tendency to increase or decrease in response to the predictor. This transforms copula regression into more than just a modeling technique; it becomes a powerful inferential framework that illuminates how the relationships between outcomes evolve in tandem with various covariates.

In conclusion, dependence is important because it indicates structural relationships that cannot be fully understood through the predictors alone. Overlooking dependence may lead to an underestimation of uncertainty, skew prediction intervals, and result in inconsistent joint forecasts. The parametric copula regression model offers a coherent approach by enabling distinct, interpretable marginal regressions while also capturing the dependence structure among responses using a well-founded copula function.

Contribution of the study: This study stands apart from earlier research by examining dependence through parametric copulas instead of relying on nonparametric or vine copula methods. This unique approach adds significant value to the work.

Recent developments in multivariate statistics have highlighted the role of copula models, yet much of the existing research has concentrated on static dependence, where the copula parameters are treated as constant. Additionally, there has been a significant focus on high-dimensional vine structures, which are primarily utilized for risk aggregation and portfolio analysis. However, there has been relatively less exploration of copula-based regression models that jointly regress two or more response variables on covariates. Notably, the modeling of the dependence parameter as a function of a predictor offers valuable insights, as it enables researchers to directly interpret how the strength or nature of dependence varies across different levels of the covariate. This approach, while promising, has not been sufficiently examined in applied fields.

In the foundational works of the late 1990s and early 2000s (e.g., Joe, 1997 [1]), scholars laid the groundwork for understanding conditional copulas and the intricate dynamics of dependence modeling. However, these early theories often faced limitations in their practical applications, typically restricted to specific types of data, including continuous time series and financial returns, which constricted their broader utility. Many authors discussed the rule copula in regression like [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]

Recent advancements in statistical modeling, particularly in copula regression, have been significantly influenced by the pioneering research conducted by Hans in 2023[3] and Klein in 2019 [21]. Their innovative work has expanded the boundaries of copula regression, integrating both distributional and boosted frameworks that adapt dynamically to varying data conditions. These contemporary methodologies not only enhance the estimation processes for both marginal and copula parameters but also demonstrate impressive flexibility across different contexts. However, despite the evident sophistication and potential of these advanced techniques, they often require the implementation of complex computational algorithms, which can be challenging to navigate. Additionally, they typically necessitate large sample sizes to yield reliable results. This reliance on robust data sets may inadvertently create barriers for practitioners who are working with smaller datasets, limiting the accessibility and applicability of these advanced statistical tools in real-world scenarios.

As a result, there is still a pressing need for more straightforward parametric copula regression models—ones that offer clarity and ease of interpretation, making them accessible and practical for those dealing with small to moderate samples while preserving the richness of the underlying statistical relationships.

In this study, the author delves into the innovative realm of parametric copula regression, exploring its application to two interrelated response variables linked through a single-dimensional predictor. This approach skillfully marries the clarity and interpretability of parametric regression with the adaptable capabilities of copula-based dependence modeling, resulting in a powerful analytical framework. The author's contribution unfolds in two significant ways:

1-The author establishes a robust modeling structure in which each response variable is intricately regressed on the predictor via a parametric quantile regression approach. Crucially, the dependence parameter of a selected copula (Gumbel) is elegantly modeled as a parametric function of the same predictor, enriching the understanding of their interrelationships.

2- The author illustrates the model's capacity to unveil varying dependence patterns across the predictor's range. Through meticulous evaluation employing bootstrap-based inference techniques and rigorous goodness-of-fit diagnostics, its performance can be assessed. Furthermore, its interpretability in practical scenarios where conventional multivariate regression methods struggle to capture crucial nuances, such as tail or asymmetric dependence is highlighted.

This comprehensive investigation not only advances methodological rigor but also enhances empirical insights, illuminating the complex tapestry of relationships that lie within our data. It also offers a constructive framework for the analysis of bivariate responses with covariate-dependent dependence by effectively bridging marginal regression models with dependence modeling through a parametric copula structure. By doing so, it enhances our theoretical understanding of regression copulas, while also providing valuable tools for practical applications in fields such as econometrics, environmental modeling, and biostatistics. This approach encourages further exploration and innovation in these areas.

This paper is structured into the following sections. Methodology and mathematical definitions are discussed in section 2. Estimation procedures and diagnostics are clarified in section 3. Real data analysis comprehending the results and the discussion is demonstrated in section 4. Conclusion is illuminated in section 5. Lastly future work is suggested in section 6.

2. Methodology and Mathematical definitions of the model

Let $\left(Y_{1i}, Y_{2i}, X_i\right)$ , $i=\mathrm{1,2},\dots n$ denote independent observations of two response variables $Y_1$, $Y_2$and a predictor $X$ . The objective is to model the joint conditional distribution of $\left(Y_1 , Y_2\right)$ given$X$ while allowing both the marginal behaviors and their dependence structure to vary systematically with the predictor.

2.1. Conditional Joint Distribution via Sklar’s Theorem

According to Sklar’s theorem, for any continuous bivariate distribution, there exists a unique copula function $C_{\theta }\left(.,.\right)$such that the joint conditional cumulative distribution function (CDF) can be written as

$$F_{Y_1,Y_2|X}\left(y_1,y_2|x\right)=C_{\theta \left(x\right)}\left(F_1\left(y_1|x\right), F_2\left(y_2|x\right)\right)$$

where $F_1\left(y_1|x\right)$and $F_2\left(y_2|x\right)$ are the conditional marginal CDFs. $C_{\theta \left(x\right)}$ is a parametric copula with parameter $\theta \left(x\right)$ possibly depending on the covariate x, and $\theta \left(x\right)$ controls the strength and the type of dependence between $Y_1$ and $Y_2$. If both marginals are continuous, differentiation of $F_{Y_1,Y_2|X}$ with respect to $y_1$ and $y_2$ yields the joint conditional density.

$$f_{Y_1,Y_2|X}\left(y_1,y_2|x\right)=c_{\theta \left(x\right)}\left(F_1\left(y_1|x\right), F_2\left(y_2|x\right)\right)f_{Y_1}\left(y_1|x\right)f_{Y_2}\left(y_2|x\right)$$

where $c_{\theta \left(x\right)}$ is the copula density given by

$$c_{\theta \left(x\right)}\left(u, v\right)={\displaystyle \frac{{\partial }^2{C}_{\theta \left(x\right)}\left(u,v\right)}{\partial u\partial v}}$$

2.2. Marginal Regression Models

Each response variable is first modeled conditionally on X using an appropriate parametric or semi-parametric regression structure. The author used parametric quantile regression (median). The marginal model can take a specific form such as : linear regression with $E\left[Y_i|X=x\right]={\beta }_{0i}+{\beta }_{1i}x$ , generalized linear model (GLM): $g_{i}E\left[Y_i|X=x\right]={\beta }_{0i}+{\beta }_{1i}x$, and Quantile regression: $Q_{Y_i|X}\left(\tau \right)={\beta }_{0i}+{\beta }_{1i}x$ , the fitted conditional CDFs ${\widehat{F}}_1\left(y_1|x\right)$ and ${\widehat{F}}_2\left(y_2|x\right)$ are then used to transform observations into pseudo-observations: $u_{1i}= {\widehat{F}}_{Y_1}\left(y_{1i}|x_i\right)$ and : $u_{2i}= {\widehat{F}}_{Y_2}\left(y_{2i}|x_i\right)$.

2.3. Dependence Structure and Copula Regression

The dependence between y1 and y2 is modeled using a parameric copula family. Typical choices include: Gaussian , Clayton, Gumbel and Frank copula. To allow the dependence to vary with the covariate, a parametric regression structure for $\theta \left(x\right)$ is defined: $\theta \left(x\right)=h\left({\gamma }_0+{\gamma }_{1}x\right)$ where $h\left(.\right)$ is a link function ensuring $\theta \left(x\right)$ lies within its permissible range such as $\theta \left(x\right)=1+exp\left({\gamma }_0+{\gamma }_{1}x\right)$ for Gumbel copula; $\theta \left(x\right)=exp\left({\gamma }_0+{\gamma }_{1}x\right)$ for Clayton copula; and $\theta \left(x\right)=tanh\left({\gamma }_0+{\gamma }_{1}x\right)$ for Gaussian copula.

2.4. Log-Likelihood and Estimation

Given $n$ independent observations, the conditional Log-Likelihhood function of the full model is

$$l\left({\beta }_1, {\beta }_2, \gamma \right)=\sum _{i=1}^n\log \left(f_{Y_1}\left(y_{1i}|x_i; {\beta }_1\right)\right) +\sum _{i=1}^n\log \left(f_2\left(y_{2i}|x_i; {\beta }_2\right)\right)+\sum _{i=1}^n\log \left(c_{\theta \left(x_i\right)}\left(u_{1i}, u_{2i}\right)\right)$$

where : $u_{1i}= F_{Y_1}\left(y_{1i}|x_i\right)$ and : $u_{2i}= F_{Y_2}\left(y_{2i}|x_i\right)$. Two common estimation strategies are:

Full Maximum Likelihhood (FML): estimate all parameters $\left({\beta }_1, {\beta }_2, \gamma \right)$ simultaneously by maximizing the joint log-likelihood . This is statistically efficient but computationally intensive.

Inference functions for the Margins (IFM): In step one; B1 and B2 are estimated from marginal regressions. And in step 2; compute the pseudo observations $\left(u_{1i}, u_{2i}\right)$ and estimate $\gamma$by maximizing $l_c\left(\gamma \right)=\sum _{i=1}^n\log \left(c_{\theta \left(x_i\right)}\left(u_{1i}, u_{2i}\right)\right)$ . IFM is simpler, computationally stable, and asymptotically consistent.

2.5. Interpretation

The estimated coefficients $\left({\beta }_1, {\beta }_2\right)$ describe how each response variable depends marginally on the predictor X, while the copula parameter $\theta \left(x\right)$ (or its implied dependence measure; Kendall’s $\tau \left(x\right)$) quantifies how the strength of dependence between $Y_1$ and$Y_2$ changes with X. A positive estimate of $\gamma$ in $\theta \left(x\right)=h\left({\gamma }_0+{\gamma }_{1}x\right)$ indicates that the dependence increases with the predictor, whereas a negative value implies weakening dependence.

3. Estimation Procedure and Model Diagnostics

This section covers how estimation is performed in practice, how uncertainty is quantified, and how to check model adequacy (both for the the marginal regression and for the copula dependence structure).

3.1. Estimation Framework

Estimation of the parametric copula regression model proceeds by either full maximum likelihood(FML) or two-step inference functions for Margins(IFM). Given that FML can be computationally demanding, the IFM approach is ofen preferred because it separates marginal and copula estimation into two manageable stages without major loss of efficiency under correct model specification.

Step1: estimation of marginal models

Each marginal regression $f_{Y_i}\left(y_i|x;{\beta }_i\right)$ is fitted separately using standard estimation techniques:

$${\widehat{\beta }}_j=\arg \underset{{\beta }_j}{\max }\sum _{i=1}^n\log \left(f_{Y_j}\left(y_{ij}|x_i;{\beta }_j\right)\right) , j=\mathrm{1,2}$$

After obtaining ${\widehat{\beta }}_1$ and ${\widehat{\beta }}_2$, the pseudo-observations ( empirical probability residuals) are computed $u_{1i}= {\widehat{F}}_{Y_1}\left(y_{1i}|x_i\right)$ and : $u_{2i}= {\widehat{F}}_{Y_2}\left(y_{2i}|x_i\right)$. These observations are approximately independent and uniformly distributed on [0,1] if the marginal models are correctly specified.

Step 2: Estimation of the copula regression

Given the pseudo observations, the dependence parameters $\gamma$ are estimated by maximizing the copula Log-Likelihood:

$$\widehat{\gamma }=\arg \underset{\gamma }{\max }\sum _{i=1}^n\log \left(c_{\theta \left(x_i;\gamma \right)}\left(u_{1i}, u_{2i}\right)\right)$$

where ${\theta }_i\left(x_i; \gamma \right)=h\left({\gamma }_0+{\gamma }_1{x}_i\right)$ is the parametric function linking the predictors to the copula parameter. You parameterize the copula PDF with thethis link function so that you can estimate the coefficients. Optimization can be performed using the derivative-free algorithms such as Nelder-Mead Simplex or gradient based methods such as BFGS, depending on smoothness of the $c_{\theta \left(x_i;\gamma \right)}$ . Standard errors for $\widehat{\gamma }$ can be obtained either from the estimated information matrix or via parametric bootstrap, which is particularly useful when asymptotic approximations may be unreliable.

Step 3: inference on dependence dynamics

After estimation, the fitted dependence parameter can be evaluated as a function of X: ${\widehat{\theta }}_i\left(x_i; \gamma \right)=h\left({\gamma }_0+{\gamma }_1{x}_i\right)$ , and the corresponding Kendall’s tau can be computed analytically according to the copula family. The function ${\theta }_i\left(x_i; \gamma \right)$ or $\widehat{\tau }\left(x\right)$quantifies how the strength of dependence changes with the predictors. A hypothesis test can be conducted for

$H_0 : {\gamma }_1=0 ( constant dependnece)$ vs. $H_1 : {\gamma }_1\ne 0 ( dependnece varies with x)$

Because the standard asymptotic tests may be inaccurate, it is common to perform a parametric bootstrap likelihood ratio test (LRT):

• Fit both the full model ( with ${\gamma }_1$) and the reduced model $\left({\gamma }_1=0\right)$.
• Compute the observed LRT statistics $T_{obs}=2\left(l_{full}-l_{reduced}\right).$
• Simulate bootstrap samples under the null model $\left({\gamma }_1=0\right)$, re-estimate both models, and compute $T_b$for each sample.
• The empirical p-value is $\widehat{p}={\displaystyle \frac{1}{B}} \sum _{b=1}^{B}I\left(T_b\ge T_{obs}\right)$.

A small p-value (for example < 0.05) provides evidence that the dependence parameter varies significantly with the predictors.

3.2. Model Diagnostics

Model checking is a critical part of copula regression analysis. Diagnostics should address both marginal adequacy and copula adequacy.

3.2.1. Marginal Model Diagnostics

To ensure that the marginal models correctly describe each response conditional on X, the following residuals are examined:

A. 

Randomized quantile residuals : $r_{ji}= {\mathsf{\Phi }}^{-1}\left(u_{ji}\right)$where ${\mathsf{\Phi }}^{-1}$is the inverse standard normal CDF. If the model fits well, these residuals should follow approximately a standard normal distribution.

B. 

Cox-Snell residuals or Pearson residuals for parametric margins

C. 

Graphical checks: scatter plots of $r_{ji}$versus Xi should show no systematic trends or hetero-sc edastic patterns. Also QQ-plots of $r_{ji}$against standard normal quantiles can help detect deviations from the assumed marginal distribution.

If marginal adequacy fails, the copula estimation step may be biased, so corrections (changing the link functions or marginal distributions) should be made first.

3.2.2. Copula Goodness-of-fit diagnostics

Once the marginal are validated, the fitted copula should be checked to ensure it correctly represents the dependence structure. Typical diagnostics include:

A. 

Empirical vs. fitted copula plots: compare the empirical copula $C_n\left(u,v\right)$with the fitted copula $C_{\widehat{\theta }\left(x\right)}\left(u,v\right)$visually or through contour plots.

B. 

Goods of fit tests: apply formal test such as Cramer Von Mises, Anderson Darling, or Kendall’s process-based tests for copulas.

C. 

Tail Dependence Checks: compare empirical upper/lower tail probabilities to those implied by the fitted copula to verify correct tail behavior (important for asymmetric copulas like Gumbel or Clayton).

D. 

Dependence vs. Predictor Plots: plot the estimated Kendall’s tau $\widehat{\tau }\left(x\right)$or $\widehat{\theta }\left(x\right)$against $x$. If the functional form is correctly specified, the trend should align with the theoretical expectations and exhibits random scatter around the fitted curve.

3.3. Uncertainty Quantification Via Bootstrap

Bootstrap methods are widely used to assess the sampling variability of parameter estimates and derived quantiles such as $\widehat{\tau }\left(x\right)$ or predicted joint probabilities. The general algorithm is:

• Generate $B$ bootstrap samples $\left(Y_{1i}^{*\left(b\right)}, Y_{2i}^{*\left(b\right)}, X_i\right)$ by resampling from the fitted model.
• Re-estimate parameters $\left({\widehat{\mathit{\beta }}}_1^{\mathit{*}\left(\mathit{b}\right)},{\widehat{\mathit{\beta }}}_2^{\mathit{*}\left(\mathit{b}\right)},{\mathit{\gamma }}^{\mathit{*}\left(\mathit{b}\right)}\right)$ for each bootstrap sample.
• Compute standard errors, confidence intervals, and empirical distributions of the estimates.
• If the bootstrap distribution of ${\gamma }_1^{*\left(b\right)}$ is centered near zero, this suggests no covariate effect on dependence.

This approach provides more accurate inference in small or moderate samples compared to asymptotic approximations.

3.4. Visualization and Interpretation

Finally, visualization helps to communicate model behavior:

• Plot fitted marginal regression lines $E\left[Y_j|X=x\right]$ alongside observed data.
• Plot estimated dependence parameter $\widehat{\theta }\left(x\right)$ or its Kendall tau $\widehat{\tau }\left(x\right)$ equivalent versus $x$ to visualize how dependence evolves with the predictor.
• Contour plots of the fitted joint density $f_{Y_1,Y_2|X}\left(y_1,y_2|x\right)$ for selected values of x can illustrate how both marginal behavior and dependence change.

Such plots provide intuitive insights into whether dependence strengthens or weakens with the covariate, and whether tail behavior aligns empirical observations.

4. Real Data Analysis: Results and Discussion

Parametric Quantile Copula regression algorithm steps

• Assess the distribution of the two response variables and their fitness.
• Assess the fitness of the type of family of copula that links the two variables.
• Apply marginal parametric quantile regression: regress each variable on the predictor using the parametric quantile regression and assess the fitness of each marginal quantile regression model then obtain the fitted theoretical CDF for each response variable.
• Use these fitted CDFs in the copula regression part of the algorithm (applying the IFM).

Organization for Economic Co-operation and Development (OECD) data platform located on the following site: https://stats.oecd.org/index.aspx?DataSetCode=BLI, delivers indicators in different segments such as economy, education, health, labor, trade and creation. The data are collected from 41 countries all over the world. In this manuscript, the author will regress two response variables or indicators; the water quality and the quality of support network on the predictor or the indicator; the life expectancy. The two response variables are distributed as median based unit Rayleigh (MBUR) introduced by Attia [33]

These illustrated flow diagrams show the steps to follow in applying parametric copula regression model.

The water quality indicator is computed as the self-reported approval with water quality and is signified by the percentage of people who report being satisfied with the quality of water as regards safety and cleanliness. The values of the water quality data are0.92, 0.92, 0.79, 0.9 ,0.62, 0.82, 0.87, 0.89, 0.93, 0.86, 0.97, 0.78, 0.91, 0.67, 0.81, 0.97, 0.8, 0.77, 0.77, 0.87, 0.82, 0.83, 0.83, 0.85, 0.75, 0.91, 0.85, 0.98, 0.82, 0.89, 0.81, 0.93, 0.76, 0.97, 0.96, 0.62, 0.82, 0.88, 0.7, 0.62, 0.72.

The quality of support network indicator measures the availability and reliability of help offered as a social support of friends, family, or community members to individuals who can rely on in time of need, and it is expressed as the percentage of people who report having someone they can count on for this support. The values of the quality support network are0.93, 0.92, 0.9, 0.93, 0.88, 0.8, 0.82, 0.96, 0.95, 0.95, 0.96, 0.94, 0.9, 0.78, 0.94, 0.98, 0.96, 0.95, 0.89, 0.89, 0.8, 0.92, 0.89, 0.91, 0.77, 0.94, 0.95, 0.96, 0.94, 0.87, 0.95, 0.95, 0.93, 0.94, 0.94, 0.85, 0.93, 0.94, 0.83, 0.89, 0.89.

The life expectancy indicator refers to the average number of years a person can expect to live based on current mortality rates and it is represented as a single number expressing the expected lifespan at birth, the values of the life expectancy data are 83, 82, 82.1, 82.1, 80.6, 76.7, 80.5, 79.3, 81.5, 78.8, 82.1, 82.9, 81.4, 81.7, 76.4, 83.2, 82.8, 82.9, 83.6, 84.4, 83.3, 75.5, 76.4, 82.7, 75.1, 82.2, 82.1, 83, 78, 81.8, 77.8, 81.6, 83.9, 83.2, 84, 78.6, 81.3, 78.9, 75.9, 73.2, 64.2.

First Step of Algorithm: The Distribution of the Response Variables.

The 2 response variable are distributed as MBUR with the following PDF , CDF and quantile functions, respectively: $f_Y\left(y;\alpha \right)={\displaystyle \frac{6}{{\alpha }^2}}\left[1-y^{{\displaystyle \frac{1}{{\alpha }^2}}}\right]y^{\left({\displaystyle \frac{2}{{\alpha }^2}}-1\right)}$ , $F_Y\left(y\right)={3y}^{{\displaystyle \frac{2}{{\alpha }^2}}}-{2y}^{{\displaystyle \frac{3}{{\alpha }^2}}}$and the quantile is

$$Q_Y\left(P\left(y\right)\right)=F^{-1}\left(y\right)={\left\{-.5\left(cos\left[{\displaystyle \frac{{cos}^{-1}\left(1-2u\right)}{3}}\right]-\sqrt{3} sin\left[{\displaystyle \frac{{cos}^{-1}\left(1-2u\right)}{3}}\right]\right)+.5\right\}}^{{\alpha }^2}$$

Figure 1 shows the boxplot of the response variables. Figure 2 shows the scatter plot between the two response variables. Figure 3 shows the scatter plot between the life expectancy water quality as well as the quality of support network. Figure 4 shows the fitted MBUR PDFs for both variables. Table 1 shows the descriptive statistics of the two response variables. Table 2 shows the Kendall Tau coefficient between the response variables. And Table 3 shows the statistical indices of both variables fitting MBUR distribution. The Kendall tau between water quality and life expectancy is 0.2859 and its statistically significant with p-value of 0.0099 while the Kendall tau between the quality of support network and the life expectancy is 0.2079 and it is statistically insignificant as the p-value is 0.0667. Although the distribution of the predictor does not affect the regression model, it is distributed as Kumaraswamy distribution as well as Beta distribution. But Kumaraswamy better fits the data rather than the Beta distribution as it has higher value of the Log-likelihood and more negative values of AIC, CAIC, BIC, and HQIC than the Beta distribution. To conduct the analysis, the predictor was transformed by dividing each value by 100 then taking the log of the result.

Table 1 depicts that the 2 response variables show left skewness that is more obvious for the quality of support network than that of the water quality. The quality of support network shows more leptokurtic shape than the shape of the water quality. Table 2 illuminates that both variables show significant positive dependence as seen by significant Kendall Tau coefficient

Table 3 highlights that the two variables fit the median based unit Rayleigh as null hypothesis fail to reject the tested MBUR distribution. Figure 4 shows the histogram and the fitted PDF curves for each response variable.

Second Step of Algorithm: The Copula Family That Links the Response Variables.

The response variables fit the Gumbel copula family. The theoretical tau is 0.3904 and the estimated dependency parameter is 1.6403. The estimated Cramer Von Mises (CVM) is 0.2154, the $AIC=-11.1298$ while $BIC= -9.4162$ . The Log-likelihood is 6.5649. Sampling from the fitted Gumbel copula using the Metropolis Hasting algorithm yielded a sampling distribution of the estimated tau with the 2.5^th and 97.5^th empirical quantiles being 0.2405 and 0.5432 respectively, a sampling distribution of the estimated dependency parameter with 2.5^th and 97.5^th empirical quantiles being 1.3166 and 2.1905 respectively, and a sampling distribution of the estimated CVM with 2.5^th and 97.5^th empirical quantiles being 0.0344 and 0.1754 respectively. Figure 5 shows the P-P plot of the theoretical vs the empirical Gumbel copula.

Figure 5. Shows the P-P plot of the empirical copula vs the theoretical copula. The theoretical copula shows near perfect alignment of the lower part and the central part with the diagonal than the upper part.

Figure 5. Shows the P-P plot of the empirical copula vs the theoretical copula. The theoretical copula shows near perfect alignment of the lower part and the central part with the diagonal than the upper part.

Figure 6. Shows the sampling distribution of the Kendall Tau coefficient. The estimated tau is 0.3904 and it lies between the 2.5^th and 97.5^th empirical quantile. The null hypothesis testing the estimated tau is failed to be rejected.

Figure 6. Shows the sampling distribution of the Kendall Tau coefficient. The estimated tau is 0.3904 and it lies between the 2.5^th and 97.5^th empirical quantile. The null hypothesis testing the estimated tau is failed to be rejected.

Figure 7, Figure 8 and Figure 9 show the sampling distributions obtained from the Metropolis Hasting algorithm for testing estimated tau, estimated dependency parameter, and the estimated CVM test. And the null hypothesis fails to reject any of the tested estimates. Hence the GoF tests for the Gumbel copula are illustrated and subsequently the two variables fit the Gumbel copula well.

Third Step of Algorithm: The Marginal Parametric Quantile Regression

As the response variables are negatively skewed, the logit link function and the log-log complementary link function can be considered as a good candidate to reparameterize the alpha parameter in the MBUR PDF. Both functions are monotone increasing functions. Hence applying the median parametric regression model to regress the water quality variable on the life expectancy variable as well as regressing the quality of support network on the life expectancy gives estimate of the alpha parameter for each observation. This is called the marginal regression model of each response variable on the predictor. The estimated alpha parameters of the marginal MBUR are utilized to acquire the fitted CDF for each marginal (or variable), these estimated CDF will be utilized in the copula regression part in which the dependency parameter of the Gumbel copula is defined in terms of each observation by the help of a link function that ensure that the dependency parameter is more than one.

Reparameterize the PDF and CDF of MBUR using the quantile function $y=c^{{\alpha }^2}$ where c is $c=-.5\left(cos\left[{\displaystyle \frac{{cos}^{-1}\left(1-2u\right)}{3}}\right]-\sqrt{3} sin\left[{\displaystyle \frac{{cos}^{-1}\left(1-2u\right)}{3}}\right]\right)+.5 , u\in \left(\mathrm{0,1}\right)$. U represents the chosen percentile, in this analysis the author uses the median so u=0.5 and ${\alpha }^2={\displaystyle \frac{ln\left(y\right)}{ln\left(c\right)}}={\displaystyle \frac{ln\left(median\right)}{ln\left(c\right)}}$

As y is the median corresponding to u=0.5 then; replacing this u into the c equation, the c will also be 0.5. Using the logit function of the median which is called $\mu$ and $i=\mathrm{1,2},..., n \& k=\mathrm{1,2}...,k$

$$Logit \left(median\right)=logit\left(\mu \right)=\phi =X^{\text{'}}B={\beta }_0+{\beta }_1{x}_{i1}+\dots +{\beta }_k{x}_{ik}$$

where n is the number of cases and k is the number of variables. Logit median is a linear combination of variables.

$$\mathrm{So} \mu ={\displaystyle \frac{e^{X^{\text{'}}B}}{1+e^{X^{\text{'}}B}}}={\displaystyle \frac{e^{\phi }}{1+e^{\phi }}} \mathrm{and} {\alpha }^2={\displaystyle \frac{ln\left(y\right)}{ln\left(c\right)}}={\displaystyle \frac{ln\left(median\right)}{ln\left(c\right)}}={\displaystyle \frac{ln\left(\frac{e^{X^{\text{'}}B}}{1+e^{X^{\text{'}}B}}\right)}{ln\left(c\right)}}={\displaystyle \frac{ln\left(\frac{e^{\phi }}{1+e^{\phi }}\right)}{ln\left(c\right)}}$$

Replace this parameter in the Log-Likelihood function and the ${\widehat{\mathit{\beta }}}_{\mathit{i}}$ coefficients can be obtained by the MLE, the author used the derivative free algorithm; Nelder Mead simplex optimizer in MATLAB.

A link function other than the logit which is also monotonically increasing is the log-log (1-median) $\log \left\{-\log \left(1-median\right)\right\}=\phi =X^{\text{'}}B= {\beta }_0+{\beta }_1{x}_{i1}+\dots +{\beta }_k{x}_{ik}$ , hence the median is $1-e^{-e^{X^{\text{'}}\beta }}$ and ${\alpha }^2={\displaystyle \frac{ln\left(y\right)}{ln\left(c\right)}}={\displaystyle \frac{ln\left(1-\mathrm{e}\mathrm{x}\mathrm{p}\left(-exp\left(X^{\text{'}}\beta \right)\right)\right)}{ln\left(c\right)}}$ . The author used the 2 link models and compare between them.

The results of the marginal regression of each of the response variable on the predictor are shown in 2 tables; each table represents the results gained from utilizing a specific link function. Table 4 shows the results of regressing the response variable on the predictor using the logit link function. While, Table 5 shows the results of regressing the response variable on the predictor using the log-log complementary link function.

Figure 9 shows the estimated curve for regressing the water quality and the quality of support network on transformed life expectancy. The graph supports the results of the significant effect of life expectancy on the water quality than its significant effect on the quality of support network. The pink curve is steeper than the blue one. Although the log-likelihood has higher value for the model incorporating the quality of support network than the model incorporating the water quality but the LRT points to the significance of the water quality over the quality of support network which is supported by Wald test statistics results and graph in Figure 9. The more the life expectancy is the more satisfaction of water quality is.

Figure 10 shows the QQ plot of the empirical quantiles versus the theoretical quantiles. The logit model regressing the water quality on the life expectancy has a curve of randomized quantile residuals that shows near perfect alignment with the diagonal all through its course, while the curve of the Cox-Snell residuals is mainly aligned with the diagonal at its lower tail and its center. Figure 11 shows the QQ plot for the logit model curve regressing the quality of support network on the life expectancy which has less efficient alignment with the diagonal as regards residuals, the randomized quantile residuals and the Cox-Snell residuals. Figure 12 and Figure 13 show the scatter plot of the predictor (transformed) and the both types of residuals (randomized quantile and Cox-Snell). There is no nonlinear relationship between the predictor and the residuals. This relationship can be expressed in Kendall tau coefficient as -0.1093 and associated p-value=0.3225 for the response variable (water quality) while it is -0.1104 and associated p-value=0.3171 for the response variable (quality of support network). These values are identical for both types of residuals for each response variable. So there is no evidence for dependency between the predictor and the residuals.

Figure 14 shows the estimated curve for regressing the water quality and the quality of support network on transformed life expectancy. The graph supports the results of the significant effect of life expectancy on the water quality than its significant effect on the quality of support network. The pink curve is steeper than the blue one. Although the log-likelihood has higher value for the model incorporating the quality of support network than the model incorporating the water quality but the LRT points to the significance of the water quality over the quality of support network which is supported by Wald test statistics results and graph in Figure 14. The more the life expectancy is the more satisfaction of water quality is.

Figure 15 shows the QQ plot of the empirical quantiles versus the theoretical quantiles. The log-log complementary model regressing the water quality on the life expectancy has a curve of randomized quantile residuals that shows near perfect alignment with the diagonal all through its course, while the curve of the Cox-Snell residuals is mainly aligned with the diagonal at its lower tail and its center. Figure 16 shows the QQ plot for the log-log complementary model curve regressing the quality of support network on the life expectancy which has less efficient alignment with the diagonal as regards the residuals, the randomized quantile residuals and the Cox-Snell residuals. Figure 17 and Figure 18 show the scatter plot of the predictor (transformed) and both types of residuals (the randomized quantile and the Cox-Snell). There is no nonlinear relationship between the predictor and the residuals. This relationship can be expressed in Kendall tau coefficient as -0.1093 and associated p-value=0.3225 for the response variable (water quality) while it is -0.1104 and associated p-value=0.3171 for the response variable (quality of support network). These values are identical for both types of residuals for each response variable. So there is no evidence for dependency between the predictor and the residuals.

Fourth Step of Algorithm: The Copula Parametric Quantile Regression

The quantile parametric copula regression is a copula regression model utilizing the fitted marginal CDFs obtained from the fitted marginal quantile parametric MBUR regression models. Regressing each copula dependency parameter on the corresponding observation with the aid of the link function ${\theta }_i\left(x_i\right)=h \left(g_0+g_1 x_i\right)=1+\exp \left(g_0+g_1 x_i\right)$ where h is the link function that ensures the dependency parameter is larger than one. The author used the Inference Function of Margin (IFM) and the Nelder Mead algorithm of MATLAB as an optimizer. The Gumbel copula PDF is

$$c_{\left(\theta ;u,v\right)}=e^{-{\left\{p+q\right\}}^{{\displaystyle \frac{1}{\theta }}}} {\left\{p+q\right\}}^{{\displaystyle \frac{1}{\theta }}-2}{\displaystyle \frac{{\left\{pq\right\}}^{\theta -1}}{uv}} \left\{{\left(p+q\right)}^{{\displaystyle \frac{1}{\theta }}}+\theta -1\right\}$$

where $={\left(-\ln u\right)}^{\theta } \& q={\left(-\ln v\right)}^{\theta }$ , reparameterize the dependency parameter ${\theta }_i\left(x_i\right)$ in this copula PDF with the ${\theta }_i\left(x_i\right)=h \left(g_0+g_1 x_i\right)=1+\exp \left(g_0+g_1 x_i\right)$, in other words, replace the right hand side of the link function in copula PDF and optimize. The u and v are the fitted CDFs obtained from the quantile marginal regression step. The relation between the theta dependence parameter and Kendall tau is $=1-{\displaystyle \frac{1}{\theta }}$ . After applying the copula regression in this step, there will be a value of theta for each $x_i$, and in Table 6, the mean of this tau is recorded as the mean of theoretical tau while the empirical tau is calculated from the u and v.

The Gumbel copula CDF is $C_{\left(\theta ;u,v\right)}=e^{-{\left\{p+q\right\}}^{{\displaystyle \frac{1}{\theta }}}}$

Table 6 enlightens the results of copula regression part using log-log complementary model and the logit model. The AIC, CAIC, BIC, HQIC for the logit model are less than those of the log-log complementary model. CVM for the logit is less than that of the log-log complementary model. The LRT is statistically significant for the logit model than the LRT for the complementary model. Figure 19, Figure 20 and Figure 21 explicate the histogram of the parametric bootstrap sampling of the estimated CVM, estimated intercept and the estimated slope, respectively obtained from the logit model depicting the descriptive statistics as shown in each figure. Figure 22 elucidates that the estimated LRT value (4.3426) lies between the 2.5^th and 97.5^th quantile of the parametric bootstrap sampling distribution. The slope is 5 % statistically significant using the logit model. The log-log complementary model shows 30% statistical significance of the estimated CVM, estimated intercept, and estimated slope using the parametric bootstrap sampling distribution as expounded in Figure 23, Figure 24 and Figure 25. Although, the estimated CVM is statistically significant at 5% significance level as shown from parametric sampling distribution in Figure 26, the estimated slope is statistically significant at 5% significant level as shown from the non-parametric bootstrap sampling distribution in Figure 27; which is less powerful than the parametric bootstrap sampling distribution. Figure 28 shows the parametric bootstrap sampling of the LRT illuminating that the LRT is significant at 30% significant level for the log-log complementary model.

From the above results, the logit model demonstrates better fit for the effect of predictor (life expectancy) on the dependency between the 2 response variables (water quality and quality of support network) than the log-log complementary model. The logit model is significant at 5% significant level while the log-log model is significant at 30 % significant level. Figure 29 and Figure 30 show that as the life expectancy increases, the dependency between water quality and quality of support network decrease for both the logit and the complementary model. Figure 31 and Figure 32 show the PP plot of the fitted theoretical copula vs. the empirical copula of the second regression part, the copula regression. Both graphs show near perfect alignment all through the course except of 6 observations which are a little bit far away from the diagonal.

5. Conclusions

The water quality and the quality of support network are linked together via Gumbel copula and they are strongly dependent to each other. They show very strong positive dependency, if one enlarges the other also enlarges, and this is logic as in the community where their residences have high expression about the good quality of the water supplies, they are more likely to support each other at the time of needs and crises. Marginally regressing the water quality on the life expectancy shows significant positive nonlinear relationship between the response and the predictor reflecting the community where their population has long life, this increase the satisfactions of the cleanliness of the water supplies. Also as the life expectancy increase, this increase the support the individuals in the community can have in times of need, but this effect of life expectancy on the quality of support network is not statistically significant as the effect of life expectancy on the water quality. But the dependency between water quality and quality of support weakens as the life expectancy increases. Using the Logit link function for the marginal median regression model shows better statistical indices than the log-log complementary link function. And this has a reflection on the second part, the parametric median copula regression part, where the indices are significant at 5% confidence level than the 30% confidence level attained while using the log-log complementary link function.

The significant effect of life expectancy on water quality can be interpreted as elderly have significant high expression of their satisfaction about the water quality and they do not show this significant high satisfaction about the quality of support network although it is present but not that significant. As the life expectancy does not have the same positive significant effect on both response variables, the dependency between the 2 response variables may weaken as the life expectancy increases. This can be interpreted as the percentage of people who can express their great satisfaction with water quality and the presence of support network are mainly the young (less life expectancy) while the elderly are the percentage who express less satisfaction with water quality and support network and this may be because they live alone so they are not receiving enough care for their water cleanliness and support so they have this impact by expressing less satisfaction for both response variables together. Also the source of water or the acquisition of the clean water all can make the expression of satisfaction vary in elderly. If they can buy the well-formed bottles of clean water or they have to check for the cleanliness for the household water supply.

This study introduces an innovative parametric copula regression framework that intricately captures the multifaceted relationship between two response variables while conditioning on a covariate. Departing from traditional regression methodologies, which often treat each response in isolation, this pioneering approach plunges into the depths of interdependencies that weave between outcomes, harnessing a dynamic parametric copula. The properties and structural nuances of this copula evolve as the predictor varies, enriching the analysis with its adaptability. Applying Sklar’s theorem, the complex joint conditional distribution is masterfully unraveled into a series of marginal regression models, seamlessly coupled with a dependence function signified by a Gumbel copula. In this framework, the copula parameter is elegantly modeled as a smooth function of the covariate, allowing for a dynamic interaction of dependence that transforms across different levels of the predictor. To estimate parameters, the inference function for margins (IFM) method was employed, ensuring robust computational stability along with impressive asymptotic efficiency. A comprehensive set of diagnostic checks based on residual analysis not only affirmed the adequacy of the fitted marginal models, but also validated the copula’s dependence structure, underscoring the reliability of the framework. Empirical findings have illuminated the model's exceptional ability to encapsulate the distinctive behaviors of marginal responses, as well as the evolving patterns of dependence between them. One striking observation revealed that the dependence parameter tended to diminish as the predictor values increased, indicating that the two outcomes became less tightly intertwined amid rising predictor levels. This flexible modeling architecture offers profound insights into how covariates influence not just individual responses but also their intricate interconnected dynamics. Such a feature proves particularly valuable in fields like finance, hydrology, and biostatistics, where grasping the co-movement between outcomes is crucial for informed decision-making and strategic planning.

6. Future Work

The proposed framework presents numerous opportunities for further development and expansion. Future research could explore the implementation of higher-dimensional copulas, which would allow for the modeling of interactions among more than two response variables concurrently. This could be particularly useful in complex systems where multiple factors are interrelated. Additionally, the use of semiparametric link functions could be investigated to enhance the flexibility in modeling the dependent relationships among variables. This approach would enable researchers to better accommodate various distributions and non-linear associations, providing a richer analysis of multivariate data. Another promising avenue for future work is the incorporation of time-varying copula parameters. This adaptation would be essential for capturing dynamic dependence structures in longitudinal or temporal datasets, where relationships between variables may evolve over time. By allowing the dependence parameters to change, the model could more accurately reflect real-world phenomena where relationships are not static. Overall, the parametric copula regression model stands out as a powerful and interpretable tool for investigating multivariate conditional dependence. It effectively bridges the traditional methods of classical regression with the advanced techniques of contemporary dependence modeling, thereby offering a comprehensive framework for analyzing the intricate relationships between multiple variables.