# Assessing the effect of genetic markers on drug immunogenicity from a mechanistic model-based approach

Julianne Duhazé, Miguel Caubet, Signe Hässler, Delphine Bachelet, Matthieu Allez, Florian Deisenhammer, Anna Fogdell-Hahn, Aude Gleizes, Salima Hacein-Bey-Abina, Xavier Mariette, et al.

► To cite this version:

Julianne Duhazé, Miguel Caubet, Signe Hässler, Delphine Bachelet, Matthieu Allez, et al. "Assessing the effect of genetic markers on drug immunogenicity from a mechanistic model-based approach." BMC Medical Research Methodology, 2020, 20 (1), pp.69. 10.1186/s12874-020-00941-z. doi:10.1186/s12874-020-00941-z

HAL Id: hal-02538677

https://hal.sorbonne-universite.fr/hal-02538677

Submitted on 9 Apr 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Assessing the effect of genetic markers on drug immunogenicity from a mechanistic model-based approach

Julianne Duhazé<fnref n="1" />,<fnref n="2" />*, Miguel Caubet<fnref n="1" />, Signe Hässler³,⁴, Delphine Bachelet⁵, Matthieu Allez⁶, Florian Deisenhammer⁷, Anna Fogdell-Hahn⁸, Aude Gleizes⁹,¹⁰, Salima Hacein-Bey-Abina¹⁰,¹¹, Xavier Mariette¹², Marc Pallardy⁹, Philippe Broët<fnref n="1" />,<fnref n="2" />,¹³* and on behalf of the ABIRISK Consortium

**Abstract**

**Background:** With the growth in use of biotherapic drugs in various medical fields, the occurrence of anti-drug antibodies represents nowadays a serious issue. This immune response against a drug can be due either to pre-existing antibodies or to the novel production of antibodies from B-cell clones by a fraction of the exposed subjects. Identifying genetic markers associated with the immunogenicity of biotherapeutic drugs may provide new opportunities for risk stratification before the introduction of the drug. However, real-world investigations should take into account that the population under study is a mixture of pre-immune, immune-reactive and immune-tolerant subjects.

**Method:** In this work, we propose a novel test for assessing the effect of genetic markers on drug immunogenicity taking into account that the population under study is a mixed one. This test statistic is derived from a novel two-part semiparametric improper survival model which relies on immunological mechanistic considerations.

**Results:** Simulation results show the good behavior of the proposed statistic as compared to a two-part logrank test. In a study on drug immunogenicity, our results highlighted findings that would have been discarded when considering classical tests.

**Conclusion:** We propose a novel test that can be used for analyzing drug immunogenicity and is easy to implement with standard softwares. This test is also applicable for situations where one wants to test the equality of improper survival distributions of semi-continuous outcomes between two or more independent groups.

**Keywords:** Genetic, Drug immunogenicity, Semi-continuous data, Two-part improper survival model, Semi-parametric

<footnote>*Correspondence: julianne.joseph-duhaze@u-psud.fr;
philippe.broet@inserm.fr</footnote>

<footnote n="1">¹CHU Ste-Justine Research Center, H3T1C5 Montréal, Canada</footnote>

<footnote n="2">²University Paris-Saclay, CESP, INSERM, 94807 Villejuif, France
Full list of author information is available at the end of the article</footnote>
**Background**

Biopharmaceuticals products (BP) such as therapeutic monoclonal antibodies are nowadays a fast-growing class of drugs whose recent use in clinic has represented a critical step forward in the treatment of many severe auto-immune diseases. Nevertheless, for some patients these BP induce an activation of the immune system, leading to the formation of antibodies against the drug. The consequences range from transient appearance of anti-drug antibodies (ADA) without any clinical significance to severe loss of efficiency by either blocking the drug or enhancing the clearance [1].

The mechanisms leading to biotherapy immunogenicity can either be patient-related (e.g: genetic background, immunological status) or treatment-related (e.g: drug characteristics and formulations) but their relative contributions to the development of ADA is currently not fully understood and still remain to be deciphered. If major achievements for minimizing product-related factors involved in immunogenicity have been recently made, thanks to the remarkable progress in biopharmaceutical engineering, there is still an urgent need for identifying non-modifiable patient-related factors that may provide a basis for stratified or personalized therapeutic approaches. However, if an extensive research has been conducted to study the immunogenic potential of the biotherapies, less has been done to identify patients who are either at high or low risk for ADA development. In this search for patient-related predictive factors of immunogenicity, the genetic diversity in immune regulatory genes, is supposed to play a major role in the development of ADA [1, 2]. If early studies about drug immunogenicity assessment have mainly relied upon response-based endpoints, time-to-event analyses are more and more often recommended for taking into account the dynamic of ADA production. For such studies, subjects that have not been previously exposed to a particular BP are followed up for a certain period of time after the first BP administration. The main outcome is the first time of ADA detection after the initial drug administration and the objective is to identify factors that are related to these time-to-events [3, 4].

The motivation behind this work is that such time-to-event analysis is not straightforward as it should take into account that the population under study is usually a mixture of pre-immune, immune-reactive and immune-tolerant subjects. Here, the so-called pre-immune subjects are those with preexisting antibodies. These preexisting ADA have been observed mainly among patients with auto-immune diseases and can originate either from the innate immune system or from the adaptive immune responses to homologous ingredients [5]. The pre-immune status of a subject can be known by doing a screening test before the first administration of the drug. In contrast, the so-called immune-tolerant subjects

are those whose immune system is in a state of unre-sponsiveness to the exposure of the drug and thus will not produce ADA. Finally, the so-called immune-reactive subjects are those who have no pre-existing ADA and who are able to produce detectable levels of antibodies. Their time to ADA detection depends on the dynamic of the B cell clones production. However, for subjects who are not pre-immune and are censored during the study, we cannot determine their status as immune-reactive or immune-tolerant subjects.

Thus, for such analysis, we have to deal with a particular kind of semicontinuous data [6, 7] with non-null proportions of zero and of infinite values coming from the pre-immune and immune-tolerant subjects, respectively. As a result, we have to use specific models that belong to the class of two-part models [6] but need to be rethought to incorporate a defective (or improper) distribution. Such defective distributions have been mainly considered in oncology for analyzing early stage of cancer. In the literature, there are two main approaches, the oldest one assumes that the population under study is composed of two subpopulations of patients: those who will not experience the event of interest and those who are likely to experience the event of interest during the follow-up. This kind of two-component mixture model incorporates the defect fraction in a parametric or semi-parametric framework (for a review see Maller and Zhou [8]). An alternative approach relies upon defining the cumulative hazard as a bounded increasing positive function [9]. This modeling which has a meaningful biological interpretation in clinical oncology belongs to the class of promotion time cure models [10].

