Shared parameter and copula models for analysis of semicontinuous longitudinal data with nonrandom dropout and informative censoring

Abstract

Analysis of longitudinal semicontinuous data characterized by subjects’ attrition triggered by nonrandom dropout is complex and requires accounting for the within-subject correlation, and modeling of the dropout process. While methods that address the within-subject correlation and missing data are available, approaches that incorporate the nonrandom dropout, also referred to informative right censoring, in the modeling step are scarce due to the computational intensity and possible intractable integration needed for its implementation. Appreciating the complexity of this problem and the need for a new methodology that is feasible for implementation, we propose to extend a framework of likelihood-based marginalized two-part models to account for informative right censoring. The censoring process is modeled using two approaches: (1) Poisson censoring for the count of visits before dropout and (2) survival time to dropout. Novel consideration was given to the proposed joint modeling approaches for the semicontinuous and censoring components of the likelihood function which included (1) shared parameter, and (2) Clayton copula. The cross-part and within-part correlations were accounted for through a complex random effect structure that models correlated random intercepts and slopes. Feasibility of implementation, and accuracy of these approaches were investigated using extensive simulation studies and clinical application.

Keywords

Copula generalized gamma informative right censoring lognormal distribution marginalized two-part model nonrandom dropout semicontinuous longitudinal data shared parameters

1. Introduction

Semicontinuous is a type of data characterized by having a mixture of high preponderance of zero values and other continuously distributed nonzero values that are typically positively skewed. These types of data are often encountered in biomedical studies, health services, and clinical research. For example, in studies interested to determine whether a gene was expressed or not, and if expressed what the level of expression is, the data are characterized by zero values for non-expressed genes and positive values for expressed genes. Semicontinuous data are viewed as being generated from two processes, one that indicates whether the values are zero and this is referred to as the binary part, and the second determines the nonzero values and is denoted as the continuous part.¹

Similar issues arise in the analysis of count data with excess zeroes. Such data are usually analyzed using the hurdle two-part (TP) conditional or zero-inflated models. The hurdle model consists of a hurdle which is the point mass at zero modeled as a binary logistic that indicates whether the value is zero or not, and a truncated or overdispersed truncated Poisson for the positive counts.^2–6 On the other hand, semicontinuous data are typically analyzed following the framework of TP models. This framework follows either the TP mixed model or the marginalized two-part (mTP) model whereby the two processes, binary and continuous, are linked by random effects structure. The TP model is appropriate when interest is in determining the conditional covariate effect on a subject-specific level. Whereas the mTP is feasible when interest is in assessing the straightforward interpretation of the covariate effects at a population-averaged marginal level, as is the case in this article.

Semicontinuous data measured at a cross-sectional level can be analyzed as a Tobit model⁷ and as a TP model.^8,9 Other approaches that are based on the general exponential dispersion family including the Tweedie family distributions¹⁰ were also presented as an alternative to TP models. In particular, the compound Poisson gamma, a special case of the Tweedie distribution, with a positive mass at zero has been used to analyze zero-inflated data.¹¹ The Tweedie-based approach is attractive since it models both the zero and nonzero parts simultaneously without resorting to mixtures of models. This translates into a simpler interpretation of the estimates and less complex estimation process compared to the TP model.¹² In terms of performance, the Tweedie-based approaches were shown to have relatively precise estimates when the proportion of zero is low (around 20% and below). However, more complex models such as the TP generalized gamma were proposed as better alternatives when the zero values are proliferated in the data¹³ as is the case with the clinical application used in our study.

Longitudinal semicontinuous data can be analyzed using the TP longitudinal approaches^1,14–18 and the mTP longitudinal models.^19–21 Other approaches include TP models that are based on estimating equations for clustered data^,22,23 hierarchical zero-inflated lognormal model for repeated skewed responses,²⁴ Bayesian-based, and hierarchical Bayesian approaches.^25–28 These methodologies assume that responses are fully recorded without any missingness, and in case of incompleteness, missing process is treated as a random non-informative missing. Recently, Biswas and Das,²⁸ proposed a TP Bayesian approach for analyzing zero-inflated continuous longitudinal response with incomplete covariates and monotone missing for the response. The binary part is modeled as probit model as in Albert and Chib (1993),²⁹ and the continuous part is modeled through varying coefficients dynamic model. Missing covariates were imputed repeatedly using their corresponding posterior predictive distributions, and missing responses were imputed using dynamic working models under different identifying restrictions^,30,31 and assuming a missing at random mechanism for the covariates and the responses.

To analyze semicontinuous longitudinal data with non-ignorable missing, Liu³² developed a joint model for the semicontinuous data and survival component. The random effects from the TP model of the longitudinal response are incorporated in the death hazard model. The logit of the probability of response being positive, the logarithm of the positive response, and the hazard for terminal event are joined by random intercepts assumed to be independent and normally distributed. This joint model is advantageous in the sense that it can be implemented using available statistical procedures; however, it has some imbedded limitations. These include the logarithmic transformation of the positive response that requires, for prediction purposes, the retransformation of the logarithmic scale back to the original scale which could complicate the analysis. In addition, the joint model includes the random intercepts only, meanwhile more complicated random effects structure that incorporates the slopes could be considered. Departure from the true underlying random effect structure that has both random intercepts and slopes was shown to introduce substantial bias in the subject-specific conditional effects of the covariates.^1,33 Moreover, in this approach, the compound symmetry was assumed for the variance–covariance structure of the random effects, but other structures could have also been considered and assessed.

A flexible Bayesian mixed-effect zero-inflated TP model for semicontinuous longitudinal data with nonrandom dropout was presented by Mahabadi.³⁴ The nonzero part of the response has a link function assumed to be the logarithmic link and is joined to the missing process through shared random effect parameters. While the logarithmic is one of the simplest link functions, other link functions might not be straightforward when back transforming to original scale. Moreover, misspecification of the random effects in the TP model may lead to bias in the estimated conditional subject-specific effects of covariates in the binary part. Meanwhile, this misspecification has a lesser effect on the estimated marginal effects under the mTP model.^1,35

Despite the great interest in the analysis of semicontinuous data, cross-sectional and longitudinal, only few methodological research have been considered in the case of informative nonrandom dropout. This could be due to the added layer of complexity that informative missing poses on the analysis of semicontinuous longitudinal data since it typically requires modeling the censoring process on its own and linking it to the two parts of the semicontinuous response.

Appreciating the paucity of research in this area, we developed a joint modeling approach under the framework of the likelihood-based mTP modeling of longitudinal semicontinuous responses. We assumed that the missing responses are nonrandom and resulting from the subject's dropout from the study. This type of subject's attrition is referred to as informative dropout or informative right censoring.³⁶ We assume that the semicontinuous response follows the lognormal as a special case of the generalized gamma distribution,³⁷ with probability of nonzero value following the logistic model. We model the informative censoring as (1) a positive count of the number of follow-up visits, which we assume to be discrete Poisson censoring and (2) a survival time-to-event or time-to-dropout censoring process. The two parts of the semicontinuous longitudinal response and the censoring process are modeled jointly in a likelihood function framework and are linked using two methods: (1) latent shared parameters random effects and (2) copula models. In the latent parameter approach, we included correlated subject-specific random intercepts and slopes that are shared between the two parts of the semicontinuous model and the censoring process. These latent random effect parameters link the different components of the likelihood function and account for the within-subject correlation that emanates from the longitudinal repeated measures on the same subject. In the copula approach, a joint distribution generated from the marginal distributions is used to link all the components of the likelihood function together. The subject-specific random effects were included to account for the within-subject correlation in the copula model.

The proposed approaches in this article have several methodological advances:

Assuming the generalized gamma with its specifications for the different emanating distributions for the semicontinuous response introduces more flexibility in the parametric part when compared to transformation function that could be complex leading to uninterpretable inferences.^{20,32,35,38,39}

Analyzing the semicontinuous longitudinal data using the mTP model instead of the TP mixed model is advantageous since mTP allows direct parametrization of the covariates in terms of the marginal mean and straightforward interpretation of the marginal effects of covariates on the population mean.¹ MTP is also robust for the misspecification of the true random effect structure in parameter estimation and to the degree of this misspecification.^1,33,35 Moreover, mTP is less sensitive to model misspecification in terms of estimating the covariates effects on the marginal population mean.³⁸

Relaxing the conditional independence assumption of the shared parameter approach by applying the copula approach for joint modeling of the two parts of the semicontinuous response and the censoring process is of added benefit. This is due to the fact that if an embedded serial correlation is present and the shared parameter approach is employed then this correlation is ignored. In this case, the shared parameter approach is more susceptible to bias. Nevertheless, the copula model relaxes the assumption of conditional independence by joining the marginal distributions into a common joint distribution which incorporates any embedded correlation leading to more accurate estimates.¹⁵

Along the same lines, compared to TP, the mTP model is less susceptible to bias when the serial correlation is present but ignored.¹⁵ Hence, modeling the semicontinuous response as mTP and joining the different components of the likelihood function using copula as proposed in two of our methodological approaches is advantageous since this methodology accounts for any embedded serial correlations between the stochastic processes and is therefore less prone to bias.

The censoring process proposed in two of our models as a positive count of visits similar to Li and Su,⁴⁰ could have an added benefit over the survival model since it overcomes the typical assumption of proportionality of hazards present in most survival models. In addition, given that in real-life situation the exact dropout time is not known, then registering the count of visits to reflect the length of stay before dropping out is more feasible than the time-to-dropout.

The proposed approaches include the Poisson and survival censoring with joint modeling based on the shared latent parameters and copula models. These approaches were assessed using extensive simulation studies and were illustrated using a longitudinal clinical study on a cohort of type 1 diabetic patients.⁴¹ Repeated measures in this cohort of patients were collected longitudinally until the diabetic subject experiences macroalbuminuria, a severe stage 3 kidney disease manifested by elevated levels of albumin excretion rate (AER) above 300 mg/24.^42,43 This complication ultimately leads to renal failure and introduces informative right censoring to this longitudinal clinical study. Clinical and sociodemographic parameters were measured longitudinally along with a fibrotic gene referred to as plasma connective tissue growth factor (CTGF) that is known to be associated with diabetic risk factors and diabetic complications. In particular, CTGF is an independent risk factor for increased levels of hemoglobin A1c (HbA1c) and hyperlipidemia,⁴² and cardiovascular events and mortality in patients with atherosclerotic disease.⁴⁴ Plasma CTGF is also an independent predictor for myocardial infarction in type 2 diabetes,⁴⁵ and renal and vascular disease in type 1 diabetic subjects.⁴² In addition, genetic variant in CTGF gene conferred a higher risk to develop nephropathy in type 1 diabetic subjects.⁴⁶

When CTGF gene is expressed its values are positive; however, when not expressed its values are zeros. This typically leads to preponderance of zeros that are mixed with positive continuous values resulting in semicontinuous data in the CTGF measured levels.

Factors that promote expression and release of CTGF gene or its inhibition are still not fully identified. The promoting factors result in the generation of positive CTGF values while the inhibitory factors attribute to the generation of zero values. Hence, given that this important clinical problem is still understudied and appreciating the crucial role CTGF displays as an independent risk factor for diabetic micro- and macro-vascular complications, we centered our clinical application on unraveling the factors that modulate expression of CTGF and regulate its levels. This clinical problem lends itself as fitting better in the context of assessing the direct covariate effect at a population-averaged level (mTP model) than subject specific (TP model). Subject-specific interpretation of gene expression is less attractive than the population-averaged level since more insight is gained as to unraveling the factors that relate to increased or decreased odds of gene marker expression and regulation. Hence, from the clinical standpoint of the genetic marker and methodological perspective, mTP models are more feasible and advantageous to be employed in our proposed framework of models. In Section 2, we present our framework of models which describes the joint modeling approaches based on shared parameters and copulas, and employing the Poisson and survival models for the missing process. In Section 3, we illustrate our framework of models and its feasibility of implementation using the real-life dataset on diabetic patients focusing on CTGF biomarker as an outcome of interest In Section 4, we provide the specifications for the simulation studies that we carried out and discuss the simulation results and performance of each model. In section 5, we present an overall discussion of our framework of models and provide concluding remarks on the future direction of research that can originate from this manuscript.

2. Models

We propose and compare approaches for analyzing longitudinal semicontinuous data with missingness that is specifically caused by nonrandom dropout leading to informative right censoring. The longitudinal semicontinuous outcome is modeled using the mTP approach with lognormal as a specification of the generalized gamma distribution.

The informative right censoring mechanism is incorporated using two options: (1) a discrete Poisson count of the number of visits for each individual subject and (2) time-to-event exponential survival model. Joint modeling was specified using the shared latent parameters, and the Clayton copula approaches.

