Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Biometrics. Author manuscript; available in PMC 2014 March 3.
Published in final edited form as:
Published online 2013 February 19. doi:  10.1111/biom.12016
PMCID: PMC3940058

Mark-specific Hazard Ratio Model with Multivariate Continuous Marks: An Application to Vaccine Efficacy


In randomized placebo-controlled preventive HIV vaccine efficacy trials, an objective is to evaluate the relationship between vaccine efficacy to prevent infection and genetic distances of the exposing HIV strains to the multiple HIV sequences included in the vaccine construct, where the set of genetic distances is considered as the continuous multivariate ‘mark’ observed in infected subjects only. This research develops a multivariate mark-specific hazard ratio model in the competing risks failure time analysis framework for the assessment of mark-specific vaccine efficacy. It allows improved efficiency of estimation by employing the semiparametric method of maximum profile likelihood estimation in the vaccine-to-placebo mark density ratio model. The model also enables the use of a more efficient estimation method for the overall log hazard ratio in the Cox model. Additionally, we propose testing procedures to evaluate two relevant hypotheses concerning mark-specific vaccine efficacy. The asymptotic properties and finite-sample performance of the inferential procedures are investigated. Finally, we apply the proposed methods to data collected in the Thai RV144 HIV vaccine efficacy trial.

Keywords: Competing risks, Competing risks, Cox regression model, Density ratio model, Failure time data, Genetic data, Mark variable, Mark-specific vaccine efficacy

1. Introduction

In studies with competing risks data, commonly a time-to-event is evaluated, and data from those who experience the event – termed “mark” data – are collected. The analysis of such data often involves the comparison of hazard rates of mark-specific failure between two treatment groups. We propose estimation and testing procedures for the continuous mark specific hazard ratio model that improves efficiency compared to alternative approaches and accommodates a multivariate mark.

The extensive genetic diversity of HIV poses a central challenge to development of an efficacious preventive HIV vaccine. Efficacy trials of such vaccines randomize HIV-uninfected volunteers to receive vaccine or placebo, and the primary objective assesses vaccine efficacy (VE), typically defined as one minus the vaccine-to-placebo hazard ratio of HIV acquisition (Halloran, Struchiner, and Longini, 1997). A secondary objective assesses VE to prevent acquisition with particular genotypes of HIV, where we expect that VE will be greatest against HIVs genetically similar to the HIV strains included in the vaccine, and will be lower against HIVs more distant from the vaccine strains (Gilbert, McKeague, and Sun, 2008) (henceforth GMS). It is important to maximize efficiency for assessing genotype-specific vaccine efficacy, to avoid discontinuation of product development owing to a false negative error.

Maximizing biological relevance and statistical power for assessing genotype-specific vaccine efficacy requires careful definition of the mark variable measuring genetic divergence of an exposing HIV to the vaccine strains. In general, viral vaccines stimulate protective immune responses to the exact protein sequences of the viruses inside the tested vaccine, but fail to stimulate protective immune responses to viruses with certain differences in their genetic sequence compared to the vaccine strains. Therefore, the genetic distance variable is defined in a way to reflect a biological model of how well the vaccine can induce immune responses that ‘cross-react’ with viruses different from the vaccine strains; under the model for cross-reactivity the vaccine is more likely to stimulate a protective immune response to HIVs with smaller distances to the vaccine strains. Thus, the statistical assessment of if and how VE varies with viral genetic distance empirically tests the fidelity of the distance for measuring cross-reactive protection. Because of challenges involved in modeling cross-reactivity, multiple immunologically meaningful genetic distances have been considered, based on (i) different conceptual definitions, (ii) different HIV sequence regions, (iii) different methods to predict which HIV amino acids stimulate antibody responses, and (iv) different reference sequences inside the vaccine.

From 2003 to 2009, a randomized Phase III preventive vaccine efficacy trial (Rerks-Ngarm et al., 2009), the RV144 ‘Thai trial’, was conducted to evaluate the combination of the recombinant canarypox vector vaccine (ALVAC-HIV [vCP1521]), expressing the HIV genes gag, pro, and the gp120 portion of env, and the recombinant gp120 subunit vaccine (AIDSVAX B/E). HIV-uninfected volunteers in Thailand, primarily at heterosexual risk for HIV infection, were randomly allocated to receive the vaccine (n1 = 8197) or placebo (n0 = 8198), and followed for 42 months for the primary endpoint of HIV infection. HIV was diagnosed in 125 subjects (51 vaccine, 74 placebo), and, based on the Cox model, overall vaccine efficacy was estimated to be 31.2% (95% CI, 1.1 to 52.1; p=0.04). This result generated great interest to understand how vaccine protection may have depended on the genetic divergence of infecting HIVs from the HIV vaccine strains. Full genome HIV sequences were measured from 121 of the 125 infected subjects, with between 2 and 13 sequences sampled per subject. In Section 5, we describe the immunologically relevant genetic distance marks that we analyze.