In this work, we propose to consider a novel two-part semiparametric improper survival model that can cope with semi-continuous time to event data and from which we can derive an efficient test statistic. Here, the time-to-event is zero for pre-immune subjects with a discrete probability mass, whereas its non-zero values (immune-reactive or immune-tolerant subjects) have an improper survival distribution. This latter improper survival distribution relies upon some biological understanding of the ADA production

The paper is organized as follows. We first present the proposed two-part survival model with its interpretation. Then, we derive a general score test for the null hypothesis of a same time-to-event distribution across the different genotypes. Next, we present the results of simulation experiments comparing the behavior of our proposed score test to a two-part logrank test. Then, the clinical relevance of using the proposed test is exemplified by the analysis of the association between genetic markers and the occurrence of ADA among a cohort of treatment-naïve patients suffering from auto-immune diseases and treated by biotherapies [11].
# Methods

## Modeling background

### Notation

Let $G$ denote the genotype for a biallelic marker. Here, being homozygous for the common variant is noted as $G=[AA]$ and corresponds to the reference group. Being heterozygous and homozygous for the alternative allele are noted as $G=[Aa]$ and $G=[aa]$, respectively.

As the underlying genetic model is unknown, two dummy variables $G_1$ and $G_2$ are created for investigating various genetic models, including: dominant, additive, recessive and overdominant. In practice, $G_1 = 0, 1, 2$ when $G=[AA]$, $[Aa]$, $[aa]$ respectively and $G_2 = 1$ when $G=[Aa]$ and 0 otherwise.

Thus, if we denote $\xi_1$ and $\xi_2$ the effects associated with $G_1$ and $G_2$ respectively (see Table 1), we have a dominant genetic model when $\xi_1 = \xi_2$, an additive genetic model when $\xi_1 \neq 0$ and $\xi_2 = 0$, a recessive genetic model when $\xi_2 = -\xi_1$, and an overdominant genetic model when $\xi_1 = 0$ and $\xi_2 \neq 0$.

In the following, we assume that the status for being a pre-immune subject can be determined before the first administration of the biotherapy. Let $Z$ be the pre-immune status variable such as $Z = 1$ if the subject is pre-immune and 0 otherwise.

Let denote $T$ the time-to-ADA detection and $C$ the censoring time. We assume that $T$ and $C$ satisfy the condition of independent and non informative censoring [12]. For each subject $i$ ($i = 1, ..., n$), $X_i = min(T_i, C_i)$ denotes the observed time of follow-up and $\delta_i = 1_{(X_i=T_i)}$ the indicator of ADA detection. We also denote $Y_i(t) = 1_{(t \le X_i)}$ the indicator of being at risk for the event at time $t$. Here, for $Z = 1$ then $T = 0$ and for $Z = 0$ then $T > 0$. For each patient $i$, the observed data consist of $(X_i, \delta_i, Z_i, G_{1i}, G_{2i})$.

### Two-part improper survival model

We define the survival distribution $S^*(t)$ for these semi-continuous data such as:

$$S^*(t; z) = \pi^z \times [1 - \pi)S(t)^{1-z}]$$

where $S(t)$ is the conditional survival function for the non-zero values of the time-to-event and $\pi$ is the probability for the zero values. Here, the survival function $S(t)$ is improper and its limiting value $S(\infty) = a$ with $0 < a < 1$ is called the tail defect and represents the probability of not experiencing the event of interest.

<table><caption>**Table 1** Effects associated to a genetic variant</caption><thead><tr><th></th><th>[AA]</th><th>[Aa]</th><th>[aa]</th></tr></thead><tbody><tr><td>G<sub>1</sub></td><td>0</td><td>1ξ<sub>1</sub></td><td>2ξ<sub>1</sub></td></tr><tr><td>G<sub>2</sub></td><td>0</td><td>1ξ<sub>2</sub></td><td>0</td></tr></tbody></table>

$G_1$ and $G_2$ are dummy variables representing a genetic bi-allelic variant. $[AA]$, $[Aa]$ and $[aa]$ are genotypes, with 'A' being the reference allele. $\xi_1$ and $\xi_2$ are the effects associated with $G_1$ and $G_2$ respectively.

In this work, we assume that $Z$ follows a logistic regression model that depends upon the genotype. For the sake of simplicity and without loss of generality, we assume a simple additive genetic model. Thus we have:

$$\pi(G_1, G_2) = \frac{e^{(\theta+\beta G_1)}}{1+e^{(\theta+\beta G_1)}}$$

where $\beta$ is the unknown regression coefficient and $\theta$ is the intercept.

For the reference group, $\pi(G_{[AA]}) = \pi_0 = \frac{e^\theta}{1+e^\theta}$. When $\beta = 0$, the proportions of pre-immune subjects are identical across the different genotypes.

For the conditional improper survival distribution $S(t)$ we introduce a new semiparametric model which is based on immunological mechanistic considerations and is presented in the next section.

### Conditional improper survival distribution

We propose to model the distribution of the time to ADA detection through a simplified mechanistic immunological model whereby each non pre-immune subject may or may not be able to produce ADA in response to the introduction of the biotherapy. From a biological perspective, it arises from the activation of unobservable BP-specific (T-dependent) B-cell clones that emerge and become immune-reactive ADA-producing clones. At the cellular level, each B-cell clone produces antibodies with a unique antigen-binding site with various levels of ADA affinities. Antibody affinity describes the strength of interaction between the antibody combining site and the relevant antigenic determinants in the drug. It represents an important qualitative parameter of the immune response. Positivity occurs as soon as any one of the B-cell clones is able to produce levels of ADA of sufficient amount and/or affinity for being detected by the assay.

Since B-cell clones are not directly observed, we assume a latent discrete probability distribution for the number of B-cell clones and a continuous distribution for their time-to-detection. From those latent distributions, we can deduce the marginal (or population) survival distribution of the time to ADA detection.

More precisely, let $K$ ($K = 0, ..., \infty$) denote a latent random variable that represents the number of (unobservable) B-cell clones. Let the random variable $T^j$, associated to the $j^{th}$ latent B-cell clones, be the time-to-detection with corresponding time-to-event survival function $A_j(t) = A_0(t)^{U_j}$. Here, $A_0(t)$ is a baseline decreasing function such as $1 \ge A_0(t) \ge 0$ and $U_j$ is a random variable that represents ADA affinity for the $j^{th}$ clone.