We employed the generalized gamma as a distribution for the mTP model since it is a flexible three-parameter distribution (e.g. shape, location, and scale parameters) that has been used to model the continuous part of the semicontinuous outcome to relax the sometimes unrealistic condition of the symmetry of the log-transformed data.³⁷ As for the censoring process, we used the Poisson distribution when we modeled the censoring part as a positive count of the individual's number of visits in the study since Poisson is a typical distribution for positive discrete data. However, other discrete distributions can be employed if they provide a better fit for the data. In our other specification of the censoring process, the survival exponential hazard was assumed since it is a traditional survival model for time to event represented by time to drop out in our approach. The aim here is to assess the feasibility of the two modeling approaches of the censoring process; the survival exponential model for time to event, versus the Poisson model for the number of visits available for the individual.

The proposed framework of methodologies has a novelty consideration that lies in the joint modeling of the different components of the two parts of the longitudinal semicontinuous measures and the missingness process achieved using the shared parameters and copula models. Copula in particular represents a methodological advance since it has not been well studied due to its computational challenge⁴⁷ especially when employed for joint modeling of the semicontinuous longitudinal outcome (zero and nonzero parts), and the censoring process as was specified here.

2.1 Models’ general specifications

Let $Y_{i j}$ be a longitudinal outcome for subject i measured at the $j^{t h}$ time point with i = 1,…,n independent subjects. The measurements were taken at time points $(t_{i 1}, t_{i 2}, \dots, t_{i m_{i}})$ that do not need to be the same for all subjects. Every individual subject i is assumed to have $m_{i}$ number of follow-up visits before dropping out and $m_{i}$ is assumed to be $1 \leq m_{i} \leq p$ where p marks the end of study period. In specific, $m_{i}$ is defined as the number of follow-up visits which reflects the length of stay for subject i in the study before either dropping out or reaching termination of the study. In addition, every subject has a survival time $T_{i}$ that reflects the observed time spent in the study before dropping out.

2.2 Longitudinal semicontinuous mTP submodel

The longitudinal semicontinuous outcome $Y_{i j}$ is expressed using the lognormal mTP model. The lognormal distribution stems from the generalized gamma family hence we start first by describing the gamma model and then we indicate the specification that results in the lognormal distribution.

The generalized gamma probability density function (PDF) is specified as

f (y; k, μ, σ) = \frac{η^{η}}{σ y Γ (η) \sqrt{η}} \exp {u \sqrt{η} - η \exp (| k | u)},

(1)

where

Γ (.)

is the standard gamma function,

u = s i g n (k) (\log y - μ) / σ

, with shape parameter k, location parameter

μ > 0

, and scale parameter

σ > 0

and

η = | k |^{- 2} > 0

. When

k \to 0

the limiting distribution of the generalized gamma reduces to lognormal distribution.

The PDF for an mTP model $(g_{m T P})$ can be written as

g_{m T P} (Y_{i j}) = {\begin{matrix} 1 - π_{i j} i f Y_{i j} = 0 \\ π_{i j} f (Y_{i j}; {X^{'}}_{i j} {\underline{β}}_{C o n t}, {\underline{β}}_{i j}) i f Y_{i j} > 0. \end{matrix}

(2)

The two parts, binary and continuous, of the mTP model expressed in terms of linear predictors can be described as follows:

\begin{aligned} \log i t (π_{i j}) = β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j} \\ and \\ E (Y_{i j}) = \exp (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{C o n t} {X^{'}}_{i}), \end{aligned}

(3)

where

π_{i j}

is the probability of nonzero value for the outcome

Y_{i j}

and is obtained from the logistic model expressed in the binary part of the mTP. Thus

π_{i j}

will take the form of

\begin{aligned} π_{i j} = & expit ((β_{0} + b_{0 i}) + (β_{1} + b_{1 i}) t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j}) \\ = & \frac{\exp ((β_{0} + b_{0 i}) + (β_{1} + b_{1 i}) t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}{1 + \exp ((β_{0} + b_{0 i}) + (β_{1} + b_{1 i}) t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})} \\ = & \frac{\exp (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}{1 + \exp (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}, \end{aligned}

(4)

{\underline{β}}_{B i n}

represents the vector of model coefficients corresponding to the vector of predictors

{\underline{θ}}_{i j}

in the binary part, and

{\underline{β}}_{C o n t}

represents the vector of marginal coefficients corresponding to the vector of predictors

X_{i j}^{'}

of the continuous part of the mTP model.

{\underline{β}}_{i j}

is the vector of individual subject random effects that can either have the intercept

β_{0 i}

and/or the slope

β_{1 i}

. The random effects are shared between the zero and nonzero components of the semicontinuous outcome

Y_{i j}

to account for the cross-parts correlations. These random effects are typically assumed to follow the bivariate normal distribution with mean dependent on the population intercept

β_{0}

and slope

β_{1}

, and variance–covariance matrix G, and can be described as follows:

{\underline{β}}_{i} = [\begin{matrix} β_{0 i} \\ β_{1 i} \end{matrix}] = [\begin{matrix} β_{0} + b_{0 i} \\ β_{1} + b_{1 i} \end{matrix}] \sim B V N (\underline{β} = [\begin{matrix} β_{0} \\ β_{1} \end{matrix}], G = [\begin{matrix} σ_{0}^{2} & σ_{01} \\ σ_{01} & σ_{1}^{2} \end{matrix}]) .

(5)

The marginal mean under the mTP generalized gamma (GG) model can be expressed as

E (Y_{i j}) = π_{i j} \exp (μ_{i j} + C (σ, k)) = \exp (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{C o n t} {X^{'}}_{i}) .

(6)

Solving for the location parameter

μ_{i j}

expressed in

E (Y_{i j})

for the GG distribution gives the following form:

μ_{i j} = β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{C o n t} X_{i}^{'} - \log (π_{i j}) - C (σ, k),

(7)

and

C (σ, k) = \frac{σ \log (k^{2})}{k} + \log [Γ (\frac{1}{k^{2}} + \frac{σ}{k})] - \log [Γ (\frac{1}{k^{2}})] .

(8)

The distribution used in the continuous nonzero part of mTP should have a finite closed-form mean that can be parameterized as in equation (7). When

C (σ) = σ^{2} / 2

and

k \to 0

the mTP lognormal distribution is obtained with

μ_{i j} = β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{C o n t} X_{i}^{'} - \log (π_{i j}) - \frac{σ^{2}}{2}

(9)

The binary part provides subject-specific estimates of the probability of having nonzero values for the outcome wherein the exponential of

{\underline{β}}_{B i n}

is interpreted as the subject-specific odds ratio for having a nonzero response attributed to a one unit increase in the respective covariate. The exponential of the parameter

{\underline{β}}_{C o n t}

in the continuous part represents the multiplicative effect on the overall mean for the whole population attributed to a one unit increase in the corresponding covariate X.

The general format of the likelihood function for mTP longitudinal models can be described as such,

\begin{aligned} L (π, μ, k, {\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, σ, G | y) = \underset{i}{Π} (1 - π_{i j} | θ_{i}, {\underline{β}}_{i})^{I (Y_{i j} = 0)} \\ {(π_{i j} | θ_{i}, {\underline{β}}_{i}) f (Y_{i j}; k, μ_{i j}, σ, {\underline{β}}_{i})}^{I (Y_{i j} > 0)} q ({\underline{β}}_{i} | \underline{β}, G), \end{aligned}

(10)

where f represents the PDF of the generalized gamma distribution or any other distribution from its family and q is the bivariate normal distribution for the random effects. Expressing

μ_{i j}

in terms of equation (7), and

π_{i j}

as denoted earlier in equation (4), the marginal likelihood function for the generalized gamma distribution can be described as

\begin{aligned} L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G | y) = \prod_{i} \int {(1 - \frac{1}{1 + \exp {- (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}})}^{I (Y_{i j} = 0)} * \\ {\frac{{| k |}^{- 2 {| k |}^{- 2}}}{(1 + \exp {- (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}) σ Y_{i j} Γ ({| k |}^{- 2}) \sqrt{{| k |}^{- 2}}} \exp {Q ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, G, k, σ)}}^{I (Y_{i j} > 0)} \\ q ({\underline{β}}_{i} | \underline{β}, G) d β_{i}, \end{aligned}

(11)

where

\begin{aligned} Q ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, G, k, σ) \\ = s i g n (k) (\log Y_{i j} - (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{C o n t} {X^{'}}_{i} - \log (\exp i t (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j}) - C (σ, k)))) / σ \sqrt{{| k |}^{- 2}} \\ - | k |^{- 2} \exp (k (\log Y_{i j} - (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{C o n t} {X^{'}}_{i} - \log (\exp i t (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})) - C (σ, k))) / σ) . \end{aligned}

(12)

In this likelihood function equation (11),

q ({\underline{β}}_{i} | \underline{β}, G)

represents the conditional bivariate normal distribution for the random effects given the mean vector of

\underline{β}

and variance–covariance matrix G. The likelihood function for the lognormal distribution is obtained by replacing the distribution of the continuous nonzero part with the corresponding PDF achieved by setting

C (σ) = σ^{2} / 2

and by having

k \to 0

2.3 Models for the missingness process

The censoring mechanism was modeled using two approaches, the Poisson censoring for the number of visits and the time-to-event survival model.

2.3.1 Poisson submodel

The censoring process is determined by the number of the last visit $m_{i}$ prior to dropout and is assumed to be dependent on the trajectories of the longitudinal outcome in addition to a set of covariates denoted as singledollarX_Cens singledollar. Hence the censoring mechanism can be described as : singledollar {m_i} \vert \underline{\beta } _isingledollar ∼ Truncated Discrete Distribution with

Pr (M = m_{i}) = F (γ_{0} + {\underline{γ}}^{'} {\underline{β}}_{i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s}) = F (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s}),

(13)

with

\underline{γ} = [\begin{matrix} γ_{1} \\ γ_{2} \end{matrix}]

and

{\underline{γ}}_{3}

is a vector of censoring parameters corresponding to the vector of covariates

{\underline{X}}_{i C e n s}

. Here we specify that conditional on the random effects

{\underline{β}}_{i}

, the number of follow-up visits

m_{i}

is assumed to follow a discrete distribution which is assumed to be the Poisson distribution with probability dependent on the individual's intercept and slope. The censoring process could also be dependent on a set of covariates

X_{i C e n s}

, but this is not a required assumption. The censoring parameters are denoted as γ₀, γ₁,

γ_{2}

, and

{\underline{γ}}_{3}

. Censoring is considered informative with respect to the longitudinal trajectories of the outcome of interest when either γ₂≠0 or

γ_{1} \neq 0

, or both are nonzero. However, censoring is non-informative when these parameters are both zero. If the

γ

parameter corresponding to a specific covariate

X_{C e n s}

is significant then this covariate has an effect on the censoring process otherwise it does not. The distribution of

m_{i}

is right-truncated since the study may terminate before observing the withdrawal of every subject and left truncated since the individual subject needs to have at least one measurement.

The Poisson probability function can be expressed as follows:

Pr [M = m_{i}] = \frac{μ_{i}^{(m_{i} - 1)} e^{- μ_{i}}}{(m_{i} - 1)! * Pr (1 \leq m_{i} \leq p)},

(14)

with

m_{i} = 1, 2, \dots, p

and mean

μ_{i} = \exp (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})

2.3.2 Survival time-to-event submodel

Let $s_{i}$ be the survival time for the $i^{t h}$ subject where i = 1,…,n independent subjects, $c_{i}$ the censoring time, and $T_{i} = min (s_{i}, c_{i})$ the observed survival time. Let $δ_{i}$ be an indicator variable such that $δ_{i} = 1$ if event is observed and in this case $T_{i} = s_{i}$ , and $δ_{i} = 0$ if censoring occurred and $T_{i} = c_{i}$ in this case. The survival function is expressed as

S_{i} (T_{i}) = \exp {- \int_{0}^{T_{i}} h_{i} (x) d x} = \exp {- Λ (T_{i})},

(15)

with

h_{i}

denoting the hazard function that is expressed as

h = h_{0} (t) \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} .

(16)

The baseline hazard function

h_{0} (t)

for simplicity purposes is assumed to follow the exponential hazard with rate

h_{0}

. Under the proportional hazard model with exponential distribution

h_{0}

reduces to 1 and in this case the hazard function can take this form,

h = \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} .

(17)

Λ (T_{i}) = \int_{0}^{T_{i}} h_{i} (x) d x

represents the cumulative hazard which turns out to be the integral of a constant under the exponential model and reduces therefore to the following formulation,

Λ (T_{i}) = h_{i} T_{i}

The survival function for the $i^{t h}$ subject can be expressed as

S_{i} (T_{i}) = \exp {- h_{i} T_{i}} .

(18)

The density function is obtained by multiplying the survival function by the hazard

g_{i} (T_{i}) = h_{i} * S_{i} = h_{i} \exp {- h_{i} T_{i}}

(19)

The conditional survival function, when expressed in terms of the censoring parameters

\underline{γ}

, will take the following formulation:

S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) = \exp {- \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} * T_{i}},

(20)

and the conditional density function is represented as such,

\begin{aligned} g_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) = \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} \\ * \exp {- \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} * T_{i}} . \end{aligned}

(21)

The likelihood function corresponding to this survival model can be expressed as follows:

\begin{aligned} L (\underline{β}, \underline{γ}, G) = & \prod_{i = 1}^{n} {\int_{- \infty}^{+ \infty} [g_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ})]}^{δ_{i}} [S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ})]^{1 - δ_{i}} q ({\underline{β}}_{i} | \underline{β}, G) d {\underline{β}}_{i} \\ L (\underline{β}, \underline{γ}, G) = & \prod_{i = 1}^{n} \int_{- \infty}^{+ \infty} h_{i}^{δ_{i}} S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) q ({\underline{β}}_{i} | \underline{β}, G) d {\underline{β}}_{i} \\ L (\underline{β}, \underline{γ}, G) = & \prod_{i = 1}^{n} {\int_{- \infty}^{+ \infty} [\exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})}]}^{δ_{i}} S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) q ({\underline{β}}_{i} | \underline{β}, G) d \underline{β} \\ L (\underline{β}, \underline{γ}, G) = & \prod_{i = 1}^{n} {\int_{- \infty}^{+ \infty} [\exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})}]}^{δ_{i}} \\ \exp {- \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} * T_{i}} q ({\underline{β}}_{i} | \underline{β}, G) d {\underline{β}}_{i} . \end{aligned}

