Modeling the Ratio of Correlated Biomarkers Using Copula Regression

Abstract

Modeling the ratio of two dependent components as a function of covariates is a frequently pursued objective in observational research. Despite the high relevance of this topic in medical studies, where biomarker ratios are often used as surrogate endpoints for specific diseases, existing models are based on oversimplified assumptions, assuming e.g. independence or strictly positive associations between the components. In this paper, we close this gap in the literature and propose a regression model where the marginal distributions of the two components are linked by Frank copula. A key feature of our model is that it allows for both positive and negative correlations between the components, with one of the model parameters being directly interpretable in terms of Kendall’s rank correlation coefficient. We study our method theoretically, evaluate finite sample properties in a simulation study and demonstrate its efficacy in an application to diagnosis of Alzheimer’s disease via ratios of amyloid-beta and total tau protein biomarkers.

Keywords: Distributional regression, Frank copula, Gamma distribution, Negative dependence, Ratio outcome

1 Introduction

A common objective in observational research is to analyze the ratio of two (possibly dependent) components $U,V\in\mathbb{R}^{+}$ (Long et al., 2016). Typical examples are, among others, (i) the LDL/HDL cholesterol ratio in cardiovascular research (Natarajan et al., 2003), defined as the ratio of the low-density lipoprotein and the high-density lipoprotein concentrations in plasma or serum, (ii) the CD4/CD8 ratio in HIV research (Caby et al., 2016), which measures the ratio of T helper cells to cytotoxic T cells in the human immune system, (iii) the testosterone over epitestosterone (T/E) ratio in antidoping research (Sottas et al., 2007), and (iv) the GEFC/REFC ratio in ophthalmic research, corresponding to the green and red emission components in fundus autofluorescence imaging (Wintergerst et al., 2022). In many of such studies, biomarker ratios are used as early indicators or even as surrogate endpoints for a specific disease. In these cases, the focus is not only on the characterization of the marginal ratio distribution, but also on modeling this distribution as a function of a set of covariates $X=(X_{1},\ldots,X_{p})^{\top}$. Usually, this amounts to specifying a regression model that includes the ratio as outcome variable (Berger et al., 2019).

When setting up a model relating the ratio outcome $R=U/V$ to the covariates $X$, a common assumption is that both components follow either log-normal or gamma distributions, thereby accounting for the nonnegativity of the component values and the skewness of their distributions (Mitchell et al., 2015; Van Domelen et al., 2021). In the former case it is easily derived that the ratio is again log-normally distributed. The latter case, which will be dealt with in this paper, is considerably less straightforward but is often preferred in practice due to its increased efficiency (Firth, 1988; Wiens, 1999; Berger et al., 2019).

In the special case where $U$ and $V$ are independently gamma distributed, the ratio $R=U/V$ follows a generalized beta distribution of the second kind, in the following abbreviated by GB2 (Kleiber and Kotz, 2003). A regression approach for the GB2 distribution has been proposed by Tulupyev et al. (2013), who studied the determinants of alcohol abuse in HIV-positive persons using the framework of vector generalized additive models (VGAMs; Yee, 2015). More recent examples include Safari-Katesari and Zaroudi (2020), Bourguignon et al. (2021), Medeiros et al. (2023) and dos Santos et al. (2023). The case of correlated gamma distributed components has earlier been studied by Lee et al. (1979) and Tubbs (1986). Based on Kibble’s bivariate gamma distribution for $(U,V)$ (Kibble, 1941), Berger et al. (2019) developed the extended GB2 (eGB2) model for the ratio of two positively correlated gamma distributed components. Their model is characterized by three parameters, of which one is directly interpretable in terms of the Pearson correlation coefficient between the two components. Conceptually, the extended GB2 model can be seen as a distributional regression model embedded in the framework of generalized additive models for location, scale and shape (GAMLSS; Rigby and Stasinopoulos, 2005; Berger and Schmid, 2020).

Despite its major importance in biostatistics, no regression modeling strategy exists (to the best of our knowledge) for ratio outcomes with two negatively correlated gamma distributed components. Negatively correlated measurements are encountered in numerous applications, for example in dementia research, where ratios of cerebrospinal fluid (CSF) biomarkers are used for the early diagnosis of Alzheimer’s disease (Koyama et al., 2012). Importantly, measurements of the widely employed amyloid-$\beta$ 42 protein and total tau protein biomarkers are known to exhibit a negative correlation (Tapiola et al., 2009). In recent publications, the Gaussian regression model has been used for modeling ratios of CSF biomarkers (e.g., Xu et al., 2020). Clearly, this model neither accounts for the characteristics of the bivariate distribution of $(U,V)$ nor for the skewness in the distribution of the ratio outcome $R$.

Motivated by these problems, and to address the current shortcomings in modeling ratio outcomes with negatively correlated gamma distributed components, we propose a regression model where the joint bivariate distribution of the two gamma distributed components is defined by Frank copula (Genest, 1987). By this choice, the model flexibly accounts for either negative or positive associations between the two components (measured by Spearman’s or Kendall’s rank correlation coefficient). It also allows for modeling different characteristics of the two marginal distributions, including possibly unequal rate and shape parameters. By relating the covariates $X$ to the parameters of the marginals, as well as to the association parameter defined by the copula, our model further allows to derive the conditional probability density function of $R\,|\,X$ as a function of covariates. This, in turn, allows for the analysis of conditional distributional parameters (like the expected value, median or quantiles), including valid inferential conclusions for these quantities.

We apply the new approach to data from a multi-center observational cohort study conducted by the German Dementia Network (DCN; Kornhuber et al., 2009). Study participants were diagnosed with either mild cognitive impairment (MCI), Alzheimer’s Disease (AD), or other dementia. The study aims at determining the diagnostic and prognostic power of clinical, laboratory and imaging methods. This task is considered to be a major challenge, as the period from the first clinical symptoms of AD to disease onset might take years to decades (Sperling et al., 2013). Consequently, as biomarker ratios like the amyloid-$\beta$ 42/total tau ratio are considered to be strong predictors of AD progression, it is of high interest to relate these measurements to patient characteristics like age, sex and educational level (Jack et al., 2015). As will be demonstrated in Section 4, the proposed copula regression model can be suitably applied to address this problem, resulting in meaningful descriptive and inferential findings regarding the associations between the biomarker ratio and individual patient characteristics.

The rest of the paper is organized as follows: Section 2 derives the distributional copula regression model, states novel theoretical results with implications for the interpretation of covariate effects, and presents estimation, prediction and inference. The efficacy of our approach is demonstrated empirically in a simulation study in Section 3 and in our main application to AD progression in Section 4. The main findings of the paper are discussed in Section 5.

2 Methods

Section 2.1 starts by deriving the distribution of the ratio of two gamma distributed components with dependence induced by a Frank copula. Details on model specification and fitting are given in Section 2.2. Section 2.3 covers the prediction of distributional parameters and inference.

2.1 Distributional concept

Let $U$ and $V$ be two gamma distributed random variables with probability density functions (PDF)

$$f_{U}(u)=\frac{\lambda_{U}^{\delta_{U}}}{\Gamma(\delta_{U})}u^{\delta_{U}-1}\exp(-\lambda_{U}u)\ \ \ \text{and}\ \ \ f_{V}(v)=\frac{\lambda_{V}^{\delta_{V}}}{\Gamma(\delta_{V})}v^{\delta_{V}-1}\exp(-\lambda_{V}v)\,,$$ (1)

where $\lambda_{U}$, $\lambda_{V}$ $>0$ denote the rate parameters and $\delta_{U}$, $\delta_{V}$ $>0$ denote the shape parameters of $f_{U}$ and $f_{V}$, respectively. We allow for dependencies between $U$ and $V$, and assume that their joint distribution can be described by Frank copula with copula function $C_{\theta}$. By Sklar’s theorem, the joint distribution of $(U,V)$ is thus given by

$$F_{U,V}(u,v) =C_{\theta}\left(F_{U}(u),F_{V}(v)\right)$$

$$=-\frac{1}{\theta}\,\log\left(1+\frac{(\exp(-\theta F_{U}(u))-1)(\exp(-\theta F_{V}(v))-1)}{\exp(-\theta)-1}\right)\,,$$ (2)