Here, we suppose that $K$ is distributed with probability mass function $\Phi$ and is supposed to be independent of $T^j$. We also supposed that the $U_j$ are independent and identically distributed with density function $\Psi$ and independent from $K$.
For patient *i* with *k* latent B-cell clone, let denote $T_i = \min_{0 \le j \le k} (T^j)$ the time-to-detection of the earliest B-cell clone. The conditional (patient-specific) survival function is expressed as:

$$S(t | K = k, (U_1 = u_1, \dots, U_k = u_k)) = \Pr(T_i > t) = \Pr(T^1 > t, \dots, T^k > t) = A_0(t)^{1_{\{k \ne 0\}}} \sum_{l=1}^{l=k} u_l = A_0(t)^{1_{\{k \ne 0\}}} w_k$$

where $w_k = \sum_{l=1}^{l=k} u_l$ is a realization of the random variable $W_k$ which is the randomly stopped sum of independent random variables [13] $U_l$ and $1_{\{k \ne 0\}} = 1$ when $k \ne 0$ and 0 otherwise. This can be interpreted as a first-activation scheme model [14].

In the following, we assume a Katz distribution for the distribution of the number of latent B-cell clones ($\Phi$) and a Gamma distribution $\Gamma(1, \tau)$ with shape parameter unity and scale parameter $\tau$ ($\tau \ge 0$) for the clonal ADA affinity ($\Psi$). As a consequence, the variable $W_k$ has a gamma distribution $\Gamma(k, \tau)$ (denoted in the following as $\Psi_{(k,\tau)}$). We also introduce $H_0(t)$ a positive non-decreasing function so that $A_0(t) = \exp\{-H_0(t)\}$.

Based upon these latter assumptions, the marginal (population) survival function is given by:

$$S(t) = \sum_{k=0}^{\infty} \left\{ \int_{0}^{\infty} A_0(t)^s \Psi_{(k, \tau)}(s) ds \right\} \Phi(k)$$

We recall that the Katz family [13] is the set of counting distributions for non-negative integers whose probability mass function satisfies the following first-order recurrence formula such as: $\Pr(x+1) = \Pr(x) \times \left(\frac{\omega+\gamma x}{1+x}\right)$; $x = 0, \dots, \infty_+$ where $\omega > 0$ and $\gamma < 1$. If $\omega + \gamma x < 0$, then $\Pr(x+j)$ is equal to zero for all $j > 0$. The mean and variance are $\mu = \frac{\omega}{1-\gamma}$ and $\sigma^2 = \frac{\omega}{(1-\gamma)^2}$. Moreover, it is worth noting that $\omega = \frac{\mu^2}{\sigma^2}$ and $\gamma$ is linked to the dispersion index such as $\frac{\sigma^2}{\mu} = (1-\gamma)^{-1}$.

The probability generating function can be written such as [13]:

$$g(s; \mu, \gamma) = \left[ (1 - \gamma)^{-1} (1 - \gamma s) \right]^{\frac{-\omega}{\gamma}} \text{ for } \gamma \neq 0, \text{ and} \\ g(s; \omega, \gamma) = \exp[-\omega(1-s)] \text{ for } \gamma = 0 \text{ with } |s| \le 1$$

This family of distribution covers a wide spectrum of discrete distributions that encompasses the binomial, the negative binomial and the Poisson distributions. If $\gamma = 0$, we have a Poisson distribution (equidispersion) with $S(t = \infty_+) = e^{-\omega}$. When $\gamma < 0$, we have underdispersion with $S(t = \infty_+) > e^{-\omega}$. When $\gamma \in [0, 1]$, we have overdispersion with $S(t = \infty_+) < e^{-\omega}$.

In the following, we assume that the distributions $\Phi$ and $\Psi$ depend upon the genotype. In practice, we assume that for the reference group [AA] the latent distribution $\Phi$ is the Poisson distribution ($\gamma = 0$) and $\Psi$ is the Gamma distribution $\Gamma(1, \tau)$. In the following, the alternative genetic variant, acting in a dominant, overdominant, additive or recessive way, is associated with changes either to the distribution of the number of clones (over/underdispersion)

or to the dispersion in ADA affinity (increase/decrease of $\tau$).

Thus, the marginal survival function for the non pre-immune subjects depends upon the genotype such as, for the reference group:

$$S(t; G_1 = 0 \cap G_2 = 0) = \exp\{-\omega[1 - (1 + H_0(t)\tau)^{-1}]\} \text{ with } \omega > 0. \text{ For the other genotypes, we have:}$$

$$S(t; G_1 \neq 0 \cup G_2 \neq 0)$$

$$= \left[ \frac{1 - (\gamma_1 G_1 + \gamma_2 G_2) (1 + H_0(t) \tau e^{(\alpha_1 G_1 + \alpha_2 G_2)})^{-1}}{1 - (\gamma_1 G_1 + \gamma_2 G_2)} \right]^{\frac{-\omega}{(\gamma_1 G_1 + \gamma_2 G_2)}}$$

where $\alpha_1, \alpha_2, \gamma_1, \gamma_2$ are the unknown regression coefficients. The parameters $\gamma_1, \gamma_2$ are linked to the distribution of the number of the B-cell clones and $\alpha_1, \alpha_2$ are linked to the dispersion in ADA affinity.

When $\gamma_1 = \gamma_2 = 0$ the improper survival distributions for each genotypes have a same defect such as $S(t = \infty_+; G_1, G_2) = e^{-\omega}$. The biological interpretation is that the alternative allele has no effect on the proportion of immune-tolerant. When $\gamma_1 \neq 0$ or $\gamma_2 \neq 0$, the improper survival distribution varies between genotypes. The biological interpretation is that the alternative allele is associated with changes in the distribution of the number of clones and consequently with changes in the proportion of immune-tolerant.

We write the improper survival model in terms of the hazard functions and take its first-order Taylor expansion of $\gamma_1 = \gamma_2 = 0$.

$$\lambda(t | G_1, G_2) = k_0(t) \times e^{\alpha G - 2 \log[1 + H_0(t) \tau e^{\alpha G}]} + \gamma G \left[1 + H_0(t) \tau e^{\alpha G}\right]^{-1} = k_0(t) e^{\Upsilon(t)}$$

where $\alpha G = \alpha_1 G_1 + \alpha_2 G_2$, $\gamma G = \gamma_1 G_1 + \gamma_2 G_2$, $k_0(t) = \omega \tau h_0(t)$ is a baseline hazard function with $h_0(t) = \frac{\partial H_0(t)}{\partial t}$ and $\Upsilon(t)$ is a function of $\gamma_1; \gamma_2; \alpha_1; \alpha_2; G_1; G_2$. This multiplicative hazard formulation will be used in the following.