(22)

The function q denotes the bivariate normal distribution of the random parameters described in equation (5).

2.4 Joint modeling approaches

The two parts of the longitudinal semicontinuous response and the censoring process were all jointly modeled in a likelihood function using the shared parameter and the copula models.

2.4.1 Shared parameters

The longitudinal semicontinuous trajectories and censoring process are jointly modeled using shared parameters. Hence, the correlation between the longitudinal and survival processes is induced by having these shared common random parameters modeled in the marginal likelihood function. Given the random effects, conditional independence is assumed between the different marginal components.³⁴ Hence, conditional marginal densities are jointly modeled in the marginal likelihood function. The two censoring processes considered for joint modeling with the longitudinal outcome include (a) Poisson censoring and (b) survival censoring model.

(a) Joint modeling of the longitudinal and Poisson censoring process using shared parameters: The number of times subject i is observed $(m_{i})$ and the semicontinuous outcome $({\underline{Y}}_{i})$ are jointly modeled using the shared parameter approach as in.^48,49 The two submodels share the latent random variables (random intercept and/or slope) with at least one of these subject-specific random effects being common to both models. The covariates are not necessarily common between the two submodels. Assuming conditional independence between ${\underline{Y}}_{i}$ and $m_{i}$ given the random effect $\underline{β_{i}}$ , the shared parameter model is obtained as,

\begin{aligned} ϕ (Y_{i}, m_{i}, {\underline{β}}_{i} | {\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3}, {\underline{X}}_{i}, {\underline{θ}}_{i}, {\underline{X}}_{i C e n s}) \\ = f (Y_{i} | \underline{β}, {\underline{β}}_{i}, {\underline{θ}}_{i}, {\underline{X}}_{i}) w (m_{i} | {\underline{β}}_{i}, \underline{β}, {\underline{X}}_{i C e n s}, γ_{0}, \underline{γ}, {\underline{γ}}_{3}) q ({\underline{β}}_{i} | \underline{β}, G) . \end{aligned}

(23)

For simplicity of notation let

ϕ (Y_{i}, m_{i}, {\underline{β}}_{i} | {\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3}, {\underline{X}}_{i}, {\underline{θ}}_{i}, {\underline{X}}_{i C e n s})

be expressed as

φ (Y_{i}, m_{i}, {\underline{β}}_{i})

, the log likelihood for the shared-parameter model can be described as,

\begin{aligned} L = L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3} | y) = \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} φ ({\underline{Y}}_{i}, m_{i}, {\underline{β}}_{i}) d {\underline{β}}_{i} . \\ Assuming conditional independence the likelihood function of the shared parameter model \\ takes the following form \\ L = L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3} | y) = \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} f ({\underline{Y}}_{i} | {\underline{β}}_{i}, m_{i}) w (m_{i} | {\underline{β}}_{i}) q ({\underline{β}}_{i} | β, G) d {\underline{β}}_{i} \\ L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G | y) = \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} {(1 - \frac{1}{1 + \exp {- (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}})}^{I (Y_{i j} = 0)} * \\ {\frac{{| k |}^{- 2 {| k |}^{- 2}}}{(1 + \exp {- (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}) σ Y_{i j} Γ ({| k |}^{- 2}) \sqrt{{| k |}^{- 2}}} \exp {Q ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, G, k, σ)}}^{I (Y_{i j} > 0)} \\ * \frac{μ_{i}^{(m_{i} - 1)} e^{- μ_{i}}}{(m_{i} - 1)! * Pr (1 \leq m_{i} \leq p)} \\ * \exp (- \frac{1}{2} [\frac{1}{1 - ρ_{β_{0 i}, β_{1 i}}^{2}} {{(\frac{β_{0 i} - β_{0}}{σ_{0}})}^{2} - 2 ρ_{β_{0 i}, β_{1 i}} (\frac{β_{0 i} - β_{0}}{σ_{0}}) (\frac{β_{1 i} - β_{1}}{σ_{1}}) + {(\frac{β_{1 i} - β_{1}}{σ_{1}})}^{2}}]) d {\underline{β}}_{i}, \end{aligned}

(24)

with

m_{i} = 1, 2, \dots, p

and mean

μ_{i} = \exp (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})

and

Q ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, G, k, σ)

is defined in equation (12).

(b) Joint modeling of the longitudinal and survival censoring processes using shared parameters assumes the following likelihood function:

\begin{aligned} L = L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3} | y) = \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} f ({\underline{Y}}_{i}, T_{i}, {\underline{β}}_{i}) d {\underline{β}}_{i} . \\ Assuming conditional independence the likelihood function take the following form \\ L = L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3} | y) = \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} f ({\underline{Y}}_{i} | {\underline{β}}_{i}, T_{i}) w_{s u r v i v a l} (T_{i} | {\underline{β}}_{i}) q ({\underline{β}}_{i} | β, G) d {\underline{β}}_{i} \\ L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G | y) = \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} {(1 - \frac{1}{1 + \exp {- (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}})}^{I (Y_{i j} = 0)} \\ {\frac{{| k |}^{- 2 {| k |}^{- 2}}}{(1 + \exp {- (β_{0 i} + β_{1 i} t_{i j} + {\underline{β}}_{B i n} {\underline{θ}}_{i j})}) σ Y_{i j} Γ ({| k |}^{- 2}) \sqrt{{| k |}^{- 2}}} \exp {Q ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, G, k, σ)}}^{I (Y_{i j} > 0)} \\ * [\exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})}]^{δ_{i}} * \exp {- \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} * T_{i}} \\ * \exp (- \frac{1}{2} [\frac{1}{1 - ρ_{β_{0 i}, β_{1 i}}^{2}} {{(\frac{β_{0 i} - β_{0}}{σ_{0}})}^{2} - 2 ρ_{β_{0 i}, β_{1 i}} (\frac{β_{0 i} - β_{0}}{σ_{0}}) (\frac{β_{1 i} - β_{1}}{σ_{1}}) + {(\frac{β_{1 i} - β_{1}}{σ_{1}})}^{2}}]) d {\underline{β}}_{i}, \end{aligned}

(25)

where

Q ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, G, k, σ)

is defined in equation (12). The survival component in this likelihood function (eq. 25) which we denote as

w_{s u r v i v a l}

in the below derivation is obtained as follows:

\begin{aligned} w_{s u r v i v a l} (T_{i} | {\underline{β}}_{i}) = [g_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ})]^{δ_{i}} [S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ})]^{1 - δ_{i}} \\ replacing g_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) and S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) by their corresponding functions \\ g_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) = \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} * \exp {- \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} * T_{i}} \\ g_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) = h_{i} * S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) \\ with h_{i} = \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} \\ and S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) = \exp {- \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} * T_{i}}, \\ then the survival component w_{s u r v i v a l} (T_{i} | {\underline{β}}_{i}) can be expressed as \\ w_{s u r v i v a l} (T_{i} | {\underline{β}}_{i}) = h_{i}^{δ_{i}} * S_{i} (T_{i} | {\underline{β}}_{i}, \underline{γ}) \\ w_{s u r v i v a l} (T_{i} | {\underline{β}}_{i}) = [\exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})}]^{δ_{i}} * \exp {- \exp {- (γ_{0} + γ_{1} β_{0 i} + γ_{2} β_{1 i} + {\underline{γ^{'}}}_{3} {\underline{X}}_{i C e n s})} * T_{i}} . \end{aligned}

2.4.2 Copula approach

Copula is another approach for joint modeling that relaxes the assumption of conditional independence by generating a joint distribution from the marginal functions. Copulas were originally derived in the theorem by Sklar⁵⁰ which we present in Appendix 1 of this manuscript. Our focus here is on the bivariate copula which represents a joint cumulative distribution function that combines the marginal distributions of two random variables into one bivariate function and incorporates the dependency between the two random variables through a correlation structure. The two random variables in our proposed framework of models are the longitudinal outcome which is highly skewed and the missing process presented by (1) the Poisson and (2) the survival exponential time to event random variables. There are different types of copulas that include but not limited to Gaussian, Clayton, and Gumbel functions among others.

In our framework of models, we employed the Clayton copula since it does not require the marginal and joint distributions to be Gaussian, but rather accommodates random variables that are skewed. This characteristic makes Clayton more applicable for situations similar to what we have in this paper where the marginal and joint distributions of the two random variables (longitudinal and missing processes) have complex non-Gaussian functions. Moreover, Clayton assumes dependence on the lower tail of the joint distribution of the two random variables. Given the proliferation of zero values and the anticipated interdependency between these small values and the decreased risk of events, we expect the correlation to be concentrated on the lower tail of the joint distribution. This assumption makes Clayton copula a reasonable function for combining the two random processes in our models, and provides another rationale to employing it in our models’ specifications. A more detailed discussion on copulas is included in Appendix 2 of this paper.

The likelihood functions corresponding to the Clayton copula joint modeling of the longitudinal semicontinuous outcome and the censoring process using (a) Poisson and (b) survival models are defined as follows:

Poisson censoring with longitudinal semicontinuous outcome using Clayton copula:

\begin{aligned} L = & L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3}, ρ | y) \\ = & \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} (1 + ρ) {[F ({\underline{Y}}_{i} | {\underline{β}}_{i}, m_{i}) W (m_{i} | {\underline{β}}_{i})]}^{- ρ - 1} {(F {({\underline{Y}}_{i} | {\underline{β}}_{i}, m_{i})}^{- ρ} W {(m_{i} | {\underline{β}}_{i})}^{- ρ} - 1)}^{\frac{- 1}{ρ - 2}} \\ * f ({\underline{Y}}_{i} | {\underline{β}}_{i}, m_{i}) w (m_{i} | {\underline{β}}_{i}) q ({\underline{β}}_{i} | β, G) d {\underline{β}}_{i} \\ L = & L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3}, ρ | y) \\ = & \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} (1 + ρ) {[\int_{0}^{{\underline{Y}}_{i}} f (y | {\underline{β}}_{i}, m_{i}) d y \sum_{0}^{m_{i}} w (m | {\underline{β}}_{i})]}^{- ρ - 1} {(\int_{0}^{{\underline{Y}}_{i}} f (y | {\underline{β}}_{i}, m_{i}) d y^{- ρ} \sum_{0}^{m_{i}} w {(m | {\underline{β}}_{i})}^{- ρ} - 1)}^{\frac{- 1}{ρ - 2}} \\ * f ({\underline{Y}}_{i} | {\underline{β}}_{i}, m_{i}) w (m_{i} | {\underline{β}}_{i}) q ({\underline{β}}_{i} | β, G) d {\underline{β}}_{i}, \end{aligned}

(26)

ρ

is the dependence parameter between the marginal functions of the semicontinuous longitudinal and Poisson censoring processes. F and W represent the cumulative distribution functions of the semicontinuous lognormal outcome and the Poisson censoring. f and w represent their respective density functions and are defined in equation (24), and q is the bivariate normal density of the random effect

{\underline{β}}_{i}

Survival censoring with longitudinal semicontinuous outcome using Clayton copula:

