Efficient screening of predictive biomarkers for individual treatment selection

TL;DR

This paper introduces an efficient statistical method for identifying predictive biomarkers in large-scale genetic studies, improving detection power and false discovery control for personalized treatment strategies.

Contribution

It proposes a novel approach combining weighted loss functions, synthetic posterior inference, and empirical Bayes estimation to enhance biomarker screening in complex clinical data.

Findings

Detected more significant biomarkers than standard methods in breast cancer data

Demonstrated higher power and false discovery rate control in simulations

Applicable to large-scale genetic and molecular studies

Abstract

The development of molecular diagnostic tools to achieve individualized medicine requires identifying predictive biomarkers associated with subgroups of individuals who might receive beneficial or harmful effects from different available treatments. However, due to the large number of candidate biomarkers in the large-scale genetic and molecular studies, and complex relationships among clinical outcome, biomarkers and treatments, the ordinary statistical tests for the interactions between treatments and covariates have difficulties from their limited statistical powers. In this paper, we propose an efficient method for detecting predictive biomarkers. We employ weighted loss functions of Chen et al. (2017) to directly estimate individual treatment scores and propose synthetic posterior inference for effect sizes of biomarkers. We develop an empirical Bayes approach, namely, we estimate unknown hyperparameters in the prior distribution based on data. We then provide efficient screening methods for the candidate biomarkers via optimal discovery procedure with adequate control of false discovery rate. The proposed method is demonstrated in simulation studies and an application to a breast cancer clinical study in which the proposed method was shown to detect the much larger numbers of significant biomarkers than existing standard methods.

Tables

Table 1. Table 1: Average numbers of true positives and false positives for the proposed methods with difference number of knots ( L 𝐿 L ) and alternative q 𝑞 q -value based methods, based on 200 simulations.

	# true positive					# false positive
FDR levels	5%	10%	15%	20%		5%	10%	15%	20%
Binary ($\pi =0.8$)
ODP-P ($L=100$)	0.83	4.91	12.59	23.44		0.13	1.06	3.76	9.56
ODP-N ($L=50$)	0.33	2.24	7.57	17.00		0.04	0.42	1.89	6.13
ODP-N ($L=100$)	0.33	2.25	7.58	17.07		0.04	0.43	1.90	6.15
ODP-N ($L=150$)	0.33	2.25	7.59	17.09		0.04	0.43	1.90	6.15
ODP-N ($L=200$)	0.33	2.25	7.59	17.10		0.04	0.43	1.90	6.15
$T_k$	0.40	1.20	2.33	3.91		0.07	0.29	0.71	1.28
$S_k$	0.35	1.11	2.08	3.44		0.07	0.27	0.56	1.05
Binary ($\pi =0.5$)
ODP-P ($L=100$)	0.36	4.80	18.47	45.64		0.07	0.82	4.25	12.67
ODP-N ($L=50$)	0.22	2.14	12.45	39.36		0.03	0.34	2.77	10.58
ODP-N ($L=100$)	0.22	2.14	12.51	39.53		0.03	0.34	2.79	10.63
ODP-N ($L=150$)	0.22	2.14	12.53	39.57		0.03	0.34	2.79	10.65
ODP-N ($L=200$)	0.22	2.14	12.54	39.59		0.03	0.34	2.79	10.66
$T_k$	0.30	0.71	1.49	2.85		0.04	0.13	0.32	0.63
$S_k$	0.29	0.71	1.31	2.33		0.03	0.13	0.26	0.51
Survival ($\pi =0.8$)
ODP-P ($L=100$)	2.86	11.45	24.35	40.81		0.34	2.20	6.68	15.90
ODP-N ($L=50$)	2.22	9.53	21.87	39.03		0.24	1.91	6.43	15.54
ODP-N ($L=100$)	2.23	9.54	21.93	39.12		0.24	1.92	6.44	15.60
ODP-N ($L=150$)	2.23	9.55	21.97	39.15		0.24	1.92	6.44	15.61
ODP-N ($L=200$)	2.23	9.55	21.98	39.17		0.24	1.92	6.44	15.62
$T_k$	1.35	3.35	6.39	10.51		0.17	0.67	1.64	3.56
$S_k$	1.16	2.80	5.66	9.74		0.15	0.50	1.29	2.93
Survival ($\pi =0.5$)
ODP-P ($L=100$)	1.34	11.84	36.93	79.33		0.13	1.92	7.91	22.08
ODP-N ($L=50$)	0.95	8.43	30.75	76.12		0.11	1.43	7.21	21.92
ODP-N ($L=100$)	0.96	8.47	30.85	76.31		0.11	1.44	7.23	22.02
ODP-N ($L=150$)	0.96	8.49	30.88	76.38		0.11	1.45	7.24	22.03
ODP-N ($L=200$)	0.96	8.50	30.90	76.40		0.11	1.45	7.24	22.04
$T_k$	0.64	1.67	3.12	6.04		0.08	0.24	0.63	1.30
$S_k$	0.57	1.46	2.92	5.07		0.06	0.24	0.54	1.10