## The proposed statistic

In this work, the general null hypothesis to be tested is defined by $H_0: \beta = \gamma_1 = \gamma_2 = \alpha_1 = \alpha_2 = 0$ and corresponds to a same survival distribution across genotypes against differences that can be associated to changes linked either to the pre-immune fraction ($\beta \neq 0$) or to the improper survival distribution for the non pre-immune subjects. From our mechanistic model, these differences are related either to the dispersion of the latent distribution for the B-cell clones ($\gamma_1 \neq 0$ or $\gamma_2 \neq 0$) or to the ADA affinity ($\alpha_1 \neq 0$ or $\alpha_2 \neq 0$).

In the following, we derive under the null hypothesis a score statistic for testing $H_0$.

The log-likelihood $LL_0(\beta, \gamma_1, \gamma_2, \alpha_1, \alpha_2; G_{1i}, G_{2i}, Z_i)$ derived under the working model can be expressed as the
sum of two terms $LL_0 = LL_1 + LL_2$ where the first term involves only the parameter $\beta$ whereas the second term involves the parameters $\gamma_1, \gamma_2, \alpha_1; \alpha_2$.

$$LL_1 = \sum_{i=1}^{n} Z_i \left\{ \theta + \beta G_{1i} - \log(1 + e^{(\theta + \beta G_{1i})}) \right\} + (1 - Z_i) \left\{ -\log(1 + e^{(\theta + \beta G_{1i})}) \right\}$$

$$LL_2 = \sum_{i=1}^{n} (1 - Z_i) \delta_i [\log(k_0(t)) + \Upsilon(t)] - (1 - Z_i) \int_{0}^{\infty} k_0(s) e^{\Upsilon(s)} ds$$

If we consider the non-parametric form for $k_0(t)$ that puts point mass $k_{0j}$ at the $j$-th ordered failure times $t_{(j)}$, then $K_0(t) = \sum_{j=1}^{J} k_{0j}(t_{(j)} \le t)$, as our improper survival model is a multiplicative model, the profile likelihood (LL2) (profile out $K_0$) leads to the classical log partial likelihood function from which a score test can be derived under the null hypothesis [15, 16]. Here the log-partial likelihood is: $LPL_2 = \sum_{i=1}^{n} (1 - Z_i) \delta_i \Upsilon(t_i) - \sum_{j=1}^{n} Y_j(t_i)(1 - Z_i)e^{\Upsilon(t_i)}$.

Thus, one can obtain two score statistics by separately computing the components of the score statistics derived from the first derivatives of $LL_1$ and the log partial likelihood $LPL_2$ with respect to $\beta$ and $\alpha_1, \alpha_2, \gamma_1, \gamma_2$ evaluated under $H_0$, respectively.

The first score statistic is as follows:

$$\hat{V}_{H_0, \beta} = \frac{\partial LL_1}{\partial \beta} = \sum_{i=1}^{n} G_{1i} \{Z_i - \pi_0\}$$

with $\pi_0 = \frac{e^\beta}{1+e^\beta}$ the probability of being a pre-immune under the null hypothesis $H_0$.

The components of the second score statistic are such as:

$$\hat{U}_{H_0, \alpha_1} = \frac{\partial LPL_2}{\partial \alpha_1} = \sum_{i=1}^{n} \delta_i W_1(t) \left\{ G_{1i} - \frac{\sum_{k=1}^{n} Y_k(t_i) G_{1k}}{\sum_{k=1}^{n} Y_k(t_i)} \right\}$$

$$\hat{U}_{H_0, \alpha_2} = \frac{\partial LPL_2}{\partial \alpha_2} = \sum_{i=1}^{n} \delta_i W_1(t) \left\{ G_{2i} - \frac{\sum_{k=1}^{n} Y_k(t_i) G_{2k}}{\sum_{k=1}^{n} Y_k(t_i)} \right\}$$

$$\hat{U}_{H_0, \gamma_1} = \frac{\partial LPL_2}{\partial \gamma_1} = \sum_{i=1}^{n} \delta_i W_2(t) \left\{ G_{1i} - \frac{\sum_{k=1}^{n} Y_k(t_i) G_{1k}}{\sum_{k=1}^{n} Y_k(t_i)} \right\}$$

$$\hat{U}_{H_0, \gamma_2} = \frac{\partial LPL_2}{\partial \gamma_2} = \sum_{i=1}^{n} \delta_i W_2(t) \left\{ G_{2i} - \frac{\sum_{k=1}^{n} Y_k(t_i) G_{2k}}{\sum_{k=1}^{n} Y_k(t_i)} \right\}$$

where $W_1(t) = 1 - 2(\frac{\Lambda(t)}{\omega})$ and $W_2(t) = 1 - (\frac{\Lambda(t)}{\omega})$ with $\Lambda(t)$ as the bounded cumulative hazard function under the null hypothesis $H_0$.

For computing these two score statistics, we should substitute $\Lambda(t), \omega$ and $\pi_0$ by efficient estimators computed under the null hypothesis. Here $\pi_0$ is the observed proportion of pre-immune subjects computed under the null hypothesis. The quantity $\hat{\Lambda}(t) = \sum_{i=1}^{n} \int_{0}^{t} \{\sum_{k=1}^{n} Y_k(s)\}^{-1} dN_i(s)$, where $N_i(t) = 1_{\{X_i \le t, \delta_i = 1\}}$ is the left-continuous version of the Nelson-Aalen estimator for the cumulative hazard [15, 17] obtained by using the pooled sample, and $\hat{\omega} = \hat{\Lambda}(t_{\max})$ is the maximum value of this estimator computed at the last observed failure time $t_{\max}$. Here, $\hat{\omega}$ is an efficient estimator of $\omega$ since we assume that the immune-reactive patients should experience the event within the one-year follow-up.

For the two scores, their corresponding information scalar and matrix, respectively $\hat{J}_{H_0}$ and $\hat{I}_{H_0}$, are obtained as presented in the Supplementary material. They are computed by using efficient estimators of $W_1(t)$ and $W_2(t)$ as given above. Then,

$$S_1 = \hat{V}_{H_0} \hat{J}_{H_0}^{-1} \hat{V}_{H_0} \text{ and}$$

$$S_2 = (\hat{U}_{H_0, \alpha_1}, \hat{U}_{H_0, \alpha_2}, \hat{U}_{H_0, \gamma_1}, \hat{U}_{H_0, \gamma_2}) \hat{I}_{H_0}^{-1} (\hat{U}_{H_0, \alpha_1}, \hat{U}_{H_0, \alpha_2}, \hat{U}_{H_0, \gamma_1}, \hat{U}_{H_0, \gamma_2})'$$

