Simulation of Survival Data With the Package rsurv (2024)

2.1 YP model and its particularcases

The YP model corresponds to a suitable choice to model survival datawith crossing survival curves. This model may also be appropriated insituations where PH and PO assumptions are not satisfied.

Following Demarqui andMayrink (2021) and references therein, the survivalfunction of the Yang and Prentice model is given by:

where$\displaystyle R_{0}(t|\boldsymbol{\theta})=\frac{1-S_{0}(t|\boldsymbol{\theta})}{S_{0}(t|\boldsymbol{\theta})}$is the baseline odds function, $\boldsymbol{\beta}$ and$\boldsymbol{\phi}$ are $p\times 1$ vectors of regressioncoefficients.

The hazard function of the YP model can be expressed as:

Let $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ two vectors of covariates.Becase

It can be show that the survival curves cross each other when$\beta_{j}\times\phi_{j}<0$, for some $j=1,\cdots,p$.Moreover, the YP model includes the PH and PO models as particular caseswhen $\boldsymbol{\beta}=\boldsymbol{\phi}$ and$\boldsymbol{\phi}=\mathbf{0}$, respectively. For this reason, theinverse of the survival function given in (2) allows thesimulation from the PH, PO and YP models.

2.2 EH model and its particularcases

The EH model, originally proposed by Etezadi-Amoli andCiampi (1987)…

According to Tseng andShu (2011), the regression coefficients$\boldsymbol{\beta}$ can be regarded as acceleration factors sincethey identify whether the acceleration or deceleration of the observedfailure time related to the baseline failure time. In contrast, theregression coefficients $\boldsymbol{\phi}$ characterize the relativehazard after transforming the time scale of the baseline hazard into anobserved time scale.

In this paper we consider a different parametrization for the EH, inline with the parametrization adopted in the R packagesurvstan (Demarqui, 2024b). Under that parametrization, the hazardand survival functions of the EH model can be expressed as:

and

It is easy to see from Equation (3) that the EH modelsincludes the AFT, PH and AH models as particular cases when$\boldsymbol{\phi}=-\boldsymbol{\beta}$,$\boldsymbol{\beta}=\mathbf{0}$, and$\boldsymbol{\phi}=\mathbf{0}$, respectively. Therefore, the inverseof the survival function given in (3) allows the simulationfrom the AFT, AH, PH and EH models.

2.3 Survival data generation procedure

The discussion presented in the previous sections suggests that thesurvival function $S(t|\boldsymbol{\Theta},\mathbf{x})$ should beexpressed as some function of $S_{0}(t|\boldsymbol{\theta})$ and$\mathbf{x}$, satisfying$S(t|\boldsymbol{\Theta},\mathbf{x})=S_{0}(t|\boldsymbol{\theta})$when $\mathbf{x}=\mathbf{0}$. The following theorem formalize thisidea, and establishes how random samples from the YP and EH regressionmodels (and its particular cases AFT, AH, PH and PO) can be draw.

Theorem 1 Let $\boldsymbol{\eta}=(\eta_{1},\eta_{2})$,where $\eta_{1}=\mathbf{x}\boldsymbol{\beta}$ and$\eta_{2}=\mathbf{x}\boldsymbol{\phi}$, and consider the survivalfunctions given in Equations (2) and (3). Then:

The proof of the theorm is given in the appendix.

In practice, the validity of Equation (5) for simulatingsurvival data from the regression models presented in the previoussection is ensured by two mild general assumptions:

It is important to notice that assumption i) does not require theexistence of a closed-form expression for$S_{0}(\cdot|\boldsymbol{\theta})$, as long as its inverse can becomputed numerically. Assumption ii), on the other hand, is necessary toavoid identifiability issues since no further assumptions/constraintsare made about the baseline survival function$S_{0}(\cdot|\boldsymbol{\theta})$.

As one can see, the result given in Equation (5) is quitegeneral, and can be applied to simulate data from a wide range ofsurvival models. The only restriction for its use regards the fact thatneither covariates nor their effects are allowed to vary with time.

