Download Bayesian Inference and Hierarchical Bayes Model in Genomics and more Study notes Biostatistics in PDF only on Docsity! Outline Intro Binomial Model Normal Model Hierarchical Bayes Intro to Genomics – Bayesian Inference (I) Grzegorz A. Rempala Department of Biostatistics Medical College of Georgia August 27-Sept 8, 2008 1/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes 1 Intro: Bayesian Inference Likelihood Function Bayes Theorem: Prior and Posterior Conjugate Prior 2 Binomial Model Example One: Beta–Binomial 3 Normal Model Example Two: Normal–Normal School Data: Hierarchical Bayesian Model 4 Hierarchical Bayes Hyper Prior Empirical Bayes 2/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes R Script ## Binomial Likelihood bin.lik <- function(p){ p^4 * (1-p)^6 } p <- seq(0,1,by=0.02) plot(p,bin.lik(p),ylab="likelihood",type="l") abline(v=0.4,lty=2) 5/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes 0.0 0.2 0.4 0.6 0.8 1.0 0. 00 00 0. 00 02 0. 00 04 0. 00 06 0. 00 08 0. 00 10 0. 00 12 p lik el ih oo d Figure: Binomial Likelihood 6/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes Bayes Theorem: Prior and Posterior In order to obtain a probability distribution of a parameter we need to perform Bayesian Inference. First we need a prior distribution. Bayes Theorem: f(θ|x) = f(x|θ)f(θ) f(x) ∝ f(x|θ)f(θ) Posterior(θ) ∝ Likelihood(θ) × Prior(θ) f(x) is normalizing constant to ensure the posterior sum/integrate to one. Prior (θ) represents the information about an parameter θ, before observing the data. Prior distribution is updated by the data to yield Posterior(θ). 7/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes Suppose we have prior belief that the coin is fair, i.e. p = 1/2. α β Mean Var 1 1 1/2 0.083 6 6 1/2 0.019 12 12 1/2 0.01 18 18 1/2 0.006 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 p P ro b. D en si ty (18,18) (12,12) (6,6) (1,1) Figure: Shapes for Beta priors. Note that Beta (1,1) is non-informative. 10/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes R Script (Beta Distribution) a <- c(1,6,12,18) b <- c(1,6,12,18) #Variance a*b/( (a+b)^2 * (a+b+1)) # plot x <- seq(0,1,by=0.01) y <- dbeta(x, a[4], b[4]) plot(x,y, type="l",xlab="p" , ylab="Prob. Density") text(0.5,dbeta(0.5,a[4],a[4]),"(18,18)") y <- dbeta(x, a[3], b[3]) lines(x, y); text(0.5,dbeta(0.5,a[3],a[3]),"(12,12)") y <- dbeta(x, a[2], b[2]) lines(x, y); text(0.5,dbeta(0.5,a[2],a[2]),"(6,6)") y <- dbeta(x, a[1], b[1]) lines(x, y); text(0.5,dbeta(0.5,a[1],a[1]),"(1,1)") 11/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes Suppose n = 5, y = 1 and y/n = 0.2. Posterior distribution f(θ|y) satisfies f(θ|y) ∝ pα+y−1(1− p)β+n−y−1 Thus it has to be Beta(α+ y, β + n− y). Hence we have Posterior mean E(θ|y) = α+ y α+ β + n . Posterior Variance var(θ|y) = E(θ|y)[1− E(θ|y)] α+ β + n+ 1 . 12/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes Question: What happens if we increase no of flips n (why) ? 0.0 0.2 0.4 0.6 0.8 1.0 0 2 4 6 8 p P os t. D en si ty m= 0.326 v= 0.003 m= 0.297 v= 0.003 m= 0.258 v= 0.003 m= 0.212 v= 0.003 Figure: Beta posteriors for n = 50 and y = 10 using four previous priors. 15/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes Normal Model Let the scores of a sample of n students in a standardized test be : y1, y2, . . . , yn. Assume they follow normal distribution, i.e. yi ∼ N(µ, σ2) where σ2 is known. Therefore ȳ ∼ N ( µ, σ2/n ) From another school district, we know that the mean score change for each school in the same district following µ ∼ N(µ0, τ2) To estimate the mean test score of the school, we can use MLE, i.e. ȳ = ∑n i=1 yi n . 16/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes We could also use the Bayesian method. The likelihood is f(ȳ|µ) ∝ n∏ i=1 exp ( −n 2 (ȳ − µ)2 σ2 ) . The prior is f(µ) ∝ exp ( −1 2 (µ− µ0)2 τ20 ) Therefore the posterior distribution is f(µ|y) ∝ exp ( −1 2 [ n σ2 (ȳ − µ)2 + 1 τ20 (µ− µ0)2 ]) ∝ exp ( −1 2 (µ− µ1)2 τ21 ) 17/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes Effect of Changes in Prior Variance 50 60 70 80 90 100 0. 00 0. 04 0. 08 tau = 5 score y. po st er io r 50 60 70 80 90 100 0. 00 0 0. 01 0 0. 02 0 0. 03 0 tau = 20 score y. po st er io r 50 60 70 80 90 100 0. 00 5 0. 01 5 0. 02 5 tau = 30 score y. po st er io r 50 60 70 80 90 100 0. 01 0 0. 01 5 0. 02 0 tau = 60 score y. po st er io r Figure: 20/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes R Script n <- 4; y.bar <- 80.5; sigmasq <- 100; mu0 <- 60; tau <- c(5,20,30,60); posterior.mean <- (mu0/tau + n*y.bar/sigmasq)/(1/tau + n/sigmasq); posterior.var <- 1/(1/tau + n/sigmasq); score <- seq(50,100,by=0.5); par(mfrow=c(2,2)); for(i in 1:4){ y.prior <- dnorm(score,mu0,tau[i]); y.posterior <- dnorm(score,posterior.mean[i], posterior.var[i]); plot(score,y.posterior,type="l", main=paste("tau =",tau[i])); abline(v=posterior.mean[i],lty=2); lines(score,y.prior,lty=4,col="blue"); } 21/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes School Example Suppose we have selected randomly 5 schools in the same school district. A sample of students’ score from each school is as follows. The variance σ2 is assumed known. School Average No of Students Variance σ2 A 80.5 (90, 75, 92, 65) 4 100 B 85 (87, 84) 2 100 C 70 (80,90,50,55,75) 5 100 D 45 ((48, 45, 42) 3 100 E 65 (55, 60 75, 70) 4 100 22/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes R Script (Part II) i = 1; posterior.mean <- (mu0/tau + n[i]*y.bar[i]/sigmasq)/ (1/tau + n[i]/sigmasq); plot(tau,posterior.mean,type="l",ylim=c(20,90), ylab="Posterior Mean"); text(tau[length(tau)],posterior.mean[length(tau)]+1, school[i]); for (i in 2:5){ posterior.mean <- (mu0/tau + n[i]*y.bar[i]/sigmasq)/ (1/tau + n[i]/sigmasq); lines(tau,posterior.mean); text(tau[length(tau)],posterior.mean[length(tau)]+1, school[i]); } 25/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes Hierarchical Bayes Model Adding complexity to the model: In previous analysis, σ2 is assumed to a known constant. We can also put a prior to the variance. The standard conjugate prior for variance is Inverse-Gamma distribution. Also the prior parameters µ0 and τ 2 are assumed to be fixed. In real situation, we are not certain about value of these parameters, therefore we can put another prior on these prior parameters we call it a hyper prior and parameters in the hyper-prior are called hyper parameters. Generally, we don’t have nice closed form solution to this hierarchical model; we need to use MCMC methods. 26/30 Greg Rempala Intro to Genomics Outline Intro Binomial Model Normal Model Hierarchical Bayes Hierarchical Bayes (II) 27/30 Greg Rempala Intro to Genomics