Thus, the statistic: $S_{H_0} = S_1 + S_2$ is asymptotically distributed under $H_0$ as a chi-square with five degrees of freedom. This test statistic is referred as TWIST, for Two-part Improper Survival Test in the following

## Statistical test comparison

We compared our proposed test to a Two-part LogRank Test, hereafter referred as TLRT. Here, the TLRT is derived in the same spirit as proposed by Lachenbruch (2001)[7] where the first part is the score test under the logistic model and the second part is a k-sample logrank test. Here, $k = 3$ for the three genotype groups.

## Results

### Data generation

Monte Carlo simulations were performed in order to evaluate the behavior of the TWIST. We reported the size of the proposed test, as well as its power properties with those obtained with a TLRT. Individual's genotype information was generated by summing the values of two Bernoulli variables independently drawn, with mean: 0.1, 0.2 and 0.3. This value corresponds to a pseudo-minor allele frequency (MAF), that is to say the proportion of [a] allele in the simulated population.

The status for being pre-immune was sampled from a Bernoulli variable, with mean of 0.05, 0.1 or 0.2. The values for $\beta$ are chosen so that the pre-immune fractions are equal across the different groups (then $\beta = 0$), or different according to the group ($\beta = \log(2), \log(\frac{1}{2})$).

For the non pre-immune and to evaluate the impact of departures from the reference latent distribution, we have considered scenarios with equi, under or overdispersed
latent distributions corresponding to Poisson, Binomial and Negative binomial distributions, respectively.

Data were generated such as, for the reference group with [AA] genotype, $S_{AA}(t) = \exp[-\omega(1 - (1 + t)^{-1})]$ and for the two other genotypes [aa] and [Aa],

$$S(t) = \left[ \frac{1 - (\gamma_1 G_1 + \gamma_2 G_2)(1 + t e^{\alpha_1 G_1 + \alpha_2 G_2})^{-1}}{1 - (\gamma_1 G_1 + \gamma_2 G_2)} \right]^{\frac{-\omega}{\gamma_1 G_1 + \gamma_2 G_2}}$$

The tail defect for the reference group was chosen such as: $S_{AA}(\infty_+) = e^{-\omega} = 0.50$. We also investigated the case where $S_{AA}(\infty_-) = e^{-\omega} = 0.30$.

The values for $\gamma_1$ and $\gamma_2$ were chosen accordingly to have over, equi or underdispersion, under an additive, dominant, recessive or overdominant genetic model. The values for $\alpha_1$ and $\alpha_2$ were chosen in such a way as to the risk $e^{\alpha_1 G_1 + \alpha_2 G_2} = 2$. Thus, when the genetic model is dominant, $\alpha_1 = 0.5 \times \log(2)$ and $\alpha_2 = 0.5 \times \log(2)$; when the genetic model is additive, $\alpha_1 = 0.5 \times \log(2)$ and $\alpha_2 = 0$; when the genetic model is recessive, $\alpha_1 = 0.5 \times \log(2)$ and $\alpha_2 = -0.5 \times \log(2)$ and when the genetic model is overdominant, $\alpha_1 = 0$ and $\alpha_2 = \log(2)$.

According to the model, $\alpha$ and $\gamma$ have the same directional effects when they are of opposite signs. Thereby, when $\gamma < 0$ and $e^{\alpha} > 1$ there is a global underdispersion, and the reverse for $\gamma \in [0 - 1]$ and $e^{\alpha} < 1$.

The censoring time $C$ was generated from a uniform distribution with parameter value computed from the chosen percentage of censoring. The percentage of censoring refers only to the percentage of censored observations without the tail defect fraction. We investigated no censoring and 20% censoring.

The total number of subjects was chosen to be 1,000 and for each configuration of parameters, 1,000 replications were performed. The levels and powers of the TWIST and the TLRT were estimated with a 0.05 significance level.

## Simulation results

The estimated level of the TWIST under its proper null hypothesis fell within the binomial range [0.036 - 0.064] (see the underlined values in Table 2b and in Tables S1 to S6 in the Supplementary material).

Just below, we comment the results of Table 2. It displays the estimated power gains in the TWIST and the TLRT for different values of the parameters $\alpha$, $\beta$ and $\gamma$, for a 10% proportion of pre-immune subjects and a minor allele frequency of 20%. In Table 2, both $\alpha$ and $\gamma$ have additive effects, hence, $\alpha_2 = 0$ and $\gamma_2 = 0$. The results for other genetic models are available in the Supplementary material.

We first comment the configurations where $\alpha$ has no effect ($e^{\alpha} = 1$, Table 2b). In this case, the TWIST's power is higher than the TLRT's when $\gamma_1 = -0.35$ (underdispersion). The TWIST's power is very close to the TLRT's when $\gamma_1 = 0$ or $\gamma_1 = 0.10$. In these cases, the difference

between the TWIST's power and the TLRT's vary from 0 to 0.02. When only $\beta$ has an effect, the TWIST has a very low power. For example, when $e^{\beta} = 0.5$, the TWIST's power is 0.04.

We will now comment the configurations where $\alpha$ has a non-null effect (Table 2a and c).

When parameters $\alpha$ and $\gamma$ have the same directional effect, the TWIST's power is always higher than the TLRT's, with a difference varying from 0.18 to 0.57. As an example, when both $\alpha$ and $\gamma$ increase the risk and $e^{\beta} = 2$, the TWIST's and TLRT's powers are 0.95 and 0.56 respectively.

When $\gamma$ has no effect, the TWIST's power is higher than the TLRT's, with a difference varying from 0.26 to 0.34.

When $\alpha$ and $\gamma$ have opposite directional effect, the TWIST's power values are close to the TLRT's. The difference between the powers of the two tests varies from -0.10 to 0.10. As an example, when $\gamma_1 = 0.10$, $e^{\alpha} > 1$ and $e^{\beta} = 2$, the TWIST's and TLRT's powers are 0.49 and 0.42 respectively.

As expected, when $e^{\beta} = 1$, the power of the TWIST is inferior to the one obtained with $e^{\beta} \neq 1$.

The results of simulations with $\alpha$ and $\gamma$ acting under a same genetic model other than additive are not displayed in this section but are joined in the Supplementary material (Table S6). The TWIST's power values obtained with $\alpha$ and $\gamma$ both dominant, overdominant or recessive are different from what we described above but follow the same global trends. Moreover, according to the genetic model, the power values are different. When both $\alpha$ and $\gamma$ are dominant or overdominant, the TWIST's power values are higher than the ones obtained under an additive genetic model. When both $\alpha$ and $\gamma$ are recessive, the TWIST's power values are lower than what is obtained under an additive genetic model.