\begin{aligned} L = & L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3}, ρ | y) \\ = & \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} (1 + ρ) {[F ({\underline{Y}}_{i} | {\underline{β}}_{i}, T_{i}) W_{s u r v i v a l} (T_{i} | {\underline{β}}_{i})]}^{- ρ - 1} {(F {({\underline{Y}}_{i} | {\underline{β}}_{i}, T_{i})}^{- ρ} W_{s u r v i v a l} {(T_{i} | {\underline{β}}_{i})}^{- ρ} - 1)}^{\frac{- 1}{ρ - 2}} \\ * f ({\underline{Y}}_{i} | {\underline{β}}_{i}, T_{i}) w_{s u r v i v a l} (T_{i} | {\underline{β}}_{i}) q ({\underline{β}}_{i} | β, G) d {\underline{β}}_{i} \\ L = & L ({\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3}, ρ | y) \\ = & \prod_{i = 1}^{n} \int_{{\underline{β}}_{i}} (1 + ρ) {[\int_{0}^{{\underline{Y}}_{i}} f (y | {\underline{β}}_{i}, T_{i}) d y \int_{0}^{T_{i}} w_{s u r v i v a l} (t | {\underline{β}}_{i}) d t]}^{- ρ - 1} {(\int_{0}^{{\underline{Y}}_{i}} f (y | {\underline{β}}_{i}, m_{i}) d y^{- ρ} \int_{0}^{T_{i}} w_{s u r v i v a l} (t | {\underline{β}}_{i}) d t^{- ρ} - 1)}^{\frac{- 1}{ρ - 2}} \\ * f ({\underline{Y}}_{i} | {\underline{β}}_{i}, T_{i}) w_{s u r v i v a l} (T_{i} | {\underline{β}}_{i}) q ({\underline{β}}_{i} | β, G) d {\underline{β}}_{i} . \end{aligned}

(27)

The survival component

w_{s u r v i v a l}

of the observed survival time

T_{i}

is defined in equation (25).

2.5 Statistical estimation and inferences

The random effects ${\underline{β}}_{i}$ are integrated out to get the marginal likelihood function. The log of the marginal likelihood function is then maximized by taking the first derivative with respect to each parameter and setting the equation to zero to obtain the maximum likelihood estimate for each of the fixed effects, ${\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G, γ_{0}, \underline{γ}, {\underline{γ}}_{3}, and ρ$ , respectively. Estimates of the fixed effect for $ρ$ pertain to the copula models only. Direct maximization of the log of the marginal likelihood function pertaining to the joint model is implemented using dual quasi-Newton optimization algorithm⁵¹ to generate consistent maximum likelihood estimates of the fixed effects. Empirical Bayes estimators using the adaptive Gaussian quadrature approach described in Pinherio and Bates⁵² was employed to obtain predicted values of the random effects ${\underline{β}}_{i}$ . The asymptotic standard errors are computed using Fisher information after substituting the maximum likelihood estimates for ${\underline{β}}_{C o n t}, {\underline{β}}_{B i n}, \underline{β}, k, σ, G,, γ_{0}, \underline{γ}, {\underline{γ}}_{3}, and ρ$ , and can be defined as follows:

V a r ({\underline{\hat{β}}}_{C o n t}, {\underline{\hat{β}}}_{B i n}, \underline{\hat{β}}, \hat{k}, \hat{σ}, \hat{G}, {\hat{γ}}_{0}, \hat{\underline{γ}}, {\underline{\hat{γ}}}_{3}, and \hat{ρ}) = d i a g {I^{- 1} ({\underline{\hat{β}}}_{C o n t}, {\underline{\hat{β}}}_{B i n}, \underline{\hat{β}}, \hat{k}, \hat{σ}, \hat{G}, {\hat{γ}}_{0}, \hat{\underline{γ}}, {\underline{\hat{γ}}}_{3}, and \hat{ρ})} .

(28)

Maximization and optimization of the likelihood function in each of the models were carried out using NLMIXED SAS procedure (SAS Institute, Cary, North Carolina, USA). All the SAS code and procedures developed to implement the proposed framework of models can be made available upon request from the authors.

3. Data example

The proposed methods were illustrated using a longitudinal cohort of type 1 diabetic patients wherein clinical parameters and genetic markers were measured repeatedly over a period of 10 years as described below. The levels of the plasma CTGF biomarker were measured on 691 subjects who participated in the Diabetes Control and Complications Trial (DCCT) cohort of type 1 diabetes.⁴¹ Subjects were enrolled in the DCCT study between 1983 and 1989, and half of them were randomized to the conventional and the other half to the intensive diabetes treatment arms. In 1993, the DCCT study was stopped when the intensive diabetes treatment was clearly shown to reduce the risks of microvascular complications.⁵³ Clinical and demographic factors were collected on this cohort of patient population and were used as covariates in our analysis to assess their differential effects on CTGF levels for the zero and nonzero parts. Plasma CTGF levels were measured longitudinally at baseline (study entry (1983–1989)), midpoint of DCCT (1988–1991), and end of DCCT (1993) by ELISA.²¹ Out of a total repeated measures of CTGF on all subjects throughout the study (n = 1985), 62% of the CTGF levels had zero values. This indicates that there was no CTGF gene expression due to an inhibition in its production or in its release into the plasma. Patients were followed until they experienced macroalbuminuria and ultimately dropped out of the study due to this severe kidney damage. This complication manifested in about 2% of the patients between the years 1983 and 1993, and therefore an incomplete panel of follow-up data and informative missing values were introduced to the recorded information on these patients.

Our proposed models were illustrated using this cohort of diabetic patients with CTGF (ng/ml) as the outcome of interest that has a semicontinuous nature. The assumption of lognormal distribution for the nonzero values of CTGF specified in the model's section (Section 2.2) was checked and compared to other distributions using the respective Q–Q plots. The distributions that were considered included normal, exponential, and gamma functions among others. Lognormal showed the best fit for this data with the majority of the observations falling on the reference line representing the expected values under this distribution.

Predictors of interest included treatment group (intensive glucose control vs. conventional group), age (in years), gender (male vs. female), smoking (binary yes/no), and duration of diabetes (in years), in addition to clinical parameters that encompassed HbA1c (in % units), systolic blood pressure (SBP) (mmHg), and lipidemia expressed in terms of high-density lipoprotein (HDL in mg/dl). Macroalbuminuria is assumed to occur when the AER levels exceed 300 mg/24. In Figure 1, we present the frequency distribution of the log scale levels of AER stratified by an indicator variable that is equal to one when CTGF is expressed (nonzero CTGF value) and zero otherwise, in addition to the corresponding summary statistics that we presented in terms of mean ± SE for the CTGF zero and nonzero groups respectively. The AER summary statistics pertaining to CTGF groups zero and nonzero were 14.66 ± 0.575 and 16.45 ± 0.904, respectively. The graphical display (Figure 1) and its corresponding summary statistics suggest a potential increase in the values of AER in the expressed nonzero group of CTGF compared to the non-expressed zero group. This could indicate that when CTGF levels are increased, AER could also increase which might put the diabetic patients at a greater risk of developing macroalbuminuria. A further conclusive assessment of the effect of CTGF on macroalbuminuria was achieved in our clinical illustration of the framework of models on this cohort of type 1 diabetic patients.

However, before proceeding with the discussion of results of our clinical application, we would like to comment on the time that each model took to converge, so that we can assess the efficiency of utilization when employed on a real-life dataset. In this regard, we noted that the shared parameters with survival censoring converged in 30.94 s, shared parameters with Poisson censoring in 55.85 s, copula model with survival censoring in 13.51 s, and copula model with Poisson censoring in 55.30 s. Accordingly, Poisson censoring needed more than double the time to converge compared to the survival models. But irrespectively, all 4 approaches successfully achieved convergence in a feasible time that was 60 s or less on average suggesting utility and efficiency of employment of the proposed framework of models on a real-life clinical application.

Now moving to the discussion of our results (Table 1), we can note here that there was a consistency in the inferences and estimates across the different models with copula except for some of the estimates under the models with shared parameters for the survival and Poisson censoring.

Table 1.

Parameter estimates for marginalized two-part (mTP^a) model assuming lognormal distribution for the nonzero component using (1) likelihood estimator assuming copula joint modeling between the longitudinal and the Poisson censoring processes; (2) likelihood estimator assuming shared parameters for joint modeling between the longitudinal and the Poisson censoring processes; (3) likelihood estimator assuming copula joint modeling between the longitudinal and the survival censoring processes; and (4) likelihood estimator assuming shared parameters for joint modeling between the longitudinal and the survival censoring processes.

Model Component	Covariate	(1) Copula Poisson censoring: Parameter estimate, (SE), P-value	(2) Shared parameters Poisson censoring: Parameter estimate, (SE), P-value	(3) Copula survival censoring: Parameter estimate, (SE), P-value	(4) Shared parameters survival censoring: Parameter estimate, (SE), P-value
Zero part	Txt group	−0.4020, (0.1299), 0.0021	−0.2987, (0.1083), 0.0060	−0.3099, (0.1198), 0.0099	−0.0759, (0.0902), 0.4005
Zero part	Smoking	0.7989, (0.1644), <0.0001	0.7941, (0.1752), <0.0001	0.7934, (0.1773), <0.0001	0.7189, (0.2269), 0.0016
Continuous nonzero part	HbA1c	0.2011, (0.0446), <0.0001	0.1974, (0.0519), 0.0002	0.2809, (0.0340), <0.0001	0.0628, (0.0313), 0.0453
	Age	0.0070, (0.0110), 0.5229	0.0049, (0.0126), 0.6965	0.0017, (0.0101), 0.8652	−0.0070, (0.0078), 0.3748
	Duration	0.0027, (0.0156), 0.8605	0.0010, (0.0203), 0.9597	0.0125, (0.0146), 0.3950	−0.0040, (0.0119), 0.7353
	SBP	0.0210, (0.0048), <0.0001	0.0155, (0.0053), 0.0042	0.0099, (0.0043), 0.0233	0.0219, (0.0034), <0.0001
	Male	−2 × 10⁻⁵, (0.1766) 0.9999	−8 × 10⁻⁴, (0.1991), 0.9967	0.0083, (0.1604), 0.9586	0.0068, (0.1243), 0.9560
	HDL	−0.0308, (0.0056), <0.0001	−0.0074, (0.0071), 0.3020	−0.0151, (0.0054), 0.0059	−0.0094, (0.0044), 0.0326
Censoring process	Gamma0	0.3978, (0.4192), 0.3430	0.2397,(0.0267),<0.0001	1.1228, (0.0286), <0.0001	2.2065, (0.0621), <0.0001
	Gamma1_Beta0i	−1.9989, (0.1042), <0.0001	−0.2325, (0.0201), <0.0001	−0.5019, (0.0021), <.0001	−1.4735, (0.0854), <0.0001
	Gamma2_Beta1i	−0.4996, (0.0649), <0.0001	−0.2129, (0.0299), <0.0001	−0.2752, (0.0262), <0.0001	−0.4484, (0.2280), 0.0496
Population shared parameters	Beta0	−0.0998, (0.2195), 0.6495	−0.2352, (0.0087), <0.0001	−0.1704, (0.2815), 0.5452	−1.1680, (0.0222), <0.0001
Population shared parameters	Beta1_Time	−0.1092, (0.0217), <0.0001	−0.0817, (0.0110), <0.0001	−0.1237, (0.0357), 0.0006	0.0196, (0.0172), 0.2541

mTP model provides estimates for the parameters in the continuous part for the entire sample (zero and nonzero values).

HbA1c: hemoglobin A1c; HDL: high-density lipoprotein.

3.1 Zero part

We start our discussion by considering the results pertaining to the zero part of the semicontinuous outcome CTGF. In this regard, our results showed a significant association between the treatment group and the probability of nonzero values for CTGF. This was demonstrated under the following three models, Poisson with copula (OR = exp(−0.402) = 0.66, p-value = 0.0021), Poisson with shared parameters for joint modeling (OR = exp(−0.2987) = 0.74, p-value = 0.006), and survival censoring with copula (OR = exp(−0.3099) = 0.73, p-value = 0.0099). From this finding, it can be deduced that patients who were in the intensive treatment group for glucose control had significantly lower odds of about 0.7 for having nonzero values for CTGF. This indicates that CTGF levels were decreasing for the intensive glucose control group compared to the conventional treatment group. This result is in line with clinical findings which showed that intensive glucose control regulates and lowers CTGF levels.^42,54 This conclusion was consistent across all the proposed models except for the survival censoring with shared parameters model that failed to show this significant association between the treatment group and the probability of nonzero values for CTGF.

The association between smoking and the probability of nonzero was consistently denoted in all four models. Poisson censoring with copula (OR = exp(0.7989) = 2.22, p-value < 0.0001) and Poisson with shared parameters (OR = exp(0.7941) = 2.21, p-value < 0.0001) for joint modeling, survival censoring with copula (OR = exp(0.7934) = 2.21. p-value < 0.0001) and survival censoring with shared parameters (OR = exp(0.7189) = 2.05, p-value = 0.0016) for joint modeling, all showed a significant association between smoking and the probability of nonzero. These results indicate that the odds of nonzero for smokers are twice that of nonsmokers. This implicates that the levels of CTGF are higher among smokers versus nonsmokers, a conclusion that is in agreement with clinical findings that showed increased levels of CTGF in pulmonary vessels among smokers.⁵⁵

3.2 Continuous nonzero part

As for the continuous nonzero part, our results showed consistency in the estimates and inferences for the predictors of the continuous part of CTGF values across the different models. Similar to the zero-part, some discrepancy was denoted in the shared parameter models.

Predictors that were consistently shown to be significantly associated with the continuous part of the CTGF included HbA1c, SBP, and HDL. In specific, HbA1c was shown to be significantly associated with the continuous nonzero part of CTGF under the Poisson censoring model with copula (effect estimate = exp(0.2011) = 1.22, p-value < 0.0001), Poisson censoring with shared parameter (effect estimate = exp(0.1974) = 1.22, p-value = 0.0002) for joint modeling, survival model for censoring with copula (effect estimate = exp(0.2809) = 1.324, p-value < 0.0001), and survival censoring with shared parameters (effect estimate = exp(0.0628) = 1.065, p-value = 0.0453) for joint modeling. Accordingly, it can be deduced that 1% increase in HbA1c contributes to an increase in the CTGF values (ng/ml) that ranges between 6.5% and 32.4%. The percent increase was obtained as (exp(0.0628)−1) × 100% = (1.065−1) × 100% = 6.5%, and (exp(0.2809)−1) × 100% = (1.324−1) × 100% = 32.4%.

Our results also showed that SBP was significantly associated with the continuous nonzero part of the CTGF values under the Poisson censoring with copula (effect estimate = exp(0.0210) = 1.02, p-value < 0.0001), Poisson with shared parameters (effect estimate = exp(0.0155) = 1.015, p-value = 0.0042) for joint modeling, survival model with copula (effect estimate = exp(0.0099) = 1.01, p-value = 0.0233), and survival model with shared parameters (effect estimate = exp(0.0219) = 1.02, p-value < 0.0001) for joint modeling. These results suggest a 2% (ng/ml) increase in CTGF levels was attributed to a 1 mmHg increase in SBP.

In addition, our results suggested that there was a consistency in the estimates and inferences pertaining to the effect of HDL on the continuous values of CTGF in the Poisson censoring with copula for joint modeling (effect estimate = exp(-0.0308) = 0.969, p-value < 0.0001), survival censoring with copula (effect estimate = exp(−0.0151) = 0.985, p-value = 0.0059), and survival censoring with shared parameters (effect estimate = exp(−0.0094) = 0.990, p-value = 0.0326) for joint modeling. These results indicate that there is an inverse relationship between HDL and CTGF continuous values. In this regard, it can be noted that as HDL increases by 1 mg/dl, CTGF levels decrease by about 1%–3% (ng/ml). These results were not consistent with the Poisson censoring with shared parameters for joint modeling which failed to show a significant association between HDL and CTGF. This is not in line with clinical findings^42,45 that showed a significant association between CTGF and hyperlipidemia; thus, do not lend support to this discrepancy in results denoted under the shared parameter models. The decreased power to capture this significant relationship could be attributed to the assumption of conditional independence between the random effects in the shared parameters structure for joint modeling.

3.3 Time effect

As for the longitudinal profile of CTGF, our results showed that CTGF values were decreasing significantly over time with an average decrease of about 11% in CTGF level every year. Percent change was obtained as exp(|Beta1_Time|)−1) × 100%). This result was demonstrated under the Poisson censoring with copula (effect estimate = exp(-0.1092) = 0.896, p-value < 0.0001), Poisson censoring with shared parameters (effect estimate = exp(−0.0817) = 0.921, p-value < 0.0001) for joint modeling, and survival model with copula for joint modeling (effect estimate = exp(−0.1237) = 0.883, p -value = 0.0006). This pattern in time for the CTGF levels was detected under all models except for the shared parameters with survival censoring.

