Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3940058

Formats

Article sections

- Summary
- 1. Introduction
- 2. Estimation method
- 3. Hypothesis testing
- 4. Simulation study
- 5. Application
- 6. Discussion
- Supplementary Material
- References

Authors

Related links

Biometrics. Author manuscript; available in PMC 2014 March 3.

Published in final edited form as:

PMCID: PMC3940058

NIHMSID: NIHMS550353

M. Juraska: ude.notgnihsaw.u@aksarujm

The publisher's final edited version of this article is available at Biometrics

See other articles in PMC that cite the published article.

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.

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 (*n*_{1} = 8197) or placebo (*n*_{0} = 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.

Let *T* denote the time to failure and *V*^{s}, a continuous, possibly multivariate, mark variable. Without loss of generality, the support of each component of ** V** is taken to be [0, 1]. Let

We assume *C* is conditionally independent of both *T* and ** V** given

We define the conditional multivariate mark-specific hazard function as

$$\lambda (t,\mathit{\upsilon}|Z=z)=\underset{{h}_{1},{h}_{21},\dots ,{h}_{2s}\searrow 0}{\text{lim}}\frac{P(T\in [t,t+{h}_{1}),\mathit{V}\in {\prod}_{i=1}^{s}[{\upsilon}_{i},{\upsilon}_{i}+{h}_{2i})|T\ge t,Z=z)}{{h}_{1}{h}_{21}\cdots {h}_{2s}},$$

(1)

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

$$VE(t,\mathit{\upsilon})=1-\lambda (t,\mathit{\upsilon}|Z=1)/\lambda (t,\mathit{\upsilon}|Z=0).$$

(2)

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.

The mark-specific hazard function factors as

$$\lambda (t,\mathit{\upsilon}|Z=z)=f(\mathit{\upsilon}|T=t,Z=z)\lambda (t|Z=z),$$

(3)

where *f*(·|*T* = *t, Z* = *z*) is the conditional density function of ** V** given

$$\frac{\lambda (t,\mathit{\upsilon}|Z=1)}{\lambda (t,\mathit{\upsilon}|Z=0)}=\frac{f(\mathit{\upsilon}|T=t,Z=1)}{f(\mathit{\upsilon}|T=t,Z=0)}\frac{\lambda (t|Z=1)}{\lambda (t|Z=0)}.$$

(4)

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)

$$f(\mathit{\upsilon}|T=t,Z=1)/f(\mathit{\upsilon}|T=t,Z=0)=g(\mathit{\upsilon},\mathit{\varphi})$$

(5)

where **ϕ** ^{d} 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* ** V**|

$$f(\mathit{\upsilon}|T=t,Z=z)=f(\mathit{\upsilon}|T=t,Z=z,\delta =1),$$

(6)

which occurs because of the assumption (*T, C*) ** V**|

$$\frac{g({\mathit{\upsilon}}_{1},\mathit{\varphi})}{g({\mathit{\upsilon}}_{2},\mathit{\varphi})}=\frac{f({\mathit{\upsilon}}_{1}|Z=1)}{f({\mathit{\upsilon}}_{1}|Z=0)}{\left\{\frac{f({\mathit{\upsilon}}_{2}|Z=1)}{f({\mathit{\upsilon}}_{2}|Z=0)}\right\}}^{-1}=\frac{P(Z=1|{\mathit{\upsilon}}_{1})}{P(Z=0|{\mathit{\upsilon}}_{1})}{\left\{\frac{P(Z=1|{\mathit{\upsilon}}_{2})}{P(Z=0|{\mathit{\upsilon}}_{2})}\right\}}^{-1}$$

