Copula-Based Regression with Mixed Covariates

Abstract

In this paper, we focused on developing copula-based modeling procedures that effectively capture the dependence between response and explanatory variables. Building upon the work of Noh et al. (J. Am. Stat. Assoc. 2013, 108, 676–688) we extended copula-based regression to accommodate both continuous and discrete covariates. Specifically, we explored the construction of copulas to estimate the conditional mean of the response variable given the covariates, elucidating the relationship between copula structures and marginal distributions. We considered various estimation methods for copulas and distribution functions, presenting a diverse array of estimators for the conditional mean function. These estimators range from non-parametric to semi-parametric and fully parametric, offering flexibility in modeling regression relationships. An adapted algorithm is applied to construct copulas and simulations are carried out to replicate datasets, estimate prediction model parameters, and compare with the OLS method. The practicality and efficacy of our proposed methodologies, grounded in the principles of copula-based regression, are substantiated through methodical simulation studies.

1. Introduction

Scientists are often led to study the relationships and dependencies between the response variable and several other covariates. However, regression analysis is the statistical tool for investigating such relationships and it is one of the most commonly used statistical methods in many scientific fields, such as medicine, biology, agriculture, economics, engineering, sociology, etc. In medical research, econometrics, and other research fields, it is very common to use regression analysis to interpret the correlation existing between different variables. However, the basic form of the regression analysis is not suitable for many cases, where the relationships are often non-linear and the probability distribution of the output variable may be an abnormal distribution.

For such dependence modeling problems, we attempt to provide a functional form that will summarize the relationship between response and explanatory variables. In several practical situations, as an example, a vector of covariates $\mathbf{X}=(X_1,\dots ,X_d)$ is used to explain, interpret, or predict the response variable Y. This is encountered in many fields, including medical fields and social science. The type of functional relationship we attempt to figure out could depend on the marginal behavior of variables or their joint behavior. In this paper, we consider the construction of dependent modeling procedures based on the separation of these two behaviors when the covariates are a mixture of continuous and discrete variables.

For this context, we consider procedures that allow the representation of a multivariate distribution as a function of its uni-variate marginals through a connection function called a copula. Copulas have been increasingly popular for modeling statistical dependence in multivariate sets of data and have been applied to various areas, including medical research, environmental science, econometrics, actuarial science, agronomy, and others. A key feature of copulas is that they provide flexible representations of the multivariate distribution by allowing for the dependence structure of the variables of interest to be modeled separately from the marginal structure and, by specifying a copula, we summarize all the dependencies between margins (see Nelsen [1] for more about this subject).

The power of this approach principally lies in the ability for a practitioner to model the dependence structure independently of the marginal behaviors. Furthermore, the advantages of using copulas in modeling are the allowance to model both linear and non-linear dependence, an arbitrary choice of a marginal distribution, and the capability of modeling extreme endpoints. However, the principal advantage of a copula regression is that there are no restrictions and no specification on the probability distributions that can be used.

It is interesting to note that copula-based regression models offer significant advantages in capturing complex dependencies between variables, making them highly useful in various fields. In finance, they allow for better portfolio risk management by modeling non-linear dependencies between asset returns and macroeconomic factors, especially during market downturns. In insurance, copula-based regression can be applied to explain pricing in terms of different dependent types of claims, such as frequency and severity. In environmental studies, regression as function of a copula is useful to establish the relationship between rainfall and river discharge, especially in the case of non-linear dependence. In healthcare, regression with a copula enables researchers to examine how lifestyle factors influence health outcomes, such as cholesterol levels, while capturing the potential interdependence among these health indicators.

In the literature, there exist many recent studies of regression based on copulas; as examples, we cite Sheikhi et al. [2] and Ali et al. [3] among others. As a new contribution to this domain, we consider in this paper the estimation problem of the mean regression function for a regression model, where $\mathbf{X}={(X_1,\dots ,X_d)}^T$ is a random vector of dimension $d\ge 1$ and Y is a random variable with cumulative distribution function (c.d.f.) $F_0$ and density function $f_0$. Y is the response variable and $\mathbf{X}$ is the set of covariates. We denote by $F_i$ the c.d.f. of the variables $X_i$ and we denote by $f_i$ its corresponding density. For a given $\mathbf{x}={(x_1,\dots ,x_d)}^T$, we will note by $\mathbf{F}\left(\mathbf{x}\right)$ the shortcut for $(F_1\left(x_1\right),\dots ,F_d\left(x_d\right))$. From the inspiring work of Sklar [4], the c.d.f. of ${(Y,{\mathbf{X}}^T)}^T$ evaluated at ${(y,{\mathbf{x}}^T)}^T$ can be expressed in terms of $C(F_0\left(y\right),\mathbf{F}\left(\mathbf{x}\right))$, where C is the copula distribution of ${(Y,{\mathbf{X}}^T)}^T$, that is, the function from ${[0,1]}^{d+1}$ to $[0,1]$ defined by

$$C(u_0,u_1,\dots ,u_d)=P\left(F_0\left(Y\right)\le u_0,F_1\left(X_1\right)\le u_1,\dots ,F_d\left(X_d\right)\le u_d\right).$$

Recently, Noh et al. [5] exploited the above decomposition to introduce a novel idea consisting of expressing the mean regression function $m\left(\mathbf{x}\right)$, in terms of the copula and margins of $\mathbf{x}$ as follows.

$$m\left(\mathbf{x}\right)={\displaystyle \frac{E\left[Yc(F_0\left(Y\right),\mathbf{F}\left(\mathbf{x}\right))\right]}{c_X\left(\mathbf{F}\left(\mathbf{x}\right)\right)}}.$$

(1)

where $c(u_0,\mathbf{v})\equiv c(u_0,v_1,\dots ,v_d)$ is the copula density corresponding to C and $c_{\mathbf{X}}\left(\mathbf{v}\right)\equiv c_{\mathbf{X}}(v_1,\dots ,v_d)$ is the copula density of $\mathbf{X}$. This shows that the mean regression function $m\left(\mathbf{x}\right)$ is the ratio of a numerator that only captures the mean dependence between Y and X and a denominator that captures the dependence within X. It is worth mentioning that the formula is only valid when the covariates are continuous. A new reformulation is needed when the covariates are not all continuous, which is the case for many real-world applications, especially in medicine.

Furthermore, Noh et al. [5] proposed a semi-parametric estimator for the regression function given in (1). Specifically, they utilized the inference function for margins (IFM) technique to estimate the copula-based regression curve. This method proceeds in two stages: first, it estimates the marginal parameters, and then it estimates the corresponding dependence parameter. These authors demonstrate, both theoretically and empirically, that the resulting estimates obtained exhibit desirable properties when the parametric copula family is adequately chosen.

Noh et al. [5] stimulated extensive research on copula-based regression. Noh et al. [6] applied the method of Noh et al. [5] to the quantile regression with i.i.d. or time series that are completely observed. De Backer et al. [7] extended the method of Noh et al. [6] to the quantile regression with censored data. Kraus and Czado [8] studied the quantile regression with complete data, using D-vine copulas. Rémillard et al. [9] discussed the asymptotic connection between the estimators of Noh et al. [6] and Kraus and Czado [8]. Chang and Joe [10] proposed an algorithm for computing the conditional distribution function via the vine copula. Furthermore, Nagler and Vatter [11] unified various copula-based regressions by formulating a general loss function which may not be continuously differentiable. Their generalized regression model includes the conditional mean regression of Noh et al. [5], the conditional quantile regression of Noh et al. [6], and the asymmetric least squares of Newey and Powell [12] as special cases. The unified framework enhances the systematic interpretation of the different existing regressions. For additional discussion into similar methods, see [13,14,15,16,17] and the literature cited therein.

As an extension of the framework by Noh et al. [5], we incorporate discrete variables into the set of covariates $X_1,\dots ,X_d$. By establishing a connection with various classes of copulas through an alternative equation to (1), we calculate the conditional mean, $m\left(\mathbf{x}\right)$, of $Y|\mathbf{X}=\mathbf{x}$. In this context, we develop the relationship between the copula and the marginals. Furthermore, we illustrate this relationship for specific families of copulas, such as Archimedean copulas and the Gaussian copula, highlighting their properties that are beneficial for our analysis.

The next step involved addressing the estimation problem. Here, we also adopt a semi-parametric approach along with the inference function for margins (IFM) method to estimate the proposed regression curve. First, we estimate the marginal distributions using their empirical distributions, and then we estimate the dependence parameter associated with the underlying copula. A simulation studies for different classes of copulas and different distributions for the output Y are considered to illustrate the usefulness of the findings.