where $F_{U,V}$, $F_{U}$ and $F_{V}$ denote the joint bivariate and marginal cumulative distribution functions (CDFs) of $U$ and $V$, respectively. The parameter $\theta\in\mathbb{R}\setminus\{0\}$ determines the association between $U$ and $V$. It can be shown that Kendall’s rank correlation coefficient $\tau\in[-1,1]$ is a monotone increasing function of $\theta$, given by

$$\tau(\theta)=1\,+\,\frac{4}{\theta}\left(\frac{1}{\theta}\int_{0}^{\theta}\frac{t}{e^{t}-1}\,dt-1\right)$$ (3)

(Harry, 2014). As a consequence, the CDF in (2.1) allows for (possibly highly) positive or negative correlations between the two components $U$ and $V$. The joint PDF of $(U,V)$ is given by

$$f_{U,V}(u,v) =\frac{\partial^{2}}{\partial u\partial v}{F_{U,V}(u,v)}=c_{\theta}\big{(}F_{U}(u),F_{V}(v)\big{)}\,f_{U}(u)\,f_{V}(v)$$

$$=\frac{-\theta\,\exp(-\theta F_{U}(u))\,\exp(-\theta F_{V}(v))\,(\exp(-\theta)-1)\,f_{U}(u)\,f_{V}(v)}{((\exp(-\theta)-1)+(\exp(-\theta F_{U}(u))-1)\,(\exp(-\theta F_{V}(v))-1))^{2}}\,,$$ (4)

where $c_{\theta}(a,b):=\partial^{2}/(\partial a\partial b)\,C_{\theta}(a,b)$ is the PDF of Frank copula.

We derive the resulting PDF of the ratio $R$, an interpretable representation thereof and the CDF in the following three propositions.

Proposition 1.

Let the PDF of $(U,V)$ be defined by (2.1). Then the PDF of the ratio $R:=U/V$ is given by

$$f_{R}(r)=\int_{0}^{1} {\left|F_{V}^{-1}(s)\right|\,c_{\theta}\left(F_{U}(r\,F_{V}^{-1}(s)),s\right)\,f_{U}(rF_{V}^{-1}(s))}\,ds$$

$$=\int_{0}^{1} \frac{\exp(-\theta\,F_{U}(r\,F_{V}^{-1}(s)))\,\exp(-\theta s)(-\theta)(\exp(-\theta)-1)}{((\exp(-\theta)-1)+(\exp(-\theta\,F_{U}(r\,F_{V}^{-1}(s)))-1)(\exp(-\theta s)-1))^{2}}$$

$$\times F_{V}^{-1}(s)\,f_{U}(r\,F_{V}^{-1}(s))\,ds\,,$$ (5)

where $|\cdot|$ denotes the absolute value function.

Proof.

Proposition 1 can be derived from Proposition 1 of Ly et al. (2019), who provided analytical results for the PDF of the quotient $U/V$ of two random variables whose dependence structure can be described by an absolutely continuous copula. ∎

Figure 1: Examples of the PDF of $R=U/V$ derived in Proposition 1 for parameters $\lambda_{U},\lambda_{V}\in\{1,2\}$, $\delta_{U},\delta_{V}\in\{2,3\}$ and $\theta\in\{-10,1,10\}$ (corresponding to rank correlation coefficients $\tau\in\{-0.67,0.11,0.67\}$). In each panel the three lines refer to $\theta=-10$ (dotted), $\theta=1$ (solid) and $\theta=10$ (dashed). Vertical lines refer to the median values of $R$.

Remark 1.

Figure 1 visualizes the PDF of $R$ for different values of the rate, shape and association parameters. The figure illustrates that the form of the PDF is strongly related to the ratio of marginal means $\mathbb{E}(U)/\mathbb{E}(V)=\lambda_{V}\delta_{U}/\lambda_{U}\delta_{V}$, which is highest in the lower left panel ($\mathbb{E}(U)/\mathbb{E}(V)=3$) where the dispersion is very large, and lowest in the upper right panel ($\mathbb{E}(U)/\mathbb{E}(V)=1/3$) where the PDFs are heavily right-skewed. Figure 1 also describes the association between the PDF and the Kendall’s rank correlation coefficient. In each of the nine cases the mode of the PDF increases as $\theta$ increases. Of note, the median of $R$ does not vary with $\theta$ (being equal for the three PDFs in each panel).

Proposition 2.

Let the PDF of $(U,V)$ be defined by (2.1). Then the PDF of the random variable $R$ in (1) can be re-written as

$$f_{R}(r)=\int_{0}^{1} {c_{\theta}\left(\frac{1}{\Gamma(\delta_{U})}\,\gamma\left(\delta_{U},r\,\Lambda\,\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})s\right)\right),s\right)}$$

$$\times\,\frac{\Lambda^{\delta_{U}}\,r^{(\delta_{U}-1)}}{\Gamma(\delta_{U})}\left(\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})s\right)\right)^{\delta_{U}}\,\exp\left(-r\,\Lambda\,\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})s\right)\right)\,ds\,,$$ (6)

where $\Lambda:=\lambda_{U}/\lambda_{V}$ denotes the ratio of the two rate parameters and $\gamma(\cdot,\cdot)$ is the lower incomplete gamma function.

Proof. The proof of Proposition 2 is given in Appendix A.

Remark 2.

By Proposition 2 the PDF of $R$ can be written as a function of the ratio of the rate parameters $\Lambda=\lambda_{U}/\lambda_{V}$. This facilitates the interpretation of the proposed regression model, as it results in a sparser representation of the covariate effects. In particular, covariate effects can be investigated using one-dimensional hypothesis tests and p-values, see our application in Section 4. Figure 2 illustrates how the median of $R$ is related to $\Lambda$. It is seen that the median decreases with $\Lambda$ independent of the values of the shape parameters and the association parameter.

Proposition 3.

Let the CDF of $(U,V)$ be defined by (2.1). Then the CDF of the random variable $R$ is given by

$$F_{R}(r)=\int_{0}^{1}\,{\frac{A\,f_{U}(r\,F_{V}^{-1}(s))\,\frac{r}{f_{V}(F_{V}^{-1}(s))}\,(\exp(-\theta s)-1)+(A-1)\exp(-\theta s)}{\exp(-\theta)+(A-1)\exp(-\theta s)-A}}\,ds\,,$$ (7)

where $A=\exp\left(-\theta\,F_{U}(rF_{V}^{-1}(s))\right)$.

Proof.

By Equation (9) of Ly et al. (2019), the CDF of $R$ is derived as

$$F_{R}(r)\, \,=\underbrace{F_{V}(0)}_{=0}\,+\,\int_{0}^{1}{\underbrace{\text{sgn}(F_{V}^{-1}(s))}_{=1}\,\frac{\partial}{\partial s}C_{\theta}\left(F_{U}(r\,F_{V}^{-1}(s)),s\right)}\,ds$$

$$\,=-\frac{1}{\theta}\,\int_{0}^{1}{\frac{\partial}{\partial s}\log\left(1+\frac{\big{(}\exp(-\theta\,F_{U}(rF_{V}^{-1}(s)))-1\big{)}\,(\exp(-\theta s)-1)}{\exp(-\theta)-1}\right)}\,ds$$

$$\,=\int_{0}^{1}\left({\frac{\exp(-\theta\,F_{U}(rF_{V}^{-1}(s)))\,f_{U}(r\,F_{V}^{-1}(s))\,\frac{r}{f_{V}(F_{V}^{-1}(s))}\,(\exp(-\theta s)-1)}{\exp(-\theta)-1+\big{(}\exp(-\theta\,F_{U}(rF_{V}^{-1}(s)))-1\big{)}\,(\exp(-\theta s)-1)}}\right.$$

$$\hskip 24.18501pt\left.+\frac{\big{(}\exp(-\theta\,F_{U}(rF_{V}^{-1}(s)))-1\big{)}\,\exp(-\theta s)}{\exp(-\theta)-1+\big{(}\exp(-\theta\,F_{U}(rF_{V}^{-1}(s)))-1\big{)}\,(\exp(-\theta s)-1)}\right)ds\,,$$

where $\text{sgn}(\cdot)$ is the sign function. Rearrangement of the last equation gives (7). ∎