for marks **υ**_{1}, **υ**_{2} [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{α + (**υ, β**)} where **ϕ** = (α, **β**^{T})^{T} and (**υ, β**) 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

$$f(\mathit{\upsilon}|T=t,Z=1)/f(\mathit{\upsilon}|T=t,Z=0)={\mathrm{e}}^{\alpha (t)+\mathrm{g\u0303}(\mathit{\upsilon},\mathit{\beta})},$$

*T* and ** V** are conditionally independent given

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

$${\mathit{U}}_{\mathit{\varphi}}(\mathit{\varphi},\lambda )=\sum _{i=1}^{n}{\delta}_{i}(-\frac{\lambda \mathit{\u0121}({\mathit{V}}_{i},\mathit{\varphi})}{1+\lambda (g({\mathit{V}}_{i},\mathit{\varphi})-1)}+{Z}_{i}\frac{\mathit{\u0121}({\mathit{V}}_{i},\mathit{\varphi})}{g({\mathit{V}}_{i},\mathit{\varphi})})\phantom{\rule{0ex}{0ex}}{U}_{\lambda}(\mathit{\varphi},\lambda )=-\sum _{i=I}^{n}{\delta}_{i}\frac{g({\mathit{V}}_{i},\mathit{\varphi})-1}{1+\lambda (g({\mathit{V}}_{i},\mathit{\varphi})-1)}$$

where ** ġ**(

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

$$\lambda (t|Z=1)/\lambda (t|Z=0)={\mathrm{e}}^{\gamma}.$$

(7)

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 ${\mathit{Z}}^{*}(t)={({Z}_{1}^{*}(t),\cdots ,{Z}_{p}^{*}(t))}^{T}$ are measured, then the Cox model (7) can be extended to λ(*t*|*Z* = 1, ** Z***(

Qin (1998) showed that (^{T}, )^{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 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 (^{T}, , )^{T}.

We restrict attention to the density ratio model (5) with weight function *g*(**υ, ϕ**) = e^{α+(υ,β)}, where **ϕ** = (α, **β**^{T})^{T} and (**υ, β**) is a polynomial function in **υ**, and the marginal Cox model adjusted for the covariates (*Z*, ** Z***

Next, define the random processes

$${\eta}_{n,\gamma}(t)={\mathrm{Z\u0304}}_{n}(t;\gamma )=\frac{\frac{1}{n}{\sum}_{k=1}^{n}{Z}_{k}I({X}_{k}\ge t){\mathrm{e}}^{\gamma {Z}_{k}}}{\frac{1}{n}{\sum}_{k=1}^{n}I({X}_{k}\ge t){\mathrm{e}}^{\gamma {Z}_{k}}}=\frac{{\xi}_{n,\gamma}(t)}{{\zeta}_{n,\gamma}(t)},t\in [0,\tau ],\phantom{\rule{0ex}{0ex}}{\eta}_{0,\gamma}(t)=\frac{E[ZI(X\ge t){\mathrm{e}}^{\gamma Z}]}{E[I(X\ge t){\mathrm{e}}^{\gamma Z}]}=\frac{{\xi}_{0,\gamma}(t)}{{\zeta}_{0,\gamma}(t)},t\in [0,\tau ],$$

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*) [0, τ] × {0, 1} × [0, 1]^{s} × {0, 1}, define the map $\mathit{\phi}={({\mathit{\phi}}_{1}^{T},{\phi}_{2},{\phi}_{3})}^{T}$ as

$$\mathit{\phi}(\mathit{\theta},{\eta}_{n,\mathit{\theta}})=\mathrm{\Delta}{\{{(-\frac{\lambda \mathit{\u0121}(\mathit{\upsilon},\mathit{\varphi})}{1+\lambda (g(\mathit{\upsilon},\mathit{\varphi})-1)}+z\frac{\mathit{\u0121}(\mathit{\upsilon},\mathit{\varphi})}{g(\mathit{\upsilon},\mathit{\varphi})})}^{T},-\frac{g(\mathit{\upsilon},\mathit{\varphi})-1}{1+\lambda (g(\mathit{\upsilon},\mathit{\varphi})-1)},z-{\eta}_{n,\mathit{\theta}}(x)\}}^{T}.$$

(8)

The dependence of on (*x*, Δ, **υ**, *z*) is suppressed in the notation. Let **Ψ** be the limit of **Ψ**_{n} as *n* → ∞ (and also *m* → ∞) and let **Ψ**(**θ**_{0}) = **0**. Let _{P} denote the zero-mean *P*–Brownian bridge process and define the class of functions *f*_{θ,t,r}(*x, z*) = *z ^{r}I*(

Theorem 1: *Let*
_{n}
*be the solution to*
**Ψ**_{n}(**θ**) = **0**
*and let*
**Ψ**(**θ**_{0}) = **0**. *Then*

$$\sqrt{n}({\mathit{\theta \u0302}}_{n}-{\mathit{\theta}}_{0})\underset{n\to \infty}{\overset{\mathcal{D}}{\to}}-{\mathbf{\Psi \u0307}}_{{\mathit{\theta}}_{0}}^{-1}\mathit{Z}$$

*where*

$$\mathit{Z}={\mathbb{G}}_{P}{\{{\mathit{\phi}}_{1}^{T}({\mathit{\theta}}_{0},{\eta}_{0,{\mathit{\theta}}_{0}}),{\phi}_{2}({\mathit{\theta}}_{0},{\eta}_{0,{\mathit{\theta}}_{0}}),{\phi}_{3}({\mathit{\theta}}_{0},{\eta}_{0,{\mathit{\theta}}_{0}})+{p}_{\delta}{l}_{{\mathit{\theta}}_{0}}\}}^{T}$$

*with the map*
$\mathit{\phi}={({\mathit{\phi}}_{1}^{T},{\phi}_{2},{\phi}_{3})}^{T}$
*defined in*
(8), *p*_{δ} = *P*(δ = 1), *and*

$${l}_{{\mathit{\theta}}_{0}}={l}_{{\mathit{\theta}}_{0}}(x,z)=\int {\zeta}_{0,{\mathit{\theta}}_{0}}^{-1}(t)({f}_{{\mathit{\theta}}_{0},t,1}(x,z)-{\eta}_{0,{\mathit{\theta}}_{0}}(t){f}_{{\mathit{\theta}}_{0},t,0}(x,z))\mathrm{d}F(t|\delta =1).$$

*Here F*(*t*|δ = 1) *is the conditional cdf of X given* δ = 1, *and*
_{θ0}
*is the continuously invertible derivative of the map*
**θ** **Ψ**(**θ**) *at*
**θ**_{0}
*and has matrix form*

$${\mathbf{\Psi \u0307}}_{{\mathit{\theta}}_{0}}=\left(\begin{array}{ccc}\hfill {\mathbf{\Psi \u0307}}_{11}\hfill & \hfill {\mathbf{\Psi \u0307}}_{12}\hfill & \hfill 0\hfill \\ \hfill {\mathbf{\Psi \u0307}}_{12}^{T}\hfill & \hfill {\mathrm{\Psi \u0307}}_{22}\hfill & \hfill 0\hfill \\ \hfill {0}^{T}\hfill & \hfill 0\hfill & \hfill {\mathrm{\Psi \u0307}}_{33}\hfill \end{array}\right)$$

*with entries*

$${\mathbf{\Psi \u0307}}_{11}=\int \mathrm{\Delta}\{-(\frac{{\rho}_{1}}{1+{\rho}_{1}(g(\mathit{\upsilon},{\mathit{\varphi}}_{0})-1)}-\frac{z}{g(\mathit{\upsilon},{\mathit{\varphi}}_{0})})\mathit{g\u0308}(\mathit{\upsilon},{\mathit{\varphi}}_{0})+(\frac{{\rho}_{1}^{2}}{{[1+{\rho}_{1}(g(\mathit{\upsilon},{\mathit{\varphi}}_{0})-1)]}^{2}}-\frac{z}{{g}^{2}(\mathit{\upsilon},{\mathit{\varphi}}_{0})})\mathit{\u0121}(\mathit{\upsilon},{\mathit{\varphi}}_{0}){\mathit{\u0121}}^{T}(\mathit{\upsilon},{\mathit{\varphi}}_{0})\}\mathrm{d}P(x,\mathrm{\Delta},\mathrm{\Delta}\mathit{\upsilon},z)\phantom{\rule{0ex}{0ex}}{\mathbf{\Psi \u0307}}_{12}=\int \frac{-\mathrm{\Delta}\mathbf{\u0121}(\mathit{\upsilon},{\mathit{\varphi}}_{0})}{{[1+{\rho}_{1}(g(\mathit{\upsilon},{\mathit{\varphi}}_{0})-1)]}^{2}}\mathrm{d}P(x,\mathrm{\Delta},\mathrm{\Delta}\mathit{\upsilon},z)\phantom{\rule{0ex}{0ex}}{\mathrm{\Psi \u0307}}_{22}=\int \frac{\mathrm{\Delta}{(g(\mathit{\upsilon},{\mathit{\varphi}}_{0})-1)}^{2}}{{[1+{\rho}_{1}(g(\mathit{\upsilon},{\mathit{\varphi}}_{0})-1)]}^{2}}\mathrm{d}P(x,\mathrm{\Delta},\mathrm{\Delta}\mathit{\upsilon},z)\phantom{\rule{0ex}{0ex}}{\mathrm{\Psi \u0307}}_{33}=\int -\mathrm{\Delta}{\eta}_{0},{\mathit{\theta}}_{0}(x)(1-{\eta}_{0,{\mathit{\theta}}_{0}}(x))\mathrm{d}P(x,\mathrm{\Delta},\mathrm{\Delta}\mathit{\upsilon},z)$$

*where*
(**υ, ϕ**) = ^{2}*g*(**υ, ϕ**)/**ϕ****ϕ**^{T}.

Denote the column vector of functions $\mathit{\phi \u0303}(\mathit{\theta},{\eta}_{0,\mathit{\theta}})={({\mathit{\phi}}_{1}^{T}(\mathit{\theta},{\eta}_{0,\mathit{\theta}}),{\phi}_{2}(\mathit{\theta},{\eta}_{0,\mathit{\theta}}),{\phi}_{3}(\mathit{\theta},{\eta}_{0,\mathit{\theta}})+{p}_{\delta}{l}_{\mathit{\theta}})}^{T}$. The following corollary describes the asymptotic variance of $\sqrt{n}({\mathit{\theta \u0302}}_{n}-{\mathit{\theta}}_{0})$.

Corollary 1: *The asymptotic random vector*
${\mathbf{\Psi \u0307}}_{{\mathit{\theta}}_{0}}^{-1}\mathit{Z}$
*in Theorem 1 is normally distributed with zero mean and covariance matrix*
$\mathbf{\Gamma}={\mathbf{\Psi \u0307}}_{{\mathit{\theta}}_{0}}^{-1}\mathbf{\Omega}{\mathbf{\Psi \u0307}}_{{\mathit{\theta}}_{0}}^{-1}$
*where*

$$\mathbf{\Omega}=P\mathit{\phi \u0303}({\mathit{\theta}}_{0},{\eta}_{0,{\mathit{\theta}}_{0}}){\mathit{\phi \u0303}}^{T}({\mathit{\theta}}_{0},{\eta}_{0,{\mathit{\theta}}_{0}})-P\mathit{\phi \u0303}({\mathit{\theta}}_{0},{\eta}_{0,{\mathit{\theta}}_{0}})P{\mathit{\phi \u0303}}^{T}({\mathit{\theta}}_{0},{\eta}_{0,{\mathit{\theta}}_{0}}).$$

Let _{n} denote the empirical estimator for **Γ** obtained by replacing *P* by the empirical probability measure _{n}, **θ**_{0} by _{n}, and η_{0,θ0} by η_{n,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^{α+(υ,β)+γ}.

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

$${H}_{0}^{0}:VE(\mathit{\upsilon})=0\text{for all}\mathit{\upsilon}\in {[0,1]}^{s},$$

(9)

which states that the vaccine provides no protection against infection with any HIV strain. If ${H}_{0}^{0}$ 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

$${H}_{0}:VE(\mathit{\upsilon})\equiv VE\text{for all}\mathit{\upsilon}\in {[0,1]}^{s}.$$

(10)

Under models (5) and (7), ${H}_{0}^{0}$ is equivalent to {${\mathrm{H\u0303}}_{0}^{0}$ : **β** = 0 and γ = 0}. The likelihood ratio test of ${\mathrm{H\u0303}}_{0}^{0}$ against {${\mathrm{H\u0303}}_{1}^{0}$ : **β** ≠ **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 ${\chi}_{s}^{2}$ and ${\chi}_{1}^{2}$ under ${\mathrm{H\u0303}}_{0}^{0}$, respectively. Simes’ procedure rejects ${H}_{0}^{0}$ if either max(*p*_{β}, *p*_{γ}) ≤ α or min(*p*_{β}, *p*_{γ}) ≤ α/2 where α is the nominal familywise level of significance. The Wald test of ${H}_{0}^{0}$ versus ${H}_{1}^{0}$ is based on the statistic $n({\mathit{\beta \u0302}}_{n}^{T},{\mathbf{\gamma \u0302}}_{n}){\mathit{\Gamma \u0302}}_{n,\mathit{\beta}\gamma}^{-1}{({\mathit{\beta \u0302}}_{n}^{T},{\mathrm{\gamma \u0302}}_{n})}^{T}$ where _{n,βγ} is the submatrix of _{n} pertaining to the components **β** and γ. Under ${H}_{0}^{0}$, the Wald test statistic is asymptotically ${\chi}_{s+1}^{2}$. We additionally propose a weighted one-sided Wald-type test of ${H}_{0}^{0}$ based on the *Z*-statistic

$$Z=(\sum _{i=1}^{s}\frac{{\mathrm{\beta \u0302}}_{n,i}}{\widehat{\text{var}}{\mathrm{\beta \u0302}}_{n,i}}-\frac{{\mathrm{\gamma \u0302}}_{n}}{\widehat{\text{var}}{\mathrm{\gamma \u0302}}_{n}})/\sqrt{\widehat{\text{var}}(\sum _{i=1}^{s}\frac{{\mathrm{\beta \u0302}}_{n,i}}{\widehat{\text{var}}{\mathrm{\beta \u0302}}_{n,i}}-\frac{{\mathrm{\gamma \u0302}}_{n}}{\widehat{\text{var}}{\mathrm{\gamma \u0302}}_{n}})},$$

(11)

which is designed to increase power to detect alternative hypotheses where both the marginal vaccine efficacy $VE=1-\frac{\lambda (t|Z=1)}{\lambda (t|Z=0)}=1-{\mathrm{e}}^{\gamma}>0$ and *V E*(**υ**) declines with all of the components of ** V**. The null hypothesis

By Proposition 1, the *T* ** V**|

$${D}_{z}=\underset{t,\mathit{\upsilon}}{\text{sup}}|{\mathrm{F\u0302}}_{T\mathit{V}}(t,\mathit{\upsilon}|z)-{\mathrm{F\u0302}}_{T}(t|z){\mathrm{F\u0302}}_{\mathit{V}}(\mathit{\upsilon}|z)|,$$

(12)

where * _{TV}* (

- Draw an independent sample (${X}_{i}^{*},{\delta}_{i}^{*}$) from the original time-on-study data (
*X*, δ_{i}_{i}) for group*Z*=*z*subjects, with replacement. - Draw a sample ${\mathit{V}}_{i}^{*}$ from the original mark data
of group**V**_{i}*Z*=*z*subjects who became infected (δ = 1), with replacement. - Compute
*D*based on the bootstrap data (${X}_{i}^{*},{\delta}_{i}^{*},{\delta}_{i}^{*}{\mathit{V}}_{i}^{*}$) for subjects_{z}*i*in group*Z*=*z*. - Estimate the critical value as the α–quantile of
*D*across the_{z}*B*replicates.

The bootstrap test has the nominal size asymptotically and is consistent against all alternatives (Romano, 1988). We test the overall hypothesis {*K*_{0} : *T* ** V** |

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^{α+βυ+γ}, υ [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 ${H}_{0}^{0}$ and *H*_{0} 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 ${\int}_{0}^{1}f(\upsilon |Z=0){\mathrm{e}}^{\alpha +\beta \upsilon}\mathrm{d}\upsilon =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 *N _{pI}* = 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

Figures 1, ,2,2, and Web Figure 2 characterize finite-sample bias of estimates and asymptotic standard error estimators of $\widehat{VE}(\upsilon )$ 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.

Finite-sample bias of estimation of *V E*(υ) = 1 − e^{α+βυ+γ} with a univariate mark. *N*_{pI} and *N*_{υI} denote the expected numbers of observed placebo and vaccine infections, respectively. This figure appears **...**

Asymptotic versus empirical (dotted) standard error estimates of $\widehat{VE}(\upsilon )$ with a univariate mark. *N*_{pI} 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 ${H}_{0}^{0}$ and *H*_{0}, 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 ${H}_{0}^{0}$, powers of the Wald and likelihood ratio test are comparable, and both tests outperform the ${\xdb}_{1}^{3}$ and ${\xdb}_{1}^{4}$ 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 ${H}_{0}^{0}$ 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%.

Powers of tests of the null hypothesis {${H}_{0}^{0}$ : V E (υ) = 0 for all υ [0,1]} against (i) the alternative hypothesis {${H}_{0}^{1}$ : non ${H}_{0}^{0}$}, and (ii) the alternative hypothesis {${H}_{2}^{0}$ : V E > and V E(υ) is a decreasing function **...**

Powers of two-sided tests of the null hypothesis {H_{0} : VE (υ) VE for all υ [0,1]}. The proposed Wald and likelihood ratio tests are compared with the alternative non- and semiparametric tests ${\xdb}_{2}^{np}$ and ${\xdb}_{}^{}$ **...**

Powers of the Wald and likelihood ratio tests of *H*_{0} are considerably higher (up to 36%) than those of GMS’s ${\xdb}_{2}^{np}$ and ${\xdb}_{2}^{sp}$ 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 {*K*_{0} : *T* *V*|*Z*}.

We additionally study simulation scenarios involving a bivariate mark ** V** = (

Assuming the model *V E*(υ) = 1 − e^{α+βυ+γ}, we examine robustness of the *V E*(υ) estimator and of the Wald and likelihood ratio tests of ${H}_{0}^{0}$ and *H*_{0} to violation of the *T* *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

$$f(\upsilon |T=t,Z=0;c)=\{ct{\mathrm{e}}^{-ct\upsilon}/(1-{\mathrm{e}}^{-ct})\}I(0\le \upsilon \le 1),c>0,$$

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 ${\int}_{0}^{1}f(\upsilon |T=t,Z=0;c){\mathrm{e}}^{\alpha (t)+\beta \upsilon}\mathrm{d}\upsilon =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* *V*|*Z*} using the diagnostic test (12) is evaluated.

Web Figure 3 depicts the Monte Carlo average deviation $\widehat{VE}(\upsilon )-\overline{VE}(\upsilon )$ where $\overline{VE}(\upsilon )={\int}_{0}^{3}VE(t,\upsilon )f(t)\mathrm{d}t/{\int}_{0}^{3}f(t)\mathrm{d}t$ 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* *V*|*Z* may frequently remain undetected. For scenarios with power to detect departures from *T* *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 ${H}_{0}^{0}$ and *H*_{0} retain correct size under violation of *T* *V*|*Z*. This robustness property holds for any violation of this assumption, because under either null hypothesis ${H}_{0}^{0}$ and *H*_{0}, α(*t*) α = 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) λ and λ(*t*|*Z* = 1) = λe^{γ0+γ1t}; thus γ_{1} governs the amount of variation in the marginal hazard ratio over time. We observed minimal differences in the patterns of average deviation $\widehat{VE}(\upsilon )-\overline{VE}(\upsilon )$, 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^{α+βυ+γ0+γ1t} does not include an interaction of *V* with *T*, and hence *V E*(*t*, υ) and *V E*(υ) have the same amount of change in υ.

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

$$d({\mathit{\text{seq}}}_{\text{subject}},{\mathit{\text{seq}}}_{\text{reference}})=\sum _{i=1}^{q}\frac{2s({a}_{i}^{s},{a}_{i}^{r})}{s({a}_{i}^{s},{a}_{i}^{s})+s({a}_{i}^{r},{a}_{i}^{r})},$$

where ${a}_{i}^{s}$ is the amino acid at site *i* in the subject’s sequence, ${a}_{i}^{r}$ 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 (*V*_{1}, *V*_{2}) measured from 110 infected subjects, with 49 and 52 unique mark values in *V*_{1} and *V*_{2}, respectively. Figure 3 shows the distributions of *V*_{1} and *V*_{2}.

Point and 95% pointwise confidence interval estimates of *V E*(υ) are plotted in Figures 4a–b. The weighted Wald test of ${H}_{0}^{0}$ : *V E*(υ) 0 yields *p*-values of 0.016 and 0.024 for the 92TH023 and A244 distances. The likelihood ratio test of *H*_{0} : *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).

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 $\widehat{VE}({\upsilon}_{1},{\upsilon}_{2})$. The weighted Wald test of ${H}_{0}^{0}$ yields *p* = 0.019. The likelihood ratio test of *H*_{0} yields *p* = 0.81. The *p*-value for the test of *T* (*V*_{1}, *V*_{2})|*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).

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*(**υ**) 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* ** V**|

$$\lambda (t,\mathit{\upsilon}|Z=z)=f(t,\mathit{\upsilon}|Z=z){S}^{-1}(t|Z=z),$$

with *f*(*t*, **υ**|*Z* = *z*) the conditional joint density function of (*T, V_{T}*)

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* *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.

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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |