Bayesian Regression in Linear Regression Analysis

Bayesian Regression xxia@uottawa.ca http://dambe.bio.uottawa.ca

Data Table 1. Data for illustrating simple linear regression with dependent variable Y (body weight loss) and independent variable X (humidity). Simple linear regression leads to estimates of = 8.704 and = -0.053, predicted Y listed as EY, and residual as Y EY Y EY Index X 1 2 3 4 5 6 7 8 9 Y 8.98 8.7040 0.2760 8.14 8.0654 0.0746 29.5 6.67 7.1340 -0.4640 43 6.08 6.4155 -0.3355 53 5.9 5.8833 0.0167 62.5 5.83 5.3776 0.4524 75.5 4.68 4.6858 -0.0058 85 4.2 4.1801 0.0199 93 3.72 3.7544 -0.0344 EY 9 0 12 Weight Loss 6 3 0 20 40 60 80 100 Humidity

Bayesian regression ? ? ? ? ?(????|?,? ? ? ? ? ?(????|?,? ???? ? ?,? ???? = We learned before three ways to deal with the integration in denominator: 1) using conjugate prior probability so that we do not need to compute the integration, using discretization or MC integration to approximate the integration, and using MCMC so that we do not need the integration 2) 3) Here we will use MCMC again, which may be visualized as a drunkard wobbling on a 2-dimensional terrain with one dimension being a (intercept) and the other being b (slope). Every step is recorded as (ai, bi), but some steps are discarded.

The prior: a Inferences: 1. Y (WtLoss) for insects with 0 humidity & no food: hould be greater than 0, i.e., f(a>0) > 0 2. Y should not be greater than the original body weight (Wo, assumed to be ~30 milligrams) The simplest prior: a uniform distribution with ? ? = 1/?? for the entire closed range of [0, Wo]. A more informative prior? Further inferences: 1. We expect ? ? = 0 = 0 and ? ? = ?? = 0, but ? 0 < ? < ?? > 0. So this is not uniform. 2. We know that beetles have about 65% of their body made of water, and that they lose only about half of their body water under extremely dry conditions (0% humidity). 3. This suggests the expected water loss is perhaps about Wo*65%/2 = 9.75 milligrams. 4. A bit of dry mass would also be lost because of starvation, so the total weight loss at 0 humidity may be ~10 milligrams. Thus, we should find a probability density distribution that peaks around ? = 10, but is effectively 0 when ? = 0 or ? = ??= 30 milligram. A gamma distribution with ? = 4,? = 2.5, which implies ? = 10,? = 5, seems to fit these suggestions. We should contrast the two priors and see their consequences

The gamma prior for a 0.09 = 4, = 2.5 mean = 10 std = 5 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0 5 10 15 20 25 30 Figure 1. Gamma distribution with ? = 4,? = 2.5, which implies ? = 10,? = 5.