GMS developed a nonparametric estimator for the univariate mark-specific hazard ratio and procedures for testing whether the mark-specific hazard ratio is (i) unity or (ii) independent of the mark. Problems (i) and (ii) are equivalent to the null hypothesis of zero vaccine efficacy against any exposing virus and the null hypothesis that vaccine efficacy does not depend on the viral divergence, respectively; both null hypotheses are also considered in Section 3. Sun, Gilbert, and McKeague (2009) developed the mark-specific proportional hazards model, which, in the case of a univariate mark, provides an alternative approach to inference about the mark-specific vaccine efficacy defined in (2).

We propose a more efficient method of estimation and hypothesis testing for the mark-specific hazard ratio in the framework of competing risks failure time analysis that accommodates multivariate marks. Our approach makes use of the maximum profile likelihood estimation method in the semiparametric density ratio/biased sampling model developed by Qin (1998). In parallel with Qin (1998), a similar method of maximum profile partial likelihood estimation was derived by Gilbert, Lele, and Vardi (1999) and Gilbert (2000) for semiparametric biased sampling models with K possibly biased samples.

1.1 Notation

Let T denote the time to failure and V [set membership] Rs, a continuous, possibly multivariate, mark variable. Without loss of generality, the support of each component of V is taken to be [0, 1]. Let C be the time to censoring. The observed right-censored failure time is X = min(T, C) with the failure indicator δ = I(TC). The mark V is always observed if δ = 1; otherwise it is unobservable. Let Z denote the indicator of assignment to the treatment group (in vaccine trials, Z = 1 indicates vaccine and Z = 0 indicates placebo). Let (Zi, Xi, δi, Vi), i = 1, …, n, be i.i.d. replicates of (Z, X, δ, V). The observed data consist of the observations (Zi, Xi, Vi) for individuals with δi = 1 and the observations (Zi, Xi) for those with δi = 0.

1.2 Assumptions

We assume C is conditionally independent of both T and V given Z, that is, C [perpendicular] T|Z and C [perpendicular] V|Z, and assume T [perpendicular] V|Z. These assumptions ensure identifiability of the Euclidean parameters in the density ratio model introduced in Section 2. In the HIV vaccine trial setting, the assumption T [perpendicular] V|Z entails that, for example, vaccine recipients with infection time T = 6 months have the same distribution of the mark V as vaccine recipients with infection time T = 2.5 years, which may approximately hold given a limited shift in the HIV sequence distribution over the period of 2 years. An analogous statement comparing vaccine recipients with infection times, say, 50 years apart would be clearly incorrect due to the genetic shift of HIV. Hence, the fact that HIV vaccine efficacy trials only last between 3–5 years is crucial for the assumption to be (approximately) met. In Section 3.1, we develop a Kolmogorov–Smirnov-type test for assessing validity of the assumption T [perpendicular] V|Z, and, in Section 4.3, we demonstrate some robustness properties of the proposed inferential procedures to violation of the assumption T [perpendicular] V|Z.

1.3 Estimands

We define the conditional multivariate mark-specific hazard function as


which is a natural generalization of the cause-specific hazard function in the presence of finitely many causes of failure (Prentice et al., 1978). GMS defined the mark-specific vaccine efficacy as


As described in GMS, V E(t, υ) has an approximate interpretation as the multiplicative vaccine effect to reduce susceptibility to HIV infection given exposure to strain υ at time t.

The remainder of this article is organized as follows. In Section 2, we introduce a semiparametric model for the mark-specific hazard ratio in (2), discuss identifiability of the Euclidean model parameters, and state the asymptotic properties of the proposed estimator for the mark-specific hazard ratio. Hypothesis testing procedures are introduced in Section 3. We summarize our findings from the simulation experiment in Section 4 and apply the proposed inferential procedures to data from the Thai trial in Section 5. Proofs of main results are available in the Web Appendix.

2. Estimation method

The mark-specific hazard function factors as


where f(·|T = t, Z = z) is the conditional density function of V given T = t and Z = z, and λ(·|Z = z) is the ordinary marginal hazard function not involving the mark. Consequently, the mark-specific hazard ratio can be written as


The factorization in (4) is advantageous because the two ratios can be estimated separately. For the mark density ratio, we consider the semiparametric density ratio model (Qin, 1998)


where ϕ [set membership] Rd is the parameter of interest and g(υ, ϕ) is a known time-independent weight function, continuously differentiable in ϕ. The nuisance parameter distribution f(υ|T = t, Z = 0) is time-independent following the assumption T [perpendicular] V|Z, and is treated nonparametrically. Due to the equality of the conditional density functions


which occurs because of the assumption (T, C) [perpendicular] V|Z, the parameter ϕ in (5) is estimable using cases with an observed failure only. Using Bayes rule, model (5) implies that