Figure 2: Median of $R=U/V$ for parameters $\Lambda\in[0.1,4]$ and $\delta_{U},\delta_{V}\in\{2,3\}$, as calculated from the formula in Proposition 2. Note that the median of $R$ does not vary with $\theta$.

2.2 Regression specification and estimation

To model the entire distribution of $R$ as a function of covariates $X=(X_{1},\ldots,X_{p})^{\top}$, we propose to relate both the logarithmic rate parameters $\lambda_{U}$ and $\lambda_{V}$ and the association parameter $\theta$ to predictor functions of the form

$$\log(\lambda_{U}|X) =\eta_{U}=\beta_{U0}+\beta_{U1}X_{1}+\ldots+\beta_{Up}X_{p}\,,$$ (8)

$$\log(\lambda_{V}|X) =\eta_{V}=\beta_{V0}+\beta_{V1}X_{1}+\ldots+\beta_{Vp}X_{p}\quad\text{and}$$ (9)

$$\theta|X =\eta_{\theta}=\beta_{\theta 0}+\beta_{\theta 1}X_{1}+\ldots+\beta_{\theta p}X_{p}\,,$$ (10)

where $\boldsymbol{\beta}_{U}=(\beta_{U0},\ldots,\beta_{Up})^{\top}$, $\boldsymbol{\beta}_{V}=(\beta_{V0},\ldots,\beta_{Vp})^{\top}$ and $\boldsymbol{\beta}_{\theta}=(\beta_{\theta 0},\ldots,\beta_{\theta p})^{\top}$ are sets of real-valued coefficients. Analogous to classical gamma regression (Chapter 5.3 of Fahrmeir et al., 2022), the use of the logarithmic transformation in (8) and (9) ensures positivity of the rate parameters. Since $\theta\in\mathbb{R}\setminus\{0\}$, no transformation is needed for the association parameter. As a result of (8) and (9) it holds that $\log(\Lambda|X)=\eta_{U}-\eta_{V}$.

Remark 3.

In principle, our approach allows to make use of the full flexibility of GAMLSS by relating all distributional parameters (including the shape parameters $\delta_{U},\delta_{V}$) to the covariates and by including nonlinear effects in the predictor functions. However, in our application we found that the specification in (8) – (10) provides a sufficient fit, thereby meeting a compromise between model fit and model complexity. Furthermore, it greatly simplifies the interpretability of the results (as we will further elaborate in Section 4). Based on these considerations, our model assumes that the shape parameters $\delta_{U}$ and $\delta_{V}$ do not depend on $X$, but can be treated as nuisance parameters.

Definition 1.

The regression model for the ratio $R=U/V$ with the distribution from Proposition 3 and with covariate-dependent parameters as specified in (8)–(10) will be termed Frank copula with Gamma Distributed Marginals (FCGAM) in the following. The FCGAM model imposes the constraint $\delta_{U}$, $\delta_{V}$ $>1$ to ensure that the two marginals both exhibit a unimodal, right-skewed distribution, which is the common form of biomarker distributions in medical applications.

Corollary 1.

For a set of i.i.d. observations $(u_{1},v_{1},\boldsymbol{x}_{1}^{\top})^{\top},\ldots,(u_{n},v_{n},\boldsymbol{x}_{n}^{\top})^{\top}$ with ratios $r_{1}=u_{1}/v_{1},\ldots,r_{n}=u_{n}/v_{n}$ and model coefficients $\boldsymbol{\gamma}=\left(\boldsymbol{\beta}_{U}^{\top},\boldsymbol{\beta}_{V}^{\top},\boldsymbol{\beta}_{\theta}^{\top},\delta_{U},\delta_{V}\right)^{\top}$, the log-likelihood function of the FCGAM model is given by

$$\ell(\boldsymbol{\beta}_{U},\boldsymbol{\beta}_{V}, \boldsymbol{\beta}_{\theta},\delta_{U},\delta_{V};u_{1},\ldots,u_{n},v_{1},\ldots,v_{n},\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})$$

$$=\sum_{i=1}^{n}\,\bigg{\{} \log\Big{(}f_{U,V}(u_{i},v_{i}|\boldsymbol{x}_{i},\boldsymbol{\beta}_{U},\boldsymbol{\beta}_{V},\boldsymbol{\beta}_{\theta},\delta_{U},\delta_{V})\Big{)}\bigg{\}}$$

$$=\sum_{i=1}^{n}\,\bigg{\{} \log\Big{(}\exp(-\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{\theta}\,F_{U}(u_{i};\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{U},\delta_{U}))\exp(-\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{\theta}\,F_{V}(v_{i};\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{V},\delta_{V}))$$

$$\ \times\ (-\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{\theta})(\exp(-\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{\theta})-1)f_{U}(u_{i};\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{U},\delta_{U})\,f_{V}(v_{i};\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{V},\delta_{V})\Big{)}$$

$$-2\, \log\,\Big{(}(\exp(-\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{\theta})-1)+(\exp(-\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{\theta}\,F_{U}(u_{i};\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{U},\delta_{U}))-1)$$

$$\ \times\ (\exp(-\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{\theta}\,F_{V}(v_{i};\boldsymbol{x}_{i}^{\top}\boldsymbol{\beta}_{V},\delta_{V}))-1)\Big{)}\bigg{\}}\,.$$ (11)

Corollary 2.

Under the usual regularity assumptions, the estimator

$$\hat{\boldsymbol{\gamma}}=(\hat{\boldsymbol{\beta}}^{\top}_{U}, \hat{\boldsymbol{\beta}}^{\top}_{V},\hat{\boldsymbol{\beta}}^{\top}_{\theta},\hat{\delta}_{U},\hat{\delta}_{V})^{\top}:=$$

$$\underset{\boldsymbol{\gamma}=\boldsymbol{\beta}_{U},\boldsymbol{\beta}_{V},\boldsymbol{\beta}_{\theta},\delta_{U},\delta_{V}}{\text{argmax}}\ell(\boldsymbol{\beta}_{U},\boldsymbol{\beta}_{V},\boldsymbol{\beta}_{\theta},\delta_{U},\delta_{V};u_{1},\ldots,u_{n},v_{1},\ldots,v_{n},\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{n})\,$$ (12)

is consistent and asymptotically normal for $n\to\infty$.

Implementational details. Maximization of the log-likelihood function in (1) can be carried out using the R function FCGAMoptim(), which is part of the supplementary material to this paper. The optimization algorithm is based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm implemented in the R function optim(), setting the additional constraint $\delta_{U}$, $\delta_{V}$ $>1$.

2.3 Prediction of distributional parameters and inference

Prediction. For a new observation with covariate values $\tilde{\boldsymbol{x}}$, predictions of the conditional PDF ${f}_{R}(r|\tilde{\boldsymbol{x}})$ can be obtained by computing the maximum likelihood estimate (MLE) and by plugging the estimated parameters $\hat{\Lambda}\,|\,\tilde{\boldsymbol{x}}=\exp(\tilde{\boldsymbol{x}}^{\top}\hat{\boldsymbol{\beta}}_{U}-\tilde{\boldsymbol{x}}^{\top}\hat{\boldsymbol{\beta}}_{V})$, $\hat{\theta}\,|\,\tilde{\boldsymbol{x}}=\tilde{\boldsymbol{x}}^{\top}\hat{\boldsymbol{\beta}}_{\theta}$ and $\hat{\delta}_{U}$, $\hat{\delta}_{V}$ in Equation (2). The predicted PDF can then be used to predict any distributional parameter of interest (like the expected value, median or quantiles). For example, denoting the predicted PDF by $\hat{f}_{R}(r|\tilde{\boldsymbol{x}})$, the predicted median can be calculated by

$$\hat{r}_{\text{med}}\,|\,\tilde{\boldsymbol{x}}=\text{min}\bigg{\{}r\in\mathbb{R}^{+}\,\Big{|}\,\int_{0}^{r}\hat{f}_{R}(s|\tilde{\boldsymbol{x}})\,ds\geq 0.5\bigg{\}}\,.$$ (13)