The internal function qsurv, reproduced below, plays a distinctrole in the algorithms implemented in the package rsurv.Specifically, this function allows the use of any baseline survivaldistribution for which there exists the implementation of its quantilefunction.

The following example shows the usefulness of the qsurv functionwhen working with customized functions or functions imported from otherpackages. There, the function qmydist mimics the functionqexp available in the package stats, whereas the functionqgengamma.orig from the package flexsurv returns thequantile function of the generalized gamma distribution, which includesthe exponential distribution as a particular case when all parametersare equal to 1.

## x1 x2 x3 x4## [1,] 0.09339179 0.09339179 0.09339179 0.09339179## [2,] 0.95202641 0.95202641 0.95202641 0.95202641## [3,] 0.17411789 0.17411789 0.17411789 0.17411789## [4,] 0.29132206 0.29132206 0.29132206 0.29132206## [5,] 0.34595405 0.34595405 0.34595405 0.34595405

The regression models implemented in the package rsurv share thesame parametrizations of those implement in the R packagesurvstan (Demarqui, 2024b). The main functions available in thepackage rsurv are: raftreg(), rahtreg(),rapheg(), raporeg(), rehreg() and rypreg().These functions permit the simulation of survival times from the AFT,AH, PH, PO, EH and YP models, respectively.

All r*reg() functions share the same basic structure.Specifically, the user must provide a numeric vector u ofquantiles used in (5), a formula containing informationregarding the structure of the linear predictors, a vector betaof regression coefficients (and additionally a vector phi for theEH and YP models), which should be compatible with the model matrixinduced by the formula argument, a dist/baseline argumentused to specify the chosen baseline distribution, along with anyadditional arguments passed to the internal function qsurvthrough the … argument.

A distinct feature of the functions r*reg() is that they return avector of generated survival times rather than a full data setcontaining additional information regarding covariates and a failureindicator variable. This brings more flexibility for the user tosimulate survival data with different characteristics and levels ofcomplexity, as we shall demonstrate throughout the examples presented inthe paper.

As a first example, suppose we want to simulate a sample of type Iright-censored survival data, assuming that the failure times aregenerated from an accelerated failure time model with loglogisticbaseline distribution. In addition, assume that we wish to consider twoexploratory variables, say age and sex, and we want to include aninteraction effect between them. Such a task can be easily accomplishedby using the function raftreg along with the functionqllogis available in the package flexsurv, as illustratedbelow:

## Rows: 1,000## Columns: 5## Rowwise:## $ age <dbl> 1.34592454, 0.99527131, 0.54622688, -1.91272392, 1.92128431, 1.~## $ sex <chr> "m", "f", "f", "m", "m", "m", "m", "f", "f", "m", "f", "m", "m"~## $ t <dbl> 15.2363453, 1.5259533, 2.1783746, 2.4354995, 58.7932958, 16.714~## $ time <dbl> 10.0000000, 1.5259533, 2.1783746, 2.4354995, 10.0000000, 10.000~## $ status <dbl> 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, ~

## age sexm age:sexm alpha gamma## 0.9494630 2.0094422 -0.4812641 1.4946497 1.0226847

In our next example we demonstrate how to generate survival data withcrossing survival curves using the YP model. In this example, we assumethat the survival times are randomly right-censored, with the censoringtimes following an exponential distribution.

## Rows: 1,000## Columns: 3## Rowwise:## $ trt <chr> "chemo", "chemo", "chemo+rad", "chemo+rad", "chemo+rad", "chemo~## $ time <dbl> 0.65743632, 1.14638933, 0.10715893, 0.09876511, 1.44704010, 0.1~## $ status <dbl> 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, ~

## short-trtchemo+rad long-trtchemo+rad alpha gamma## 1.920795 -1.446502 1.404515 1.043123

The two examples presented above clearly demonstrate the flexibilityprovided by the functions r*reg() in the simulation of survivaldata with different types of censoring. In the next section we look atthe simulation of survival data with more sophisticated structures.

3.1 Cure rate models

