Bayesian Variable Selection Method for Modeling Dose-Response Microarray Data Under Simple Order Restrictions Bayes2013, Rotterdam Martin Otava I-Biostat, Universiteit Hasselt 22.05.2013 Martin Otava (I-Biostat, UHasselt) BVS 1 / 20
Research team Hasselt University (Belgium): Martin Otava Ziv Shkedy Durham University (UK): Adetayo Kasim Imperial College London (UK): Bernet Kato Zoetis (Belgium): Dan Lin Janssen Pharmaceutica (Belgium): Luc Bijnens Hinrich W.H. G¨ ohlmann Willem Talloen Martin Otava (I-Biostat, UHasselt) BVS 2 / 20
Model formulation Dose-response modeling Increasing dose of therapeutical compound. Variety of possible responses: Toxicity. Inhibition or stimulation. Gene expression level. Goal: Determine if there is any relationship. If so, what is the shape of the profile. Select threshold doses (e.g. MED). 6.5 ● ● Gene expression 7.8 Gene expression 6.4 ● ● ● 7.6 6.3 ● ● ● ● ● ● 6.2 7.4 ● ● ● ● ● 6.1 7.2 ● ● ● ● ● 6.0 7.0 ● ● ● 0 1 2 3 0 1 2 3 Dose Dose Martin Otava (I-Biostat, UHasselt) BVS 3 / 20
Model formulation Order constraints Compound effect becomes stronger when dose is increased. Monotone restriction 2.0 (non-decreasing or 1.8 non-increasing). 1.6 Response 1.4 Zero effect is meaningful. 1.2 1.0 No parametrical assumptions 0.8 about dose-response curve 0 1 2 3 shape. Dose Martin Otava (I-Biostat, UHasselt) BVS 4 / 20
Model formulation Basic Model One-way ANOVA model formulation: Y ij = µ i + ε ij i = 0 , . . . , K − 1 ε ij ∼ N (0 , σ 2 ) j = 0 , 1 , 2 , . . . , n i = ⇒ necessary to incorporate order constraints. Testing the hypothesis H 0 : µ 0 = µ 1 = µ 2 = . . . = µ K − 1 against ordered alternative (one inequality strict) H up : µ 0 ≤ µ 1 ≤ µ 2 ≤ . . . ≤ µ K − 1 H dn : µ 0 ≥ µ 1 ≥ µ 2 ≥ . . . ≥ µ K − 1 Martin Otava (I-Biostat, UHasselt) BVS 5 / 20
Bayesian variable selection model Reformulation of model New notation (non-decreasing trend): µ 0 , i = 0 , i E ( Y ij ) = µ i = � µ 0 + δ ℓ , i = 1 , . . . , K − 1 ℓ =1 with priors: N ( η µ , σ 2 µ 0 ∼ µ ) , N ( η δ i , σ 2 ∼ δ i ) I (0 , A ) , i = 1 , . . . , K − 1 . δ i ⇒ δ i ≥ 0 . Martin Otava (I-Biostat, UHasselt) BVS 6 / 20
Bayesian variable selection model Set of all models g_0 g_1 g_2 1.5 1.5 1.5 Response Response Response 1.0 1.0 1.0 0.5 0.5 0.5 0.0 0.0 0.0 0 1 2 3 0 1 2 3 0 1 2 3 Dose Dose Dose g_3 g_4 g_5 1.5 1.5 1.5 Response Response Response 1.0 1.0 1.0 0.5 0.5 0.5 0.0 0.0 0.0 0 1 2 3 0 1 2 3 0 1 2 3 Dose Dose Dose g_6 g_7 1.5 1.5 Response Response 1.0 1.0 0.5 0.5 0.0 0.0 0 1 2 3 0 1 2 3 Dose Dose Martin Otava (I-Biostat, UHasselt) BVS 7 / 20
Bayesian variable selection model Sub-hypotheses H up : µ 0 ≤ µ 1 ≤ µ 2 ≤ . . . ≤ µ K − 1 Model Up: Mean Structure z µ 0 = µ 1 = µ 2 = µ 3 (0,0,0) g 0 g 1 µ 0 < µ 1 = µ 2 = µ 3 (1,0,0) µ 0 = µ 1 < µ 2 = µ 3 (0,1,0) g 2 g 3 µ 0 < µ 1 < µ 2 = µ 3 (1,1,0) g 4 µ 0 = µ 1 = µ 2 < µ 3 (0,0,1) µ 0 < µ 1 = µ 2 < µ 3 (1,0,1) g 5 g 6 µ 0 = µ 1 < µ 2 < µ 3 (0,1,1) (1,1,1) g 7 µ 0 < µ 1 < µ 2 < µ 3 Martin Otava (I-Biostat, UHasselt) BVS 8 / 20
Bayesian variable selection model Modification to BVS The distribution of δ is continuous. = ⇒ probability of all models except one equals zero! Instead of only sampling δ i we need to select which δ i occurs in model. Let be z i indicator of δ i occurring in the model. � 1 , is included in the model , δ i z i = 0 , δ i is not included in the model . i � ⇒ E ( Y ij ) = µ 0 + z ℓ δ ℓ . ℓ =1 Martin Otava (I-Biostat, UHasselt) BVS 9 / 20
Bayesian variable selection model BVS model formulation Basic model: Y ij ∼ N ( µ i , σ 2 ) Hyper Priors: Modeling of mean: σ − 2 Γ(10 − 3 , 10 − 3 ) , ∼ i N (0 , 10 6 ) , η µ ∼ � E ( Y ij ) = µ i = µ 0 + z ℓ δ ℓ . σ − 2 Γ(10 − 3 , 10 − 3 ) , ∼ µ ℓ =1 N (0 , 10 6 ) , η δ i ∼ Priors: σ − 2 Γ(10 − 3 , 10 − 3 ) . ∼ δ i N ( η µ , σ 2 π i ∼ U(0 , 1) . µ 0 ∼ µ ) , N ( η δ i , σ 2 ∼ δ i ) I (0 , A ) , δ i z i ∼ Bernoulli( π i ) , Martin Otava (I-Biostat, UHasselt) BVS 10 / 20
Results interpretation Posterior mean of µ i Posterior distribution for all dose-specific means. Use posterior mean of such distribution as our estimation. Connection of Bayesian model averaging. = ⇒ posterior model probabilities are weights. R � µ BVS = ˆ w r ˆ µ r r =0 Martin Otava (I-Biostat, UHasselt) BVS 11 / 20
Results interpretation Posterior probability of model Vector z = ( z 1 , . . . , z K − 1 ) uniquely defines the model. i =1 z i 2 i − 1 = Transformation G ( z ) = 1 + � K − 1 ⇒ unique value for each model. In each MCMC iteration we sample one vector z = ( z 1 , . . . , z K − 1 ). Posterior mean of indicator G ( z ) = r + 1 translates into posterior probability of the model g r . = ⇒ For posterior probabilities holds: P [ G ( z ) = r + 1 | data] = P ( g r | data). Martin Otava (I-Biostat, UHasselt) BVS 12 / 20
Results interpretation Example: BVS model Incorporating models with equal means results into less decreasing profile. ● ● ● ● ● Posterior means are averages of ● ● ● ● 35 ● ● ● ● ● ● means of particular models at ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● each MCMC iteration. ● ● ● ● Weight ● ● 30 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● 25 ● ● ● R Order restricted ● � ¯ BVS µ BVS = ˆ P ( g r | data)ˆ ● ● µ r 20 ● ● ● r =0 0 1 2 3 Dose Connection to model averaging. Martin Otava (I-Biostat, UHasselt) BVS 13 / 20
Results interpretation Example: Posterior probabilities Posterior probabilities of 0.6 particular models. 0.5 Posterior probability 0.4 Model g 0 represents H 0 . 0.3 Model g 1 is strongly supported 0.2 by the data. 0.1 0.0 Connection to model selection. g_0 g_1 g_2 g_3 g_4 g_5 g_6 g_7 Model Martin Otava (I-Biostat, UHasselt) BVS 14 / 20
Hypothesis testing Hypothesis testing Depends on: data on hand, prior distributions, set of alternative hypotheses. We use objective priors and consider the set of all possible alternative hypotheses. Use ¯ P ( g 0 | data), estimation of P ( H 0 | data), to reject H 0 . Questions: How to select threshold for deciding if H 0 is rejected? There is no straightforward control mechanism like Type I error. Simulation study can give us insight in the properties of BVS. Martin Otava (I-Biostat, UHasselt) BVS 15 / 20
Hypothesis testing Simulation study Under the H 0 and under model g 7 . P ( H 0 | data) < τ used as criterion for rejecting H 0 by BVS. P H 0 (data ∗ ) < τ used as criterion for rejecting H 0 by LRT and MCTs. What happens to false rejections and false non-rejections while varying threshold τ ? When maintaining approximately same empirical Type I error as MCTs or LRT, BVS seems to achieve similar power. How to select threshold for BVS in practice? = ⇒ future research. Martin Otava (I-Biostat, UHasselt) BVS 16 / 20
Hypothesis testing Simulation study - Results 1.0 1.0 0.8 0.8 BVS LRT MCT−W Type I error 0.6 0.6 MCT−M Power BVS 0.4 0.4 LRT MCT−W MCT−M 0.2 0.2 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Threshold Threshold Martin Otava (I-Biostat, UHasselt) BVS 17 / 20
Discussion Conclusion Model uncertainty taken into account! Model selection: ¯ P ( g r | data). µ = � R r =0 ¯ Estimation of means: ˆ P ( g r | data)ˆ µ r . Inference: ¯ P ( g 0 | data). BVS framework address all perspectives simultaneously. According to simulations seems to perform comparably with LRT and MCTs. Martin Otava (I-Biostat, UHasselt) BVS 18 / 20
Discussion Future research How to select threshold for rejecting H 0 using P ( H 0 | data)? How to fit BVS models with different types of restrictions (e.g. umbrella profiles)? How do BVS models behave when used for multiplicity adjustment? Martin Otava (I-Biostat, UHasselt) BVS 19 / 20
Discussion Thank you for your attention! Martin Otava (I-Biostat, UHasselt) BVS 20 / 20
Recommend
More recommend