Inference. Despite the asymptotic results from Corollary 2, more reliable finite-sample confidence intervals have been established in additive models (Wood, 2017). This is particularly the case for the quantities of interest here (such as the median of $R$ above). The reason is that these are nonlinear transformations of the original model coefficients such that confidence intervals would show an additional finite-sample bias due to the application of the $\Delta$-rule. Following Wood (2017), we thus propose to construct confidence intervals of $\boldsymbol{\gamma}$ using a Bayesian approach, which we accordingly refer to as credible intervals. Assuming flat priors on $\boldsymbol{\gamma}$, the posterior distribution of $\boldsymbol{\gamma}$ is given by

$$\boldsymbol{\gamma}\,|\,u_{1},\ldots,u_{n},v_{1},\ldots,v_{n}\,\sim\,N\big{(}\hat{\boldsymbol{\gamma}},\hat{J}^{-1}(\hat{\boldsymbol{\gamma}})\big{)}\,,$$ (14)

where $\hat{J}(\hat{\boldsymbol{\gamma}})$ is the Hessian of the negative log-likelihood evaluated at $\hat{\boldsymbol{\gamma}}$ (Equation (6.26) of Wood, 2017). Consequently, approximate $(1-\alpha)\%$ credible intervals for the coefficients $\boldsymbol{\gamma}$ can be obtained by drawing a large sample from the posterior distribution (14) and by calculating the $\alpha/2$ and $(1-\alpha/2)$ percentiles from this sample (Wood, 2017, p. 293). In our simulations (Section 3) and in the analysis of the DCN study data (Section 4) we used samples of size 10,000 throughout.

3 Simulations

We conducted three simulation studies to investigate the performance of the FCGAM model. Our main aims were (a) to analyze the model fit and the coverage of the credible intervals, (b) to evaluate how the performance of the FCGAM approach is affected by the sample size and the choice of the association parameter $\theta$, and (c) to benchmark our method against alternative ones, in particular against the extended GB2 model by Berger et al. (2019) which assumes the correlation between $U$ and $V$ to be positive.

3.1 Experimental design

In all simulations the ratio outcome was generated according to the PDF of the FCGAM model derived in Proposition 1. Similar to the application data in Section 4, we considered two standard normally distributed covariates $X_{1},X_{2}\sim N(0,1)$ and two binary covariates $X_{3},X_{4}\sim B(1,0.5)$, which were equi-correlated with Pearson correlation coefficient 0.4. For each $n\in\{200,500,1000\}$ we simulated $1000$ independent data sets.

In Simulation Study 1, we considered scenarios with fixed negative correlation (the case which motivated our development of the FCGAM model), setting $\beta_{\theta 0}\in\{-1,-5,-10\}$ and $\beta_{\theta 1}=\ldots=\beta_{\theta 4}=0$. This resulted in the respective rank correlation coefficients $\tau\in\{-0.11,-0.46,-0.67\}$. The rate parameters were related to the four covariates through the coefficients $\beta_{U}=(0,0.4,-0.4,0.2,-0.2)^{\top}$ and $\beta_{V}=(0,-0.2,0.2,-0.4,0.4)^{\top}$. The shape parameters were set to $\delta_{U}=2$ and $\delta_{V}=6$ in all scenarios.

In Simulation Study 2, we considered scenarios with fixed positive correlation (the case which has already been covered by the eGB2 model but also applies to the FCGAM model), setting $\beta_{\theta 0}\in\{1,5,10\}$ and $\beta_{\theta 1}=\ldots=\beta_{\theta 4}=0$. This resulted in the respective rank correlation coefficients $\tau\in\{0.11,0.46,0.67\}$. To ensure that the outcome values were in a meaningful range (comparable to Simulation Study 1) we set the regression coefficients to $\beta_{U}=(0,0.4,-0.4,0.2,-0.2)^{\top}$ and $\beta_{V}=(0,0.2,-0.2,0.4,-0.4)^{\top}$, and the shape parameters to $\delta_{U}=2$ and $\delta_{V}=2$.

In Simulation Study 3, we evaluated how the model fit of the FCGAM model was affected when falsely assuming a dependence of $\theta$ on $X_{1},\ldots,X_{4}$, or when ignoring a present dependence of $\theta$ on $X_{1},\ldots,X_{4}$. For this we reconsidered the data sets from Simulation Study 1 with $\tau=-0.11$, and additionally considered scenarios where the association parameter $\theta$ was related to the four covariates through the coefficient vector $\beta_{\theta}=(0,1,-1,0.5,-0.5)^{\top}$ (resulting in covariate-dependent rank correlation coefficients $\tau_{i}$, with the remaining parameters as in Simulation Study 1). In both cases we fitted the FCGAM model with covariate-dependent $\theta$ (according to (10)) and with constant $\theta=\beta_{\theta 0}$.

Benchmark methods. We evaluated the fits of the 1000 FCGAM models by computing the predictive log-likelihood values on 1000 independent test data sets. The test data sets (of size $n$ each) were also used to compare the FCGAM model to alternative models. To this purpose, we evaluated the predictive log-likelihood values of the following benchmark methods, where (ii), (iii) and (vi) are univariate regression models for $R$, (iv) is a univariate regression model for $\log(R)$, and (v) and (vii) are distributional regression models:

• The copula-based FCGAM model introduced in Section 2.2.

• The extended GB2 model (Berger et al., 2019) (eGB2) assuming a positive correlation between $U$ and $V$.

• The simple GB2 model (GB2) assuming zero Pearson correlation between $U$ and $V$.

• A Gaussian regression model with log-transformed outcome values (LN).

• A Gaussian GAMLSS with log-transformed outcome values, where both the mean and the standard deviation were related to the covariates (LN.LSS). The standard deviation was modeled using the log link.

• A Gamma regression model with the original outcome values (GA). The mean parameter was related to the covariates and was modeled using the log link.

• A Gamma GAMLSS with the original outcome values, where both the mean and the scale parameters were related to the covariates (GA.LSS) using the log link.

3.2 Results

Point estimates of the FCGAM coefficients. Figure 3 presents the coefficient estimates $\hat{\beta}_{U}$ in Simulation Study 1 with negative (but covariate-independent) correlation between $U$ and $V$. The boxplots show that on average the estimated coefficients are very close to the true ones, regardless of the association parameter $\theta$. Accordingly, the finite-sample bias of the MLEs is small in all scenarios (with varying $n$ and $\theta$). From Figure 3 it can also be seen that, as expected, the variance of the estimates decreases with increasing sample size, in particular for the two binary covariates $X_{3}$ and $X_{4}$. In contrast, the correlation (determined by the value of $\theta$) has only a small impact on the variance of the estimates. The coefficient estimates $\hat{\beta}_{V}$ (presented in Supplementary Figure S1) exhibit even smaller variances in the scenarios with $n=500$ and $n=1000$.

The coefficient estimates $\hat{\beta}_{U}$ and $\hat{\beta}_{V}$ from Simulation Study 2 with positive correlation between $U$ and $V$ are shown in Supplementary Figures S2 and S3, respectively. In both cases the bias is small throughout all scenarios. Regarding the variance of the estimates, the results are largely the same as in Figure 3.

Figure 3: Point estimates of the FCGAM coefficients in Simulation Study 1. The boxplots visualize the MLEs of the coefficients $\beta_{U1}=0.4$, $\beta_{U2}=-0.4$, $\beta_{U3}=0.2$ and $\beta_{U4}-0.2$ that were obtained from fitting the FCGAM model to 1000 data sets of size $n$ each. The red lines refer to the true values of the coefficients.

Coverage proportions of the credible intervals. The coverage proportions of the 95% credible intervals obtained from the FCGAM fits are presented in Table 1. They range between $0.928-0.958$ (Simulation Study 1) and between $0.928-0.962$ (Simulation Study 2), which is close to the nominal coverage of 95%. There were only minor differences with regard to sample size and the correlation coefficient. This result demonstrates that not only point estimation but also inference works well for highly positive or negative correlations and fairly small samples.

Table 1: Coverage proportions of the FCGAM credible intervals. For each coefficient $\beta_{Uj}$, $j=1,\ldots,4,$ and $\beta_{Vj}$, $j=1,\ldots,4,$ the table contains the coverage proportion of the 95% credible interval, as obtained from fitting the FCGAM model to 1000 independent data sets of size $n$ each.