Rodrigues etal. (2009) presented an unified formulation for cure ratemodels, through the use of probability-generating functions, thatincludes several cure rate models available in the literature. In theirframework, cure rate models are formulated based on a two-stage processthat possesses an appealing biological motivation in terms ofincidence-latency of disease.

Suppose we are interested in modeling the time to relapse to cancer andlet $N$ denote the number of carcinogenic cells (often calledclonogens) left active after a patient has gone through initialtreatment. Then, for an individual with $N$ clonogens left active, let$\xi_{j}$ be the non-observable random time (i.e., thepromotion time) for the $j$-th clonogen to produce a detectable cancermass, $j=0,1,...,N$.

The time to relapse to cancer, which is in fact the observable quantity,is defined as:

Then, under the assumptions that $N$ is independent of$\xi_{j},~{}\forall~{}j\geq 0$, and$\xi_{j}\overset{\text{i.i.d.}}{\sim}S_{\xi}(\cdot|\boldsymbol{\Theta},\mathbf{x})$,$j\geq 0$, it follows that the population survival function can beexpressed as:

The next corollary of Theorems 1 and 2 provides a general expression tosimulate survival data with a cure fraction.

Corollary 1 Given a cure fraction probability $\pi=P(N=0)$and a survival function $S_{\xi}(t|\boldsymbol{\Theta},\mathbf{x})$satisfying (4), then time to relapse to cancer $T$ can begenerated as follows:

where $U\sim U(0,1)$.

The proof of the corollary in show in the Appendix.

It is important to note from Equation (9) that, given$V=A_{p}^{-1}(U)$, all the time generation procedures describedpreviously can be straightforwardly applied in the generation ofsurvival times in the presence of a cure fraction, provided that theinverse of the probability generating function, $A_{p}^{-1}(u)$, isknown.

Common distributions assumed for the incidence sub-model (i.e.,the distribution of the latent r.v. $N$) are the Bernoulli, Poisson,negative binomial and the Bell distributions. Covariates are usuallyincluded in the incidence sub-model. If $\mathbf{z}$ is a$1\times q$ vector of covariates, and $\boldsymbol{\kappa}$ is thecorresponding $q\times 1$ vector of regression coefficients, than$\pi=P(N=0)$ can be related to the linear predictor$\mathbf{z}\boldsymbol{\kappa}$ through appropriate choices of linkfunctions.

The inverse of the probability generating functions presented in Table1 are implemented in the function inv_pgf().Similarly to the r*reg() functions, the functioninv_pgf() requires a formula argument specifying thestructure of the linear predictor, along with a kappa argumentpassing the vector of regression coefficients, that should be compatiblewith the model matrix induced by the formula argument, and anincidence argument determining the incidence sub-model, and itscorresponding link function. Specifically, the functionrsurv::bernoulli() allows the specification of logit, probit,cloglog, and cauchit link functions, whereas the log, identity andsquare root links can be used with the functionsrsurv::poisson(), rsurv::bell() andrsurv::negbin(). When the negative binomial distribution isassumed for the incidence sub-model, the function inv_pgf()further requires an additional argument zeta, representing theextra negative-binomial parameter.

The next example illustrates how survival data with cure fraction can beeasily generated using the package rsurv. In that example weconsider a regression structure only for the incidence sub-model. Morespecifcally, we simulate survival data from a mixture cure rate modelwith a probit link function. We further assume that the survival timesare randomly right-censored.

## Rows: 1,000## Columns: 4## Rowwise:## $ trt <chr> "A", "A", "B", "B", "B", "A", "B", "A", "B", "A", "A", "A", "A"~## $ age <dbl> 0.2193297, 0.9952571, -0.4572095, 0.4864704, 1.9512081, 0.49597~## $ time <dbl> 0.17789062, 0.17142535, 0.50301782, 0.64774117, 0.40046281, 0.4~## $ status <dbl> 0, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, ~

3.2 Fraily models

Frailty models are often employed in survival analysis to accommodateunknown and/or unobserved risk factors through the inclusion of randomeffects, known as frailties (a term coined by Vaupeletal. (1979) torepresent individual’s heterogeneity). These type of models have becomeincreasingly popular in multivariate survival analysis since they allowone to take into account the association among individual survival timeswithin subgroups (or clusters) of subjects.