3.4 Censoring effect

With respect to the censoring process, our results indicated the presence of informative right censoring that is dependent on both the baseline values of CTGF and the fluctuation in its values over time captured under all models. In this regard, an inverse association was detected between the baseline values of CTGF and the censoring process as evidenced by the parameter estimate Gamma1_Beta0i and its corresponding significance under all models. The corresponding parameter estimates for the effect of CTGF levels at baseline on the censoring process were as follows; Poisson with copula (effect estimate = exp(−1.9989) = 0.135, p-value < 0.0001), Poisson with shared parameters (effect estimate = exp(−0.2325) = 0.793, p-value < 0.0001) for joint modeling, survival censoring with copula (effect estimate = exp(−0.5019) = 0.605, p-value < 0.0001), and survival censoring shared parameters (effect estimate = exp(−1.4735) = 0.229, p-value < 0.0001) for joint modeling. This suggests that an increase in CTGF at baseline contributes to a decrease in the survival life of the diabetic patient before developing macroalbuminuria. Similarly an increase in the CTGF values over time contributes to a significant decrease in the survival time before experiencing macroalbuminuria. This was evidenced by the estimate of the effect of the CTGF slopes over time on the censoring process determined by Gamma2_Beta1i which was significant under the Poisson censoring with copula (effect estimate = exp(−0.4996) = 0.606, p-value < 0.0001), Poisson with shared parameters (effect estimate = exp(−0.2129) = 0.808, p-value < 0.0001) for joint modeling, survival censoring with copula (effect estimate = exp(−0.2752) = 0.759, p-value < 0.0001), and survival censoring with shared parameters (effect estimate = exp(−0.4484) = 0.638, p-value = 0.0496). Hence, these results indicated that there is an informative right censoring that leads to informative dropout due to macroalbuminuria which is contributed to the longitudinal pattern of CTGF and its baseline measures. The higher the increase in CTGF values over time and the higher the baseline measures of CTGF the shorter the time to develop macroalbuminuria. Accordingly, it can be inferred that CTGF baseline measures and longitudinal trajectories are consistently associated with macroalbuminuria indicating that the missing mechanism was dependent on both the intercept and slope of CTGF.

4. Simulation studies

The proposed approaches were investigated under extensive simulation studies that aimed at comparing the models’ performances and illustrating the feasibility of implementation and successful convergence of the framework of models under various conditions. The performance indicators for all four models were determined under different censoring levels and preponderance of zero values. These included low (10%), middle (30%), and high (50%) proportions of zero values introduced to the outcome and examined under medium to high levels of censoring ranging respectively between 15% and 20% and 48% and 53%. For feasibility purposes, this aim was thoroughly addressed by modeling the semicontinuous outcome as a function of time. Conducting a simulation study with all the covariates simulated for the zero and nonzero parts studied under the different settings is overwhelming, and its complexity is compounded by the choice of the different values for the parameters for all the covariates generated from the clinical application. Hence, for practicality and feasibility purposes we resorted to expressing the outcome as a function of time while achieving the aim of illustrating feasibility of implementation and models’ performance assessment. However, to determine the models’ performance when covariates were simulated and incorporated in the semicontinuous outcome, we extended our simulation study to include covariates in the binary and continuous parts and we executed this simulation study under the extreme condition of high-level censoring and proliferation of zero values. Details of the simulation studies are presented in the following section.

4.1 Simulation study design with time effect

The value of the parameter for the time effect β₁ we used in our simulation study was chosen to be the closest to the parameter value estimates we got under the four different models in our clinical application. On the other hand, the parameters for the censoring process were chosen in a way to give medium to high levels of censoring in the observations. This was reflected by having censoring levels ranging from 15 to more than 50% generated by the choice of $γ_{0} = 6$ , $γ_{1} = - 5.2,$ $γ_{2} = 6.2$ , and $γ_{2} = 15.2$ , respectively.

A random number generator which was a function of the censoring parameters was used to generate the number of observations for each sampled data and the respective observed survival time. The continuous part of the outcome $Y_{i j}$ was simulated from a lognormal distribution modeled as a function of time $t_{i j}$ and individual random intercepts $β_{0 i}$ and slopes $β_{1 i}$ as indicated in this equation $G^{- 1} (Y_{i j}) = β_{0 i} + β_{1 i} t_{i j} + e_{i j}$ and equivalently $G^{- 1} (E (Y_{i j})) = β_{0 i} + β_{1 i} t_{i j}$ . $G^{- 1}$ is the canonical link function pertaining to the natural parameter of the family of the outcome distribution which is in this case the logarithmic function for the lognormal distribution, and $t_{i j} = 0, 1, 3, 5, 8$ represents the time of the j^th repeated measure of the outcome for the i^th subject. $β_{0 i}$ and $β_{1 i}$ were simulated from a bivariate normal distribution with parameters $β_{0} = 0.2$ , $β_{1} = - 0.2$ , $σ_{β_{0}}^{2} = 0 .1445,$ $σ_{β_{1}}^{2} = 0 .0187,$ $σ_{β_{0} β_{1}} = - 0.0166$ , $σ_{ε}^{2} = 1.9.$ As indicated earlier, the censoring parameters were $γ_{0} = 6$ , $γ_{1} = - 5.2$ , and $γ_{2} = 6.2$ for the medium level of censoring, and $γ_{2} = 15.2$ for the high level of censoring.

Performance is determined by the accuracy and precision of the estimates assessed by bias, relative bias, mean squared error for the population parameters denoted as MSE(a)_β0 for the intercepts, and MSE(a)_β1 for the slopes, and mean squared errors for the individual subject-specific parameters denoted as MSE(b)_ β0i for the intercepts, and MSE(b)_β1i for the slopes. Monte Carlo simulated data were generated using 5000 replications, sample size of 200, and repeated measures ranging between 0 and 5 observations.

The performance indicators were compared for all four models (1) Poisson censoring with copula for joint modeling, (2) Poisson censoring with shared parameters for joint modeling, (3) survival time-to-event model for the censoring process and copula for joint modeling, and (4) survival time-to-event model for the censoring process and shared parameters for joint modeling. The objective is to assess the robustness of each of the proposed models to the increase (or decrease) in the proportion of zeros in the semicontinuous outcome and the respective sensitivity of each of the models to the severity of the censoring in the missing process. Results of the simulation studies were outlined in Tables 2 and 3, along with the specifications of the parameters and the simulation conditions. We also present below a detailed discussion of the results of the simulations studies with time effect.

Table 2.
Simulation with time effect assuming medium-level censoring of 15%–20%. Summary^a of performance of four estimators: (1) likelihood estimator assuming copula joint modeling between the longitudinal and the Poisson censoring processes; (2) likelihood estimator assuming shared parameters for joint modeling between the longitudinal and the Poisson censoring processes; (3) likelihood estimator assuming copula joint modeling between the longitudinal and the survival censoring processes; and (4) likelihood estimator assuming shared parameters for joint modeling between the longitudinal and the survival censoring processes.

Estimator

Simulations conditions Performance indicator (1) Copula Poisson censoring (2) Shared parameters Poisson censoring (3) Copula survival censoring (4) Shared parameters survival censoring

10% zeros Bias_ β0 $\times$ 10 0.023 0.046 −0.019 0.123

15%–20% censoring Bias_β1 $\times$ 10 0.026 −0.034 −0.001 0.003

Relative Bias_ β0 0.011 0.023 −0.009 0.061

Relative Bias_β1 −0.012 0.017 0.000 −0.001

MSE(a)_ β0 $\times$ 100 0.002 0.004 0.011 0.056

MSE(a)_β1 $\times$ 100 0.005 0.002 0.126 0.001

MSE(b)_ β0i $\times$ 10 0.118 0.225 0.296 0.749

MSE(b)_β1i $\times$ 10 0.031 0.061 0.188 0.046

30% zeros Bias_ β0 $\times$ 10 0.077 0.185 −0.001 −0.029

15%–20% censoring Bias_β1 $\times$ 10 −0.129 −0.079 −0.010 0.056

Relative Bias_ β0 0.038 0.092 −0.000 −0.014

Relative Bias_β1 0.061 0.039 0.005 −0.028

MSE(a)_ β0 $\times$ 100 0.015 0.082 0.091 3.477

MSE(a)_β1 $\times$ 100 0.032 0.020 0.166 0.125

MSE(b)_ β0i $\times$ 10 0.117 0.172 0.100 0.389

MSE(b)_β1i $\times$ 10 0.040 0.045 0.068 0.038

50% zeros Bias_ β0 $\times$ 10 0.018 0.002 0.033 1.351

15%–20% censoring Bias_β1 $\times$ 10 0.028 0.097 −0.087 −0.157

Relative Bias_ β0 0.009 0.001 0.016 0.675

Relative Bias_β1 −0.014 −0.048 0.043 0.078

MSE(a)_ β0 $\times$ 100 0.033 0.126 0.001 15.609

MSE(a)_β1 $\times$ 100 0.050 0.053 0.007 0.555

MSE(b)_ β0i $\times$ 10 0.127 0.145 0.217 0.239

MSE(b)_β1i $\times$ 10 0.047 0.039 0.049 0.035

		Estimator
10% zeros	Bias_ β0 $\times$ 10	0.023	0.046	−0.019	0.123
15%–20% censoring	Bias_β1 $\times$ 10	0.026	−0.034	−0.001	0.003
	Relative Bias_ β0	0.011	0.023	−0.009	0.061
Relative Bias_β1	−0.012	0.017	0.000	−0.001
	MSE(a)_ β0 $\times$ 100	0.002	0.004	0.011	0.056
MSE(a)_β1 $\times$ 100	0.005	0.002	0.126	0.001
	MSE(b)_ β0i $\times$ 10	0.118	0.225	0.296	0.749
	MSE(b)_β1i $\times$ 10	0.031	0.061	0.188	0.046
30% zeros	Bias_ β0 $\times$ 10	0.077	0.185	−0.001	−0.029
15%–20% censoring	Bias_β1 $\times$ 10	−0.129	−0.079	−0.010	0.056
	Relative Bias_ β0	0.038	0.092	−0.000	−0.014
Relative Bias_β1	0.061	0.039	0.005	−0.028
	MSE(a)_ β0 $\times$ 100	0.015	0.082	0.091	3.477
MSE(a)_β1 $\times$ 100	0.032	0.020	0.166	0.125
	MSE(b)_ β0i $\times$ 10	0.117	0.172	0.100	0.389
MSE(b)_β1i $\times$ 10	0.040	0.045	0.068	0.038
50% zeros	Bias_ β0 $\times$ 10	0.018	0.002	0.033	1.351
15%–20% censoring	Bias_β1 $\times$ 10	0.028	0.097	−0.087	−0.157
	Relative Bias_ β0	0.009	0.001	0.016	0.675
Relative Bias_β1	−0.014	−0.048	0.043	0.078
	MSE(a)_ β0 $\times$ 100	0.033	0.126	0.001	15.609
MSE(a)_β1 $\times$ 100	0.050	0.053	0.007	0.555
	MSE(b)_ β0i $\times$ 10	0.127	0.145	0.217	0.239
	MSE(b)_β1i $\times$ 10	0.047	0.039	0.049	0.035

Datasets are simulated using the following parameters: β0 = 0.2, β1 =−-0.2, $γ_{0} = 6$ , $γ_{1} = - 5.2$ , $γ_{2} = 6.2$ , $σ_{β_{0}}^{2} = 0.1445$ , $σ_{β_{1}}^{2} = 0.0187$ , $σ_{β_{0} β_{1}} = - 0.0166$ , and $σ_{ε}^{2} = 1.9$ , with 5000 simulations.

Table 3.

Simulation with time effect assuming high-level censoring of 48%–53%. Summary^a of performance of four estimators: (1) likelihood estimator assuming copula joint modeling between the longitudinal and the Poisson censoring processes; (2) likelihood estimator assuming shared parameters for joint modeling between the longitudinal and the Poisson censoring processes; (3) likelihood estimator assuming copula joint modeling between the longitudinal and the survival censoring processes; and (4) likelihood estimator assuming shared parameters for joint modeling between the longitudinal and the survival censoring processes.

		Estimator
Simulations conditions	Performance indicator	(1) Copula Poisson censoring	(2) Shared parameters Poisson censoring	(3) Copula survival censoring	(4) Shared parameters survival censoring
10% zeros	Bias_ β0 $\times$ 10	0.038	0.031	−0.044	−0.023
48%–53% censoring	Bias_β1 $\times$ 10	−0.160	0.188	0.007	−0.084
	Relative Bias_ β0	0.019	0.015	−0.022	−0.011
	Relative Bias_β1	0.080	−0.094	−0.003	0.042
	MSE(a)_ β0 $\times$ 100	0.143	0.007	0.005	2.500
	MSE(a)_β1 $\times$ 100	0.862	0.047	0.028	0.198
	MSE(b)_ β0i $\times$ 10	0.106	0.279	0.411	0.484
	MSE(b)_β1i $\times$ 10	0.026	0.051	0.319	0.066
30% zeros	Bias_ β0 $\times$ 10	−0.010	0.653	0.000	0.666
48%–53% censoring	Bias_β1 $\times$ 10	−0.005	−0.002	0.012	0.170
	Relative Bias_ β0	−0.005	0.326	0.000	0.333
	Relative Bias_β1	0.002	0.001	−0.006	−0.085
	MSE(a)_ β0 $\times$ 100	0.000	6.315	0.062	8.293
	MSE(a)_β1 $\times$ 100	0.000	0.424	0.062	0.090
	MSE(b)_ β0i $\times$ 10	0.100	0.156	0.719	0.355
	MSE(b)_β1i $\times$ 10	0.029	0.036	0.069	0.055
50% zeros	Bias_ β0 $\times$ 10	−0.018	0.074	0.000	0.502
48%–53% censoring	Bias_β1 $\times$ 10	−0.125	0.628	−0.004	−0.015
	Relative Bias_ β0	−0.009	0.037	0.000	0.251
	Relative Bias_β1	0.062	−0.314	0.002	0.007
	MSE(a)_ β0 $\times$ 100	0.093	14.011	0.000	6.771
	MSE(a)_β1 $\times$ 100	0.254	2.144	0.003	0.446
	MSE(b)_ β0i $\times$ 10	0.123	0.129	0.252	0.178
	MSE(b)_β1i $\times$ 10	0.031	0.042	0.043	0.018

4.1.1. Effect of the proportion of zeros on models’ performance

We start our discussion by assessing the performance of each of the models under the different proportions of zeros. As indicated in Table 2, under the levels of censoring that were in the range of 15%–20%, our results suggested that the shared parameters models were more sensitive to the variation of the proportion of zeros compared to the copula models for joint modeling. In this regard, if we consider the shared parameters model with Poisson censoring (model 2 in Table 2), our results showed an increase in the mean square errors for $β_{0}$ and $β_{1}$ by 32-fold and 27-fold when the proportion of zeros increased from 10% to 50%, and by respective increase of 21-fold and 10-fold when the proportion of zeros increased from 10% to 30%. The Poisson censoring with copula model (model 1 in Table 2) showed an increase in the mean square errors for $β_{0}$ and $β_{1}$ by 17-fold and 10–fold, respectively, when the proportion of zeros increased from 10% to 50%, and by a respective increase of 8-fold and 6-fold when the proportion of zeros increased from 10% to 30%. Hence, as the proportion of zeros increases, both the Poisson shared parameter and Poisson copula models exhibited an increase in bias and MSE. However, Poisson copula had lesser increase in its performance indicators compared to the shared Poisson model.

As for the survival censoring with shared parameters (model 4 in Table 2), it is evident that the bias and relative bias increased with the increase of the zero proportion from 10% to 50%, by 10-fold for the estimates of the intercept $β_{0}$ and 52-fold for the estimates of the slope $β_{1}$ . Similarly, and still in model 4 (Table 2), the mean square errors for the intercept $β_{0}$ increased by about 3 × 10² fold and that of the slope $β_{1}$ increased by about 5 × 10² when the proportion of zeros increased from 10% to 50%. In addition, respective increase of 62-fold and 125-fold was denoted for the mean square errors for $β_{0}$ and $β_{1}$ when the proportion of zeros increased from 10% to 30% (model 4 in Table 2). The survival censoring with copula model (model 3 in Table 2) did not exhibit an evident sensitivity in its performance to the increase in the proportion of zeros.

The effect of the proportion of zeros on the performance of each of the models was also verified in the context of high level of censoring that had a range of 48%–53%. As shown in the corresponding results presented in Table 3, the shared parameter model with Poisson censoring was sensitive to the variation in proportions of zeros. In this regard, the Poisson censoring with shared parameters (model 2 in Table 3) exhibited an evident increase in the mean square errors of $β_{0}$ from 0 to 6.315 × 10⁻², and to14.011 × 10⁻² when the proportion of zeros increased from 10% to 30% and from 10% to 50%, respectively. The same can be noted for the mean square errors for $β_{1}$ that increased by 46-fold and by 9-fold when the proportion of zeros increased from 10% to 30% and from 10% to 50% respectively (model 4 in Table 3).

The shared parameter model with survival censoring (model 4 in Table 3) was also sensitive to the variation in the zero proportion. This was evident in the 22-fold increase in bias and relative bias, and 3-fold increase in the mean square errors in the estimate of the intercept $β_{0}$ that corresponded to the increase in the zero proportion from 10% to 50%. The mean square errors for $β_{0}$ also increased by about 3 fold when the proportion of zeros increased from 10% to 30% (model 4 in Table 3). The copula models with Poisson and survival censoring (models 1 and 3 in Table 3) appeared to be more robust to the variation in the proportion of zeros whereby the performance indicators did not exhibit any significant change with the variation in the proportions of zeros.

Hence, it can be drawn that joint modeling of the longitudinal measures and censoring process through copula was less sensitive to the proportion of zeros compared to the shared parameters approach. Moreover, it can also be noted that all four models exhibited a good overall performance but the copula approach was in general more accurate and precise than the shared parameters and this was apparent in the various levels of the zero proportion including the severe case of high proportion of zero (50%) in both Tables 2 and 3.

4.1.2. Effect of the censoring level on models’ performance

The effect of the censoring level on the performance of the models and their sensitivity to the increase in the severity of the missing process were assessed in these simulation studies. This can be determined by comparing the performance of each model under the censoring levels of 15%–20% (Table 2) to that under the censoring levels of 48%-53% (Table 3). Our results showed that the shared parameter models exhibited some sensitivity to the level of censoring. In this regard, under the 10% and 50% proportion of zeros, and when the level of censoring increased from between 15% and 20% (model 2 in Table 2) to 48%and 53% (model 2 in Table 3), the Poisson censoring with shared parameters had an increase in bias for the estimate of $β_{1}$ by about 6-old. Similarly, in the same context but under 30% zeros, when the censoring level increased from 15% to 20% (model 2 in Table 2) to 48% to 53% (model 2 in Table 3), the bias in $β_{0}$ increased by about 4-fold and its mean square errors by 77-fold, and the mean square errors of $β_{1}$ increased by 21-fold as well. Along the same lines, the mean square errors for $β_{0}$ and $β_{1}$ , under the 50% zeros increased from 0.126 × 10⁻² to 14.011 × 10⁻², and from 0.053 × 10⁻² to 2.144 × 10⁻² (model 2 Table 2 vs. model 2 Table 3).

The results of the survival censoring with shared parameters also suggested some sensitivity of this model to the level of censoring. Under the proportion of 10% zeros, if we compare the results of model 4 in Table 2 to the results of model 4 in Table 3 we notice an increase in bias of the estimate of $β_{1}$ by 28-fold and in its mean square errors by 2 × 10² fold, and in the mean square error of $β_{0}$ by 45-fold. In the same context of comparing model 4 performance listed in Tables 2 and 3, but with the increase of the zero proportion to 30%, we can notice that $β_{0}$ had an increase bias of 23-fold and mean square error of 2-fold. The estimate of $β_{1}$ also had an increase in bias by 3-fold.

As for the copula joint modeling, the results of the Poisson model (model 1 compared between Table 2 and Table 3) showed an increase in bias for $β_{1}$ by 6-fold and 4-fold, and in its mean square errors by 1.7 × 10² fold and 5-fold under the 10% and 50% zeros, respectively. Still under the Poisson censoring with copula model the mean square errors for $β_{0}$ under the 10% zeros also exhibited an increase of 72-fold when censoring increased from the range of 15%–20% (model 1 in Table 2) to 48%–52% (model 1 in Table 3). However, the survival censoring with copula was more robust to the level of censoring with an increase of 7 fold in only the mean square errors of the individual subjects $β_{0 i}$ (model 3 in Table 2 vs. model 3 in Table 3).

4.1.3. Comparison of the Poisson and survival censoring models

Our simulation results did not show a clear distinction in the performance of one censoring model compared to the other. This is attributed to the fact that no censoring model outweighed the other in terms of consistently having lower bias and MSEs. Hence, a definitive conclusion cannot be drawn on the approach that is better for modeling the missing mechanism.

4.1.4 Summary of the simulation studies results

Our simulation results showed that all four models exhibited a good performance with minimum bias and mean errors, suggesting an overall accuracy and precision in its estimates. Some of the models were more sensitive to the proportion of zeros in the semicontinuous than others. In particular the shared parameters models were more sensitive in their performance to the zero proportion compared to joint modeling with copula. Similar conclusion can be drawn for the sensitivity of the models to the severity of the censoring, whereby the models with shared parameters for joint modeling were less robust compared to the models with copula. In specific, the survival censoring with copula appeared to be the most robust to the proportion of zeros. Overall, it appears that a slight gain in accuracy and precision was attained under the copula joint modeling compared to the shared parameters approaches. This could be attributed to the relaxation of the conditional independence under the copula structure and the incorporation of the dependence between the longitudinal and censoring processes through the joint distribution of the two corresponding marginal functions.

5. Simulation study with time and covariates effects

To assess the feasibility of implementation and performance of the different models in the presence of covariates we carried out additional simulation studies with covariates included in the two parts of the models executed under the severe condition of high proliferation of zeros (50%) and censoring levels (48–53%). We simulated a normally distributed predictor x₁ ∼ N(2.5,1.5) which was included in the continuous part of the semicontinuous outcome with true beta of 0.06, and binary predictor x₂ with probability of success of 50% and true beta of 0.07 which was included in the binary part. Detailed specifications and results of this simulation study are displayed in Appendix 3 Table A1. The simulation results with covariates effects were in line with the results of the simulation results with time effects whereby the copula approach showed better performance than the shared parameter. In this regard, the bias under the copula approach was reduced by at least two fold compared to the shared parameter for both Poisson and exponential hazard and for all the parameters including the ones for the covariates. The mean squared errors were also reduced significantly under the copula models compared to shared parameter approach. Comparing the copula models with Poisson and exponential hazard as specifications for the censoring process, we noted that the copula with exponential hazard model outperformed the copula with Poisson censoring in majority of the parameters indicating overall better accuracy and precision. In this regard, copula exponential hazard had at least 4 fold decrease in bias compared to the Poisson copula in the estimates pertaining to the covariates x₁ and x₂ and its associated mean square errors were negligible. However, overall, all 4 approaches presented in this framework of models had minimum attributed bias and mean squared errors despite that the simulation study was carried out under sever conditions of censoring and proliferation of zeros.

6. Discussion

Semicontinuous data are characterized by two distinct features: (1) preponderance of zeros that form a point mass at this value and (2) nonzero positive values that are right skewed. Analysis of this data are typically done using TP approaches and gets intricate when measured in a longitudinal context. Analysis becomes even more complex when compounded by informative right censoring and missingness due to subject's attrition and dropout that is dependent on the rate of change in the longitudinal trajectories of the outcome. We propose here a framework of MTP models whereby the binary part of the semicontinuous outcome is modeled as logistic regression and the continuous nonzero part is modeled using the lognormal as a special case of mTP generalized gamma. The censoring process is modeled using two approaches: (1) survival time-to-drop out and (2) Poisson censoring for the number of follow-up visits preceding the subject's dropout. Poisson censoring is flexible in the sense that it relaxes the proportional hazard assumption imbedded in most survival models, and is also feasible since the exact dropout time is typically unknown.⁴⁰ This framework of models is likelihood based and its different components that include the binary and continuous parts of the semicontinuous outcome and the censoring process are all linked together using two different approaches: (1) shared latent random effects and (2) Clayton copula model. Shared random effects approach assumes conditional independence between the responses and their corresponding missingness given the random effects. In addition, it assumes conditional independence between the vectors of responses measured over time on each subject given the random effects.³⁴ This assumption of conditional independence leads to bias in the estimates of the overall conditional means when a serial correlation is present and embedded in the realization of the stochastic processes of the two parts of the semicontinuous responses and is being ignored by employing the shared random effect approach.¹⁵ Nevertheless, the degree of bias under this same setting is attenuated when the mTP model is employed since the marginal mean is more robust to ignoring this serial correlation and thus less susceptible to bias compared to conditional mean.³⁵ To relax the assumption of conditional independence imposed by the shared parameter approach we used the Clayton copula for joint modeling of all the components of the likelihood function. In general, a copula is a multivariate joint distribution that couples multiple marginal distributions of different processes using a correlation structure. Linking the two parts of the semicontinuous responses and the censoring process using copula has not been well studied, probably due to computational and programing challenges in implementing these models compared to the shared parameters approach. This poses itself as one of the novelties of our proposed corresponding models. Clayton copula was previously used with bivariate frailty models⁵⁶ and is more suitable for cases where the marginal distributions are non-normal as in the semicontinuous responses and censoring process. As for the latent variables, the vector of random effects was shared between the two parts of the semicontinuous response and the censoring process to incorporate the cross-part correlation between the unobserved responses and dropout mechanism.³⁴ In addition, the longitudinal nature of the responses was accommodated by having this vector of random effect shared between the two parts of the semicontinuous outcome similar to the approaches in^14,17,27 to account for the within-part correlations. Failing to incorporate the within-part and cross-part correlations will, respectively, result in inaccurate inferences and biased estimates.³⁹ Hence, this vector of random effects will account for the cross-part and within-part correlations and will therefore result in consistent and efficient estimation. The structure of the random effects employed in our proposed framework of models includes both random intercept and slope. This is considered as a methodological advance since majority of available approaches limit the random effects to just the intercept. Having a complex random structure along with an estimation process that requires simultaneous fitting of both parts of the semicontinuous response leads to challenges in the computational process. This could be the reason why majority of approaches omit the slope from the random effect structure.²⁰ Simultaneous fitting of the two parts of the semicontinuous response is needed to account for the cross-parts correlation because, if not performed as such, bias will be introduced to the estimates.^18,57,58 Moreover, other approaches resorted to dropping the random effects (intercept and slope) completely from the model²² to avoid computational burden triggered by the complex numerical integration needed in the estimation process for the covariate effects on the marginal overall mean.^57,59 Omitting the random effects fully from the model induces the assumption of independence between the binary and continuous parts of the semicontinuous response. Hence, these approaches account only for the within-cluster correlation in the repeated measures in each part separately by adjusting the standard errors using the robust sandwich estimate of the covariance matrix. This assumption of independence between the two parts of the semicontinuous response leads to ignoring the cross-parts correlation and introduces bias in the estimation process of the parameters.^18,57,58 In addition, when the random structure should be entailed of both random effects (intercept and slope), but only intercept was included, then this misspecification has a greater impact on the induced bias in the estimation of the conditional covariate effects in the binary part of the TP model compared to the marginal effect in the binary part of the mTP model.¹ Alternatively, when the random structure requires only having the intercept, and in the model we include both the random intercept and slope, then this misspecification should not introduce a significant bias in the estimate since the model is considered over-specified in this case.⁶⁰ Thus, under-specification is the one that should raise concern in the estimation process. An added advantage of the mTP model is that it directly parametrizes the effect of covariates on the marginal mean while keeping the original untransformed scale of the semicontinuous response, and allowing incorporation of any closed-form distributions for the marginal mean. However, the TP models commonly use the natural logarithmic transformation of the positive part of the semicontinuous response resulting in the lognormal distribution. This transformation might not be sufficient to normalize the data especially in the presence of strong skewness and could impose unrealistic condition of the symmetry on the logarithmic scale. Hence, a more flexible parametric model is needed to relax this assumption and to also avoid retransformation to the original scale especially when inferences become uninterpretable due to the complexity of the transformation function.³⁵ In this regard, generalized gamma was proposed as an alternative to transformation of the response to a different scale to achieve symmetry. Different specifications of the scale and shape parameters of the generalized gamma lead to different distributions including the lognormal as a special case.^37,57,61,62 Hence, our framework of models presents a methodological advance since it assumes the generalized gamma with the lognormal as a special case for this flexible parametric distribution, and employs the mTP approach for modeling the two parts of the semicontinuous data. MTP model has embedded advantages due to its robustness to misspecifications in the model and in the structure of the random effects,^20,33 and to ignoring the serial correlation in the stochastic processes of the two parts of the semicontinuous response.¹⁵ The lognormal distribution was chosen as a special case of generalized gamma in our framework of models since in general it leads to smaller bias and optimal estimation of marginal treatment effect compared to the other distributions that are also special cases of the generalized gamma including gamma and Weibull. In addition, unlike the lognormal, the more generalized gamma distribution can exhibit a convergence issue, a challenge that is expected given its complexity.⁶³

Accuracy and precision of our proposed approaches were examined using rigorous simulation studies and sensitivity analysis. The effect of the degree of preponderance of zeros and the severity of the censoring process on the performance of each model was considered. Our simulation results showed that all four models exhibited accuracy in the estimation but were not fully robust to the predominance of zeros and severity of censoring. This is in line with what has been denoted in a previous cross-sectional study with balanced observations³⁸ whereby an increase in bias was captured with the proliferation of zeros. Our sensitivity analysis showed that joined modeling using copula is more robust to both the proportion of zero values and the severity of censoring compared to the shared parameters approach. This could probably be due to relaxing the assumption of conditional independence by coupling the different marginal distributions into a joint multivariate distribution through a correlation structure. The two modeling approaches for the censoring process, Poisson and survival time-to-drop out, did not exhibit a distinctive outweigh in performance of one approach over the other. Accordingly it cannot be confirmed that one censoring model is consistently more accurate in the estimates than the other. This framework of models was also illustrated using a longitudinal cohort of type 1 diabetic subjects with CTGF fibrotic gene as the semicontinuous response. Typically CTGF is analyzed by omitting the zero values and considering the continuous part only. This would definitely result in a possible loss of important information especially on the factors that could regulate expression and release of this gene and therefore modulation of its excreted levels in blood. Parallel to the simulation results, the illustration using the diabetes data example showed a better concordance between inferences drawn from the copula model and clinical findings than the conclusions deduced from the shared random effect approach.

The different components of this framework of models were expressed in terms of covariates that are not necessarily the same. This includes the zero and nonzero parts of the longitudinal semicontinuous outcome, the mean of the Poisson censoring process, and the hazard function. The parameters $β_{0}$ and $β_{1}$ are assumed to be scalar values since the continuous and binary parts of the outcome are modeled and are expressed as functions of covariates. However, if there is a rationale to incorporate covariate effects on $β_{0}$ and $β_{1}$ then the proposed framework of models can be extended to accommodate expressing these parameters in terms of covariates while avoiding non-identifiability issues. The proposed mTP model is formulated in a way that the location parameter $μ_{i j}$ and the probability of nonzero value $π_{i j}$ are both expressed in terms of the covariates. As for the scale parameter $σ$ , it is assumed to be constant for all the positive values. However, our models can be extended to allow this scale parameter $σ$ to be expressed in terms of covariates such that $σ = σ_{i} = g^{- 1} (s w_{i})$ where $g^{- 1}$ is the inverse link function and $w_{i}$ is a set of covariates that could be the same or different from the covariates included in the different parts of the model. The link function could be specified as an exponential function whereby $g^{- 1} (s w_{i}) = \exp (s w_{i})$ as in the work of Liu et al.³⁹ To determine whether this extension is feasible, we attempted to modify the framework of models to allow this heterogeneity. We were successful in achieving convergence and attaining estimates of the parameters for all the four models (results are not shown). This extension of the scale parameter will allow heterogeneity in the variance of the positive values of the semicontinuous outcome and could be the focus of future research work.

Another parameter that can also be expressed in terms of covariates is the scalar-valued dependence parameter ρ in the copula model. Typically, ρ is assumed to be a scalar value,⁶⁴ the reason why we followed the same specification. However, our models can be also extended to allow this dependence parameter to be expressed in terms of specific covariates while keeping the range of its values between −1 and 1. In this regard, and to assess if this extension is feasible, we carried out a preliminary attempt and we modified the four models so that the copula dependence parameter is expressed as a function of covariates. Successful convergence was achieved and parameter estimates were attained indicating that this framework of models can be modified to allow the copula dependence parameter to be expressed as a function of the covariates (results are not shown). This specification will allow a more flexible assumption on the interdependency between the components of the copula when expressed as a function of specific covariates. This can be studied thoroughly in a follow-up future research paper.

To conclude, our work presents a methodological advance to the analysis of semicontinuous longitudinal data with informative right censoring manifested by subject's attrition and dropout that is a function of the longitudinal trajectories. The dearth of available models that tackle this problem could be due to its inherent complexity rendering its implementation difficult and possibly infeasible. The novelty consideration here is in the joint modeling of the censoring process and the semincontinuous response using the shared parameters and copula models with a complex random effect structure. In particular, Clayton copula represents a new method in linking the different components of the likelihood function. Successful convergence of our models using the adaptive Gaussian quadrature approach and dual quasi-Newton optimization algorithm for numerical integration was consistently achieved. This confirms the feasibility of implementation of this new framework of methodological approaches and its utility in real-life applications.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

This work was supported by the National Institutes of Health Grants HL077192 (AAJ) and 5 P01 HL055782.

Research ethics and patient consent

The DCCT was approved by the Institutional Review Boards of all participating DCCT centers and all participants provided written informed consent.

Data availability

Clinical illustration used in this study was adopted from the DCCT/EDIC study investigators who should be contacted directly to request data acquisition.

ORCID iD

Miran A. Jaffa

Appendix 1

Sklar's theorem for bivariate copula: Let H be a joint distribution function with margins $F_{1}$ and $F_{2}$ . A two-dimensional copula $C_{ρ} : [0, 1]^{2} \to [0, 1]$ exists such that $H (x_{1}, x_{2}) = C_{ρ} (F_{1} (x_{1}), F_{2} (x_{2}))$ where $ρ$ being the dependence parameter between the marginal functions $F_{1}$ and $F_{2}$ . The joint density function of the random variables $X_{1}$ and $X_{2}$ can be expressed as follows:

h (x_{1}, x_{2}) = \frac{\partial^{2} H (x_{1}, x_{2})}{\partial x_{1} \partial x_{2}} = \frac{\partial^{2} C_{ρ} (F_{1} (x_{1}), F_{2} (x_{2}))}{\partial x_{1} \partial x_{2}} = c_{ρ} (u, v) f_{1} (x_{1}) f_{2} (x_{2})

where

u = F_{1} (x_{1}), v = F_{2} (x_{2}), c_{ρ} (u, v) = \frac{\partial^{2} C_{ρ} (u, v)}{\partial_{u} \partial_{v}}

f_{1} (x_{1})

and

f_{2} (x_{2})

are the corresponding marginal densities.