Simulation Study 1		$\beta_{U1}$	$\beta_{U2}$	$\beta_{U3}$	$\beta_{U4}$	$\beta_{V1}$	$\beta_{V2}$	$\beta_{V3}$	$\beta_{V4}$
$n$=200	$\theta=-1$	0.938	0.932	0.936	0.938	0.949	0.949	0.954	0.950
	$\theta=-5$	0.958	0.942	0.955	0.941	0.943	0.935	0.946	0.952
	$\theta=-10$	0.937	0.938	0.932	0.947	0.952	0.946	0.940	0.949
$n$=500	$\theta=-1$	0.949	0.937	0.940	0.938	0.952	0.953	0.946	0.930
	$\theta=-5$	0.954	0.938	0.928	0.935	0.952	0.932	0.935	0.949
	$\theta=-10$	0.934	0.948	0.939	0.951	0.955	0.947	0.958	0.943
$n$=1000	$\theta=-1$	0.949	0.937	0.955	0.947	0.944	0.946	0.948	0.937
	$\theta=-5$	0.947	0.942	0.936	0.936	0.957	0.950	0.950	0.950
	$\theta=-10$	0.933	0.945	0.934	0.953	0.943	0.948	0.949	0.943
Simulation Study 2		$\beta_{U1}$	$\beta_{U2}$	$\beta_{U3}$	$\beta_{U4}$	$\beta_{V1}$	$\beta_{V2}$	$\beta_{V3}$	$\beta_{V4}$
$n$=200	$\theta=1$	0.935	0.929	0.934	0.935	0.947	0.933	0.959	0.954
	$\theta=5$	0.952	0.940	0.954	0.941	0.949	0.951	0.954	0.947
	$\theta=10$	0.938	0.939	0.939	0.947	0.951	0.944	0.948	0.953
$n$=500	$\theta=1$	0.941	0.948	0.936	0.928	0.940	0.942	0.946	0.940
	$\theta=5$	0.948	0.941	0.929	0.929	0.952	0.944	0.941	0.962
	$\theta=10$	0.931	0.947	0.937	0.951	0.948	0.944	0.949	0.937
$n$=1000	$\theta=1$	0.950	0.944	0.954	0.944	0.952	0.939	0.953	0.935
	$\theta=5$	0.955	0.942	0.934	0.934	0.956	0.952	0.949	0.947
	$\theta=10$	0.928	0.946	0.934	0.955	0.935	0.956	0.948	0.945

Distributional prediction. The root mean squared error (RMSE) of the estimated conditional median values computed from (13) are given in Table 2. In Simulation Study 1 the performance is quite similar for all three values of $\theta$, whereas in Simulation Study 2 the RMSE considerably decreases with increasing value of $\theta$. This indicates that estimating the median value works best for highly positive correlations where the PDF of $R$ is rather diffuse with a large mode value (compare Figure 1). It is also seen from Table 2 that the means and the standard deviations of the RMSE decrease with increasing sample size.

Table 2: RMSE of the estimated conditional median values of $R$. The table presents the mean RMSE of the estimated conditional median of $R$, as obtained from fitting the FCGAM model to 1000 independent data sets of size $n$ each. Standard deviations of the RMSE values (across the 1000 data sets) are given in brackets.

Simulation Study 1	$\theta=-1$	$\theta=-5$	$\theta=-10$
$n$ = 200	0.076 (0.035)	0.084 (0.039)	0.081 (0.037)
$n$ = 500	0.049 (0.021)	0.053 (0.023)	0.051 (0.021)
$n$ = 1000	0.034 (0.015)	0.038 (0.016)	0.036 (0.016)
Simulation Study 2	$\theta=1$	$\theta=5$	$\theta=10$
$n$ = 200	0.154 (0.056)	0.101 (0.035)	0.063 (0.021)
$n$ = 500	0.097 (0.033)	0.064 (0.022)	0.040 (0.014)
$n$ = 1000	0.069 (0.023)	0.046 (0.015)	0.029 (0.010)

Figure 4: Comparison of the FCGAM model to alternative methods in Simulation Study 1. The boxplots visualize the predictive log-likelihood values obtained from the FCGAM model and from the benchmark methods (ii) to (vii). All models were fitted to 1000 independent data sets and evaluated on independently generated test data sets of the same size. In each panel, the dashed horizontal line indicates the median predictive log-likelihood value of the best performing method.

Comparison to alternative models. Figure 4 and Supplementary Figure S4 show the prediction accuracy (i.e. the predicted log-likelihood values on the test sets) of the FCGAM model and the benchmark methods (ii) to (vi). In Simulation Study 1 with negative correlation, it can be observed that the FCGAM model achieves the highest accuracy in all scenarios (Figure 4). The superiority is even more evident when the sample size and the value of the correlation coefficient are increased. The extended GB2 and simple GB2 methods yield similar performances as the Gaussian models with log-transformed outcome (LN and LN.LSS), whereas the Gamma regression models (GA and GA.LSS) result in the lowest accuracy. For both LN and GA the GAMLSS models are not superior to their simple counterparts.

In Simulation Study 2 with positive correlation, the results change considerably (Supplementary Figure S4). As expected, the performance of the FCGAM and eGB2 models is largely the same, as the eGB2 model also assumes gamma distributed components with positive correlation. The simple GB2 model (assuming uncorrelated components) and the Gaussian models with log-transformed outcomes (LN and LN.LSS) perform comparably well in the scenarios with $\theta=1$, but deteriorated with increasing correlation ($\theta=5$ and $\theta=10$). Again, the Gamma regression models (GA and GA.LSS) exhibit the worst performance by far.

Table 3: Analysis of misspecified models for the association parameter in Simulation Study 3. The table presents the mean RMSE of the estimated conditional median values (upper part) and the mean of the predictive log-likelihood values (lower part), as obtained from fitting the FCGAM model to 1000 independent data sets and evaluating the fits on 1000 independently generated test data sets. Standard deviations (across the 1000 data sets) are given in brackets. The left part of the table refers to the scenarios with covariate-dependent correlation coefficients (in the observed range $\tau_{i}\in[-0.483,\ldots,0.464]$), whereas the right part refers to the scenarios with fixed negative correlation $\tau=-0.11$. The terms “modeled $\theta$” and “constant $\theta$” refer to the FCGAM models with a covariate-dependent predictor function for $\theta$ (as in (10)) and an intercept-only predictor function for $\theta$ ($\beta_{\theta 1}=\ldots=\beta_{\theta p}=0$ in (10)), respectively.

RMSE	Covariate-dependent correlation		Fixed correlation
	modeled $\theta$	constant $\theta$	modeled $\theta$	constant $\theta$
$n$ = 200	0.069 (0.031)	0.082 (0.042)	0.079 (0.032)	0.077 (0.034)
$n$ = 500	0.042 (0.018)	0.058 (0.026)	0.050 (0.020)	0.047 (0.018)
$n$ = 1000	0.031 (0.013)	0.049 (0.021)	0.036 (0.015)	0.034 (0.014)
Predictive	Covariate-dependent correlation		Fixed correlation
log-likelihood	modeled $\theta$	constant $\theta$	modeled $\theta$	constant $\theta$
$n$ = 200	-6.769 (19.176)	-7.116 (19.133)	-20.165 (19.163)	-19.356 (19.056)
$n$ = 500	-10.736 (29.031)	-12.853 (29.087)	-44.585 (29.092)	-43.979 (29.050)
$n$ = 1000	-18.237 (41.177)	-23.330 (40.965)	-85.656 (40.989)	-85.071 (40.967)

Misspecified models for the association parameter in Simulation Study 3. The RMSE of the estimated conditional median values and the predictive log-likelihood values of the FCGAM fits are summarized in Table 3. It is seen that ignoring the dependence of $\theta$ on the covariates (left part of Table 3) decreases both the predictive ability and the model fit. In the scenario with $n=1000$ (large sample size), the difference in predictive log-likelihood values of 5.093 suggests “considerably less” empirical support for the model with constant $\theta$ (according to the rules of thumb provided in Burnham and Anderson, 2002). On the other hand, when unnecessarily modeling the dependence of $\theta$ on $X_{1},\ldots,X_{4}$ (right part of Table 3) the predictive ability and the model fit are mostly unaffected (showing only negligible differences in the RMSE and the predictive log-likelihood values).

Overall summary. Taken together, we make the following key empirical observations:

1. Point estimates from the FCGAM model are reliable and nearly unbiased even for small sample sizes.