for marks υ1, υ2 [set membership] [0, 1]s. Thus g(υ1, ϕ)/g(υ2, ϕ) is the odds of being assigned vaccine for an individual infected with strain υ1, relative to that for an individual infected with strain υ2.

A common choice of the weight function g(υ, ϕ) is exp{α + g(υ, β)} where ϕ = (α, βT)T and g(υ, β) is a polynomial function. In this case, the following proposition characterizes the necessary and sufficient condition for identifiability of ϕ (see Web Appendix A for a proof).

Proposition 1: In the semiparametric density ratio model


T and V are conditionally independent given Z if and only if α(t) is constant.

Following Qin (1998), the profile score functions for the parameter of interest ϕ and the Lagrange multiplier λ [set membership] [0, 1] (that arises by profiling out the nuisance parameter) are


where ġ(υ, ϕ) = [partial differential]g(υ, ϕ)/[partial differential]ϕ. The maximum profile likelihood estimator ([phi with hat]T, [lambda with circumflex])T for (ϕT, λ)T is defined as the solution to the system Uϕ(ϕ, λ) = 0 and Uλ(ϕ, λ) = 0.

For the marginal hazard ratio in (4), we propose to use the Cox regression model


The parameter γ may be estimated with the standard maximum partial likelihood estimator (MPLE), or, alternatively, with the more efficient method of Lu and Tsiatis (2008) that leverages auxiliary data predictive of the failure time (we implemented the Lu and Tsiatis method in the R speff2trial package). If additional, possibly time-dependent, covariates Z*(t)=(Z1*(t),,Zp*(t))T are measured, then the Cox model (7) can be extended to λ(t|Z = 1, Z*(t) = z*(t)) = λ(t|Z = 0, Z*(t) = 0)eγ+γ*Tz*(t), and any estimation method for (γ, γ*T)T available for the Cox model can be employed. For example, using Z*(t) = (t, Zt)T specifies λ(t|Z=1,t)/λ(t|Z=0,t)=eγ+γ2*t, allowing the overall treatment effect to change over time.

2.1 Asymptotic properties of the proposed estimator

Qin (1998) showed that ([phi with hat]T, [lambda with circumflex])T is a consistent and asymptotically normal (CAN) estimator if the proportions of failures in the placebo and treatment group converge to positive constants. The MPLE [gamma with circumflex] and the Lu and Tsiatis (2008) estimator have also been shown to be CAN. The present goal is to characterize the joint asymptotic distribution of ([phi with hat]T, [lambda with circumflex], [gamma with circumflex])T.

We restrict attention to the density ratio model (5) with weight function g(υ, ϕ) = eα+g(υ,β), where ϕ = (α, βT)T and g(υ, β) is a polynomial function in υ, and the marginal Cox model adjusted for the covariates (Z, Z*T (t))T. Although Theorem 1 applies in this setting, for expositional simplicity the results are presented for g (υ, β) a quadratic form in υ and the marginal Cox model (7) with sole covariate the treatment group indicator Z. In this special case, the mark-specific hazard ratio model in (4) takes the form eα+g(υ,β)+γ. Let θ = (ϕT, λ, γ)T denote the vector of Euclidean parameters, and let the random map Ψn(θ) denote the set of all estimating functions for θ. The estimator [theta w/ hat]n for θ is obtained as the solution to Ψn(θ) = 0.

Next, define the random processes


where τ is selected such that P(T > τ) > 0. We will use the notation η·,γ or η·,θ to denote the processes indexed by the component γ or by the full parameter vector θ. For (x, Δ, υ, z) [set membership] [0, τ] × {0, 1} × [0, 1]s × {0, 1}, define the map φ=(φ1T,φ2,φ3)T as


The dependence of [var phi] on (x, Δ, υ, z) is suppressed in the notation. Let Ψ be the limit of Ψn as n → ∞ (and also m → ∞) and let Ψ(θ0) = 0. Let GP denote the zero-mean P–Brownian bridge process and define the class of functions fθ,t,r(x, z) = zrI(xt)eγz for r = 0, 1. Theorem 1 characterizes the asymptotic distribution of n(θ̂nθ0) (see Web Appendix B for a proof).

Theorem 1: Let [theta w/ hat]n be the solution to Ψn(θ) = 0 and let Ψ(θ0) = 0. Then




with the map φ=(φ1T,φ2,φ3)T defined in (8), pδ = P(δ = 1), and


Here F(t|δ = 1) is the conditional cdf of X given δ = 1, and [Psi with dot above]θ0 is the continuously invertible derivative of the map θ [mapsto] Ψ(θ) at θ0 and has matrix form


with entries


where g(υ, ϕ) = [partial differential]2g(υ, ϕ)/[partial differential]ϕ[partial differential]ϕT.

Denote the column vector of functions φ̃(θ,η0,θ)=(φ1T(θ,η0,θ),φ2(θ,η0,θ),φ3(θ,η0,θ)+pδlθ)T. The following corollary describes the asymptotic variance of n(θ̂nθ0).