The rest of the paper is organized as follows. Section 2, discusses different copula concepts in the multivariate setting. Section 3 outlines the copula-based regression model proposed for case where the set of covariates includes both discrete and continuous variables. Section 4 covers the estimation procedure of the proposed regression model. Section 5 is dedicated to a simulation study that assesses the performance of the suggested copula-based regression. Conclusion and remarks come in Section 6.

2. Preliminaries

Copulas are a mathematical concept used in multivariate analysis to describe the dependency structure between components of a multivariate random vector. They play a central role in various fields employing multivariate statistical analysis, such as risk management and finance. Therefore, copulas provide a framework for modeling the relationships between variables by describing their joint distribution independently of their marginal distributions.

This section provides a brief overview of the copula concept, which will be utilized in the development of the proposed model. According to Nelsen [1], a multivariate copula is defined as follows.

Definition 1.

A d-dimensional copula is a function from ${[0,1]}^d$ to $[0,1]$ with following properties:

1. 

For every $\mathbf{u}=(u_1,\dots ,u_d)\in {[0,1]}^d$,

(i) 

$C\left(\mathbf{u}\right)=0$ if at least one of coordinate of u is 0.

(ii) 

$C\left(\mathbf{u}\right)=u_k$ if all coordinates of u are 1 except $u_k$.

2. 

For every $\mathbf{a}=(a_1,\dots ,a_d)$ and $\mathbf{b}=(b_1,\dots ,b_d)$ in ${[0,1]}^d$ such that $a_i\le b_i$, $i=1,\dots ,d$,

$${\Delta }_{a_d}^{b_d}{\Delta }_{a_{d-1}}^{b_{d-1}}\cdots {\Delta }_{a_1}^{b_1}C\left(\mathbf{u}\right)\ge 0,\mathit{for}\mathit{all}\mathbf{u}\in {[0,1]}^d,$$

where,

$${\Delta }_{a_k}^{b_k}C\left(\mathbf{u}\right)=C(u_1,\dots ,u_{k-1},b_k,u_{k+1},\dots ,u_d)-C(u_1,\dots ,u_{k-1},a_k,u_{k+1},\dots ,u_d)$$

Sklar’s Theorem is a fundamental result in copula theory. It enables us to express the joint distribution of a multivariate random vector in terms of their marginal distributions and a copula function. It can be stated as follows (see, Nelsen [1]).

Theorem 1.

Let H be a d-dimensional distribution function with marginal distributions $F_1,F_2,\dots ,F_d$. Then, there exists a d-copula C such that for all $(x_1,\dots ,x_d)\in {\mathbb{R}}^n$,

$$H(x_1,\dots ,x_d)=C(F_1\left(x_1\right),F_2\left(x_2\right),\dots ,F_d\left(x_d\right)).$$

If $F_1,F_2,\dots ,F_d$ are all continuous, then C is unique; otherwise, C is uniquely determined on $\mathit{Ran}\left(F_1\right)\times \cdots \times \mathit{Ran}\left(F_d\right)$. Conversely, if C is a d-copula and $F_1,F_2,\dots ,F_d$ are distribution functions, then the function H defined by is a d-dimensional distribution function with margins $F_1,F_2,\dots ,F_d$.

It is well-known that Sklar’s Theorem has numerous practical applications in various fields involving multivariate data analysis. For instance, Sklar’s Theorem is commonly employed to analyze dependencies among different financial assets. It enables us to understand how the dependence structure between the prices of different assets might affect the overall risk of a portfolio.

2.1. Archimedean Copulas

Archimedean copulas constitute an important class of parametric copulas. This type of copula describes the dependence structure between random variables with greater flexibility through a single function called the generator. The latter is often expressed in terms of dependence parameters that control the strength of dependence among the components of a given random vector.

The generator of a d-dimensional Archimedean copula is an increasing and continuous function $\varphi$ defined from $[0,\infty )$ to $[0,1]$ such that $\varphi \left(0\right)=0$ and $\varphi (\infty )=1$. Suppose that $\varphi$ is differentiable up to the order $d-2$ with derivatives noted by ${\varphi }^{\left(i\right)}$ for $i=1,\dots ,d-2$ and let $\psi$ be its inverse, that is, $\psi ={\varphi }^{-1}$. Hereafter is the definition of a multivariate Archimedean copula. For details on this subject, see McNeil and Nešlehová [18].

Definition 2.

The d-dimensional Archimedean copula $C_{\varphi }$ is defined through its generator ϕ as follows:

$$C_{\varphi }(u_1,\dots ,u_d)=\varphi (\psi \left(u_1\right)+\cdots +\psi \left(u_d\right)),\mathit{for}\mathit{all}(u_1,\dots ,u_d)\in {[0,1]}^d,$$

where the generator ϕ is subject to the conditions that ${(-1)}^i{\varphi }^{\left(i\right)}\left(x\right)\ge 0$ for $i=1,\dots ,d-2$, and ${(-1)}^{d-2}{\varphi }^{(d-2)}$ is non-increasing and convex.

Hereafter, we present the Clayton copula and the Frank copula, both considered among the most popular multivariate Archimedean copulas. The d-dimensional Clayton copula is defined as follows:

$$C_{\varphi }(u_1,\dots ,u_d)={\left(u_1^{-\theta }+\cdots +u_d^{-\theta }-d+1\right)}^{-\frac{1}{\theta }},\mathrm{for}\mathrm{all}(u_1,\dots ,u_d)\in {[0,1]}^d$$

(2)

It is an Archimedean copula whose inverse generator function is defined, for all $t>0$, by

$${\psi }_{\theta }\left(t\right)={\left(1+\theta t\right)}^{-\frac{1}{\theta }},\mathrm{for}\mathrm{all}\theta \in \left(0,\infty \right).$$

(3)

Likewise, the d-dimensional Frank family copula is expressed, for all $(u_1,\dots ,u_d)\in {[0,1]}^d$, by:

$$C_{\varphi }(u_1,\dots ,u_d)=-{\displaystyle \frac{1}{\theta }}ln\left(1+{\displaystyle \frac{(e^{-\theta u_1}-1)(e^{-\theta u_2}-1)\cdots (e^{-\theta u_d}-1)}{{(e^{-\theta }-1)}^{d-1}}}\right),\theta >0.$$

(4)

Its inverse generator function is given, for all $t>0$, by

$${\psi }_{\theta }\left(t\right)=-ln\left({\displaystyle \frac{e^{-\theta t}-1}{e^{-\theta }-1}}\right)\mathrm{for}\mathrm{all}(u_1,\dots ,u_d)\in {[0,1]}^d),\theta >0.$$

(5)

For $d\ge 3$, the Frank copula describes only the positive dependence, whereas in the two-dimensional case, this copula models both positive and negative association.

2.2. Gaussian Copulas

The d-Gaussian copula $C_{\Sigma }$ is defined through the standardized d-multivariate normal distribution ${\Phi }_{\Sigma }$. The correlation matrix $\Sigma$ represents the dependence parameters of this copula. Specifically, $C_{\Sigma }$ is expressed by

$$C_{\Sigma }(u_1,\dots ,u_d)={\Phi }_{\Sigma }\left({\Phi }^{-1}\left(u_1\right),\dots ,{\Phi }^{-1}\left(u_d\right)\right)\mathrm{for}\mathrm{all}(u_1,\dots ,u_d)\in {[0,1]}^d),$$

where $\Phi$ denotes the standard normal distribution. In other words, the multivariate Gaussian copula is explicitly given by,

$$C_{\Sigma }(u_1,\dots ,u_d)={\int }_{-\infty }^{{\Phi }^{-1}\left(u_1\right)}\dots {\int }_{-\infty }^{{\Phi }^{-1}\left(u_d\right)}{\displaystyle \frac{1}{2\pi \det (\Sigma )}}exp\left\{-{\displaystyle \frac{1}{2}}\left({\mathbf{x}}^{\top }({\Sigma }^{-1}-I)\mathbf{x}\right)\right\}d\mathbf{x},$$

where I is the identical matrix. The bivariate Gaussian copula is reduced to

$$C_{\rho }(u_1,u_2)={\int }_{-\infty }^{{\Phi }^{-1}\left(u_1\right)}{\int }_{-\infty }^{{\Phi }^{-1}\left(u_2\right)}{\displaystyle \frac{1}{2\pi \sqrt{1-{\rho }^2}}}exp\left\{-{\displaystyle \frac{x^2-2\rho x+y^2}{2(1-{\rho }^2)}}\right\}dxdy,$$

where $\rho$ represents the Pearson correlation coefficient, a parameter within the range $[-1,1]$, serving as the dependence parameter for this copula.

3. Model Description

Starting from a random vector $(Y,\mathbf{X})$, where Y is a continuous random variable with cumulative distribution function $F_0$, assume that $\mathbf{X}=\left({\mathbf{X}}_q,{\mathbf{X}}_{d-q}\right)$, where the random vectors ${\mathbf{X}}_q=(X_1,\dots ,X_q)$ and ${\mathbf{X}}_{d-q}=(X_{q+1},\dots ,X_d)$ are continuous and discrete, respectively, and suppose without loss of generality that, for any $i=q+1,\dots ,d$, $\mathrm{Ran}\left(X_i\right)\subset \mathbb{N}=\{0,1,2,\dots \}$. Denote by $F_1,\dots ,F_d$ the distribution functions of $X_1,\dots ,X_d$, respectively, and, for all ${\mathbf{x}}_d\in {\mathbb{R}}^q\times {\Pi }_{q+1}^d\mathrm{Ran}\left(X_i\right)$, let ${\mathbf{F}}_d\left(\mathbf{x}\right)=\left(F_1\left(x_1\right),\dots ,F_d\left(x_d\right)\right)$. Let C and $C_{\mathbf{X}}$ be the copulas of $(Y,\mathbf{X})$ and $\mathbf{X}$, respectively. Moreover, for $(u,{\mathbf{v}}_d)=(u,v_1,\dots ,v_d)\in {(0,1)}^{d+1}$, set,

$$\begin{array}{c}\hfill {\partial }_{q+1}C(u,{\mathbf{v}}_d)={\displaystyle \frac{{\partial }^{q+1}}{\partial u\partial v_1\cdots \partial v_q}}C(u,{\mathbf{v}}_d)\mathrm{and}{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{v}}_d\right)={\displaystyle \frac{{\partial }^q}{\partial v_1\cdots \partial v_q}}{C}_{\mathbf{X}}\left({\mathbf{v}}_d\right).\end{array}$$

(6)

For $i=q+1,\dots ,d$, set ${\delta }_i=({\delta }_{i,1},\dots ,{\delta }_{i,d})$ such that ${\delta }_{i,i}=1$ and ${\delta }_{i,j}=0$, for $i\ne j$. For $i=q+1,\dots ,d$, let ${\Delta }_i$ be the forward difference operator defined by

$$\begin{array}{ccc}\hfill {\Delta }_i{\partial }_{q+1}C(F_0\left(k\right),{\mathbf{F}}_d\left({\mathbf{x}}_d\right))& =& {\partial }_{q+1}C(F_0\left(k\right),{\mathbf{F}}_d\left({\mathbf{x}}_d\right))-{\partial }_{q+1}C(F_0\left(k\right),{\mathbf{F}}_d({\mathbf{x}}_d-{\delta }_i)),\hfill \\ \hfill {\Delta }_i{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left({\mathbf{x}}_d\right)\right)& =& {\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left({\mathbf{x}}_d\right)\right)-{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d({\mathbf{x}}_d-{\delta }_i)\right).\hfill \end{array}$$

and set ${\Delta }_q^d={\prod }_{i=q+1}^d{\Delta }_i$.

Proposition 1.

For all $\mathbf{x}\in {\mathbb{R}}^q\times {\Pi }_{q+1}^d\mathit{Ran}\left(X_i\right)$, the conditional mean, of Y given $\mathbf{X}=\mathbf{x}$ is expressed by,

$$\begin{array}{c}\hfill m\left(\mathbf{x}\right)={\displaystyle \frac{E\left\{Y{\Delta }_q^d{\partial }_{q+1}C(F_0\left(Y\right),{\mathbf{F}}_d\left(\mathbf{x}\right))\right\}}{{\Delta }_q^d{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left(\mathbf{x}\right)\right)}}.\end{array}$$

(7)

Proof. 

For all $\mathbf{x}\in {\mathbb{R}}^q\times {\Pi }_{q+1}^d\mathrm{Ran}\left(X_i\right)$, let $f\left(.|\mathbf{x}\right)$ be the conditional density function of Y given $\mathbf{X}=\mathbf{x}$. Clearly, one has