2. The FCGAM model outperforms the eGB2 model in case of negative correlation and is en par with the eGB2 model when the correlation is positive.

3. In all scenarios, the Gamma regression models perform worst. In particular, they perform worse than the Gaussian models with log-transformed outcome.

4. Falsely modeling the association parameter does not deteriorate predictive performance to a large degree, whereas the FCGAM model with covariate-dependent $\theta$ improves the fit when the true association depends on the covariates.

4 Cohort Study of the German Dementia Competence Network

Background. The multi-center cohort study conducted by the German Dementia Competence Network (DCN; Kornhuber et al., 2009) enrolled patients aged older than 50 years that were diagnosed with either mild cognitive impairment (MCI), Alzheimer’s disease (AD) or other dementia. Recruitment took place between 2003 and 2007. The main objective of the original study was to establish biomarkers for the diagnosis and prognosis of AD using clinical, laboratory and imaging measurements. Here, we investigate covariates that are potentially associated with amyloid-$\beta$ 42, amyloid-$\beta$ 40 and total tau protein concentrations measured in cerebrospinal fluid samples. These analyses are of high relevance for clinical routine in the neurosciences, since biomarkers enable the detection of AD pathology long before the occurrence of the first clinically obvious symptoms (Sperling et al., 2013). Thus, relating covariates to biomarker values provides insight into disease pathology and prevention at the individual patient level. In the neurosciences, amyloid-$\beta$ 42, amyloid-$\beta$ 40 and total tau protein concentrations are usually not analyzed separately but in terms of their ratios. More specifically, the amyloid-$\beta$ 42/40 ratio and amyloid-$\beta$ 42/total tau ratio are considered to be strong predictors of AD progression (Koyama et al., 2012). Therefore we focus on the group of MCI patients and relate their ratios to patient-related risk factors for dementia.

(a)

(b)

Figure 5: Analysis of the DCN study data. Distribution of the amyloid-$\beta$ 42/40 ratios (a) and the amyloid-$\beta$ 42/total tau ratios (b) in patients with MCI ($n=330$).

Description of the data. In the DCN study, amyloid-$\beta$ and total tau baseline concentrations were measured in 374 patients diagnosed with MCI. In all other MCI patients, CSF biosamples were not collected due to either logistic reasons or lack of consent to the invasive procedure of lumbar puncture. Exclusion of patients that did not meet the inclusion criterion (age $\leq$ 50 years; $7$ patients) and of patients with missing values in at least one of the considered risk factors ($37$ patients) resulted in an analysis data set of $n=330$ patients. For details on the handling of missing values we refer to Berger et al. (2019). The unconditional distributions of the amyloid-$\beta$ 42/40 ratio and the amyloid-$\beta$ 42/total tau ratio are visualized in Figure 5. While the values of the amyloid-$\beta$ 42/40 ratios are all smaller than 0.3, the amyloid-$\beta$ 42/total tau ratios range between 0.2 and 13, exhibiting a heavily right-skewed distribution. Kendall’s rank correlation coefficient between the two components is measured to be $\tau=0.307$ (amyloid-$\beta$ 42/40) and $\tau=-0.269$ (amyloid-$\beta$ 42/total tau). Thus, our analysis had to deal with both positive and negative correlations between the ratio components. As mentioned before, this problem was our main motivation for the development of the FCGAM model.

The risk factors included in the analysis are summarized in Table 4. These were: (i) sex, (ii) age in years, (iii) educational level (measured by the number of years of education), and (iv) a binary variable indicating whether a patient was a carrier of the apolipoprotein E$\epsilon$4 (ApoE $\epsilon$4) allele, which is a strong genetic predictor of AD.

Table 4: Description and summary statistics of the two ratio outcomes and the covariates used for the analysis of the DCN study data (Q1 = first quartile, Q3 = third quartile). All numbers refer to a subset of patients diagnosed with MCI ($n=330$). For details on the collection of the data, see Kornhuber et al. (2009).

Variable	Summary statistics
	min	Q1	median	Q3	max	mean	sd
Amyloid-$\beta$ 42/40	0.03	0.08	0.10	0.14	0.26	0.11	0.04
Amyloid-$\beta$ 42/total tau	0.19	0.91	2.13	3.72	12.95	2.70	2.34
Age (years)	51	60	66	73	89	66.51	8.11
Education (years)	2	11	11	13	19	12.18	2.96
Sex	male:	194 (58.8%)			female:	136 (41.2%)
ApoE $\epsilon$4	no:	182 (55.2%)			yes:	148 (44.8%)

Model fitting I. In a preliminary analysis, we fitted GA and GA.LSS models for the components amyloid-$\beta$ 42, amyloid-$\beta$ 40 and total tau, where either the rate parameters only (GA) or both the rate and the shape parameters (GA.LSS) were related to the four covariates. According to the Bayesian information criterion (BIC) the simple GA models (BIC $=4899.988$ for amyloid-$\beta$ 42, BIC $=6248.201$ for amyloid-$\beta$ 40 andBIC $=4574.974$ for total tau) showed better fits than the respective GA.LSS models (BIC $=4914.396$ for amyloid-$\beta$ 42, BIC $=6266.375$ for amyloid-$\beta$ 40 and BIC $=4581.902$ for total tau). This result indicates that it is sufficient to relate the two rate parameters to the covariates. Furthermore, it supports the assumptions of the proposed FCGAM model, which treats the shape parameters $\delta_{U}$ and $\delta_{V}$ as nuisance parameters.

The fits of the FCGAM model with covariate-dependent association parameter are presented in Supplementary Table S1. According to the credible intervals given in columns 4 and 6, none of the risk factors is found to affect the association parameter $\theta$. Applying Equation (3) yielded the mean estimated rank correlations $\hat{\tau}(\hat{\theta})=0.35$ (range: 0.21 to 0.49) for amyloid-$\beta$ 42/40 and $\hat{\tau}(\hat{\theta})=-0.23$ (range: -0.41 to 0.03) for amyloid-$\beta$ 42/total tau. Both estimates are close to the respective unconditional rank correlations.

Model fitting II. Based on the above findings and to further reduce model complexity, we fitted FCGAM models with constant association parameter $\theta$ (setting the coefficients $\beta_{\theta,\text{Age}},\ldots,\beta_{\theta,\text{ApoE $\epsilon 4$}}$ to zero). We then calculated the BIC from these reduced models along with their counterparts obtained from the models with covariate-dependent $\theta$. For amyloid-$\beta$ 42/40, the BIC values were $11063.74$ (constant $\theta$) and $11084.75$ (modeled $\theta$). For amyloid-$\beta$ 42/total tau, the BIC values were $9249.99$ (constant $\theta$) and $9267.485$ (modeled $\theta$). This result suggests that the reduced models with constant $\theta$ meet a better compromise between model fit and model complexity than the respective full models with covariate-dependent $\theta$. This observation is confirmed by randomized quantile residuals of the reduced FCGAM models (Figure 6) indicating only slight deviations from normality.

(a)

(b)

Figure 6: Analysis of the amyloid-$\beta$ 42/40 ratios (a) and the amyloid-$\beta$ 42/total tau ratios (b) in the DCN study data. The figure presents normal quantile-quantile plots of the quantile residuals obtained from the FCGAM model fits.