Corollary 1: The asymptotic random vector Ψ̇θ01Z in Theorem 1 is normally distributed with zero mean and covariance matrix Γ=Ψ̇θ01ΩΨ̇θ01 where


Let [Gamma]n denote the empirical estimator for Γ obtained by replacing P by the empirical probability measure Pn, θ0 by [theta w/ hat]n, and η0,θ0 by ηn,[theta w/ hat]n in the definition of Γ. Corollary 1 leads to the construction of Wald confidence intervals (pointwise in υ) for the components of θ, and, subsequently, for the parameter V E(υ) = 1 − eα+g(υ,β)+γ.

3. Hypothesis testing

For illustrating the testing procedures, we consider the simplified weight function g(υ, ϕ) = eα+βTυ, ϕ = (α, βT)T, which leads to the mark-specific vaccine efficacy function V E(υ) = 1 − eα+βTυ. We develop likelihood ratio and Wald tests to evaluate the null hypothesis

H00:VE(υ)=0for allυ[0,1]s,

which states that the vaccine provides no protection against infection with any HIV strain. If H00 is rejected, the question arises as to whether vaccine efficacy depends on the viral divergence; thus, we develop the likelihood ratio and Wald tests for the null hypothesis

H0:VE(υ)VEfor allυ[0,1]s.

Under models (5) and (7), H00 is equivalent to {00 : β = 0 and γ = 0}. The likelihood ratio test of 00 against {10 : β0 or γ ≠ 0} uses Simes (1986) procedure, in which the profile likelihood ratio test statistic for β in model (5) and the partial likelihood ratio test statistic for γ in model (7) are evaluated separately. P-values pβ and pγ are obtained based on the fact that the likelihood ratio statistics are asymptotically χs2 and χ12 under 00, respectively. Simes’ procedure rejects H00 if either max(pβ, pγ) ≤ α or min(pβ, pγ) ≤ α/2 where α is the nominal familywise level of significance. The Wald test of H00 versus H10 is based on the statistic n(β̂nT,γ̂n)Γ̂n,βγ1(β̂nT,γ̂n)T where [Gamma]n,βγ is the submatrix of [Gamma]n pertaining to the components β and γ. Under H00, the Wald test statistic is asymptotically χs+12. We additionally propose a weighted one-sided Wald-type test of H00 based on the Z-statistic


which is designed to increase power to detect alternative hypotheses where both the marginal vaccine efficacy VE=1λ(t|Z=1)λ(t|Z=0)=1eγ>0 and V E(υ) declines with all of the components of V. The null hypothesis H0 is equivalent to {H0 : β = 0}, and thus the density ratio model (5) alone serves to construct the likelihood ratio and Wald test statistics for H0.

3.1 Diagnostic test for T [perpendicular] V|Z

By Proposition 1, the T [perpendicular] V|Z assumption is necessary for parameter identifiability in the time-independent density ratio model. For a given treatment group Z = z, we propose a diagnostic test of the null hypothesis {K0 : T [perpendicular] V|Z = z} based on the statistic


where FTV (t, υ|z) is the nonparametric maximum likelihood estimator of the joint cdf of (T, V) for group z developed by Huang and Louis (1998), FT (t|z) is one minus the Kaplan-Meier estimator of the conditional survival function of T, and FV (υ|z) is the empirical distribution function of the observed values of V. The critical values for the distribution of (12) under K0 can be calculated using a bootstrap algorithm as follows:

  1. Draw an independent sample (Xi*,δi*) from the original time-on-study data (Xi, δi) for group Z = z subjects, with replacement.
  2. Draw a sample Vi* from the original mark data Vi of group Z = z subjects who became infected (δ = 1), with replacement.
  3. Compute Dz based on the bootstrap data (Xi*,δi*,δi*Vi*) for subjects i in group Z = z.
  4. Repeat Steps (1)(3) B times.
  5. Estimate the critical value as the α–quantile of Dz across the B replicates.

The bootstrap test has the nominal size asymptotically and is consistent against all alternatives (Romano, 1988). We test the overall hypothesis {K0 : T [perpendicular] V |Z} using Simes’ procedure applied to the tests of K0 performed separately for the two groups Z = 1 and Z = 0.

4. Simulation study

4.1 Assessment of the proposed methods under model validity

The simulation setup mimics typical placebo-controlled HIV vaccine efficacy trials. Initially, we investigate finite-sample performance of the V E(υ) estimation and testing procedure in the model V E(υ) = 1 − eα+βυ+γ, υ [set membership] [0, 1], i.e., the data are generated following models (5) with g(υ, α, β) = eα+βυ and (7). We compare the performance of the estimator to that of the semiparametric estimator of Sun, Gilbert, and McKeague (SGM) (2009), and compare the performance of the Wald and likelihood ratio tests of H00 and H0 to those of GMS.