Table 2. Table 2: The numbers of detected genes in breast cancer data.

FDR levels	5%	10%	15%	20%
ODP-P	11	45	104	197
ODP-N	10	40	81	140
$T_k$	1	2	2	2
$S_k$	0	0	0	0

Equations

$$L_n(f)=\sum _{i=1}^n\frac{M\{Y_i,T_{i}f({\text{boldmathX}}_i,{\text{boldmathZ}}_i)\}}{T_i\pi ({\text{boldmathX}}_i,{\text{boldmathZ}}_i)+(1-T_i)/2},$$

$$\sum _{i=1}^n\frac{M(Y_i,T_i{\alpha }_k+T_i{X}_{ik}{\beta }_k+T_i{\text{boldmathZ}}_i^t{\text{boldmathω}}_k)}{T_i\pi ({\text{boldmathX}}_i,{\text{boldmathZ}}_i)+(1-T_i)/2}.$$

$$L_k({\alpha }_k,{\beta }_k,{\text{boldmathω}}_k)=\exp \left[-\frac{1}{A}\sum _{i=1}^n\frac{M(Y_i,T_i{\alpha }_k+T_i{X}_{ik}{\beta }_k+T_i{\text{boldmathZ}}_i^t{\text{boldmathω}}_k)}{T_i\pi ({\text{boldmathX}}_i,{\text{boldmathZ}}_i)+(1-T_i)/2}\right],$$

$$G({\beta }_k)=\pi g_0({\beta }_k)+(1-\pi )g_1({\beta }_k),$$

$$\begin{array}{rl}G({\beta }_k|\mathrm{Data})& =\frac{{\text{PL}}_k({\beta }_k)G({\beta }_k;\pi ,g_1)}{\int {\text{PL}}_k({\beta }_k)G({\beta }_k;\pi ,g_1)\mathrm{d}{\beta }_k}\\ & =\frac{\pi {\text{PL}}_k({\beta }_k)g_0({\beta }_k)+(1-\pi ){\text{PL}}_k({\beta }_k)g_1({\beta }_k)}{\pi {\text{PL}}_k(0)+(1-\pi )\int {\text{PL}}_k({\beta }_k)g_1({\beta }_k)\mathrm{d}{\beta }_k},\end{array}$$

$$\begin{array}{rl}\prod _{k=1}^p& \int {\text{PL}}_k({\beta }_k)G({\beta }_k;\pi ,g_1)\mathrm{d}{\beta }_k\\ & =\prod _{k=1}^p\left\{\pi {\text{PL}}_k(0)+(1-\pi )\int {\text{PL}}_k({\beta }_k)g_1({\beta }_k)\mathrm{d}{\beta }_k\right\}.\end{array}$$

$$P({\gamma }_k=1|\text{Data})=\frac{(1-\pi )\int {\text{PL}}_k({\beta }_k)g_1({\beta }_k)\mathrm{d}{\beta }_k}{\pi {\text{PL}}_k(0)+(1-\pi )\int {\text{PL}}_k({\beta }_k)g_1({\beta }_k)\mathrm{d}{\beta }_k},$$

$${\text{ODS}}_k=\frac{\int {\text{PL}}_k({\beta }_k){\hat{g}}_1({\beta }_k)\mathrm{d}{\beta }_k}{{\text{PL}}_k(0)},$$

$$\text{FDR}(\lambda )=\frac{1}{|\Theta (\lambda )|}\sum _{k\in \Theta (\lambda )}\hat{P}({\gamma }_k=0|\text{Data}),$$

$$Y=I\left(\sum _{k=1}^p{\gamma }_k{X}_k+\sum _{k=1}^p{\delta }_k{X}_k^2+T\sum _{k=1}^p{\beta }_k{X}_k+{\xi }_{1}T{Z}_1+{\xi }_{2}T{Z}_2+\epsilon >0\right),$$

$$f(x)=\pi {\delta }_0(x)+(1-\pi )0.3\varphi (x;0.2,(0.1{)}^2)+(1-\pi )0.7\varphi (x;-0.5,(0.1{)}^2),$$

$$T_k=\frac{{\hat{\beta }}_k}{s_k},S_k=\frac{{\hat{\beta }}_{k1}-{\hat{\beta }}_{k0}}{\sqrt{s_{k1}^2+s_{k0}^2}},$$

$$\tilde{Y}=\exp \left(\sum _{k=1}^p{\gamma }_k{X}_k+\sum _{k=1}^p{\delta }_k{X}_k^2+T\sum _{k=1}^p{\beta }_k{X}_k+{\xi }_{1}T{Z}_1+{\xi }_{2}T{Z}_2+\epsilon \right),$$