Appendix 2

Appendix 3

Table A1.

Simulation with time and covariates effects assuming high-level censoring of 48–53%. Summary^a of performance of four estimators: (1) likelihood estimator assuming copula joint modeling between the longitudinal and the Poisson censoring processes; (2) likelihood estimator assuming shared parameters for joint modeling between the longitudinal and the Poisson censoring processes; (3) likelihood estimator assuming copula joint modeling between the longitudinal and the survival censoring processes; and (4) likelihood estimator assuming shared parameters for joint modeling between the longitudinal and the survival censoring processes.

		Estimator
Simulations Conditions	Performance indicator	(1) Copula Poisson censoring	(2) Shared Parameters Poisson censoring	(3) Copula survival censoring	(4) Shared parameters survival censoring
50% zeros 48%–53% censoring	Bias_ β0 $\times$ 10	0.003	−0.492	0.002	0.933
	Bias_β1 $\times$ 10	−0.103	0.449	0.005	−0.043
	Bias_βx1 $\times$ 10	0.282	0.476	−0.052	0.368
	Bias_βx2 $\times$ 10	−0.274	−0.096	−0.038	0.050
	Relative Bias_ β0	0.001	−0.246	0.001	0.466
	Relative Bias_β1	−0.051	0.224	0.002	−0.021
	Relative Bias_βx1	0.483	0.793	−0.090	0.630
	Relative Bias_βx2	−0.392	−0.138	−0.055	0.714
	MSE(a)_ β0 $\times$ 100	0.048	6.291	0.000	8.107
	MSE(a)_β1 $\times$ 100	0.171	1.271	0.000	0.279
	MSE(b)_ β0i $\times$ 10	0.123	1.369	0.142	2.299
	MSE(b)_β1i $\times$ 10	0.030	0.439	0.042	0.213
	MSE_β˟1 $\times$ 100	0.131	6.010	0.000	1.585
	MSE_β˟2 $\times$ 100	0.102	2.827	0.000	5.802

Datasets are simulated using the following parameters: β0 = 0.2, β1 = −0.2, $γ_{0} = 6$ , $γ_{1} = - 5.2$ , $γ_{2} = 15.2$ , $σ_{β_{0}}^{2} = 0.1445$ , $σ_{β_{1}}^{2} = 0.0187$ , $σ_{β_{0} β_{1}} = - 0.0166$ , and $σ_{ε}^{2} = 1.9$ , with 5000 simulations. Predictor x₁∼N(2.5,1.5) was included in the continuous part with true beta of 0.06, and binary predictor x₂ with probability of success of 50% and true beta of 0.07 was included in the binary part.

References

Tom

Farewell

. A likelihood-based two-part marginal model for longitudinal semicontinuous data. Stat Methods Med Res 2015; 24: 194–205.

Hasan

Sneddon

. Zero-inflated Poisson regression for longitudinal data. Communications in Statistics, Simulation and Computation 2009; 38: 638–653.

Tong

Sun

. Semiparametric analysis of panel count data with correlated observation and follow-up times. Lifetime Data Anal 2009; 15: 177–196.

Greene

. Models for count data with endogenous participation. Empir Econ 2009; 36: 133–173.

Todem

Zhang

Ismail

, et al. Random effects regression models for count data with excess zeros in caries research. J Appl Stat 2010; 37: 1661–1679.

Maruotti

. A two-part mixed-effects pattern-mixture model to handle zero-inflation and incompleteness in a longitudinal setting. Biom J 2011; 53: 716–734.

Tobin

. Estimation of relationships for limited dependent variables. Econometrica 1958; 26: 24–36.

Duan

Manning

Morris

, et al. A comparison of alternative models for the demand for medical care. Journal of Business and Economic Statistics 1983; 1: 115–126.

Manning

Morris

Newhouse

, et al. A two-part model of the demand for medical care: preliminary results from the health insurance study. In: van der Gaag

Perlman

(eds) Health, economics, and health economics. North-Holland: Amsterdam, 1981, pp.103–123.

10.

Jorgensen

Exponential dispersion models. Journal of the Royal Statistical Society. Series B (Methodological) 1987; 49: 127–162.

11.

Jorgensen

. The theory of dispersion models. 1st ed. Boca Raton, Florida: Chapman and Hall/CRC, 1997.

12.

Hass

Levine

Sands

, et al. The modeling of medical expenditure data from a longitudinal survey using the generalized method of moments (GMM) approach. Stat Med 2016; 35: 2652–2664.

13.

Kurz

. Tweedie distributions for fitting semicontinuous health care utilization cost data. BMC Med Res Methodol 2017; 17: 71.

14.

Olsen

Schafer

. A two-part random effects model for semicontinuous longitudinal data. J Am Stat Assoc 2001; 96: 730–745.

15.

Albert

Shen

. Modeling longitudinal semicontinuous emesis volume data with serial correlation in an acupuncture clinical trial. Journal of the Royal Statistical Society: Series C 2005; 54: 707–720.

16.

Berk

Lachenbruch

. Repeated measures with zeros. SMMR 2002; 11: 303–316.

17.

Tooze

Grunwald

Jones

. Analysis of repeated measures data with clumping at zero. SMMR 2002; 11: 341–355.

18.

Tom

Farewell

. Bias in 2-part mixed models for longitudinal semicontinuous data. Biostatistics 2009; 10: 374–389.

19.

Hall

Zhang

. Marginal models for zero inflated clustered data. Statistical Modeling 2004; 4: 161–180.

20.

Smith

Neelon

Preisser

, et al. A marginalized two-part model for longitudinal semicontinuous data. SMMR 2017; 26: 1949–1968.

21.

Jaffa

Gebregziabher

Garrett

, et al. Analysis of longitudinal semicontinuous data using marginalized two-part model. J Transl Med 2018; 16: 01.

22.

Lin

Shih

. Analyzing excessive No changes in clinical trials with clustered data. Biometrics 2005; 60: 257–267.

23.

Albert

. On the interpretation of marginal inference with a mixture model for clustered semicontinuous data. Biometrics 2005; 61: 879–880.

24.

Elashoff

Robbins

, et al. A hierarchical zero-inflated log-normal model for skewed responses. Stat Methods Med Res 2011; 20: 175–189.

25.

Ghosh

Albert

. A Bayesian analysis for longitudinal semicontinuous data with an application to an acupuncture clinical trial. Comput Stat Data Anal 2009;. 53: 699–706.

26.

Robinson

Zeger

Forrest

. A hierarchical multivariate two-part model for profiling providers’ effects on health care charges. J Am Stat Assoc 2006; 101: 911–923.

27.

Zhang

Strawderman

Cowen

, et al. Bayesian Inference for a two-part hierarchical model: an application to profiling providers in managed health care. J Am Stat Assoc 2006; 101: 934–945.

28.

Biswas

Das

. A Bayesian approach of analysing semi-continuous longitudinal data with monotone missingness. Statistical Modeling 2020; 20: 148–170.

29.

Albert

Chib

. Bayesian Analysis of binary and polychotomous response data. J Am Stat Assoc 1993; 88: 669–679.

30.

Little

. Modeling the drop-out mechanism in repeated-measures studies. J Am Stat Assoc 1995; 90: 1112–1121.

31.

Molenberghs

Michiels

Kenward

, et al. Missing data mechanisms and pattern mixture models. Stat Neerl 1998; 52: 153–161.

32.

Liu

. Joint modeling longitudinal semi-continuous data and survival, with application to longitudinal medical cost data. Statist Med. 2009; 28: 972–986.

33.

Heagerty

Kurland

. Misspecified maximum likelihood estimates and generalized linear mixed models. Biometrika 2001; 88: 973–985.

34.

Mahabadi

. A Bayesian shared parameter model for incomplete semicontinuous longitudinal data: an application To toenail dermatophyte onychomycosis study. Journal of Statistical Theory and Applications 2014; 13: 317–332.

35.

Yiu

Tom

. Two-part models with stochastic processes for modeling longitudinal semicontinuous data: computationally efficient inference and modeling the overall marginal mean. Stat Methods Med Res. 2018; 27: 3679–3695.

36.

Bailey

. Analysing changes in the presence of informative right censoring caused by death and withdrawal. Stat Med 1988; 7: 337–346.

37.

Manning

Basu

Mullahy

. Generalized modeling approaches to risk adjustment of skewed outcomes data. J Health Econ 2005; 20: 465–488.

38.

Smith

Neelon

Maciejewski

, et al. Two parts are better than one: modeling marginal means of semicontinuous data. Health Serv Outcomes Res Method 2017; 17: 198–218.

39.

Liu

Strawderman

Johnson

and M et al. Analyzing repeated measures semi-continuous data, with application to an alcohol dependence study. Stat Methods Med Res 2016; 25: 133–152.

40.

. Accommodating informative dropout and death: a joint modeling approach for longitudinal and semi-competing risks data. J R Stat Soc Ser C Appl Stat 2018; 67: 145–163.

41.

The DCCT Research Group. The diabetes control and complications trial (DCCT): design and methodologic considerations for the feasibility phase. Diabetes 1986; 35: 530–545.

42.

Jaffa

Usinger

McHenry

, et al. Connective tissue growth factor and susceptibility to renal and vascular disease risk in type 1 diabetes. J Clin Endocrinol Metab 2008; 93: 1893–1900.

43.

Molitch

Steffes

Sun W

, et al. Epidemiology of diabetes interventions and complications study group. Development and progression of renal insufficiency with and without albuminuria in adults with type 1 diabetes in the diabetes control and complications trial and the epidemiology of diabetes interventions and complications study. Diabetes Care 2010; 33: 1536–1543.

44.

Gerritsen

Falke

van Vuuren

, et al. Plasma CTGF is independently related to an increased risk of cardiovascular events and mortality in patients with atherosclerotic disease: the SMART study. Growth Factors 2016; 4: 149–158.

45.

Hunt

Jaffa

Garrett

, et al. Plasma connective growth factor (CTGF/CCN2) levels predict myocardial infraction in the veterans affairs diabetes trial (VADT) cohort. Diabetes Care 2018; 41: 840–846.

46.

Wang

Carter

Jaffa

, et al. Genetic variant in the promoter of connective tissue growth factor gene confers susceptibility to nephropathy in type 1 diabetes. J Med Genet. 2010; 47: 391–397.

47.

Liu

. A likelihood reformulation method in non-normal random effects models. Stat Med 2008; 27: 3105–3124.

48.

Molenberghs

Kenward

Lesaffre

. The analysis of longitudinal ordinal data with nonrandom dropout. Biometrika 1997; 84: 33–34.

49.

Daniels

Hogan

. Missing data in longitudinal studies. New York: CRC Press, 2008.

50.

Sklar

. Fonctions de répartition à n dimensions et leurs marges. Publ. Inst Statist Univ. Paris 1959; 8: 229–231.

51.

Fletcher

. Practical methods of optimization. 2nd ed. New York: Willey, 1987, p.49.

52.

Pinheiro

Bates

. Approximations to the Log-likelihood function in the nonlinear mixed effects model. J Comput Graph Stat 1995; 4: 12–35.

53.

The Diabetes Control and Complications Trial Research Group. The effect of intensive treatment of diabetes on the development and progression of long term complications in insulin-dependent diabetes mellitus. N Engl J Med 1993; 329: 977–986.

54.

Roestenberg

van Nieuwenhoven

Wieten

, et al. Connective tissue growth factor Is increased in plasma of type 1 diabetic patients with nephropathy. Diabetes Care 2004; 27: 1164–1170.

55.

Zhou

Zeng

, et al. Expression variations of connective tissue growth factor in pulmonary arteries from smokers with and without chronic obstructive pulmonary disease. Sci Rep 2015; 5: 8564.

56.

Oakes

. Bivariate survival models induced by frailties. J Am Stat Assoc 1989; 84: 487–493.

57.

Liu

Cowen

Strawderman

, et al. A flexible two-part random effects model for correlated medical costs. J Health Econ 2010; 29: 110–123.

58.

Albert

. Letter to the editor. Biometrics 2005; 47: 879–881.

59.

Tom

Farewell

. A correlated formulation for marginal inference derived from two-part mixed models for longitudinal semi-continuous data. SMMR 2016; 25: 2014–2020.

60.

Preisser

Das

Long

, et al. Marginalized zero-inflated negative binomial regression with application dental carries. Stat Med 2016; 35: 1722–1735.

61.

Azzalini

. A class of distributions which includes the normal ones. Scand. J. Stat 1985; 12: 171–178.

62.

Chai

Bailey

. Use of log-skew-normal distribution in analysis of continuous data with a discrete component at zero. Stat. Med. 2008; 27: 3643–3655.

63.

Voronca

Gebregziabher

Durkalski

, et al. Marginalized Two part models for generalized gamma family of distributions. Cornell University New York: Cornell University Library. 2015; arXiv: 1511.05629[stat.ME].

64.

Zimmer

D M

. Analyzing comovements in housing prices using vine copulas. Econ Inq 2015; 53: 1156–1169