We performed simulations with the same parameters as described above but with 20% censoring (see Table S1 in the Supplementary material). As expected, the TWIST's power values when there is 20% censoring are lower than when there is no censoring. Moreover, the difference between the TWIST's power and the TLRT's is lower with 20% censoring compared to no censoring.

Tables S2 to S6 in the Supplementary material display the results of simulations performed with all parameters equal to those of Table 2, unless otherwise specified. For all these configurations, the trends are similar to what we detailed above about Table 2, with changes in power values and gains. We describe briefly just below the results displayed in those additional tables.

With a minor allele frequency of 10%, as expected, the TWIST's power is lower than when the MAF is 20%, and with a 30% MAF, it is higher (see Tables S4 and S5).

With a proportion of pre-immune subjects in the reference group of 5% or 20%, the TWIST's power is similar
<table><caption>**Table 2** Power of the TWIST compared to the TLRT when α and γ are additive</caption><thead><tr><th rowspan="2">e<sup>β</sup></th><th colspan="2">Underdispersion<br>(γ<sub>1</sub> = -0.35 and γ<sub>2</sub> = 0)</th><th colspan="2">Equidispersion<br>(γ<sub>1</sub> = 0 and γ<sub>2</sub> = 0)</th><th colspan="2">Overdispersion<br>(γ<sub>1</sub> = 0.10 and γ<sub>2</sub> = 0)</th></tr><tr><th>TWIST</th><th>TLRT</th><th>TWIST</th><th>TLRT</th><th>TWIST</th><th>TLRT</th></tr></thead><tbody><tr><td></td><td colspan="6" style="text-align: left;">(a) α has an additive effect and e<sup>α</sup> > 1</td></tr><tr><td>2</td><td>0.95</td><td>0.56</td><td>0.42</td><td>0.08</td><td>0.49</td><td>0.42</td></tr><tr><td>1</td><td>0.93</td><td>0.36</td><td>0.29</td><td>0.02</td><td>0.37</td><td>0.28</td></tr><tr><td>0.50</td><td>0.95</td><td>0.47</td><td>0.32</td><td>0.03</td><td>0.42</td><td>0.32</td></tr><tr><td></td><td colspan="6" style="text-align: left;">(b) α has no effect, e<sup>α</sup> = 1</td></tr><tr><td>2</td><td>0.50</td><td>0.44</td><td>0.09</td><td>0.10</td><td>0.59</td><td>0.60</td></tr><tr><td>1</td><td>0.36</td><td>0.25</td><td><u>0.04</u></td><td><u>0.02</u></td><td>0.50</td><td>0.48</td></tr><tr><td>0.50</td><td>0.45</td><td>0.34</td><td>0.04</td><td>0.04</td><td>0.53</td><td>0.53</td></tr><tr><td></td><td colspan="6" style="text-align: left;">(c) α has an additive effect and e<sup>α</sup> < 1</td></tr><tr><td>2</td><td>0.22</td><td>0.32</td><td>0.41</td><td>0.12</td><td>0.95</td><td>0.77</td></tr><tr><td>1</td><td>0.11</td><td>0.14</td><td>0.30</td><td>0.04</td><td>0.91</td><td>0.65</td></tr><tr><td>0.50</td><td>0.14</td><td>0.18</td><td>0.35</td><td>0.05</td><td>0.96</td><td>0.74</td></tr></tbody></table>

(b) The underlined values show the estimated level of the type I error.

(c) The tail defects for the [AA], [Aa] and [aa] groups are respectively: 0.50, 0.55 and 0.59 for γ₁ = -0.35, and 0.50, 0.44 and 0.37 for γ₁ = 0.10. The proportion of pre-immune subjects is 10% and the MAF is 20%.

to what we obtain when this proportion is 10%, with slight variations according to the configuration of the parameters (see Additional file 1: Tables S4 and S5).

When the tail defect is 30%, the TWIST's power is higher than when the tail defect is 50% (see Table S6).

We also performed simulations with α having a different genetic model than γ, those results are displayed in Tables S7 and S8 of the Supplementary material and display similar trends than what is described above, with different power values according to the genetic model.

## Biotherapeutic immunogenicity study

We used our novel test to look for association between genetic variants and occurrence of ADA. The population study consists of 469 patients from the ABIRISK consortium cohort [11] suffering from auto-immune diseases and treated by biotherapies. These patients were naive for the biotherapies before the study and were followed during twelve months. The outcome was the time between the date of the first dose of biotherapy and the first detection of anti-drug antibodies (ADA). Patients without ADA occurrence were censored at the date of their last clinic visit. Among the 469 subjects, 17 (3.6%) had pre-existing ADA and 129 (27.5%) developed ADA during the one-year follow-up.

Exploratory analyses were performed on a list of 1,697 genes linked to immunity. We selected the genetic variants with a minor allele frequency greater than 0.15, for a total of 19,745 variants. At a 10% false discovery rate level [18], five genetic variants were significantly associated with ADA occurrence with the TWIST, and none with the TLRT. Among these five variants, we highlight

one from the fibroblast growth factor signaling pathway which test statistic is 28.96 (q-value = 0.093) with the TWIST and 21.87 with the TLRT (q-value = 0.34). Taken separately, the first and second components of the TWIST statistic have values of 8.02 and 20.94 respectively. It's worth noting that if we restricted the analyses to the 452 non pre-immune subjects with the same 10% false discovery rate, this genetic variant was not selected (q-value = 0.39).

This latter finding emphasizes the use of the TWIST and the inclusion of the pre-immune subjects for the analyses. Figure 1 shows the survival curves for the three genotypic groups of the highlighted variant. We found that the rare variant is associated with a higher occurrence of ADA in both pre-immune and non pre-immune subjects.

## Discussion

In this paper, we proposed a test for assessing the effect of genetic markers on drug immunogenicity. This test is based on a score statistic derived from a two-part semi-parametric improper survival model which allows to test the effect of genetic markers on the time-to-detection of ADA. It extends classical two-part models for semi-continuous data to the case where we can expect a mass at infinity. This model has a mechanistic interpretation of the ADA production by the B-cell clones. It allows different types of departures from equality of the two-part improper survival distributions. The biological interpretation of these departures are linked either to the proportion of pre-immune subjects or to the dynamic of ADA production. In this latter case, genetic variants can be
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Survival curve of the fibroblast growth factor signaling pathway variant. We named 'A' the reference allele and 'a' the alternative allele</figcaption></figure>