$$P(T=1|X)=\text{logit}^{-1}\Big{(}\gamma_{0}+\sum_{k=1}^{p}\gamma_{k}X_{k}\Big{)}.$$

$$\sum_{k=1}^{p}\left[z_{k}\Big{\{}\text{PL}_{k}(0)+\log\pi\Big{\}}+(1-z_{k})\sum_{\ell=1}^{L}I(w_{k}=\ell)\Big{\{}\text{PL}_{k}(a_{\ell})+\log p_{\ell}\Big{\}}\right].$$

$$Q(\varphi |{\varphi }^{(t)})$$

$${\xi }_k^{(t)}=\frac{{\pi }^{(t)}{\text{PL}}_k(0)}{{\pi }^{(t)}{\text{PL}}_k(0)+(1-{\pi }^{(t)}){\sum }_{\ell =1}^L{p}_{\ell }^{(t)}{\text{PL}}_k(a_{\ell })}\mathrm{and}{\eta }_{k\ell }^{(t)}=\frac{p_{\ell }^{(t)}{\text{PL}}_k(a_{\ell })}{{\sum }_{j=1}^L{p}_j^{(t)}{\text{PL}}_k(a_j)}.$$

$${\pi }^{(t+1)}=\frac{1}{p}\sum _{k=1}^p{\xi }_k^{(t)}\mathrm{and}{p}_{\ell }^{(t+1)}=\frac{1}{p}\sum _{k=1}^p{\eta }_{k\ell }^{(t)},\ell =1,\dots ,L.$$

Full text

**Efficient Screening of Predictive Biomarkers for Individual Treatment Selection

**

Shonosuke Sugasawa1,2∗, Hisashi Noma2,3

1Center for Spatial Information Science, The University of Tokyo, Chiba, Japan

2Research Center for Medical and Health Data Science, The Institute of Statistical Mathematics, Tokyo, Japan

3Department of Data Science, The Institute of Statistical Mathematics, Tokyo, Japan

*Correspondence: Shonosuke Sugasawa

         Center for Spatial Information Science

         The University of Tokyo

         5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8568, Japan

         E-mail: sugasawa@csis.u-tokyo.ac.jp

Summary

The development of molecular diagnostic tools to achieve individualized medicine requires identifying predictive biomarkers associated with subgroups of individuals who might receive beneficial or harmful effects from different available treatments. However, due to the large number of candidate biomarkers in the large-scale genetic and molecular studies, and complex relationships among clinical outcome, biomarkers and treatments, the ordinary statistical tests for the interactions between treatments and covariates have difficulties from their limited statistical powers. In this paper, we propose an efficient method for detecting predictive biomarkers. We employ weighted loss functions of Chen et al. (2017) to directly estimate individual treatment scores and propose synthetic posterior inference for effect sizes of biomarkers. We develop an empirical Bayes approach, namely, we estimate unknown hyperparameters in the prior distribution based on data. We then provide efficient screening methods for the candidate biomarkers via optimal discovery procedure with adequate control of false discovery rate. The proposed method is demonstrated in simulation studies and an application to a breast cancer clinical study in which the proposed method was shown to detect the much larger numbers of significant biomarkers than existing standard methods.

Key words: Empirical Bayes; False discovery rate; Optimal discovering procedure; Propensity score

Introduction

Due to the advances of disease biology, it has been revealed that there is substantial molecular heterogeneity among individual patients in many diseases. This implies that the benefits and harms of many treatments might also be heterogeneous, and accurate molecular diagnostic methods could maximize the treatment benefits for individual patients (Gabriel and Normand, 2012). To develop the molecular diagnostic tools, a key task is the identification of predictive biomarkers associated with subgroups of individuals who might receive beneficial or harmful effects from different available treatments (Matsui et al., 2015; Lipkovich et al., 2017). It is typically explained by interactions between the treatment and candidate biomarkers, but the conventional interaction tests have substantial difficulties for these analyses because of their serious limitations of statistical powers (Matsui et al., 2018).

On the other hand, efficient estimating methods of individual treatment effects (ITE) have been also developed for large-scale genetic and molecular studies. Specifically, there is a growing number of literatures regarding the efficient estimation of ITE (e.g. Chen et al., 2017; Kehl and Ulm, 2006; Lu et al., 2018; Sugasawa and Noma, 2019; Tian et al., 2014; Wager and Athey, 2018; Zhang et al., 2017) among many others. Although these methods can effectively estimate ITE, the estimated model is typically too complicated, and it would be difficult to understand which biomarkers are actually associate with ITE. Besides, the regularization methods often provide us not only the estimates of model parameters but also selection results of biomarkers included in the estimated model (Lu et al., 2013; Zhang and Zhang, 2018), but the regularization method do not necessarily guarantee the control of the degree of false discoveries, i.e. false discovery rate (FDR; Benjamini and Hochberg, 1995). Since the primary purpose of large-scale genetic and molecular studies is screening relevant candidate markers, the prioritization and assessment of the accuracy of selections are relevant tasks.