We specify the failure times T to be exponential with rates λ and λeγ in the placebo and vaccine groups. The rate λ = log(0.85)/(−3) is chosen so that 15% of placebo recipients are infected by 3 years. We specify the censoring times C to be Uniform(0, 15) in each group which implies a 20% chance of censoring by 3 years. The observed time on study X is defined as min(T, C, 3). The vaccine-to-placebo assignment ratio is 1:1. A univariate mark V for placebo and vaccine recipients is generated from distributions with density functions f(υ|Z = 0) = {2e−2υ/(1 − e−2)} I(0 ≤ υ ≤ 1) and f(υ|Z = 1) = f(υ|Z = 0)eα+βυ, where, for a given value of β, the value of the parameter α = α(β) is defined as the solution to 01f(υ|Z=0)eα+βυdυ=1. Only the mark values for subjects with T ≤ 3 and δ = 1 would be observed in a real study and hence are used in the analysis.

We consider sample sizes n = 556, 741, and 1481 per arm so that the expected numbers of observed placebo infections by year 3 are NpI = 75, 100, and 200, respectively. The simulation scenarios are characterized by the following values of (β, γ): (0, 0), (0.5, −0.8), (1.2, −0.2), and (2.1, −1.3). The true V E(υ) functions are depicted in Web Figure 1; in subsequent tables and figures, the simulation scenarios are labeled by the pertaining V E(0) and V E(1) values. For all simulations the nominal significance level is taken to be 5% for two-sided tests and 2.5% for one-sided tests. The results are based on 104 replicated data sets.

Figures 1, ,2,2, and Web Figure 2 characterize finite-sample bias of estimates and asymptotic standard error estimators of VE^(υ) as well as coverage probabilities of 95% pointwise Wald confidence intervals for V E(υ). In each simulation scenario, bias is the smallest in the central support region of the mark and slightly increases in the right tail. As expected, bias decreases with an increasing number of events. Web Figure 4 compares bias of the proposed and the SGM estimator. Asymptotic and empirical standard error estimates are consistently in good accordance with one another, and Web Figure 5 illustrates a substantial efficiency gain relative to the SGM method in this setting. Coverage probabilities for all mark values approach the nominal confidence level with an increasing number of failures.

Figure 1
Finite-sample bias of estimation of V E(υ) = 1 − eα+βυ+γ with a univariate mark. NpI and NυI denote the expected numbers of observed placebo and vaccine infections, respectively. This figure appears ...
Figure 2
Asymptotic versus empirical (dotted) standard error estimates of VE^(υ) with a univariate mark. NpI and NυI denote the expected numbers of observed placebo and vaccine infections, respectively. This figure appears in color in the electronic ...

Tables 1 and and22 summarize the observed sizes and powers of the Wald and likelihood ratio tests of H00 and H0, and include a comparison with the alternative non- and semiparametric tests of GMS. The sizes of all considered tests are in good agreement with the nominal significance level. As for two-sided testing of H00, powers of the Wald and likelihood ratio test are comparable, and both tests outperform the Û13 and Û14 tests of GMS, with power gains up to 19%, for all specified values of the model parameters. The one-sided weighted Wald-type test (11) of H00 is designed to be sensitive to alternative hypotheses with positive marginal vaccine efficacy and V E(υ) declining in υ, for which it is conjectured to be more powerful than the standard log-rank test that assesses the marginal vaccine efficacy only; this conjecture is verified in the results with power gain reaching up to 27%.