associated with changes in the number of clones or to their ADA affinities. As seen from our simulations, under the null hypothesis, our proposed test maintains a correct type I error in all the investigated configurations. Looking at the simulation results under the alternative hypotheses, we see that power gains of the proposed test are better than that of the two-part logrank test in most situations. In particular, we see that striking gains in power can be achieved by the TWIST as compared to the two-part logrank test when the parameters have the same directional effect, which is the case that we could expect the most in applications on real data. Since for complex diseases the underlying model is unknown, we have not restricted our test to one of the common genetic models. For various combinations of the classical genetic models, we obtain good performances with the proposed test. Moreover, it is worth noting that a straightforward extension of this test statistic can be done to investigate specific genetic models. Power gains decrease when the plateau or the percentage of censoring increase, which is not surprising since there are less events. Likewise, when the groups are too unbalanced, for example when the minor allele frequency is low or when the underlying genetic model is recessive, power gains decrease. The simulation results also show that our test performs well when there is a non negligible proportion of pre-immune subjects.

Notwithstanding the good behavior of the proposed test, some limitations should be mentioned. In this work, we have considered a proportional effect of the studied

variable on the time-to-detection associated with the B-cell clones. However, if this hypothesis is not true (e.g. in case of accelerated failure time), the power of the test will be reduced. Here, we have considered a Katz family for the distribution of the number of latent B-cell clones. This modelisation provides a flexible probabilistic modeling approach but at the expense of increasing the number of parameters to be tested. Moreover, it is worth keeping in mind that this test performs best in situations where a sufficient length of follow-up is provided. It means that we need a window of observation long enough for the last potential event among immune reactive to occur within it. In other words, we should provide a follow-up adequate for detecting the presence of immune-tolerant subjects in the study population. If this condition is unmet, the power of the test is reduced.

We used the proposed test for analyzing the relationship between some genetics loci and drug immunogenicity among a cohort of patient with auto-immune diseases and naive for the tested drugs. In this study, 3.6% of the patients were pre-immune subjects. We found that the alternative alleles of five variants on immunity genes are significantly associated with an increased risk of developing ADA. The result for one variant from the fibroblast growth factor signaling pathway is interesting since it has an effect on both the fraction of pre-immune subjects and the time of ADA occurrence. For this case, the inclusion of the pre-immune subjects in the analyses allows to detect the association of the variant with ADA occurrence.
Thus, we think that our test provides a new way of
assessing the effect of genetic markers on drug immuno-
genicity, while incorporating some biological understand-
ing of the ADA production. Our model-based approach
can be considered as a convenient representation of com-
plex patterns of genetic effects which allows detection of
departure from the null hypothesis model.

It is worth noting that since the proposed test statis-
tic gives equal weights to its two parts, its global power
can be reduced when one of them has a smaller effect
size (e.g. pre-immune subjects). Following the idea intro-
duced in the paper of Hu and Proschan [19], it could be
interesting to consider a linear combination of the two
components of our proposed test statistic. As in the work
of Hu and Proschan, we could assign a lesser weight to
the component with the smaller a priori effect than to
the other. However, it requires additional work for deriv-
ing the probability density function of the optimal linear
combination of the two components of the test statistic.
Moreover, further works should also be done for deriving
simple test statistics for testing separately each compo-
nent for the non pre-immune subjects.

In view of the gains in power obtained by the proposed
test, its use can be recommended in many clinical situa-
tions where unwanted drug immunogenicity occurs such
as in oncology and clinical immunology. It can also be
used to evaluate response to vaccination for populations
where pre-immune and non pre-immune patients exist.
Moreover, the proposed test can also be extended to take
other factors into account such as different treatments, by
developing a stratified version with strata defined by the
levels of the factors.

It can also be used to evaluate response to vaccination for populations where pre-immune and non pre-immune patients exist. Moreover, the proposed test can also be extended to take other factors into account such as different treatments, by developing a stratified version with strata defined by the levels of the factors.

Conclusion

In this study, we proposed a novel test statistic for assess-
ing the effect of genetic markers on drug immunogenicity
taking into account that the population under study is a
mixed one. This test statistic, which is easy to implement
with standard softwares, is also applicable in situations
where one wants to test the equality of improper survival
distributions of semi-continuous outcomes between two
or more independent groups.

**Supplementary information**

**Supplementary information** accompanies this paper at
https://doi.org/10.1186/s12874-020-00941-z.

Additional file 1: Supplementary material.

Abbreviations

ADA: Anti-drug antibodies; BP: Biopharmaceutical products; MAF: Minor allele frequency; TLRT: Two-part logrank test; TWIST: Two-part improper survival test;

Acknowledgments

This research work has been motivated by questions raised by studies conducted as part of the ABIRISK consortium (Anti-Biopharmaceutical Immunization: Prediction and analysis of clinical relevance to minimize the risk, Innovative Medicines Initiative). The methodological research was conducted as part of sustainability projects of the ABIRISK consortium. We thank the reviewers for their valuable and supportive comments.

Authors' contributions

JD and PB developed the original test statistic. PB coordinated the project and is JD's PhD thesis advisor. JD, MC and PB participated in writing the original draft. JD, SH, DB and PB analyzed the data. SH, DB, MA, FD, AFH, AG, SHBA, XM, MP and PB participated in the collection of the data. MP coordinated the ABIRISK project. All authors read and approved the final manuscript.

Funding

The research of the first author was supported by a grant from University Paris-Saclay (France). The funder had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

The data analyzed in this study were collected in the context of the ABIRISK project by ABIRISK partners. Access to the minimal dataset underlying the findings can be obtained by interested researchers upon request to the ABIRISK Sustainability Scientific Committee by submission of an analysis plan. The analysis plan should explain the purpose of the use of the data and confirm the intention to use the data only for replication studies concerning anti-drug inhibitors, since this is the limitation of the ethical permission on how this data can be used. The contact person of the ABIRISK Sustainability Scientific Committee to whom the requests should be sent is Marc Pallardy (marc.pallardy@inserm.fr).

Ethics approval and consent to participate