Main results. The results obtained from the reduced FCGAM models are shown in Table 5 and Supplementary Figures S5 and S6. The upper part of Table 5 refers to the parameter $\Lambda=\lambda_{U}/\lambda_{V}$, reporting the differences $\hat{\beta}_{\Lambda j}:=\hat{\beta}_{Uj}-\hat{\beta}_{Vj}$. Note that the coefficient estimates are very similar to the respective estimates of the more complex model in Table S1. For example, for amyloid-$\beta$ 42/total tau one obtains $\hat{\beta}_{\Lambda,\text{ApoE $\epsilon 4$}}=0.3786$ (Table 5) and $\hat{\beta}_{\Lambda,\text{ApoE $\epsilon 4$}}=0.2411+0.1406=0.3817$ (Table S1). The credible intervals in Table 5 were obtained by drawing a sample of size 10,000 from the posterior distribution in (14) and by calculating the $2.5\%$ and $97.5\%$ percentiles from the sampled differences ${\beta}_{Uj}-{\beta}_{Vj}$. According to the results of the FCGAM model, there is strong evidence for an effect of the risk factors age and ApoE $\epsilon$4 on theamyloid-$\beta$ 42/40 and amyloid-$\beta$ 42/total tau ratios. As depicted in Supplementary Figures S5(a) and S6(a), both the expected amyloid-$\beta$ 42/40 ratio and the expected amyloid-$\beta$ 42/total tau ratio decrease with increasing age, implying a higher risk of progression to AD in older patients. Similarly, the expected ratios of ApoE $\epsilon$4 carriers are strongly reduced compared to patients not carrying the allele (Supplementary Figures S5(d) and S6(d), confirming the important role of this genetic risk factor in AD progression). The figures also illustrate how the estimated median values as well as the $10\%$ and $90\%$ percentiles of the distributions change with the covariates. In contrast to age and ApoE $\epsilon$4, Table 5 shows no evidence for an effect of sex and educational level on the two ratio outcomes. These results are in full agreement with the findings by Berger et al. (2019), who fitted an eGB2 model with amyloid-$\beta$ 42/40 outcome to the DCN study data.

Table 5: Analysis of the amyloid-$\beta$ 42/40 ratios (left) and the amyloid-$\beta$ 42/total tau ratios (right) in the DCN study data. The table presents the coefficient estimates with $95\%$ credible intervals (calculated by the procedure described in Section 2.3), as obtained from fitting FCGAM models with constant association parameter $\theta$.

		amyloid-$\boldsymbol{\beta}$ 42/40		amyloid-$\boldsymbol{\beta}$ 42/total tau
Parameter	Covariate	$\hat{\beta}$	$95\%$ CI	$\hat{\beta}$	$95\%$ CI
$\Lambda$	Age (years)	$0.0089$	$[0.0041;0.0137]$	$0.0256$	$[0.0156;0.0377]$
	Education (years)	$-0.0017$	$[-0.0152;0.0117]$	$-0.0150$	$[-0.0461;0.0168]$
	Sex (male)	.	.	.	.
	Sex (female)	$0.0724$	$[-0.0073;0.1547]$	$0.0283$	$[-0.1544;0.2146]$
	ApoE $\epsilon$4 (no)	.	.	.	.
	ApoE $\epsilon$4 (yes)	$0.1967$	$[0.1157;0.2766]$	$0.3786$	$[0.1965;0.5584]$
$\theta$		$3.5325$	$[2.7671;4.2888]$	$-2.0611$	$[-2.8064;-1.3073]$
$\delta_{U}$		$6.0586$	$[5.1169;7.0218]$	$5.8090$	$[4.9062;6.7092]$
$\delta_{V}$		$10.0151$	$[8.5182;11.5039]$	$2.6718$	$[2.2900;3.0542]$

5 Discussion

The main contribution of this work is a copula-based regression model that relates the ratio of two gamma distributed components to a set of covariates. Our model is primarily designed for the analysis of ratio outcomes in medical research, which is an important task, for instance, in neurology (Novellino et al., 2021), infectiology (Caby et al., 2016) and pharmacology (Cawley et al., 2022). Importantly, when biomarker ratios are used as clinical metrics or indicators of clinical outcomes, our model may be used to relate the respective ratio values to a set of risk factors and/or confounding variables. A prototypical example is given by the prognosis of AD progression considering ratios of amyloid-$\beta$ and total tau protein biomarkers, as presented in Section 4 of this paper.

Conceptually, the FCGAM model developed in this paper has the following advantages: First, by assuming the ratio components to follow marginal gamma distributions, the FCGAM model represents the two biomarkers by real-valued random variables with positive support and right-skewed (marginal) distributions. These distributional characteristics, which are common to most biomarkers encountered in medical research, are directly incorporated in the definition of the proposed copula model. As a consequence, the resulting ratio density incorporates the full information contained in the marginal densities of the components of the ratio. We emphasize that this property does not apply to simpler modeling approaches approximating the ratio by a single log-normal or gamma-distributed variable. In fact, without consideration of the paired components themselves, these approximations inevitably bear the risk of a loss of information (“neglected companion”; Kerkhof et al., 2019). This issue has been demonstrated by the results of our simulation study (Section 3), which resulted in an increased estimation accuracy of the proposed copula-based approach in all data-generating scenarios. We also stress that linking the two marginal distributions by a copula does in general not restrict our model to the use of two gamma distributions for the ratio components. In fact, although our model can be seen as the most relevant use case in many medical applications, the marginal distributions can in principle be replaced by arbitrary parametric distributions. For instance, our model can in a straightforward manner be extended to situations where one biomarker is discrete or ordinal.

Second, the proposed FCGAM model has a high flexibility regarding the direction of the association between the two ratio components. Importantly, by the choice of Frank copula, the FCGAM model allows for both positive and negative values of the (rank) correlation between the components $U$ and $V$, thereby improving previous modeling approaches that restricted this correlation to be zero (Yee, 2015) or positive (Berger et al., 2019). As demonstrated in the simulation study in Section 3, the FCGAM model indeed performs better in terms of estimation accuracy when the association between $U$ and $V$ is negative. On the other hand, it does not perform worse than the aforementioned approach when the association between $U$ and $V$ is positive.

Third, although the proposed model incorporates the full information contained in the marginal densities $f_{U}$ and $f_{V}$, it provides a rather simple interpretation of the associations between the ratio $U/V$ and the covariates. This is because the FCGAM model reduces the original five-parameter set $(\lambda_{U},\delta_{U},\lambda_{V},\delta_{V},\theta)^{\top}$ (including all parameters of the marginal densities and the association parameter $\theta$) to the restricted set $(\Lambda,\delta_{U},\delta_{V},\theta)^{\top}$ with $\Lambda=\lambda_{U}/\lambda_{V}$. As a consequence, when treating $\delta_{U}$, $\delta_{V}$ (and possibly also $\theta$) as nuisance parameters, the associations between $U/V$ and each of the covariates can be investigated using one-dimensional coefficient estimates and single-parameter hypothesis tests. Similarly, the association between the components $U$ and $V$ has a natural interpretation in terms of Kendall’s rank correlation, being related to $\theta$ by the one-to-one relationship given in Equation (3).

Finally, beside the flexibility in specifying other marginal distributions than the gamma distribution, the FCGAM model may be extended in many other ways. For example, Frank copula could be replaced by other copulas (noting that the results on ratio densities are also valid for other absolutely continuous copulas; see Ly et al., 2019). When there is particular interest in the tail dependencies of $U$ and $V$, benchmark experiments to identify the best fitting copula and/or marginal distributions could be performed using resampling techniques (e.g. bootstrapping or subsampling). It should be noted, however, that other copulas from the literature might be less flexible regarding the range of $\theta$ (and thus also the range of possible associations between the components $U$ and $V$; see e.g. Ghosh et al., 2022, for a recent overview of copulas allowing for modeling negative dependence). For example, it is not possible to model negative associations between $U$ and $V$ using non-rotated Gumbel or Joe copulas.

Despite our biostatistical focus, the proposed FCGAM methodology is a general statistical modeling approach that can in principle be used in any research discipline. Interesting areas are e.g. environmental research (Perri and Porporato, 2022) and information engineering (Mekić et al., 2012).

6 Appendices

It contains the proof of Proposition 2 (Appendix A), additional simulation results (Appendix B) as well as further results of the analysis of the DCN study data (Appendix C).

7 Acknowledgments

Moritz Berger acknowledges support by the grant BE 7543/1-1 of the German research foundation (DFG). Nadja Klein acknowledges support by the Emmy Noether grant KL 3037/1-1 of the DFG. The analysis of the DCN study data was supported by the German Federal Ministry of Education and Research (Kompetenznetz Demenzen, grant 01GI0420)

Appendix Appendix A Proof of Proposition 2

By Proposition 1, the PDF of the ratio $R=U/V$ is given by

$$f_{R}(r;\lambda_{U},\lambda_{V},\delta_{U},\delta_{V},\theta)=\int_{0}^{1} \,c_{\theta}\left(F_{U}(r\,F_{V}^{-1}(s;\lambda_{V},\delta_{V});\lambda_{U},\delta_{U}),s\right)$$