Table 1
Powers of tests of the null hypothesis {H00 : V E (υ) = 0 for all υ [set membership] [0,1]} against (i) the alternative hypothesis {H01 : non H00}, and (ii) the alternative hypothesis {H20 : V E > and V E(υ) is a decreasing function ...
Table 2
Powers of two-sided tests of the null hypothesis {H0 : VE (υ) [equivalent] VE for all υ [set membership] [0,1]}. The proposed Wald and likelihood ratio tests are compared with the alternative non- and semiparametric tests Û2np and Û ...

Powers of the Wald and likelihood ratio tests of H0 are considerably higher (up to 36%) than those of GMS’s Û2np and Û2sp tests. We conjecture that the power gain can be ascribed to the additional assumptions made by our methods compared to those of GMS. Web Table 1 supports nominal size of the diagnostic test (12) of {K0 : T [perpendicular] V|Z}.

4.2 Simulation scenarios with a bivariate mark variable

We additionally study simulation scenarios involving a bivariate mark V = (V1, V2), where V1 and V2 are generated as independent univariate marks from distributions described in Section 4.1. The point and interval estimator for bivariate mark-specific vaccine efficacy exhibits similar properties as in scenarios with a univariate mark (Web Appendix D). For bivariate marks, the likelihood ratio tests of H00 and H0 are slightly more powerful than the respective Wald tests, and the diagnostic test (12) of the null hypothesis {T [perpendicular] (V1, V2)|Z} has nominal size.

4.3 Robustness analysis of the proposed methods under model mis-specification

Assuming the model V E(υ) = 1 − eα+βυ+γ, we examine robustness of the V E(υ) estimator and of the Wald and likelihood ratio tests of H00 and H0 to violation of the T [perpendicular] V|Z assumption. To this end, we generate samples (T, V) with correlation ρTV. For an infected placebo recipient, we first generate T = t, and second generate V from the conditional density function


where c governs the magnitude of ρTV. Using Proposition 1, for an infected vaccine recipient with T = t, we generate V from the conditional density function f(υ|T = t, Z = 1; c) = f(υ|T = t, Z = 0; c)eα(t)+βυ, where α(t) = α(t, β, c) is the solution to 01f(υ|T=t,Z=0;c)eα(t)+βυdυ=1. The data generating mechanism implies that the true V E(t, υ) = 1 − eα(t)+βυ+γ. The values of c are chosen such that ρTV varies over 0.1–0.5. In each scenario, power to reject the null hypothesis {T [perpendicular] V|Z} using the diagnostic test (12) is evaluated.

Web Figure 3 depicts the Monte Carlo average deviation VE^(υ)VE¯(υ) where VE¯(υ)=03VE(t,υ)f(t)dt/03f(t)dt and f(t) is the density function of the failure time T. For scenarios where the magnitude of correlation is detected with moderate power, the deviation is minimal. Such scenarios are of particular interest because, in these cases, violation of the assumption T [perpendicular] V|Z may frequently remain undetected. For scenarios with power to detect departures from T [perpendicular] V|Z greater than 90%, the deviation tends to increase for mark values in the right tail of the distribution, although this is fairly innocuous because these are the scenarios where the method would fail the diagnostics and hence would not be used.

Web Tables 2 and 3 indicate that the likelihood ratio and Wald tests of H00 and H0 retain correct size under violation of T [perpendicular] V|Z. This robustness property holds for any violation of this assumption, because under either null hypothesis H00 and H0, α(t) [equivalent] α = 0 for all t > 0, i.e., the estim and V E(υ) in the time-independent model coincides with the true V E(t, υ) characterizing the data generating mechanism regardless of the value of ρTV. Web Tables 2 and 3 show that power of the likelihood ratio and Wald tests is minimally sensitive to correlation between T and V and tends to be higher than that of GMS’s tests.

Finally, we performed a robustness analysis with respect to violation of the proportional marginal hazards assumption in which we considered distributions of T|Z = 0 and T|Z = 1 characterized by the hazard functions λ(t|Z = 0) [equivalent] λ and λ(t|Z = 1) = λeγ01t; thus γ1 governs the amount of variation in the marginal hazard ratio over time. We observed minimal differences in the patterns of average deviation VE^(υ)VE¯(υ), standard error estimates, and power relative to those presented for correctly specified models (results not reported). We conjecture that the similar performance occurs because V E(t, υ) = 1 − eα+βυ+γ01t does not include an interaction of V with T, and hence V E(t, υ) and V E(υ) have the same amount of change in υ.

5. Application

5.1 Genetic distance definition

In the RV144 Thai trial, full length HIV genomes were measured from 121 of the 125 infected subjects; 3 are missing data because their viral load was too low for the Sanger sequencing technology to work, and 1 dropped out. The 4 subjects with missing data are excluded from the analysis. The recent RV144 ‘immune correlates’ study showed that vaccine recipients with high levels of IgG antibodies that bind to the V1/V2 region of the HIV Envelope gp120protein had a significantly lower rate of HIV infection (Haynes et al., 2012). To help discern the potential role of V1/V2 antibody epitopes as an immunologic cause of the observed partial vaccine efficacy, we focus the genetic distance definitions on the V1/V2 region (85 amino acids long). Three HIV gp120 sequences were included in the vaccine: 92TH023 in ALVAC canarypox vector, and A244, MN in AIDSVAX gp120. 92TH023 and A244 are subtype E-viruses whereas MN is subtype B, and 110 (91%) of the 121 subjects were infected with subtype E viruses. The E viruses are closer genetically to the infecting sequences than MN, and thus are more likely to stimulate protective immune responses. Accordingly, our analysis focuses on the 92TH023 and A244 insert sequences, and excludes the 11 subjects with non-E infecting sequences.

To maximize biological relevance and statistical power of the analysis, V1/V2 distances were computed only including the subset of the 85 amino acids that were predicted to be part of antibody epitopes that were stimulating vaccine-induced responses. Web Appendix E describes this process, which yielded q = 9 important amino acid sites that were used for defining the distances.

For each of the two insert reference sequences and each sequence from a subject, the genetic distance was computed using the Blosum90 scoring matrix, adjusted such that the addition or subtraction of an N-linked glycosylation site is scored 3-fold or 1.5-fold higher than the largest Blosum90 score depending on whether the mismatch is in position 2 of the aminoacid triple comparing the N-linked glycosylation site. First, a distance was computed as


where ais is the amino acid at site i in the subject’s sequence, air is the amino acid at site i in the reference sequence, and s(x, y) is the score from the modified Blosum90 matrix. Second, the distances were re-scaled to values between 0 and 1 (2 identical sequences have distance0). Between 2 and 13 sequences (total 1030 sequences) were measured per infected subject, and the final distance was defined as the subject-specific median distance. This yields a data set of a bivariate mark (V1, V2) measured from 110 infected subjects, with 49 and 52 unique mark values in V1 and V2, respectively. Figure 3 shows the distributions of V1 and V2.

Figure 3
RV144 trial: distribution of the distances to the 92TH023 and A244 vaccine insert sequence by vaccine/placebo group. This figure appears in color in the electronic version of this article.

5.2 Inference about mark-specific vaccine efficacy

Point and 95% pointwise confidence interval estimates of V E(υ) are plotted in Figures 4a–b. The weighted Wald test of H00 : V E(υ) [equivalent] 0 yields p-values of 0.016 and 0.024 for the 92TH023 and A244 distances. The likelihood ratio test of H0 : V E(υ) = V E yields p = 0.52 for the 92TH023 and p = 0.83 for the A244 distance. These p-values are consistent with the broad confidence bands in Figures 4a and 4b. The Thai trial had low power to detect a change in V E(υ) with viral divergence. The assumption of conditional independence between T and V given Z was verified using the Kolmogorov–Smirnov-type test, with p-values of 0.88and 0.98 for the placebo and vaccine groups, based on 1000 bootstrap iterations. The Qin and Zhang (1997) goodness-of-fit test does not reject the validity of the mark density ratio model for either the 92TH023 (p = 0.49) or the A244 distance (p = 0.51). Moreover, a test of proportional hazards assumed in (7) (Grambsch and Therneau, 1994) does not reject(p = 0.21).

Figure 4
Estimated mark-specific vaccine efficacy for (a) the univariate distance to the 92TH023 insert sequence, (b) the univariate distance to the A244 insert sequence, and (c) the bivariate distance with components defined in (a) and (b). 95% pointwise confidence ...

Figure 4c displays VE^(υ1,υ2). The weighted Wald test of H00 yields p = 0.019. The likelihood ratio test of H0 yields p = 0.81. The p-value for the test of T [perpendicular] (V1, V2)|Z is 0.64 and 0.89 in the placebo and vaccine groups again supporting the appropriateness of the method. The goodness-of-fit test does not reject the validity of the bivariate-mark density ratio model (p = 0.46).

6. Discussion

This article develops inferential procedures for the mark-specific hazard ratio model with a multivariate continuous mark variable. We show that the mark-specific hazard ratio factors as the mark density ratio and the ordinary marginal hazard ratio. This factorization is appealing because it enables the two components to be estimated separately. The semiparametric method of maximum profile likelihood estimation in the mark density ratio model yields improved efficiency of mark-specific hazard ratio estimation relative to current alternative approaches. We derive the joint asymptotic distribution of the estimators for the mark density ratio and the marginal hazard ratio. This result allows estimation of mark-specific vaccine efficacy (V E(υ)) with pointwise confidence intervals and testing of hypotheses of interest about V E(υ). The proposed weighted Wald-type test is recommended for testing V E(υ) [equivalent] 0 for all υ, and the likelihood ratio test for testing V E(υ) invariant in υ.

Due to growing complexity in HIV vaccine design, it becomes increasingly desirable to evaluate V E(υ) as a function of a multivariate mark representing a set of genetic distances. Our approach allows flexible semiparametric modeling of a multivariate mark variable, and the estimation method does not impose practical constraints on the mark dimension.

A limitation of the methods is that they do not apply if the T [perpendicular] V|Z assumption fails; in this scenario, an alternative approach may be used that is based on a different factorization of the mark-specific hazard function,


with f(t, υ|Z = z) the conditional joint density function of (T, VT)T given Z = z, and S(t|Z = z) the conditional marginal survival function of T given Z = z. This approache liminates the T [perpendicular] V|Z assumption as well as the proportional marginal hazards assumption in exchange for flexible parametric modeling of the failure time in the density ratio model. It would be of interest to develop estimation and testing methods for this approach.

For a univariate mark variable, the present method is an alternative to the SGM method, which estimates the same V E(υ) estimand and tests the same null hypotheses. The following issues are relevant for the choice of method. First, the present method makes two assumptions not made by SGM: T [perpendicular] V|Z and a parametric form for the mark density ratio. Therefore, where the goodness-of-fit tests reject either of these assumptions, the SGM method is more appropriate. However, if the data support both assumptions, then the present method maybe preferred for its substantially improved efficiency, a key issue given the concern for false negative errors. Second, SGM assumes a mark-specific proportional hazards model, such that V E(t, υ) is independent of t. The present method partially relaxes this assumption via the factorization 1 − V E(t, υ) = g(v, ϕ) {1 − V E(t)}, which allows flexible modeling of V E(t), for example via the Cox model with time-dependent covariates, which allows V E(t, υ) to vary with time. Third, the SGM method uses kernel smoothing to estimate V E(υ) over the range of mark values, thereby providing more flexible estimation. However, issues of band width selection and complicated behavior of the estimator in the tails of the mark variable make the SGM method more complicated to implement.

Although the presented methods are motivated by a specific scientific application, they provide a general approach to the analysis of survival data in the presence of a continuous, possibly multivariate, mark variable. The R code implementing the proposed methods is available upon request and is posted at the second author’s webpage.

Supplementary Material

Supp Web Appendix


The authors thank the reviewers and Jon Wellner for helpful suggestions, and thank the participants, investigators, and sponsors of the RV144 Thai trial, including the U.S. Military HIV Research Program (MHRP); U.S. Army Medical Research and Materiel Command; NIAID; U.S. and Thai Components, Armed Forces Research Institute of Medical Science Ministry of Public Health, Thailand; Mahidol University; SanofiPasteur; and Global Solutions for Infectious Diseases. This research was supported by NIH NIAID grant 2 R37 AI054165-10.


Supplementary Materials

Web Appendices referenced in Sections 2–5 are available with this paper at the Biometrics website on Wiley Online Library.


  • Gilbert P. Large sample theory of maximum likelihood estimates in semiparametric biased sampling models. Annals of Statistics. 2000;28:151–194.
  • Gilbert P, Lele S, Vardi Y. Maximum likelihood estimation in semiparametric selection bias models with application to AIDS vaccine trials. Biometrika. 1999;86:27–43.
  • Gilbert P, McKeague I, Sun Y. The 2-sample problem for failure rates depending on a continuous mark: an application to vaccine efficacy. Biostatistics. 2008;9:263–276. [PubMed]
  • Grambsch P, Therneau T. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika. 1994;81:515–526.
  • Halloran M, Struchiner C, Longini I. Study designs for evaluating different efficacy and effectiveness aspects of vaccines. American Journal of Epidemiology. 1997;146:789–803. [PubMed]
  • Haynes B, Gilbert P, McElrath M, Zolla-Pazner S, Tomaras G, Alam S, Evans D, Montefiori D, Karnasuta D, Sutthent R, Liao H, DeVico A, Lewis G, Williams C, Fong Y, Janes H, DeCamp A, Huang Y, Rao M, Billings E, Karasavvas N, Robb M, Ngauy V, de Souza M, Paris R, Ferrari G, Bailer R, Soderberg K, Andrews C, Berman P, Frahm N, De Rosa S, Alpert M, Yates N, Shen X, Koup R, Pitisuttithum P, Kaewkungwal J, Nitayaphan S, Rerks-Ngarm S, Michael N, Kim J. Immune-correlates analysis of an HIV-1 vaccine efficacy trial. New England Journal of Medicine. 2012;366:1275–1286. [PMC free article] [PubMed]
  • Huang Y, Louis T. Nonparametric estimation of the joint distribution of survival time and mark variables. Biometrika. 1998;85:785–798.
  • Lu X, Tsiatis A. Improving the efficiency of the log-rank test using auxiliary covariates. Biometrika. 2008;95:679–694.
  • Prentice R, Kalbfleisch J, Peterson A, Flournoy N, Farewell V, Breslow N. The analysis of failure times in the presence of competing risks. Biometrics. 1978;34:541–554. [PubMed]
  • Qin J. Inferences for case-control and semiparametric two-sample density ratio models. Biometrika. 1998;85:619–630.
  • Qin J, Zhang B. A goodness-of-fit test for logistic regression models based on case-control data. Biometrika. 1997;84:609–618.
  • Rerks-Ngarm S, Pitisuttithum P, Nitayaphan S, Kaewkungwal J, Chiu J, Paris R, Premsri N, Namwat C, de Souza M, Adams E, Benenson M, Gurunathan S, Tartaglia J, McNeil J, Francis D, Stablein D, Birx D, Chunsuttiwat S, Khamboonruang C, Thongcharoen P, Robb M, Michael N, Kunasol P, Kim J. for the MOPH-TAVEG Investigators. Vaccination with ALVAC and AIDSVAX to Prevent HIV-1 Infection in Thailand. New England Journal of Medicine. 2009;361:2209–2220. [PubMed]
  • Romano J. A bootstrap revival of some nonparametric distance tests. Journal of the Americal Statistical Association. 1988;83:698–708.
  • Simes R. An improved Bonferroni procedure for multiple tests of significance. Biometrika. 1986;73:751–754.
  • Sun Y, Gilbert P, McKeague I. Proportional hazards models with continuous marks. Ann. Stat. 2009;37:394–426. [PMC free article] [PubMed]