$$\begin{array}{ccc}\hfill m\left(\mathbf{x}\right)& =& E\left(Y|\mathbf{X}=\mathbf{x}\right)\hfill \\ & =& \int yf\left(y|\mathbf{x}\right)dy\hfill \\ & =& \int {\displaystyle \frac{y{\Delta }_q^d{\partial }_{q+1}C\left(F_0\left(y\right),{\mathbf{F}}_d\left(\mathbf{x}\right)\right)f_0\left(y\right)f_1\left(x_1\right)\cdots f_q\left(x_q\right)}{{\Delta }_q^d{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left(\mathbf{x}\right)\right)f_1\left(x_1\right)\cdots f_q\left(x_q\right)}}dy\hfill \\ & =& \int {\displaystyle \frac{y{\Delta }_q^d{\partial }_{q+1}C\left(F_0\left(y\right),{\mathbf{F}}_d\left(\mathbf{x}\right)\right)f_0\left(y\right)}{{\Delta }_q^d{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left(\mathbf{x}\right)\right)}}dy\hfill \\ & =& {\displaystyle \frac{E\left\{Y{\Delta }_q^d{\partial }_{q+1}C(F_0\left(Y\right),{\mathbf{F}}_d\left(\mathbf{x}\right)\right\}}{{\Delta }_q^d{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left(\mathbf{x}\right)\right)}}.\hfill \end{array}$$

□

Remark 1.

A developed expression of the conditional mean can be obtained by using the expanded expression of ${\Delta }_q^d$ described by

$$\begin{array}{c}\hfill {\Delta }_q^d=1+\sum _{S\subset \{1,\dots ,d-q\}}{(-1)}^{\left|S\right|}\prod _{i\in S}{T}_i,\end{array}$$

(8)

where $T_i=1-{\Delta }_i$ and $\left|S\right|$ denotes the cardinality number of any nonempty subset S of $\{1,\dots ,d-q\}$. Therefore, one sees from (7) and (8) that $m\left(\mathbf{x}\right)$ can expressed by,

$$\begin{array}{c}\hfill {\displaystyle \frac{E\left\{Y{\partial }_{q+1}C(F_0\left(Y\right),{\mathbf{F}}_d\left({\mathbf{x}}_d\right)\right\}+{\sum }_{S\subset \{1,\dots ,d-q\}}{(-1)}^{\left|S\right|}E\left\{Y{\prod }_{i\in S}{T}_i{\partial }_{q+1}C(F_0\left(Y\right),{\mathbf{F}}_d\left({\mathbf{x}}_d\right)\right\}}{{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left(\mathbf{x}\right)\right)+{\sum }_{S\subset \{1,\dots ,d-q\}}{(-1)}^{\left|S\right|}{\prod }_{i\in S}{T}_i{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left(\mathbf{x}\right)\right)}},\end{array}$$

(9)

where

$$\begin{array}{ccc}\hfill T_i{\partial }_{q+1}C(F_0\left(k\right),{\mathbf{F}}_d\left({\mathbf{x}}_d\right))& =& {\partial }_{q+1}C(F_0\left(k\right),{\mathbf{F}}_d({\mathbf{x}}_d-{\delta }_i)),\hfill \\ \hfill T_i{\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d\left({\mathbf{x}}_d\right)\right)& =& {\partial }_q{C}_{\mathbf{X}}\left({\mathbf{F}}_d({\mathbf{x}}_d-{\delta }_i)\right).\hfill \end{array}$$

3.1. Archimedean Copula-Based Predicted Mean

Suppose that the dependence structure of $(Y,\mathbf{X})$ is described by an Archimedean class of copulas C with generator $\varphi$. This means the copulas C and $C_{\mathbf{X}}$ are expressed, for all $(u,{\mathbf{v}}_d)\in {(0,1)}^{d+1}$, by

$$C(u,{\mathbf{v}}_d)=\psi \left(\varphi \left(u\right)+\sum _{j=1}^d\varphi \left(v_j\right)\right)\mathrm{and}{C}_{\mathbf{X}}\left({\mathbf{v}}_d\right)=\psi \left(\sum _{j=1}^d\varphi \left(v_j\right)\right),$$

where the function $\psi$ represents the inverse of the generator $\varphi$. Therefore, the partial derivative of the copulas C and $C_{\mathbf{X}}$ are given by

$${\partial }_{q+1}C(u,{\mathbf{v}}_d)={\psi }^{(q+1)}\left(\varphi \left(u\right)+\sum _{j=1}^d\varphi \left(v_j\right)\right){\varphi }^{\prime }\left(u\right)\prod _{j=1}^q{\varphi }^{\prime }\left(v_j\right),$$

and

$${\partial }_q{C}_{\mathbf{X}}\left({\mathbf{v}}_d\right)={\psi }^{\left(q\right)}\left(\sum _{j=1}^d\varphi \left(v_j\right)\right)\prod _{j=1}^q{\varphi }^{\prime }\left(v_j\right).$$

Hence, the regression curve is given by

$$\begin{array}{c}\hfill m\left(\mathbf{x}\right)={\displaystyle \frac{E\left\{Y{\varphi }^{\prime }\left(F_0\left(Y\right)\right){\Delta }_q^d{\psi }^{(q+1)}\left(\varphi \left(F_0\left(Y\right)\right)+{\sum }_{j=1}^d\varphi \left(F_j(x_j\right)\right)\right\}}{{\Delta }_q^d{\psi }^{\left(q\right)}\left({\sum }_{j=1}^d\varphi \left(F_j(x_j\right)\right)}}.\end{array}$$

(10)

To exemplify the above conditional mean, let us examine the scenario where $d=2$ and $q=1$, implying that covariate $X_1$ is continuous, while covariate $X_2$ is discrete. In such instances, we have:

$$\begin{array}{ccc}\hfill m\left(\mathbf{x}\right)& =& {\displaystyle \frac{E\left\{Y{\varphi }^{\prime }\left(F_0\left(Y\right)\right){\Delta }_1^2{\psi }^{″}\left\{\varphi \left(F_0\left(Y\right)\right)+\varphi \left(F_1\left(x_1\right)\right)+\varphi \left(F_2\left(x_2\right)\right)\right\}\right\}}{{\Delta }_1^2{\psi }^{\prime }\left\{\varphi \left(F_1\left(x_1\right)\right)+\varphi \left(F_2\left(x_2\right)\right)\right\}}},\hfill \end{array}$$

(11)

where

$$\begin{array}{ccc}\hfill {\Delta }_1^2{\psi }^{″}\left\{\varphi (F_0\left(Y\right)+\varphi \left(F_1\left(x_1\right)\right)+\varphi \left(F_2\left(x_2\right)\right)\right\}& =& {\psi }^{″}\left\{\varphi \left(F_0\left(Y\right)\right)+\varphi \left(F_1\left(x_1\right)\right)+\varphi \left(F_2\left(x_2\right)\right)\right\}\hfill \\ & & -{\psi }^{″}\left\{\varphi \left(F_0\left(Y\right)\right)+\varphi \left(F_1\left(x_1\right)\right)+\varphi \left(F_2(x_2-1)\right)\right\},\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill {\Delta }_1^2{\psi }^{\prime }\left\{\varphi \left(F_1\left(x_1\right)\right)+\varphi \left(F_2\left(x_2\right)\right)\right\}& =& {\psi }^{\prime }\left\{\varphi \left(F_1\left(x_1\right)\right)+\varphi \left(F_2\left(x_2\right)\right)\right\}-{\psi }^{\prime }\left\{\varphi \left(F_1\left(x_1\right)\right)+\varphi \left(F_2(x_2-1)\right)\right\}.\hfill \end{array}$$

Example 1.

Illustrating Equation (11), let us assume that $(Y,X_1,X_2)$ follows the Clayton copula $C_{\theta }$ as described in (2). Specifically, for all $\theta \in (0,\infty )$,

$$C_{\theta }(u,v_1,v_2)={\left(u^{-\theta }+v_1^{-\theta }+v_2^{-\theta }-2\right)}^{-\frac{1}{\theta }},\mathit{for}\mathit{all}(u,v_1,v_2)\in {[0,1]}^3.$$

The generator ${\varphi }_{\theta }$ of this copula and its inverse ${\psi }_{\theta }$ given in (3) satisfy

$${\psi }_{\theta }\left(t\right)={\displaystyle \frac{t^{-\theta }-1}{\theta }},{\psi }_{\theta }^{\prime }\left(t\right)=-{\left(1+\theta t\right)}^{-\frac{1}{\theta }-1},{\psi }_{\theta }^{″}\left(t\right)=(\theta +1){\left(1+\theta t\right)}^{-\frac{1}{\theta }-2}.$$

Hence, standard calculations show that (11) reduces to

$$\begin{array}{ccc}\hfill m\left(\mathbf{x}\right)& =& {\displaystyle \frac{(1+\theta )E\left\{Y{F}_0{\left(Y\right)}^{-\theta -1}{\Delta }_1^2{\left(F_0{\left(Y\right)}^{-\theta }+F_1{\left(x_1\right)}^{-\theta }+F_2{\left(x_2\right)}^{-\theta }-2\right)}^{-\frac{1}{\theta }-2}\right\}}{{\Delta }_1^2{\left(F_1{\left(x_1\right)}^{-\theta }+F_2{\left(x_2\right)}^{-\theta }-1\right)}^{-\frac{1}{\theta }-1}}}.\hfill \end{array}$$

(12)

Likewise, let us express Equation (11) when $(Y,X_1,X_2)$ follows the Frank copula $C_{\theta }$ expressed in (4), namely,

$$C_{\varphi }(u,v_1,v_2)=-{\displaystyle \frac{1}{\theta }}ln\left(1+{\displaystyle \frac{(e^{-\theta u}-1)(e^{-\theta v_1}-1)(e^{-\theta v_2}-1)}{{(e^{-\theta }-1)}^2}}\right),\theta \in \mathbb{R}\setminus \left\{0\right\}.$$

Calculations similar to those used previously lead to

$$\begin{array}{ccc}\hfill m\left(\mathbf{x}\right)& =& {\displaystyle \frac{E\left\{Y{\Delta }_1^2{K}_{1,\theta }(F_0\left(Y\right),F_1\left(x_1\right),F_2\left(x_2\right))\right\}}{{\Delta }_1^2{K}_{2,\theta }(F_1\left(x_1\right),F_2\left(x_2\right))}},\hfill \end{array}$$

(13)

where

$$\begin{array}{c}\hfill K_{1,\theta }(u,v_1,v_2)={\displaystyle \frac{e^{-\theta u}{e}^{-\theta v_1}(e^{-\theta v_2}-1)\left[{(1-e^{-\theta })}^2+e^{-\theta u}(e^{-\theta v_1}-1)(e^{-\theta v_2}-1)\right]}{{\left[{(1-e^{-\theta })}^2+(e^{-\theta u}+1)(e^{-\theta v_1}-1)(e^{-\theta v_2}-1)\right]}^2}},\end{array}$$

(14)

and

$$\begin{array}{c}\hfill K_{2,\theta }(v_1,v_2)={\displaystyle \frac{e^{-\theta v_1}(e^{-\theta v_2}-1)}{1-e^{-\theta }+(e^{-\theta v_1}-1)(e^{-\theta v_2}-1)}}.\end{array}$$

(15)

3.2. Gaussian Copula-Based Predicted Mean

This section presents the expression of the regression curve when the copula C of $(Y,\mathbf{X})$ is Gaussian. This means that the copula C is expressed in terms of the standardized $(d+1)$-multivariate normal distribution ${\Phi }_{\Sigma }$ and the correlation matrix $\Sigma$, which is assumed to be non-singular, as follows.

$$C(u,\mathbf{v};\Sigma )={\Phi }_{\Sigma }\left({\Phi }^{-1}\left(u\right),{\Phi }^{-1}\left(\mathbf{v}\right)\right)\forall (u,\mathbf{v})\in {[0,1]}^{d+1},$$

where $\Phi$ denotes the standard normal distribution and where we note that

$${\Phi }^{-1}\left(\mathbf{v}\right)=({\Phi }^{-1}\left(v_1\right),{\Phi }^{-1}\left(v_2\right),\dots ,{\Phi }^{-1}\left(v_d\right)).$$

To derive ${\partial }_{q+1}C(u,{\mathbf{v}}_d)$ and ${\partial }_q{C}_{\mathbf{X}}\left({\mathbf{v}}_d\right)$, let us first decompose the correlation matrix $\Sigma$ as follows.

$$\Sigma =\left(\begin{array}{cc}{\Sigma }_{11}& {\Sigma }_{12}\\ {\Sigma }_{21}& {\Sigma }_{22}\end{array}\right),\mathrm{with}\mathrm{sizes}\left(\begin{array}{cc}(q+1)\times (q+1)& (q+1)\times (d-q)\\ (d-q)\times (q+1)& (d-q)\times (d-q)\end{array}\right),$$

where ${\Sigma }_{11}$ and ${\Sigma }_{22}$ represent the correlation matrices of the $(q+1)$-continuous random vector $(Y,{\mathbf{X}}_q)$ and the $(d-q)$-discrete random vector ${\mathbf{X}}_{d-q}$, respectively. Furthermore, ${\Sigma }_{12}$ denotes the correlation matrix between the random vectors ${\mathbf{X}}_{d-q}$ and $(Y,{\mathbf{X}}_q)$ and ${\Sigma }_{21}={{\Sigma }_{12}}^{\top }$.

Consider the $(d+1)$-uniform random vector $(U,V_1,\dots ,V_d)$ with distribution C, and set ${\mathbf{V}}_q=(V_1,\dots ,V_q)$ and ${\mathbf{V}}_{d-q}=(V_{q+1},\dots ,V_d)$. Therefore, one observes, for all $(u,{\mathbf{v}}_q,{\mathbf{v}}_{d-q})\in (0,1)\times {(0,1)}^q\times {(0,1)}^{d-q}$,

$$\begin{array}{c}\hfill {\partial }_{q+1}C(u,{\mathbf{v}}_d)=\mathrm{P}\left({\mathbf{V}}_{d-q}\le {\mathbf{v}}_{d-q}\mid U=u,{\mathbf{V}}_q={\mathbf{v}}_q\right)c_{q+1}(u,{\mathbf{v}}_q),\end{array}$$

where $c_{q+1}$ is the copula density of $(U,{\mathbf{V}}_q)$. Let ${\Phi }^{-1}\left({\mathbf{V}}_{d-q}\right)=({\Phi }^{-1}\left(V_1\right),\dots ,{\Phi }^{-1}\left(V_{d-q}\right))$ and ${\Phi }^{-1}(u,{\mathbf{v}}_q)=({\Phi }^{-1}\left(u\right),{\Phi }^{-1}\left(v_1\right),\dots ,{\Phi }^{-1}\left(q\right))$. Since $({\Phi }^{-1}\left(U\right),{\Phi }^{-1}\left(V_1\right),\dots ,{\Phi }^{-1}\left(V_d\right))$ is normally distributed, then the conditional random vector ${\Phi }^{-1}\left({\mathbf{V}}_{d-q}\right)\mid U=u,{\mathbf{V}}_q={\mathbf{v}}_q$ is distributed as

$$N\left({\Sigma }_{21}{\Sigma }_{11}^{-1}{\Phi }^{-1}(u,{\mathbf{v}}_q),{\Sigma }_{22}-{\Sigma }_{21}{\Sigma }_{11}^{-1}{\Sigma }_{21}^{\top }\right),$$

with distribution function $G_{q+1}(·\mid {\Phi }^{-1}(u,{\mathbf{v}}_q))$. It follows that

$$\begin{array}{ccc}\hfill {\partial }_{q+1}C(u,{\mathbf{v}}_d)& =& G_{q+1}({\Phi }^{-1}\left({\mathbf{V}}_{d-q}\right)\mid {\Phi }^{-1}(u,{\mathbf{v}}_q))c_{q+1}(u,{\mathbf{v}}_q),\hfill \end{array}$$

where $c_{q+1}(u,{\mathbf{v}}_q)$ is the copula density of the random vector $(U,{\mathbf{V}}_q)$. It remains to derive ${\partial }_q{C}_{\mathbf{X}}\left({\mathbf{v}}_d\right)$, where $C_{\mathbf{X}}$ is the Gaussian copula $\mathbf{X}$ with correlation matrix

$$\tilde{\Sigma }=\left(\begin{array}{cc}{\tilde{\Sigma }}_{11}& {\tilde{\Sigma }}_{12}\\ {\tilde{\Sigma }}_{21}& {\Sigma }_{22}\end{array}\right),\mathrm{with}\mathrm{sizes}\left(\begin{array}{cc}q\times q& q\times (d-q)\\ (d-q)\times q& (d-q)\times (d-q)\end{array}\right),$$

where ${\tilde{\Sigma }}_{11}$ and ${\Sigma }_{22}$ represent the correlation matrices of the q-continuous random vector ${\mathbf{X}}_q$ and the $(d-q)$-discrete random vector ${\mathbf{X}}_{d-q}$, respectively. Similarly, ${\tilde{\Sigma }}_{12}$ denotes the correlation matrix between the random vectors ${\mathbf{X}}_{d-q}$ and ${\mathbf{X}}_q$ and ${\tilde{\Sigma }}_{21}={\tilde{\Sigma }}_{12}^{\top }$. It follows that, for all $({\mathbf{v}}_q,{\mathbf{v}}_{d-q})\in \times {(0,1)}^q\times {(0,1)}^{d-q}$,

$$\begin{array}{c}\hfill {\partial }_q{C}_{\mathbf{X}}\left({\mathbf{v}}_d\right)=\mathrm{P}\left({\mathbf{V}}_{d-q}\le {\mathbf{v}}_{d-q}\mid ,{\mathbf{V}}_q={\mathbf{v}}_q\right)c_q\left({\mathbf{v}}_q\right),\end{array}$$

where $c_q$ is the copula density of the random vector ${\mathbf{V}}_q$. Using the fact that the conditional random vector ${\Phi }^{-1}\left({\mathbf{V}}_{d-q}\right)\mid {\mathbf{V}}_q={\mathbf{v}}_q$ is distributed as

$$N\left({\tilde{\Sigma }}_{21}{\tilde{\Sigma }}_{11}^{-1}{\Phi }^{-1}\left({\mathbf{v}}_q\right),{\Sigma }_{22}-{\tilde{\Sigma }}_{21}{\tilde{\Sigma }}_{11}^{-1}{\tilde{\Sigma }}_{21}^{\top }\right),$$

with distribution function $G_q(·\mid {\Phi }^{-1}\left({\mathbf{v}}_q\right))$. It follows that

$$\begin{array}{c}\hfill {\partial }_q{C}_{\mathbf{X}}\left({\mathbf{v}}_d\right)=G_q({\Phi }^{-1}\left({\mathbf{V}}_{d-q}\right)\mid {\Phi }^{-1}\left({\mathbf{v}}_q\right))c_q\left({\mathbf{v}}_q\right),\end{array}$$

Therefore, the predicted mean is given by

$$\begin{array}{c}\hfill m\left(\mathbf{x}\right)={\displaystyle \frac{E\left\{Y{c}_{q+1}(F_0\left(Y\right),{\mathbf{F}}_q\left({\mathbf{x}}_q\right)){\Delta }_q^d{G}_{q+1}({\Phi }^{-1}\left({\mathbf{F}}_{d-q}\left({\mathbf{x}}_{d-q}\right)\right)\mid {\Phi }^{-1}(F_0\left(Y\right),{\mathbf{F}}_q\left({\mathbf{x}}_q\right))\right\}}{c_q\left({\mathbf{F}}_q\left({\mathbf{x}}_q\right)\right){\Delta }_q^d{G}_q({\Phi }^{-1}\left({\mathbf{F}}_{d-q}\left({\mathbf{x}}_{d-q}\right)\right)\mid {\Phi }^{-1}\left({\mathbf{F}}_q\left({\mathbf{x}}_q\right)\right)}}.\end{array}$$

Example 2.

Consider the case $d=2$ and $q=1$. This means that the covariate $X_1$ is continuous and the covariate $X_2$ is discrete. Assume further that the copula of $(Y,X_1,X_2)$ is Gaussian with the correlation matrix

$$\Sigma =\left(\begin{array}{ccc}1& {\rho }_{12}& {\rho }_{13}\\ {\rho }_{12}& 1& {\rho }_{23}\\ {\rho }_{13}& {\rho }_{23}& 1\end{array}\right).$$

In such a case, we have,

$$\begin{array}{c}\hfill m\left(\mathbf{x}\right)={\displaystyle \frac{E\left\{Y{c}_2(F_0\left(Y\right),F_1\left(x_1\right)){\Delta }_1^2{G}_2({\Phi }^{-1}\left(F_2\left(x_2\right)\right)\mid {\Phi }^{-1}(F_0\left(Y\right),F_1\left(x_1\right))\right\}}{c_1\left(F_1\left(x_1\right)\right){\Delta }_1^2{G}_1({\Phi }^{-1}\left(F_2\left(x_2\right)\right)\mid {\Phi }^{-1}\left(F_1\left(x_1\right)\right))}},\end{array}$$

where $c_1\left(F_1\left(x_1\right)\right)=1$. It remains to calculate the copula density $c_2$ of $(U,V_1)$ and the conditional normal distribution $G_1(·\mid {\Phi }^{-1}\left(F_1\left(x_1\right)\right)$ and $G_2(·\mid {\Phi }^{-1}(F_0\left(Y\right),F_1\left(x_1\right))$. Since the copula of $(U,V_1)$ is Gaussian with correlation matrix, ${\Sigma }_{11}=\left(\begin{array}{cc}1& {\rho }_{12}\\ {\rho }_{12}& 1\end{array}\right).$ Then, the copula density $c_2(F_0\left(Y\right),F_1\left(x_1\right))$ is given by

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\sqrt{1-{\rho }_{12}^2}}}exp\left\{{\displaystyle \frac{2{\rho }_{12}{\Phi }^{-1}\left(F_0\left(Y\right)\right){\Phi }^{-1}\left(F_1\left(x_1\right)\right)-{\rho }_{12}^2\left[{\Phi }^{-1}{\left(F_0\left(Y\right)\right)}^2+{\Phi }^{-1}{\left(F_1\left(x_1\right)\right)}^2\right]}{2(1-{\rho }_{12}^2)}}\right\}.\end{array}$$

Also, we have

$${\Sigma }_{12}={\Sigma }_{21}^{\top }=({\rho }_{13},{\rho }_{23}),{\Sigma }_{22}=1,{\tilde{\Sigma }}_{11}=1,{\tilde{\Sigma }}_{12}={\tilde{\Sigma }}_{21}={\rho }_{23}.$$

Standard calculations show that, from this, $G_2(·\mid {\Phi }^{-1}(F_0\left(Y\right),F_1\left(x_1\right)))$ is the distribution of $N(\mu (Y,x_1),{\sigma }^2)$, such that

$$\begin{array}{c}\hfill \mu (Y,x_1)={\displaystyle \frac{({\rho }_{13}-{\rho }_{12}{\rho }_{23}){\Phi }^{-1}\left(F_0\left(Y\right)\right)+({\rho }_{23}-{\rho }_{12}{\rho }_{13}){\Phi }^{-1}\left(F_1\left(x_1\right)\right)}{1-{\rho }_{12}^2}},\end{array}$$

and

$$\begin{array}{c}\hfill {\sigma }^2={\displaystyle \frac{1-{\rho }_{12}^2-{\rho }_{13}^2-{\rho }_{23}^2+2{\rho }_{12}{\rho }_{13}{\rho }_{23}}{1-{\rho }_{12}^2}}.\end{array}$$

Likewise, $G_1(·\mid {\Phi }^{-1}\left(F_1\left(x_1\right)\right)$ is the distribution of $N(\tilde{\mu }\left(x_1\right),{\tilde{\sigma }}^2)$, such that

$$\begin{array}{c}\hfill \tilde{\mu }\left(x_1\right)={\rho }_{23}{\Phi }^{-1}\left(F_1\left(x_1\right)\right)\mathit{and}{\tilde{\sigma }}^2=1-{\rho }_{23}^2.\end{array}$$

$$\begin{array}{c}\hfill m\left(\mathbf{x}\right)={\displaystyle \frac{E\left\{Y{c}_2(F_0\left(Y\right),F_1\left(x_1\right)){\Delta }_1^2\Phi \left({\sigma }^{-1}[{\Phi }^{-1}(F_2\left(x_2\right)-\mu (Y,x_1)]\right)\right\}}{{\Delta }_1^2\Phi \left({\tilde{\sigma }}^{-1}\left[{\Phi }^{-1}\right(F_2\left(x_2\right)-\tilde{\mu }\left(x_1\right)\right)\left]\right)}},\end{array}$$

where

$$\begin{array}{ccc}\hfill {\Delta }_1^2\Phi \left({\sigma }^{-1}[{\Phi }^{-1}(F_2\left(x_2\right)-\mu (Y,x_1)]\right)& =& \Phi \left({\sigma }^{-1}[{\Phi }^{-1}(F_2\left(x_2\right)-\mu (Y,x_1)]\right)\hfill \\ & & -\Phi \left({\sigma }^{-1}[{\Phi }^{-1}(F_2(x_2-1)-\mu (Y,x_1)]\right),\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill {\Delta }_1^2\Phi \left({\tilde{\sigma }}^{-1}[{\Phi }^{-1}(F_2\left(x_2\right)-\tilde{\mu }\left(x_1\right)]\right)& =& \Phi \left({\tilde{\sigma }}^{-1}[{\Phi }^{-1}(F_2\left(x_2\right)-\tilde{\mu }\left(x_1\right)]\right)\hfill \\ & & -\Phi \left({\tilde{\sigma }}^{-1}[{\Phi }^{-1}(F_2(x_2-1)-\tilde{\mu }\left(x_1\right)]\right).\hfill \end{array}$$

Example 3.

As a continuation of Example 2, in order to give a closed form for the conditional mean $m\left(y\right|x)$, we consider the case where the variables Y and $X_1$ are distributed as standard normal and where the correlation matrix of the Gaussian copula is determined by ${\rho }_{12}=0$ and ${\rho }_{13}={\rho }_{23}=0.6$. Thus, the conditional expectation is given by:

$$\begin{array}{c}\hfill m\left(\mathbf{x}\right)={\displaystyle \frac{E\left[Y\left\{\Phi \left(1.89[{\Phi }^{-1}(F_2\left(x_2\right)-0.6Y-0.6x_1]\right)-\Phi \left(1.89[{\Phi }^{-1}(F_2(x_2-1)-0.6Y-0.6x_1]\right)\right\}\right]}{\Phi \left(1.25\left[{\Phi }^{-1}(F_2\left(x_2\right)-0.6x_1)\right]\right)-\Phi \left(1.25\left[{\Phi }^{-1}(F_2(x_2-1)-0.6x_1)\right]\right)}},\end{array}$$

for any discrete random variables $X_2$ and for any $(y,x_1,x_2)\in {\mathbb{R}}^2\times \mathit{Ran}\left(X_2\right)$.

4. Estimation

Consider a sample of n observations $(Y_1,{\mathbf{X}}_1),\dots ,(Y_n,{\mathbf{X}}_n)$ from the random vector $(Y,\mathbf{X})$. For $i=1,\dots ,n$, denote ${\mathbf{X}}_i=(X_{i1},\dots ,X_{id})$. To estimate the conditional mean described in (7), we need to estimate the marginal distributions $F_0,F_1,\dots ,F_d$ as well as the partial derivatives of the copulas C and $C_{\mathbf{X}}$, namely, ${\partial }_{q+1}C$ and ${\partial }_q{C}_{\mathbf{X}}$, provided in (6). In this paper, we use a semi-parametric methodology that first consists of estimating margins $F_0,F_1,\dots ,F_d$ through their rescaled empirical distributions given by

$$F_{n0}\left(y\right)={\displaystyle \frac{1}{n+1}}\sum _{i=1}^n\mathrm{I}(Y\le y)\mathrm{and}{F}_{nj}\left(x\right)={\displaystyle \frac{1}{n+1}}\sum _{i=1}^n\mathrm{I}(X_{ij}\le x),j=1,\dots ,d,$$

respectively, where ${\mathrm{I}}_A$ stands for the indicator function for given event A. An alternative method for estimating these quantities is through the use of the kernel smoothing technique. Typically, this approach yields more accurate estimations compared to the method based on empirical distributions. The idea behind this method is to estimate the distributions $F_0,F_1,\dots ,F_d$ using

$${\widehat{F}}_0\left(y\right)={\displaystyle \frac{1}{n+1}}\sum _{i=1}^{n}K\left({\displaystyle \frac{y-Y_i}{h}}\right)\mathrm{and}{\widehat{F}}_j\left(x\right)={\displaystyle \frac{1}{n+1}}\sum _{i=1}^{n}K\left({\displaystyle \frac{x-X_{ij}}{h}}\right),j=1,\dots ,d.$$

$K(·)$ represents a non-negative function known as the kernel, while h signifies the bandwidth. It is well known that the selection of the bandwidth is crucial and significantly influences the accuracy of the estimation.

The second step is to estimate parametrically the copula of $(Y,\mathbf{X})$. To this end, assume that the copula C is identified as a member of some parametric family, $\mathcal{F}=\left\{C_{\theta },\mathit{\theta }\in \Theta \right\}$, where $\Theta \subset \mathbb{R}.$ This means that there exists ${\mathit{\theta }}_0\in \Theta$ such that $C=C_{{\theta }_0}$ and $C_{\mathbf{X}}=C_{\mathbf{X},{\mathit{\theta }}_{\mathbf{0}}}$. Therefore, the copula C is then estimated by $C_{\widehat{\theta }}$, where $\widehat{\theta }$ is an estimator of ${\theta }_0$. This estimator is typically obtained by maximizing, in terms of $\theta$, the expressed pseudo-likelihood function,

$$L\left(\theta \right)=\sum _{i=1}^{n}ln\left\{c_{\theta }\left({\widehat{F}}_0\left(Y_i\right),{\widehat{\mathbf{F}}}_d\left({\mathbf{X}}_i\right)\right)\right\}$$

where ${\widehat{\mathbf{F}}}_d\left({\mathbf{X}}_i\right)=({\widehat{F}}_1\left(X_{i1}\right),\dots ,{\widehat{F}}_d\left(X_{id}\right))$. In other words, the estimator of the $\theta$ is given by

$$\widehat{\theta }{=}_{\theta \in \Theta }L\left(\theta \right).$$

Finally, the conditional mean $m\left(\mathbf{x}\right)$ is estimated using (7) as follows:

$$\begin{array}{c}\hfill \widehat{m}\left(\mathbf{x}\right)={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{Y_i{\Delta }_q^d{\partial }_{q+1}{C}_{\widehat{\theta }}({\widehat{F}}_0\left(Y_i\right),{\widehat{\mathbf{F}}}_d\left(\mathbf{x}\right))}{{\Delta }_q^d{\partial }_q{C}_{\mathbf{X},\widehat{\mathit{\theta }}}\left({\widehat{\mathbf{F}}}_d\left(\mathbf{x}\right)\right)}}.\end{array}$$

(16)

Example 4.

Let us examine the above estimation procedure in the scenario where $d=2$ and $q=1$, the situation entails $X_1$ being a continuous covariate, while $X_2$ is discrete. Additionally, let us assume that the copula governing $(Y,X_1,X_2)$ is the Clayton copula with parameter $\theta \in (0,\infty )$. The theoretical conditional mean $\widehat{m}\left(\mathbf{x}\right)$ is provided in (12). Its estimated couterpart can be derived from (19) as follows:

$$\begin{array}{ccc}\hfill \widehat{m}\left(\mathbf{x}\right)& =& {\displaystyle \frac{(1+\widehat{\theta }){\sum }_{i=1}^n{Y}_i{\widehat{F}}_0{\left(Y_i\right)}^{-\widehat{\theta }-1}{\Delta }_1^2{\left({\widehat{F}}_0{\left(Y_i\right)}^{-\widehat{\theta }}+{\widehat{F}}_1{\left(x_1\right)}^{-\widehat{\theta }}+{\widehat{F}}_2{\left(x_2\right)}^{-\widehat{\theta }}-2\right)}^{-\frac{1}{\widehat{\theta }}-2}}{n{\Delta }_1^2{\left({\widehat{F}}_1{\left(x_1\right)}^{-\widehat{\theta }}+{\widehat{F}}_2{\left(x_2\right)}^{-\widehat{\theta }}-1\right)}^{-\frac{1}{\widehat{\theta }}-1}}}.\hfill \end{array}$$

(17)

The estimators $\widehat{\theta }$, ${\widehat{F}}_0$, ${\widehat{F}}_1$, and ${\widehat{F}}_2$ can be computed using a sample $(Y_i,X_{1i},X_{2i})$, $i=1,\dots ,n$, selected from the distribution of $(Y,X_1,X_2)$. Similarly, in the case where the dependence structure of $(Y,X_1,X_2)$ is modeled by a Frank copula, the estimated conditional mean $\widehat{m}\left(\mathbf{x}\right)$ can be derived from (13) and (19) as follows:

$$\begin{array}{ccc}\hfill \widehat{m}\left(\mathbf{x}\right)& =& {\displaystyle \frac{1}{{\Delta }_1^2{K}_{2,\widehat{\theta }}({\widehat{F}}_1\left(x_1\right),{\widehat{F}}_2\left(x_2\right))}}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{Y}_i{\Delta }_1^2{K}_{1,\widehat{\theta }}({\widehat{F}}_0\left(Y_i\right),{\widehat{F}}_1\left(x_1\right),{\widehat{F}}_2\left(x_2\right)),\hfill \end{array}$$

(18)

$K_{1,\theta }$ and $K_{2,\theta }$ are given in (14) and (15), respectively.

5. Simulation Study

The objective of this section is to conduct simulations to compare the proposed conditional mean estimator with some competitors. To achieve this, we focus on the case where $d=2$ with mixed covariates; specifically, $X_1$ is continuous, and $X_2$ is discrete. In this case, the proposed estimator is deduced from its general form expressed in (19) as follows,

$$\begin{array}{c}\hfill \widehat{m}(x_1,x_2)={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\displaystyle \frac{Y_i\left[{\partial }_2{C}_{\widehat{\theta }}({\widehat{F}}_0\left(Y_i\right),{\widehat{F}}_1\left(x_1\right),{\widehat{F}}_2\left(x_2\right))-{\partial }_2{C}_{\widehat{\theta }}({\widehat{F}}_0\left(Y_i\right),{\widehat{F}}_1\left(x_1\right),{\widehat{F}}_2(x_2-1))\right]}{{\partial }_1{C}_{\mathbf{X},\widehat{\theta }}\left({\widehat{F}}_1\left(x_1\right),{\widehat{F}}_2\left(x_2\right)\right)-{\partial }_1{C}_{\mathbf{X},\widehat{\theta }}\left({\widehat{F}}_1\left(x_1\right),{\widehat{F}}_2(x_2-1)\right)}},\end{array}$$

(19)

where

$${\partial }_2{C}_{\widehat{\theta }}(u,v_1,v_2)={\displaystyle \frac{{\partial }^2{C}_{\widehat{\theta }}(u,v_1,v_2)}{\partial u\partial v_1}}\mathrm{and}{\partial }_1{C}_{\mathbf{X},\widehat{\theta }}(v_1,v_2)={\displaystyle \frac{\partial C_{\mathbf{X},\widehat{\theta }}(v_1,v_2)}{\partial v_1}}.$$

As scenarios, we consider the most common cases to show the improvement of our estimator over the OLS estimator. However, for the copula of $(Y,X_1,X_2)$, we consider Clayton, Frank, and Gumbel with parameter $\theta \in (0,\infty )$ and for the variables $Y\sim N(\mu ,{\sigma }^2)$ or $Y\sim t\left(df\right)$, while $X_1\sim {\mathcal{U}}_{[a,b]}$ and $X_2\in \{1,2,3\}$ with distribution $F_2\left(1\right)=p_1$, $F_2\left(2\right)=p_1+p_2$ and $F_2\left(3\right)=1$. The generalized inverse of $F_2$ is

$$F_2^{-1}\left(t\right)=inf\left\{v\in \left\{1,2,3\right\}:F_2\left(v\right)>t\right\},$$

or equivalently,

$$F_2^{-1}\left(t\right)=\mathbb{I}(0\le t<p_1)+2\mathbb{I}(p_1\le t<p_1+p_2)+3\mathbb{I}(p_1+p_2\le t\le 1).$$

Simulation algorithm:

• Given $n,a,b,p_1,p_2,\mu ,{\sigma }^2$, and $\theta$.
• For $i=1,\dots ,n$.
• Generate $u_i,v_{i,1},v_{i,2}$ from a copula $C_{\theta }$.
• Set $y_i=F_0^{-1}\left(u_i\right)$, $x_{i,1}=F_1^{-1}\left(v_{i,1}\right)$ and $x_{i,2}=F_2^{-1}\left(v_{i,2}\right)$.
• Use the generated sample $(y_i,x_{i,1},x_{i,2})$, $i=1,\dots ,n$ to estimate $\theta$ and define the empirical distributions of $F_0$, $F_1$, and $F_2$.
• Evaluate the estimator $\widehat{m}(x_1,x_2)$ for $(x_1,x_2)$ belonging to the grid defined by

  $$F=E\times \{1,2,3\}\mathrm{where}E=\{i/K,i=1,\dots ,K\}.$$

For fixed $(x_1,x_2)\in F$, we first compute the theoretical value $m(x_1,x_2)$ and then evaluate $\widehat{m}(x_1,x_2)$ using J random samples of size n. We denote the corresponding estimates by ${\widehat{m}}^{\left(j\right)}(x_1,x_2)$, where $j=1,\dots ,J$. To assess the performance, we employ the empirical integrated mean squared error (IMSE), which is formulated as follows:

$$\mathrm{IMSE}={\displaystyle \frac{1}{J}}\sum _{j=1}^J\left[{\displaystyle \frac{1}{\left|F\right|}}\sum _{(x_1,x_2)\in F}{({\widehat{m}}^{\left(j\right)}(x_1,x_2)-m(x_1,x_2))}^2\right],$$

(20)

where $\left|F\right|$ denotes the cardinality number of the grid F. Notably, $\mathrm{IMSE}$ can be decomposed into the square of empirical bias, ${\mathrm{IBIAS}}^2$, and the empirical variance, $\mathrm{IVAR}$, as follows:

$$\mathrm{IMSE}={\mathrm{IBIAS}}^2+\mathrm{IVAR}.$$

In this simulation study, different values of the parameters are considered, which represent different dependence scenarios ranging from weak to strong, with Kendall’s tau, $\tau$, values lying in the interval $(0.2,0.75)$. With a sample size $n=200$, the response, Y, is generated from $N(0,1)$ distribution and Student’s t-distribution with 3 degrees of freedom. Also, $X_1$ is generated from a Uniform(0, 1) and $X_2\in \{1,2,3\}$ with distribution $F_2\left(1\right)=p_1$, $F_2\left(2\right)=p_1+p_2$ and $F_2\left(3\right)=1$, where $p_1=0.3$ and $p_2=0.5$. In this context, we report and compare the integrated mean square error (IMSE) and the integrated mean absolute error (IMAE) with the respective errors derived from the least squares (ls) regression method. This comprehensive approach ensured the reliability of the comparison by accounting for variability in outcomes across multiple realizations. The reported values in

Table 1. Simulation results for normal distribution with $n=200$.

		IMSE		IMAE
Copula	Parameters	${\displaystyle \widehat{{\mathit{m}}_{\mathit{c}}}}$	${\displaystyle \widehat{{\mathit{m}}_{\mathit{ls}}}}$	${\displaystyle \widehat{{\mathit{m}}_{\mathit{c}}}}$	${\displaystyle \widehat{{\mathit{m}}_{\mathit{ls}}}}$
Clayton	$\theta =0.50(\tau =0.20)$	0.005	0.015	0.084	0.071
	$\theta =2.00(\tau =0.50)$	0.012	0.033	0.116	0.104
	$\theta =6.00(\tau =0.75)$	0.032	0.095	0.147	0.136
Frank	$\theta =2.37(\tau =0.25)$	0.00026	0.120	0.011	0.283
	$\theta =5.73(\tau =0.50)$	0.00004	0.305	0.004	0.459
	$\theta =14.14(\tau =0.75)$	0.00001	0.613	0.002	0.670
Gumbel	$\theta =1.25(\tau =0.20)$	0.020	0.050	0.105	0.151
	$\theta =2.00(\tau =0.50)$	0.036	0.095	0.134	0.208
	$\theta =4.00(\tau =0.75)$	0.062	0.232	0.144	0.285
Gausian	${\rho }_{12}=0.4,{\rho }_{13}=0.4,{\rho }_{23}=-0.4$	0.079	0.080	0.236	0.234
	$(\tau =0.26)$
	${\rho }_{12}=0.9,{\rho }_{13}=0.9,{\rho }_{23}=0.85$	0.051	0.071	0.172	0.210
	$(\tau =0.69)$

and

Table 2. Simulation results for Student’s t with $df=3$ and $n=200$.

		IMSE		IMAE
Copula	Parameters	${\displaystyle \widehat{{\mathit{m}}_{\mathit{c}}}}$	${\displaystyle \widehat{{\mathit{m}}_{\mathit{ls}}}}$	${\displaystyle \widehat{{\mathit{m}}_{\mathit{c}}}}$	${\displaystyle \widehat{{\mathit{m}}_{\mathit{ls}}}}$
Clayton	$\theta =0.50(\tau =0.20)$	0.011	0.044	0.121	0.121
	$\theta =2.00(\tau =0.50)$	0.022	0.079	0.157	0.158
	$\theta =6.00(\tau =0.75)$	0.061	0.291	0.206	0.247
Frank	$\theta =2.37(\tau =0.25)$	0.00112	0.270	0.023	0.420
	$\theta =5.73(\tau =0.50)$	0.00024	0.673	0.010	0.679
	$\theta =14.14(\tau =0.75)$	0.00016	1.556	0.008	1.059
Gumbel	$\theta =1.25(\tau =0.20)$	0.051	0.138	0.157	0.246
	$\theta =2.00(\tau =0.50)$	0.076	0.264	0.199	0.352
	$\theta =4.00(\tau =0.75)$	0.120	0.709	0.194	0.576
Gausian	${\rho }_{12}=0.4,{\rho }_{13}=0.4,{\rho }_{23}=-0.4$	0.400	0.525	0.431	0.448
	$(\tau =0.26)$
	${\rho }_{12}=0.9,{\rho }_{13}=0.9,{\rho }_{23}=0.85$	0.087	0.301	0.220	0.435
	$(\tau =0.69)$

, corresponding to normal distribution and Student’s t-distribution, respectively, represent the averages calculated from a total of 100 realizations. The results show that the proposed method consistently outperformed the least squares regression method across all the scenarios. This dominance was evident in both metrics, IMSE and IMAE, and across all varieties of Kendall’s tau values and sample sizes. We also analyzed the evolution of MSE with sample size, confirming a clear reduction as n grows, improving estimator accuracy and stability (see

Table 3. Simulation results for normal distribution with $n=50$ and $n=100$.

IMSE		$\mathit{n}\mathbf{=}\mathbf{50}$		$\mathit{n}\mathbf{=}\mathbf{100}$
Copula	Parameters	$\widehat{{\mathit{m}}_{\mathit{c}}}$	$\widehat{{\mathit{m}}_{\mathit{ls}}}$	$\widehat{{\mathit{m}}_{\mathit{c}}}$	$\widehat{{\mathit{m}}_{\mathit{ls}}}$
Clayton	$\theta =0.50(\tau =0.20)$	0.031	0.079	0.014	0.042
	$\theta =2.00(\tau =0.50)$	0.047	0.103	0.031	0.080
	$\theta =6.00(\tau =0.75)$	0.078	0.257	0.070	0.199
Frank	$\theta =2.37(\tau =0.25)$	0.001	0.159	0.0006	0.138
	$\theta =5.73(\tau =0.50)$	0.00017	0.371	0.00009	0.306
	$\theta =14.14(\tau =0.75)$	0.00006	0.661	0.00002	0.598
Gumbel	$\theta =1.25(\tau =0.20)$	0.027	0.088	0.026	0.067
	$\theta =2.00(\tau =0.50)$	0.059	0.127	0.043	0.103
	$\theta =4.00(\tau =0.75)$	0.085	0.269	0.064	0.237
Gausian	${\rho }_{12}=0.4,{\rho }_{13}=0.4,{\rho }_{23}=-0.4$	0.127	0.113	0.101	0.095
	$(\tau =0.26)$
	${\rho }_{12}=0.9,{\rho }_{13}=0.9,{\rho }_{23}=0.85$	0.362	0.091	0.101	0.078
	$(\tau =0.69)$

). Specifically, we considered $n=50$ and $n=100$, which are relatively small. As n increases, the estimator improves significantly in terms of IMSE.

Particularly, the proposed method revealed a more accurate and robust performance, indicating a lower IMSE and IMAE between the estimated and actual values than the least squares method. This enhanced performance can be attributed to the proposed method’s ability to more effectively capture and account for the underlying correlation structure represented by Kendall’s tau in the data. Unlike the least squares method, which assumes a specific form of relationship (linear), the proposed method offers a more flexible and robust approach to analyzing data with varying degrees of correlation and complexity.

6. Conclusions

This paper extends the copula-based regression model introduced by Noh et al. [5] by addressing the scenario where covariates are mixed, encompassing both continuous and discrete explanatory variables. Unlike the original model, which dealt exclusively with continuous covariates, the proposed approach broadens the applicability of the copula-based regression framework. The parameter estimation has been performed using the inference function for margins (IFM), which first estimates the marginal parameters and then estimates the corresponding dependence parameter. Through detailed examples, we demonstrated the estimation of the proposed regression equation and conducted a comprehensive simulation study under various scenarios involving different types of copulas. The results of the simulation study indicate that the suggested model performs favorably compared to classical regression approaches, showcasing its potential to handle mixed-covariate data effectively. This extension provides a valuable contribution to the field of regression analysis, offering a new regression tool for researchers and practitioners dealing with diverse explanatory data types. An interesting potential research direction involves extending this concept to regression with multivariate responses using the same mixed covariates. This extension is particularly relevant in various practical applications where multiple outcomes need to be modeled simultaneously. For instance, in environmental studies, multivariate regression can be used to assess how industrial emissions simultaneously impact both air and water quality, accounting for the complex interactions between pollutants. In healthcare, it enables researchers to examine how lifestyle factors influence multiple health outcomes, such as blood pressure, cholesterol levels, and blood sugar levels.