The package rsurv permits the inclusion of frailty terms throughthe linear predictors, i.e.,$\eta_{1}=\mathbf{x}\boldsymbol{\beta}+w$ and$\eta_{2}=\mathbf{x}\boldsymbol{\phi}+w$, where $w$ is thefrailty term. This convention adopted by the package rsruv allowsthe inclusion of frailties into the linear predictors using thestats::offset() function.

The function rfrailty() implemented in the package rsurvfacilitates the process of simulating survival data with frailties. Thisfunction permits random generation of frailties from the most commonlyfrailty distributions used in practice, namely the gamma distribution,the gaussian distribution, and the positive stable (ps) distribution,which are briefly described in Table 2. Since thegamma and positive stable distributions have support on $\Re^{+}$, bydefault, the function rfrailty() returns samples from positiverandom frailties in the log-scale, that is, $w=\log(z)$.

Our next example shows how the functions r*reg() andrfrailty() can be combined to generate survival data from sharedfrailty models. For this type of survival frailty models, elementsbelonging to the same cluster share the same random effect (frailty). Togenerate frailties using the function rfrailty() we need to passa vector of length $n$ containing information regarding the clusterconfiguration to the cluster argument, the assumed distributionto the frailty argument, and finally, the parameter of thefrailty distribution (sigma or alpha, depending on thechoice of passed to frailty argument).

## Rows: 1,000## Columns: 6## Rowwise:## $ id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3,~## $ age <dbl> 1.34592454, 0.99527131, 0.54622688, -1.91272392, 1.92128431, 1~## $ sex <chr> "m", "f", "f", "m", "m", "m", "m", "f", "f", "m", "f", "m", "m~## $ frailty <dbl> -0.15917714, -0.15917714, -0.15917714, -0.15917714, -0.1591771~## $ time <dbl> 0.066662943, 0.370877454, 3.527734703, 0.109155541, 0.11269257~## $ status <dbl> 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,~

## Call:## emfrail(formula = Surv(time, status) ~ age * sex + cluster(id),## data = simdata, distribution = emfrail_dist(dist = "gamma"))#### Regression coefficients:## coef exp(coef) se(coef) adj. se z p## age 0.9711 2.6409 0.0675 0.0682 14.2469 0## sexm 2.1159 8.2973 0.0949 0.0976 21.6777 0## age:sexm -0.4746 0.6221 0.0811 0.0813 -5.8405 0## Estimated distribution: gamma / left truncation: FALSE#### Fit summary:## Commenges-Andersen test for heterogeneity: p-val 1.25e-19## no-frailty Log-likelihood: -4679.899## Log-likelihood: -4636.705## LRT: 1/2 * pchisq(86.4), p-val 7.39e-21#### Frailty summary:## estimate lower 95% upper 95%## Var[Z] 0.233 0.152 0.350## Kendall’s tau 0.105 0.071 0.149## Median concordance 0.102 0.068 0.145## E[logZ] -0.121 -0.185 -0.078## Var[logZ] 0.263 0.164 0.418## theta 4.284 2.859 6.570## Confidence intervals based on the likelihood function

3.3 Interval-censored survivaldata

Consider a nonnegative random variable $T$ representing the time untilthe occorrence of an event of interest. According to Sun (2006), anobservation on $T$ is interval-censored when its exact value isunknow, and only an interval of the form $(L,R]$ can be observed, sothat

It is easy to see from (10) that right and left-censoredobservations on $T$ occur when $R=\infty$ and $L=0$,respectively. Naturally, when $L=R$, $T$ can be observed exactly.

Two important types of interval-censored survival data that often arisein practice are the so-called type I (also known as current statussurvival data in the literature) and type II interval-censored survivaldata. Type I interval-censored survival data occurs when all observedintervals result in either left-censored ($L=0$) or right-censored($R=\infty$) observations. Type II interval-censored survival data,on the other hand, emerges when we have at least one finite intervalsatisfying $L>0$.