In this article, we propose an efficient method for screening predictive biomarkers associated with ITE under control of FDR. We employ the weighted loss function proposed by Chen et al. (2017) designed to directly estimate ITE score without explicitly specifying the main effect, and construct synthetic likelihood for effect sizes of candidate biomarkers. The synthetic likelihood is then combined with a latent semiparametric distribution of true effect sizes, which enables us derive the synthetic posterior distribution for effect sizes. The underlying semiparametric distribution of effect sizes can be estimated by the EM algorithm (Dempster et al., 1977), and we propose the optimal discovery procedure (ODP) to detect predictive biomarkers with adequate control of FDR. Although most of multiple testing procedures with FDR control are based on summary statistics such as $z$-values or $p$-values for effect sizes of biomarkers (e.g. Storey and Tibshirani, 2003; Sun and Cai, 2007), the proposed method more directly use data via synthetic likelihood. Also the proposed method would be able to increase statistical power to detect relevant biomarkers by borrowing information from related biomarkers via an empirical Bayes approach. In fact, through simulation studies, we numerically show that the proposed method can detect much more true predictive biomarkers than standard methods. Moreover, we apply the proposed method to an observational study of breast cancer, and we found that the proposed methods were able to detect a large number of predictive biomarkers whereas the standard methods detected less than 2 even when FDR is allowed to be $20\%$.