$$\times\,\left|F_{V}^{-1}(s;\lambda_{V},\delta_{V})\right|\,f_{U}(rF_{V}^{-1}(s;\lambda_{V},\delta_{V}),\lambda_{U},\delta_{U})\,ds$$

$$=\int_{0}^{1} c_{\theta}\left(\frac{1}{\Gamma(\delta_{U})}\,\gamma\left(\delta_{U},\lambda_{U}\,r\,F_{V}^{-1}(s;\lambda_{V},\delta_{V})\right),s\right)$$

$$\times\,F_{V}^{-1}(s;\lambda_{V},\delta_{V})\ \frac{\lambda_{U}^{\delta_{U}}}{\Gamma(\delta_{U})}\left(r\,F_{V}^{-1}(s;\lambda_{V},\delta_{V})\right)^{\delta_{U}-1}$$

$$\times\,\exp\left(-\lambda_{U}\,r\,F_{V}^{-1}(s;\lambda_{V},\delta_{V})\right)\,ds$$

$$=\int_{0}^{1} c_{\theta}\left(\frac{1}{\Gamma(\delta_{U})}\,\gamma\left(\delta_{U},\lambda_{U}\,r\ \frac{\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})\,s\right)}{\lambda_{V}}\right),s\right)$$

$$\times\,\frac{\lambda_{U}^{\delta_{U}}}{\Gamma(\delta_{U})}\left(\frac{\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})\,s\right)}{\lambda_{V}}\right)^{\delta_{U}}r^{\delta_{U}-1}$$

$$\times\,\exp\left(-\lambda_{U}\,r\ \frac{\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})\,s\right)}{\lambda_{V}}\right)\,ds$$

$$\overset{\Lambda=\lambda_{U}/\lambda_{V}}{=}\int_{0}^{1} {c_{\theta}\left(\frac{1}{\Gamma(\delta_{U})}\,\gamma\left(\delta_{U},r\,\Lambda\,\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})s\right)\right),s\right)}$$

$$\times\,\frac{\Lambda^{\delta_{U}}\,r^{(\delta_{U}-1)}}{\Gamma(\delta_{U})}\,\left(\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})s\right)\right)^{\delta_{U}}$$

$$\times\,\exp\left(-r\,\Lambda\,\gamma^{-1}\left(\delta_{V},\Gamma(\delta_{V})s\right)\right)\,ds$$

$$=f_{R} (r;\Lambda,\delta_{U},\delta_{V},\theta)\,,$$

where $\gamma(\cdot,\cdot)$ denotes the lower incomplete gamma function. $\square$

Appendix Appendix B Further Simulation Results

Figure S1: Point estimates of the FCGAM coefficients in Simulation Study 1. The boxplots visualize the MLEs of the coefficients $\beta_{V1}=-0.2$, $\beta_{V2}=0.2$, $\beta_{V3}=-0.4$ and $\beta_{V4}=0.4$ that were obtained from fitting the FCGAM model to 1000 data sets of size $n$ each. The red lines refer to the true values of the coefficients.

Figure S2: Point estimates of the FCGAM coefficients in Simulation Study 2. The boxplots visualize the MLEs of the coefficients $\beta_{U1}=0.4$, $\beta_{U2}=-0.4$, $\beta_{U3}=0.2$ and $\beta_{U4}=-0.2$ that were obtained from fitting the FCGAM model to 1000 data sets of size $n$ each. The red lines refer to the true values of the coefficients.

Figure S3: Point estimates of the FCGAM coefficients in Simulation Study 2. The boxplots visualize the MLEs of the coefficients $\beta_{V1}=0.2$, $\beta_{V2}=-0.2$, $\beta_{V3}=0.4$ and $\beta_{V4}=-0.4$ that were obtained from fitting the FCGAM model to 1000 data sets of size $n$ each. The red lines refer to the true values of the coefficients.

Figure S4: Comparison of the FCGAM model to alternative methods in Simulation Study 2. The boxplots visualize the predictive log-likelihood values obtained from the FCGAM model and from the benchmark methods (ii) to (vii). All models were fitted to 1000 independent data sets and evaluated on independently generated test data sets of the same size. In each panel, the dashed horizontal line indicates the median predictive log-likelihood value of the best performing method.

Appendix Appendix C Further Results of the Analysis of the DCN Study Data

Table S1: Analysis of the amyloid-$\beta$ 42/40 ratios (left) and the amyloid-$\beta$ 42/total tau ratios (right) in the DCN study data. The table presents the coefficient estimates with $95\%$ credible intervals (calculated by the procedure described in Section 2.3), as obtained from fitting FCGAM models with covariate-dependent association parameter $\theta$.

		amyloid-$\boldsymbol{\beta}$ 42/40		amyloid-$\boldsymbol{\beta}$ 42/total tau
Parameter	Covariate	$\hat{\beta}$	$95\%$ CI	$\hat{\beta}$	$95\%$ CI
$\lambda_{U}$	Age (years)	$0.0073$	$[0.0019;0.0127]$	$0.0092$	$[0.0037;0.0148]$
	Education (years)	$0.0079$	$[-0.0071;0.0224]$	$0.0035$	$[-0.0123;0.0192]$
	Sex (male)	.	.	.	.
	Sex (female)	$0.0017$	$[-0.0905;0.0933]$	$-0.0247$	$[-0.1178;0.0680]$
	ApoE $\epsilon$4 (no)	.	.	.	.
	ApoE $\epsilon$4 (yes)	$0.1774$	$[0.0871;0.2681]$	$0.2411$	$[0.1503;0.3325]$
$\lambda_{V}$	Age (years)	$-0.0015$	$[-0.0056;0.0025]$	$-0.0160$	$[-0.0242;-0.0079]$
	Education (years)	$0.0100$	$[-0.0013;0.0212]$	$0.0166$	$[-0.0071;0.0406]$
	Sex (male)	.	.	.	.
	Sex (female)	$-0.0751$	$[-0.1460;-0.0061]$	$-0.0746$	$[-0.2127;0.0663]$
	ApoE $\epsilon$4 (no)	.	.	.	.
	ApoE $\epsilon$4 (yes)	$-0.0141$	$[-0.0811;0.0526]$	$-0.1406$	$[-0.2745;-0.0065]$
$\theta$	Age (years)	$0.0369$	$[-0.0645;0.1347]$	$0.0637$	$[-0.0356;0.1640]$
	Education (years)	$0.1247$	$[-0.1453;0.3928]$	$0.0165$	$[-0.2416;0.2750]$
	Sex (male)	.	.	.	.
	Sex (female)	$-0.5679$	$[-2.1432;1.0406]$	$-0.9324$	$[-2.4944;0.6044]$
	ApoE $\epsilon$4 (no)	.	.	.	.
	ApoE $\epsilon$4 (yes)	$0.5610$	$[-0.9635;2.0978]$	$-1.2897$	$[-2.8033;0.1928]$
$\delta_{U}$		$5.7822$	$[4.9088;5.6759]$	$5.6900$	$[4.8084;6.5605]$
$\delta_{V}$		$10.0740$	$[8.4946;11.6237]$	$2.6023$	$[2.2232;2.9831]$

(a)

(b)

(c)

(d)

Figure S5: Analysis of the amyloid-$\beta$ 42/40 ratios in the DCN study data. The black lines refer to the estimated PDFs for a covariate profile of a randomly selected study participant (60 years of age, Sex = male, Education = 8 years, ApoE $\epsilon$4 = no). The gray lines refer to a situation where the participant would have been 70 years of age (a), would have been female (b), would have had 18 years of education (c), and would have been a carrier of the ApoE $\epsilon$4 allele (d). The vertical lines correspond to the estimated mean values (solid), median values (dashed) and $10\%$ and $90\%$ percentiles (dotted).

(a)

(b)

(c)

(d)

Figure S6: Analysis of the amyloid-$\beta$ 42/total tau ratios in the DCN study data. The black lines refer to the estimated PDFs for a covariate profile of a randomly selected study participant (60 years of age, Sex = male, Education = 8 years, ApoE $\epsilon$4 = no). The gray lines refer to a situation where the participant would have been 70 years of age (a), would have been female (b), would have had 18 years of education (c), and would have been a carrier of the ApoE $\epsilon$4 allele (d). The vertical lines correspond to the estimated mean values (solid), median values (dashed) and $10\%$ and $90\%$ percentiles (dotted).