The generation of type I and type II interval-censored survival data isfacilitated by the function rinterval available in the packagersurv. For the random generation of type II interval-censoredsurvival data, the algorithm described in Kiani andArasan (2012), whichemulates a situation in which patients are expected to attendpre-scheduled visits to the physician with attendance probability $p$,is implemented. The function rinterval requires an argumenttime passing a sample of failure times, an argument tauthat can be either a vector of censoring times (for type Iinterval-censored survival data) or a time-grid of scheduled visits (fortype II interval-censored survival data), and an argument typeindicating the type of interval-censored survival data. Finally, theprob argument, only evaluated when type = II, can beregarded as the attendance probability of a scheduled visit.

The code presented below examplifies how to simulate type Iinterval-censored survival data using the function functionrinterval(). There, it is assumed that the exact failure timesare generated from a PH model with an exponential baseline distribution,whereas the censoring times are generated from a Weibull distribution.

## Rows: 300## Columns: 6## $ age <dbl> 1.34592454, 0.99527131, 0.54622688, -1.91272392, 1.92128431, 1.3~## $ sex <chr> "m", "f", "m", "m", "m", "f", "m", "f", "f", "m", "m", "f", "f",~## $ time <dbl> 0.05387373, 0.05326879, 0.16546309, 2.17932784, 0.04415598, 0.01~## $ tau <dbl> 0.8276881, 1.4807894, 3.2631929, 3.3893988, 2.7803610, 0.4764430~## $ left <dbl> 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000~## $ right <dbl> 0.8276881, 1.4807894, 3.2631929, 3.3893988, 2.7803610, 0.4764430~

#### Model: Cox PH## Dependency structure assumed: Independence## Baseline: exponential## Call: ic_par(formula = cbind(left, right) ~ age + sex, data = simdata1,## model = "ph", dist = "exponential")#### Estimate Exp(Est) Std.Error z-value p## log_scale -0.3877 0.6786 0.1038 -3.736 1.871e-04## age 0.9186 2.5060 0.1241 7.402 1.341e-13## sexm 0.2245 1.2520 0.1965 1.142 2.534e-01#### final llk = -94.91423## Iterations = 9

The next example illustrates how type II interval-censored survival datacan be generated using the function rinterval(). In that example,we assume that patients are expected to return to the physician atcertain pre-specified times, with attendancy probability equals to 0.7.

## Rows: 300## Columns: 5## $ age <dbl> 1.34592454, 0.99527131, 0.54622688, -1.91272392, 1.92128431, 1.3~## $ sex <chr> "m", "f", "m", "m", "m", "f", "m", "f", "f", "m", "m", "f", "f",~## $ time <dbl> 2.583140967, 2.152683999, 0.350258378, 0.045283358, 20.802515190~## $ left <dbl> 2, 1, 0, 0, 5, 2, 0, 1, 0, 1, 0, 0, 1, 2, 2, 0, 0, 1, 3, 5, 2, 3~## $ right <dbl> 3, Inf, 1, 1, Inf, 3, 1, 2, 1, 3, 1, 1, 2, 4, 5, 2, 1, 2, Inf, I~

#### Model: Accelerated Failure Time## Dependency structure assumed: Independence## Baseline: exponential## Call: ic_par(formula = cbind(left, right) ~ age + sex, data = simdata2,## model = "aft", dist = "exponential")#### Estimate Exp(Est) Std.Error z-value p## log_scale 0.3638 1.439 0.0716 5.081 3.747e-07## age 1.0450 2.842 0.0908 11.500 0.000e+00## sexm 0.4921 1.636 0.1480 3.326 8.825e-04#### final llk = -266.4422## Iterations = 7

3.4 Multivariate survival data, dependent censoring andcompetingrisks

Copulas provide a way to construct multivariate distributions based uponindependent marginal distributions (Hofertetal., 2018). A multivariatesurvival function can be constructed by the specification of independentmarginal survival models and an appropriated copula function.