This paper is organized as follows. In Section 2, we describe the proposed method and demonstrate the optimal discovery procedure. In Section 3, we conducted simulation studies to confirm the effectiveness of the proposed method compared with conventional testing approaches. We then apply the proposed method to a dataset of a breast cancer study in Section 4. We give some discussions in Section 5. Finally, R-code implementing the proposed method is available at GitHub (https://github.com/sshonosuke/SB-ITS).

Proposed Method

Notations and assumptions

We adopt the notation based on the potential outcome framework in causal inference (Rubin, 2005). Let $T\in\{-1,1\}$ be the treatment indicator with $T=1$ indicating a patient being treated and $T=-1$ indicating the opposite. Let $Y^{(T)}$ denote the potential outcome indicating the response of a patient (e.g., survival time, binary disease status) if the patient receives treatment $T$. In practice, only one of the potential outcomes $Y^{(1)}$ and $Y^{(-1)}$ can be observed for each patient, that is, $Y=I(T=1)Y^{(1)}+I(T=-1)Y^{(-1)}$, where $I(\cdot)$ is the indicator function. We employ strongly ignorable assumption (Rosenbaum and Rubin, 1983; Rubin, 2005), that is, $T$ is independent of $(Y^{(1)},Y^{(-1)})$ given the $p$-dimensional covariates {\text{\boldmathX}}=(X_{1},\ldots,X_{p})^{t} of potential predictive biomarkers of the treatment (e.g., genotypes, gene expressions). Moreover, let $Z$ be a vector of potentially confounding covariates with $X$. For the treatment assignment, we assume that probability of treatment assignment is a function of $X$ and $Z$, that is, P(T=1|{\text{\boldmathX}},{\text{\boldmathZ}})=\pi({\text{\boldmathX}},{\text{\boldmathZ}}), where \pi({\text{\boldmathX}},{\text{\boldmathZ}}) is known as propensity score, and \pi({\text{\boldmathX}},{\text{\boldmathZ}})=1/2 under randomized clinical trial or \pi({\text{\boldmathX}},{\text{\boldmathZ}}) needs to be estimated (e.g. via regression modeling) in observational studies. We first assume that \pi({\text{\boldmathX}},{\text{\boldmathZ}}) is known, and provide discussions on estimating \pi({\text{\boldmathX}},{\text{\boldmathZ}}) in Section 4. The observed data \{(Y_{i},T_{i},{\text{\boldmathX}}_{i},{\text{\boldmathZ}}_{i}),i=1,\ldots,n\} consists of $n$ independent identically distributed copies of (Y,T,{\text{\boldmathX}},{\text{\boldmathZ}}).

Weighted loss function and synthetic posterior

Our goal is to detect biomarkers among $X$ that are associated with individual treatment effect (ITE) denoted by $\Delta(X)$, e.g. \Delta({\text{\boldmathX}},{\text{\boldmathZ}})={\rm E}[Y|{\text{\boldmathX}},{\text{\boldmathZ}},T=1]-{\rm E}[Y|{\text{\boldmathX}},{\text{\boldmathZ}},T=-1], with adjustment of confounding with $Z$. Traditional approaches for estimating \Delta({\text{\boldmathX}},{\text{\boldmathZ}}) use parametric models including main effect (function of $X$ and $Z$) and interaction effect (function of $X$, $Z$ and $T$), thereby the estimation of main effect (nuisance part) may lead to inefficient estimation of the interaction effect. To avoid the problem, we here employ the weighted loss function (Chen et al., 2017) given by

[TABLE]

where \pi({\text{\boldmathX}}_{i},{\text{\boldmathZ}}_{i})=P(T_{i}=1|{\text{\boldmathX}}_{i},{\text{\boldmathZ}}_{i}) is the propensity score, and $M(\cdot,\cdot)$ is a loss function. As noted in Chen et al. (2017), the choice of $M(u,v)$ determines the measure of treatment effect. For example, under $M(u,v)=(u-v)^{2}$, the minimizer of the true risk ${\rm E}[L_{n}(f)]$ is \{{\rm E}[Y_{i}|{\text{\boldmathX}}_{i},{\text{\boldmathZ}}_{i},T_{i}=1]-{\rm E}[Y_{i}|{\text{\boldmathX}}_{i},{\text{\boldmathZ}}_{i},T_{i}=-1]\}/2, which is difference of expectations between treatment and opposite groups.

Now we consider a linear predictor for the $k$th biomarker: f(X_{k},{\text{\boldmathZ}};\alpha_{k},\beta_{k},{\text{\boldmath\omega}}_{k})=\alpha_{k}+\beta_{k}X_{k}+{\text{\boldmathZ}}^{t}{\text{\boldmath\omega}}_{k}, where $\alpha_{k}$ is the biomarker-specific intercept, $\beta_{k}$ is the effect size of the $k$th biomarker and {\text{\boldmath\omega}}_{k} is a vector of regression coefficients of $Z$. Then, the weighted loss function (1) is reduced to

[TABLE]

Since our primary interest is the effect size $\beta_{k}$, the other parameters $\alpha_{k}$ and {\text{\boldmath\omega}}_{k} are nuisance parameters. We define synthetic likelihood for the effect sizes $\beta_{k}$’s and the nuisance parameters based on the loss function (2), which is given by

[TABLE]

where $A$ is a scaling constant such that A=n^{-1}\sum_{i=1}^{n}\{T_{i}\pi({\text{\boldmathX}}_{i},{\text{\boldmathZ}}_{i})+(1-T_{i})/2\}^{-1}, and the maximizer of $L_{k}$ is the same as the minimizer of (2). Note that w_{i}=A^{-1}\{T_{i}\pi({\text{\boldmathX}}_{i},{\text{\boldmathZ}}_{i})+(1-T_{i})/2\}^{-1} can be interpreted as the weight for the $i$th observation such that $\sum_{i=1}^{n}w_{i}=n$. Regarding the concrete form of the loss function $M(u,v)$, we simply adopt negative log-likelihood functions to make our empirical Bayes approaches valid. Therefore, the logarithm of the synthetic likelihood (3) can be regarded as a weighted log-likelihood.

From (3), $\beta_{k}$ can be interpreted as interaction effect of $T_{i}$ and $X_{ik}$, so that $\beta_{k}=0$ means that the $k$th biomarker is irrelevant to ITE; otherwise, the $k$th biomarker is a predictive biomarker. Hence, we consider statistical testing whether $\beta_{k}$ is zero or not for each $k$. To this end, we first note that $\alpha_{k}$ and {\text{\boldmath\omega}}_{k} are nuisance parameters in (2). For constructing (profile) synthetic likelihood of $\beta_{k}$ from (2), we consider the following two methods:

(Plug-in method).   Compute the point estimator ${\widehat{\alpha}}_{k}$ and {\widehat{\text{\boldmath\omega}}}_{k} of $\alpha_{k}$ and {\text{\boldmath\omega}}_{k}, respectively, obtained by maximizing the synthetic likelihood function (2), and replace $\alpha_{k}$ and {\text{\boldmath\omega}}_{k} in (2) with ${\widehat{\alpha}}_{k}$ and {\widehat{\text{\boldmath\omega}}}_{k}.

• 2.

(Normal-approximation method).   Compute the mode ${\widehat{\beta}}_{k}$ and the inverse value of Hessian at the mode, $s_{k}$, of the synthetic likelihood function (2) with respect to $\beta_{k}$, and define the synthetic likelihood function of $\beta_{k}$ as $\phi({\widehat{\beta}}_{k};\beta_{k},s_{k})$.

The synthetic likelihood of $\beta_{k}$ will be denoted by $\text{PL}_{k}(\beta_{k})$. Since information from $n$ individuals is summarized in ${\widehat{\beta}}_{k}$ and $s_{k}$, the normal approximation method have computational advantage compared to the plug-in method which needs to compute summation over $n$ individuals. However, the normal approximation might be poor when $n$ is not large, which may lead to less power than the plug-in method.

We consider the multiple hypothesis tests, H${}_{0}:\beta_{k}=0$ vs H${}_{1}:\beta_{k}\neq 0$. In order to express the null and non-null biomarkers, we introduce the following latent structure for the effect sizes:

[TABLE]

where $\pi$ is the prior probability of being null, that is, $\pi=P(\beta_{k}=0)$, and the functions $g_{0}$ and $g_{1}$ represent the distributions of the null and non-null biomarkers, respectively. For null biomarkers, we use the one-point distribution on [math], that is, $g_{0}(\cdot)=\delta_{0}(\cdot)$, where $\delta_{c}(x)$ represent the one-point distribution on $x=c$. The form of $g_{1}$ is not specified, so that the latent structure of $\beta_{k}$ can be seen as a semiparametric model where a nonparametric distribution is assumed for non-null biomarkers. Combined with the profile likelihood ${\rm PL}_{k}$ and the prior (4), we define the following synthetic posterior distribution of $\beta_{k}$:

[TABLE]

which are independent for $k=1,\ldots,p$.

The synthetic posterior distribution (5) have hyperparameters $\pi$ and $g_{1}$. We consider an empirical Bayes approach, that is, we estimate these parameters from the synthetic marginal likelihood:

[TABLE]

We employ the smoothing-by-roughening approach (Shen and Louis, 1999) in which the nonparametric estimate of $g_{1}$ is supported by fixed discrete mass points, that is, we approximate $g_{1}(x)$ as $\sum_{\ell=1}^{L}p_{\ell}\delta_{a_{\ell}}(x)$ with mixing probabilities $p_{\ell}$’s such that $\sum_{\ell=1}^{L}p_{\ell}=1$, and fixed knots $a_{1},\ldots,a_{L}$. Then, the above synthetic likelihood can be efficiently maximized by an EM algorithm (Dempster et al., 1977) whose details are provided in Appendix. Note that the mass probabilities on $a_{1},\ldots,a_{L}$ and [math] are $(1-\pi)p_{1},\ldots,(1-\pi)p_{L}$ and $\pi$, respectively, which could be directly estimated via a simpler EM algorithm. However, biomarkers-specific indices directly extracted from the proposed EM algorithm are useful for the optimal discovery procedure described in the subsequent section.

Regarding the selection of knots, some guidance is given in Shen and Louis (1999), which recommends using sufficiently large range of knots to cover the main part of the true underlying distribution. We here suggest a simple data-dependent method. Let $\underline{\beta}$ and $\overline{\beta}$ be the minimum and maximum values of the point estimates of effect sizes $\{{\widehat{\beta}}_{1},\ldots,{\widehat{\beta}}_{p}\}$. Then, we set knots $a_{1},\ldots,a_{L}$ as equally spaced points from $\underline{\beta}$ to $\overline{\beta}$, that is, $a_{\ell}=\underline{\beta}+(\ell-1)(\overline{\beta}-\underline{\beta})/(L-1)$ with $\ell=1,\ldots,L$, where $L$ is the number of knots. The choice of $L$ would be less sensitive to the results as long as $L$ is large to a certain extent. We set $L=100$ as a default choice and we cary out sensitivity analysis of $L$ in our simulation studies in Section 3.

Biomarker-specific indices

Some biomarker-specific indices are useful for screening biomarkers. Let $\gamma_{k}$ be the indicator variable for null/non-null status for the $k$th biomarkers, such that $\gamma_{k}=1$ if the $k$th biomarker is non-null and $\gamma_{k}=0$ otherwise. The value of $\gamma_{k}$ is unknown and has the prior probability $P(\gamma_{k}=1)=1-\pi$. From the synthetic posterior (5), we can compute the posterior probability of being non-null as

[TABLE]

where the integral can be expressed as summation over the discrete mass points used to estimate $g_{1}$. By plugging in the hyperparameter estimates of $\pi$ and $g_{1}$, the estimated posterior probability of being non-null can be obtained. Other biomarker-specific indices (e.g. posterior mean of $\beta_{k}$) can also be estimated based on the synthetic posterior distribution (5).

Screening biomarkers based on the optimal discovering procedure

For ranking biomarkers, we apply the optimal discovery procedure (Storey, 2007; Storey et al., 2007), and use the following statistic:

[TABLE]

which is the model-based version of the optimal discovery statistic (Cao et al., 2009; Noma and Matsui, 2012, 2013). Since the optimal discovery procedure is known to provide maximum expected number of significant biomarkers under given a certain false discovery rate (Noma and Matsui, 2012), we may efficiently detect predictive biomarkers by using the statistic $\text{ODS}_{k}$. We select a set of $\Theta(\lambda)$ of biomarkers whose ODS values are equal to or greater than $\lambda$. To estimate the number of false positives in the selected set, we oblate a false discovery rate (FDR) based on the fitted model (e.g. McLachlan et al., 2006).

[TABLE]

where $|\Theta(\lambda)|$ represents the size of $\Theta(\lambda)$, the number of selected biomarkers. Note that the summation of $\widehat{P}(\gamma_{k}=0|\text{Data})$ (posterior probability of being null) represents the expected number of null (false-positive) biomarkers in $\Theta(\lambda)$. For specified value of $\lambda$, we may compute FDR, so that we may compute the reasonable value of $\lambda$ under user-specified FDR (e.g. 5%).

Simulation study

We here assess the performance the proposed methods together with some alternative methods. We considered two types of responses: binary and censored survival responses. We set $n=1000$ (the number of samples) and $p=3000$ (the number of candidate biomarkers) throughout this study. We generated the covariates $(X_{1},\ldots,X_{p})$ from a mean zero multivariate normal distribution with covariance matrix whose $(i,j)$-element is $(0.1)^{|i-j|}$. We also generated two confounding covariates $(Z_{1},Z_{2})$ from a bivariate normal distribution with $E[Z_{1}]=E[Z_{2}]=0.1X_{1}-0.1X_{10}$, ${\rm Var}(Z_{1})={\rm Var}(Z_{2})=1$ and ${\rm Cov}(Z_{1},Z_{2})=0.2$. The treatment indicator $T$ was indecently generated from a Bernoulli distribution with probability $\text{logit}^{-1}(0.2X_{1}+0.1X_{2}+0.1Z_{1})$, where $\text{logit}^{-1}(x)=\exp(x)/(1+\exp(x))$.

We first consider the case with binary response. We generated $n$ independent samples from the following model:

[TABLE]

where ${\varepsilon}\sim N(0,5^{2})$. For the true parameter values of the main effect, we set $\gamma_{1}=\gamma_{3}=\gamma_{5}=0.2$, $\gamma_{2}=\gamma_{4}=\gamma_{6}=-0.2$, $\delta_{1}=\delta_{3}=\delta_{5}=0.2$, $\delta_{2}=\delta_{4}=\delta_{6}=-0.2$, $\xi_{1}=\xi_{2}=0.1$, and the other values of $\gamma_{k}$’s and $\delta_{k}$’s were set to [math]. For effect sizes $\beta_{k}$ in the interaction term, we consider the following generating distribution

[TABLE]

so that the distribution of $\beta_{k}$ consists of three components, one-point distribution on [math] representing null effect, two normal distributions with positive and negative means representing biomarkers that have positive and negative values on the individual treatment effects, respectively. We considered two cases of null probability $\pi$ of $0.5$ and $0.8$. For the simulated dataset, we applied the proposed ODP methods with the use of negative log-binomial likelihood function for the loss function in which $Z_{1}$ and $Z_{2}$ are adopted as adjustment covariates. As explained in Section 2, we considered two approaches for constructing the synthetic profile likelihood: the plug-in method and the normal-approximation method, which are denoted by ODP-P and ODP-N, respectively. For estimating unknown non-null distribution, we employed the data-dependent method described in the end of Section 2.2 and set $L=100$ as the number of knots. We also applied ODP-N with $L=50,150,200$ to check sensitivity of the choice of $L$. As conventional screening methods, we used the following standardized test statistics for the $k$th biomarkers:

[TABLE]

where ${\widehat{\beta}}_{kj}$ and $s_{kj}$ are estimate and standard error based on the logistic regression with $k$th biomarker and two confounding covariates in treatment $(j=1)$ and control $(j=0)$ groups. The FDR for these tests were estimated using the method by Storey and Tibshirani (2003).

We also consider the case with survival responses. We generated $n$ independent survival times from the following regression model:

[TABLE]

where ${\varepsilon}\sim N(0,5^{2})$ and all the model parameters were the same as those in the previous parts. The censoring time was generated from the uniform distribution $U(20,60)$, which induced censoring survival time $Y$. The censoring rate was around $30\%$ in each simulated dataset. For the simulated dataset, we applied the proposed two ODP methods (ODP-P and ODP-N) with the use of negative log Cox partial likelihood function for the loss function, and the same numbers of knots. As competitors, we adopted the same form of test statistics as (7), where ${\widehat{\beta}}_{kj}$ and $s_{kj}$ are obtained from Cox regression with only $k$th biomarker in treatment $(j=1)$ and control $(j=0)$ groups.

Based on 200 simulations, we calculated the average number of significant biomarkers and true positives at FDR$=$5, 10, 15, or 20% for three types of responses, which are summarized in Table 1. Overall, the proposed methods detected the larger numbers of biomarkers than the standard methods, which clearly shows the efficiency of the propose methods. Moreover, the numbers of true positives of the proposed methods are much larger than the standard methods when we allow larger FDR. Also the efficiency gain of the proposed method emerged as $\pi$ (null probability) gets smaller. Comparing the two proposed methods, ODP-P tends slightly more efficient than ODP-N possibly because the approximation error of the normal approximation used in ODP-N may not be negligible and it may reduce the power of ODP-N. It is also observed that the number of knots has little effect on the results.

Application to a breast cancer clinical study

In this section, we illustrate the practical application of the proposed method in a breast cancer clinical study (Loi et al., 2007). The dataset consisted of 414 estrogen receptor (ER)-positive breast cancer patients with microarray gene expression profiling data. We applied the proposed ODP method to detect the significant gene expression associated with the individual treatment effects. The data are available from the NCBI GEO database (GSE 6532). Here we adopted the time to relapse-free survival as the outcome. After excluding patients with incomplete information, there were 268 and 125 patients receiving tamoxifen and alternative treatments, respectively. For each patient, $p=44928$ gene expression measurements were available, dented by $X_{1},\ldots,X_{p}$, and there is no covariates for adjustment. For estimating the propensity score $\pi(X)$, we employed the following logistic linear regression:

[TABLE]

We estimated the regression coefficients penalized likelihood with lasso (Tibshirani, 1996), where the tuning parameter was selected by 10-fold cross validation. With use of the estimated propensity score, we applied the two types of proposed ODP methods, ODP-P and ODP-N, with negative log-binomial likelihood and $L=200$ as the number of knots. The estimated null probabilities were $0.79$ and $0.72$ for ODP-N and ODP-P, respectively, and estimated underlying distributions of effect sizes are shown in the right panel of Figure 1 with histogram of $t$-values of each biomarker. From Figure 1, both methods produced similar estimates of the underlying distribution, and the true underlying distribution of effect sizes would not have a simple form and seem to have two modes. The numbers of detected microarrays with different FDR values based on the proposed methods and the two standard methods describe in Section 3 are reported in Table 2. The proposed methods were able to detected much larger numbers of microarrays than the direct use of $q$-values with statistics $T_{k}$ and $S_{k}$. Also, ODP-P found more microarrays than ODP-N possibly because the normal approximation is not accurate in this situation where $n$ is around 400.

Discussions

We developed an efficient method for detecting predictive biomarkers associated with individual treatment effects based on the optimal discovery procedure with the synthetic posterior for effect sizes. The key feature of our approach is to employ the idea of direct estimation of individual treatment effect via the weighted loss function (Chen et al., 2017) and combine the semi-parametric prior for interaction effects, which enables flexible and efficient estimation of the underlying signal components. Employing the empirical Bayes method, the signal components can be accurately estimated from the data. Then, the estimated model is combined with the optimal discovery procedure to make up a effective screening procedure under which the treatment effects are heterogenous, i.e. there exist predictive biomarkers. The advantage of the proposed method compared with existing approaches was confirmed through simulation study and the application to breast cancer data.

Although we adopted the independence structure for prior latent distribution for effect sizes of candidate biomarkers, they are typically correlated in practice. However, the optimality of the optimal discovering procedure is still valid even under potential correlations among effect sizes or candidate biomarkers (e.g. Storey, 2007; Noma and Matsui, 2012) and FDR is adequately controlled as long as underlying (marginal) distribution of $\beta_{k}$ is correctly specified. This condition would be naturally satisfied since we employed semiparametric hierarchical mixture modeling. Moreover, the underlying distribution of $\beta_{k}$ could be unbiasedly estimated even under potential correlation while the standard error could be larger than that of the model with taking account of potential correlations, as investigated in Matsui and Noma (2011). In fact, our simulation study in Section 3 adopted situations where the biomarkers are correlated, thereby the independence assumption for the effect sizes is violated, but the proposed method performed reasonably well.

Although we employed the weighted loss function Chen et al. (2017) for constructing synthetic posterior distribution (5), other types of loss functions can be readily used in the proposed method by simply changing the formula of synthetic likelihood (2). For instance, A-learning method (Murphy, 2003; Lu et al., 2013; Ciarleglio et al., 2015), which also can directly estimate individual treatment effect, would be useful alternatives. Since the detailed investigation and comparison among several types of loss functions would extend the scope of this paper, it is left to a future study.

Acknowledgements

This work was supported by JST CREST Grant Number JPMJCR1412, JST AIP-PRISM Grant Number JPMJCR18Z9 and JSPS KAKENHI Grant Numbers JP17K19808, JP15K15954, JP16H07406, Japan.

Appendix (EM algorithm)

Let $\phi=(\pi,p_{1},\ldots,p_{L})$ be the vector of unknown hyperparameters. We introduce two types of latent variables, $z_{k}$ and $w_{k}$ $(k=1,\ldots,p)$, such that $P(z_{k}=1)=1-P(z_{k}=0)=\pi$ and $P(w_{k}=\ell)=p_{\ell}$. Then, logarithm of the synthetic likelihood (6) can be augmented as

[TABLE]

With current parameter values $p_{\ell}^{(t)}$ and $\pi^{(t)}$, we have the following objective function in the M-step:

[TABLE]

where

[TABLE]

The updating steps are given by

[TABLE]

Hence, each iteration in the EM algorithm requires computing $\xi_{k}^{(t)}$ and $\eta_{k\ell}^{(t)}$ as the E-step and updating $\pi$ and $p_{\ell}$’s as the M-step, which is continued until numerical convergence.