The ethics approvals for the ABI-MS-P01 study (EudraCT Number 2012-005450-30) were obtained from Medical Ethics Committee of the General University Hospital in Prague (reference 125/12, Evropský grant 1.LF UK-CAGEKID) for Czech Republic, Institutional committee of Heinrich Heine University, Düsseldorf, Germany (protocol reference 4451) for Düsseldorf, from Ethikkommission der Fakultät für Medizin der Technischen Universität München, München, Germany (reference 335/13) for München, from Ethikkommission Nordwest- und Zentralschweiz, Basel (reference 305/13) for Switzerland, from Ethikkommission der Medizinischen Universität Innsbruck, Innsbruck (reference AN2013-0040 331/2.1) for Austria, from Comité Ético de Investigación Clínica de l'Hospital Universitari Vall d'Hebrón, Barcelona (reference EPA(AG)66/2013(3866)) for Spain, from Stockholm Regional Ethics Committee, Stockholm (reference Dnr. 2013/1034-31/3 and Dnr. 2015/749-32) for Sweden. The ethics approvals for the ABI-RA-P01 study (ClinicalTrials.gov, NCT02116504) were obtained from Comité de Protection des Personnes Ile de France VII (reference 13-048) for France, from the Medical Ethical Committee of the Academisch Medisch Centrum, Amsterdam (reference 2013-304#B20131074) for the Netherlands, from the Local Ethics Committee of AOU Careggi (reference Protocol N° 2012/0035982) for Italy, from NRES Committee London, City & East (reference 14/LO/0506) for the United Kingdom. The ethics approvals for the ABI-IBD-P01 study were obtained from Comité de Protection des Personnes Ile de France IV (reference 2013/24) for France, from Comité d'éthique hospitalo-facultaire universitaire de Liège (reference 2015/55) for Belgium. All participants gave their written informed consent.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

<footnote n="1">¹CHU Ste-Justine Research Center, H3T1C5 Montréal, Canada. ²University Paris-Saclay, CESP, INSERM, 94807 Villejuif, France.</footnote>
³Immunology-Immunopathology-Immunotherapy (i3) Laboratory, UMR-S 959, INSERM, Sorbonne Université and Biotherapy Unit, Inflammation-Immunopathology-Biotherapy Department, Hôpital Pitié-Salpetrière, AP-HP, 75000 Paris, France. ⁴CESP, INSERM UMR 1018, Faculty of Medicine, Paris-Sud University, UVSQ, Paris-Saclay University, 94807 Villejuif, France. ⁵Inserm CIC-EC 1425, Centre d'Investigation Clinique and Département d'Epidémiologie Biostatistique et Recherche Clinique, Hôpital Bichat, AP-HP, 75000 Paris, France. ⁶Department of Gastroenterology, Hôpital Saint-Louis, AP-HP, Université Paris-Diderot, 75000 Paris, France. ⁷Innsbruck Medical, A-6010 - A-6080 Innsbruck, Austria. ⁸Department of Clinical Neuroscience, Karolinska Institutet, SE-171 76 Stockholm, Sweden. ⁹INSERM UMR996, Faculté Pharmacie, Université Paris Sud, 92290 Châtenay-Malabry, France. ¹⁰Clinical Immunology Laboratory, AP-HP, Le Kremlin-Bicêtre Hospital, Paris-Sud University Hospitals, 94270 Le Kremlin-Bicêtre, France. ¹¹UTCBS, CNRS UMR 8258, INSERM U1022, Faculty of Pharmacy, Paris-Descartes-Sorbonne-Cité University, 75000 Paris, France. ¹²Centre for Immunology of Viral Infections and Autoimmune Diseases, INSERM UMR1184, AP-HP, Université Paris-Sud, Hôpitaux Universitaires Paris-Sud, 94270 Le Kremlin-Bicêtre, France. ¹³Abirisk consortium WP4, 94807 Villejuif, France.

Received: 19 June 2019 Accepted: 25 February 2020

Published online: 20 March 2020

## References

1. Van de Weert M, Horn Moller E. Immunogenicity of biopharmaceuticals. New York: Springer; 2008.

2. Pineda C, Castañeda Hernández G, Jacobs I, Alvarez D, Carini C. Assessing the immunogenicity of biopharmaceuticals. BioDrugs. 2016;30(3):195-206.

3. Shankar G, Arkin S, Cocea L, Devanarayan V, Kirshner S, Kromminga A, et al. Assessment and reporting of the clinical immunogenicity of therapeutic proteins and peptides-harmonized terminology and tactical recommendations. AAPS J. 2014;16(4):658-73.

4. Rup B, Pallardy M, Sikkema D, Albert T, Allez M, Broet P, et al. Standardizing terms, definitions and concepts for describing and interpreting unwanted immunogenicity of biopharmaceuticals: recommendations of the Innovative Medicines Initiative ABIRISK consortium. Clin Exp Immunol. 2015;181(3):385-400.

5. Xue L, Rup B. Evaluation of pre-existing antibody presence as a risk factor for posttreatment pnti-drug antibody induction: analysis of human clinical study data for multiple biotherapeutics. AAPS J. 2013;15(3):893-6.

6. Farewell V, Long D, Tom B, Yiu S, Su L. Two-Part and related regression models for longitudinal data. Annu Rev Stat Appl. 2017;4:283-315.

7. Lachenbruch P. Comparisons of two-part models with competitors. Stat Med. 2001;20(8):1215-34.

8. Maller R, Zhou X. Survival analysis with long-term survivors. New York: Wiley; 1997.

9. Yakovlev A, Tsodikov A. Stochastic models of tumor latency and their biostatistical applications. Singapore: World Scientific; 1996.

10. Tsodikov A, Ibrahim J, Yakovlev A. Estimating cure rates from survival data: an alternative to two-component mixture models. J Am Stat Assoc. 2003;764(98):1063-78.

11. ABIRISK: Anti-Biopharmaceutical Immunization: prediction and analysis of clinical relevance to minimize the RISK. http://www.abirisk.eu/. Accessed 14 May 2019.

12. Fleming T, Harrington D. Counting processes and survival analysis. New York: Wiley; 1991.

13. Johnson N, Kemp A, Kotz S. Continuous univariate distributions. New York: Wiley; 2005.

14. Cooner F, Banerjee S, Carlin BP, Sinha D. Flexible cure rate modeling under latent activation schemes. J Am Stat Assoc. 2007;102(478):560-72.

15. Andersen P, Borgan O, Gill R, Keiding N. Statistical models based on counting processes. New York: Wiley; 1993.

16. Cox D. Partial likelihood. Biometrika. 2001;62(2):269-76.

17. Nelson W. Theory and applications of hazard plotting for censored failure data. Technometrics. 1972;14:945-65.

18. Dalmasso C, Broet P, Moreau T. A simple procedure for estimating the false discovery rate. Bioinformatics. 2005;21:660-8.

19. Hu Z, Proschan M. Two-part test of vaccine effect: Z. HU AND M. PROSCHAN. Stat Med. 2015;34(11):1904-11. http://doi.wiley.com/10.1002/sim.6412.

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_4.png" /></figure>