Let $C_{\theta}(u_{1},\cdots,u_{d})$ a d-variated copula. Then, itfollows that the joint survival function of $(T_{1},\cdots,T_{d})$can be expressed as:

To generate a random sample of $\mathbf{T}=(T_{1},\cdots,T_{d})$it is sufficient to generate a random sample from$(U_{1},\cdots,U_{d})\sim C_{\theta}(u_{1},\cdots,u_{d})$ andthen take$T_{i}=S_{i}^{-1}(U_{i}|\boldsymbol{\Theta}_{i},\mathbf{x}_{i})$.

## Rows: 1,000## Columns: 4## $ age <dbl> -0.17469076, -0.15529194, -1.03781319, -1.56364665, 1.03961867, -0~## $ sex <chr> "m", "m", "f", "f", "f", "m", "m", "m", "m", "f", "f", "f", "f", "~## $ t1 <dbl> 0.016749142, 0.167458832, 0.245770610, 0.695707186, 0.061162858, 0~## $ t2 <dbl> 0.084588793, 0.583402263, 0.175592119, 0.288876906, 0.298716691, 0~

## t1 t2## t1 1.0000000 0.5306627## t2 0.5306627 1.0000000

## Rows: 1,000## Columns: 7## Rowwise:## $ age <dbl> -0.17469076, -0.15529194, -1.03781319, -1.56364665, 1.03961867~## $ sex <chr> "m", "m", "f", "f", "f", "m", "m", "m", "m", "f", "f", "f", "f~## $ t1 <dbl> 0.016749142, 0.167458832, 0.245770610, 0.695707186, 0.06116285~## $ t2 <dbl> 0.084588793, 0.583402263, 0.175592119, 0.288876906, 0.29871669~## $ c <dbl> 2.26609325, 4.57257662, 3.63517003, 2.02494996, 2.66336179, 2.~## $ status1 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,~## $ status2 <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0,~

## t1 t2## t1 1.0000000 0.5121882## t2 0.5121882 1.0000000

In the next example we show how to simulate survival data with dependentcensoring. In that case, the observable survival time is given by$Y=\min\{T,C\,A\}$, where $T$, $C$, and $A$ represents thetime to failure, the dependent censoring time, and the censoring timedue to random censoring.

## Rows: 1,000## Columns: 6## Rowwise:## $ age <dbl> -0.17469076, -0.15529194, -1.03781319, -1.56364665, 1.03961867~## $ sex <chr> "m", "m", "f", "f", "f", "m", "m", "m", "m", "f", "f", "f", "f~## $ a <dbl> 3.1886887, 3.6066255, 3.0834830, 0.8921447, 0.8049795, 0.96584~## $ y <dbl> 0.016749142, 0.167458832, 0.175592119, 0.288876906, 0.06116285~## $ status1 <dbl> 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1,~## $ status2 <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0,~

Notice that the example presented above can be viewed as a particularcase of competing risks survival data. In the next example, we considera scenario where elements are subjected to three competing causes offailure.

## Rows: 1,000## Columns: 7## Rowwise:## $ age <dbl> -0.53451004, -1.48095715, -0.76892438, 0.43319025, 0.76357712,~## $ sex <chr> "m", "m", "m", "m", "m", "m", "m", "f", "f", "f", "f", "m", "f~## $ a <dbl> 2.17617891, 1.49961178, 4.65348449, 1.04816011, 1.12082428, 0.~## $ y <dbl> 0.18034466, 0.51001897, 0.07337358, 0.12327447, 0.13908667, 0.~## $ status1 <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,~## $ status2 <dbl> 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0,~## $ status3 <dbl> 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1,~

Proof of Theorem 1

We start the proof of Theorem 1 by the YP model given in Equation(2). For this step of the proof, define$\kappa_{1}=e^{\eta_{1}}$ and $\kappa_{2}=e^{\eta_{2}}$.

Now we present the proof of Theorem 1 for the EH model given in Equation(3).

Proof of Corollary 1

It follows immediately from Theorem 2 that:

The rest of the proof follows from (9) and the directapplication of Theorem 1.