The prior for b Inferences: 1. b should be negative, i.e., increased humidity will result in reduced body weight loss, so the maximum of b is 0. 2. We have previously reasoned that the body weight loss would be around 10 milligram when humidity is 0, but it could be as much as 20 or even 25. If the maximum possible a is 25, and the most extreme case with no body weight loss at the highest humidity (93%), then we would have a minimum slope equal to ????= ? 0 93= 0.2688 3. This would suggest a uniform distribution of ? ? = between 0 and 0.2688. Do you want to have a more informative prior? You have to support it with reason. ?=25 0 1 0 ????= 3.72 in the range lnPrior<-function(B) { a = B[1] b = B[2] a.prior = dgamma(a, shape=4, scale=2.5, log = T) # a.prior = dunif(a, min=0, max = 30, log = T) # alternative a.prior b.prior = dunif(b, min = -0.2688, max = 0, log = T) return(a.prior+b.prior) }

Likelihood The probability density?(??|? = 0,?2 , where ??= ?? ??, from a normal distribution with mean equal to 0 and variance (s2) equal to the variance of the residuals (r). lnL <- function(B, X, Y) { a <- B[1] b <- B[2] EY <- a+b*X Res <- Y-EY lnL <- dnorm(Res,0,sd(Res),log=TRUE) return(sum(lnL)) } This lnL function returns?(????|?,? , i.e., the probability density of observing 'data' given the linear model with parameters a and b.

N <- 50000; max.a <- 30;min.a <- 0; max.b <- 0; min.b <- -0.2688; X<-c(0,12,29.5,43,53,62.5,75.5,85,93); Y<-c(8.98,8.14,6.67,6.08,5.9,5.83,4.68,4.2,3.72); sd.a <- 0.25; sd.b <- 0.004 # use 1/20 of sd.a.prior and 1/20 of sd.b.prior rnd <- sample(1:10000,N,replace=T)/10000 # between 0 and 1 a <- rep(0,N); b<-rep(0,N) a[1] <- 10; b[1] = -0.1 for (i in seq(1:(N-1))) { a[i+1] <- rnorm(1,mean=a[i],sd=sd.a) b[i+1] <- rnorm(1,mean=b[i],sd=sd.b) if(a[i+1] < min.a) { a[i+1] <- min.a; } else if(a[i+1] > max.a) { a[i+1] <- max.a; } if(b[i+1] < min.b) { b[i+1] <- min.b; } else if(b[i+1] > max.b) { b[i+1] <- max.b; } jointP.new <- lnPrior(c(a[i+1],b[i+1])) + lnL(c(a[i+1],b[i+1]),X=X,Y=Y) jointP.old <- lnPrior(c(a[i],b[i])) + lnL(c(a[i],b[i]),X=X,Y=Y) if(jointP.new < jointP.old) { Alpha <- exp(jointP.new -jointP.old ) if(rnd[i] > Alpha) { a[i+1] <- a[i]; b[i+1] <- b[i] a[i] <- 0; b[i] <- 0 } } } MCMC Check rejection rate: length(a[a==0])/N

Sample the last 10000 steps par(mfrow=c(1,2)); posta <- a[(N-9999):N]; posta <- posta[posta>0] freq <- hist(posta) meana<-mean(posta); sda<-sd(posta); meana; sda postb <- b[(N-9999):N]; postb <- postb[postb<0] freq <- hist(postb) meanb<-mean(postb); sdb<-sd(postb); meanb;sdb # simple way of generating credible interval quantile(posta,c(0.025,0.975)); quantile(postb,c(0.025,0.975))

MCMC Results 70 (C) (A) a=8.68504 sa=0.23171 150 a= 8.65036 sa=0.21332 60 50 Frequency 100 Frequency 40 30 50 20 10 0 0 8.0 8.5 9.0 9.5 8.0 8.2 8.4 8.6 8.8 9.0 9.2 posta posta (D) (B) 100 b=-0.05288 sb= 0.00422 b =-0.05217 sb= 0.00399 60 80 Frequency Frequency 60 40 40 20 20 0 0 -0.065 -0.060 -0.055 -0.050 -0.045 -0.040 -0.065 -0.060 -0.055 -0.050 -0.045 -0.040 postb postb Figure 2. R output from running Bayesian MCMC to obtain regression coefficients (a and b) and their standard errors (sa and sb). (A) and (B): output with the gamma prior for a as depicted in Figure 1. (C) and (D): output with a uniform prior for a.

Credible & high density intervals res.a<-hist(posta,breaks=100); res.b<-hist(postb,breaks=100); #1 a<-res.a$mids; da<-res.a$density; b<-res.b$mids; db<-res.b$density spl.a<-smooth.spline(x=a,y=da,df=10); spl.b<-smooth.spline(x=b,y=db,df=10); # create equal-interval (ei) a and b, named ei.a and ei.b ei.a<-seq(min(a),max(a),by=0.001); ei.b<-seq(min(b),max(b),by=0.00001); ei.da<-pmax(0,predict(spl.a,x=ei.a)$y); #6 ei.db<-pmax(0,predict(spl.b,x=ei.b)$y); samp.ei.a<-sample(ei.a,1e5,replace=TRUE,prob=ei.da) samp.ei.b<-sample(ei.b,1e5,replace=TRUE,prob=ei.db) quantile(samp.ei.a,c(0.025,0.975)); #10 produce credible interval quantile(samp.ei.b,c(0.025,0.975)); #11 # Get high-density interval samp.ei.a <- sample(ei.da, 1e5, replace = TRUE, prob = ei.da) samp.ei.b <- sample(ei.db, 1e5, replace = TRUE, prob = ei.db) critDensity.a <- quantile(samp.ei.a, 0.05) # horizontal line to cut the peak critDensity.b <- quantile(samp.ei.b, 0.05) #16 ei.a.high <-ei.a[ei.da >= critDensity.a]; ei.b.high <-ei.b[ei.db >= critDensity.b]; lo.a<-ei.a.high[1]; hi.a<-ei.a.high[length(ei.a.high)] lo.b<-ei.b.high[1]; hi.b<-ei.b.high[length(ei.b.high)] lo.a; hi.a; lo.b; hi.b; #21 produce high-density interval

Credible & high density intervals par(mfrow=c(1,2)); plot(da~a,type="l"); lines(spl.a, col="blue"); abline(v=lo.a); abline(v=hi.a); abline(h=critDensity.a); plot(db~b,type="l"); lines(spl.b, col="blue"); abline(v=lo.b); abline(v=hi.b); abline(h=critDensity.b); #30 (B) (A) 2.0 100 80 1.5 60 db da 1.0 40 0.5 20 0.0 0 8.0 8.5 9.0 9.5 -0.065 -0.060 -0.055 -0.050 -0.045 -0.040 -0.035 a b