#************************# # R code for chapter 9 # # Adam Branscum # # 1-22-2009 # #************************# ############ Bank salary data Bankdata <- read.table(file="H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData\\BankSalaryData.txt",header=T, sep="") attach(Bankdata) options(contrasts=c("contr.treatment","contr.poly")) Bankdata # Figure 9.1. Scaterplot matrix pairs(~BegSal+Educ+Exper+Time,labels=c("Beginning Salary","Education","Experience","Time"),data=Bankdata) male <- Sex-1 # Sex 1=female; 2=male n <- length(BegSal) r <- 5 fit <- lm(BegSal~Educ+Exper+Time+male) summary(fit) SSE <- t(BegSal-fit\$fitted.values)%*%(BegSal-fit\$fitted.values) MSE <- SSE/(n-r) tauhat <- 1/MSE sdtau <- tauhat*sqrt(2/(n-r)) new <- data.frame(Educ=16, Exper=54.5, Time=33, male=1) predict.lm(fit, newdata=new,interval="confidence",se.fit=T) predict.lm(fit, newdata=new,interval="prediction",se.fit=T) a <- qchisq(0.025,88) b <- qchisq(0.975,88) a/SSE b/SSE sqrt(SSE/b) sqrt(SSE/a) ############ FEV data FEVdata <- read.table(file="H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData\\FEVdataAge10to19.txt",header=T, sep="\t") attach(FEVdata) options(contrasts=c("contr.treatment","contr.poly")) library(R2WinBUGS) Ytilde <- c(2.8,3,4,3.3) D <- diag(c(.04,.04,.04,.09)) Xtilde <- matrix(c(1,11,0,0, 1,13,1,13, 1,16,0,0, 1,18,1,18),4,4,byrow=T) Xtildeinv <- solve(Xtilde) beta0 <- c(t(Xtildeinv %*% Ytilde)) C0 <- Xtildeinv %*% D %*% t(Xtildeinv) C0inv <- solve(C0) a <- 1.73 b <- 0.78 n <- length(FEV) r <- dim(Xtildeinv)[1] FEVdataBUGS <- list("n","r","beta0","C0inv","a","b","FEV","Age","Smoke") parameters <- c("beta","tau","meanFEVs","meanFEVns","RM","MD","FEV20s","FEV20ns") inits <- list(list(tau=1,beta=c(0,0,0,0),FEV20s=2, FEV20ns=2)) FEV.fit <- bugs(FEVdataBUGS, inits, parameters,"FEVWBModel.txt", working.directory = "H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=60000, n.thin=1, n.burnin=10000) print(FEV.fit,digits=3) attach.bugs(FEV.fit) # Code for proper reference prior library(R2WinBUGS) Xtilde <- matrix(c(1,11,0,0, 1,13,1,13, 1,16,0,0, 1,18,1,18),4,4,byrow=T) Xtildeinv <- solve(Xtilde) r <- dim(Xtildeinv)[1] beta0 <- rep(0,r) c1 <- 0.001 C0inv <- diag(rep(c1,r)) a <- 0.001 b <- 0.001 n <- length(FEV) FEVdataBUGS <- list("n","r","beta0","C0inv","a","b","FEV","Age","Smoke") parameters <- c("beta","tau","meanFEVs","meanFEVns","RM","MD","FEV20s","FEV20ns") inits <- list(list(tau=1,beta=c(0,0,0,0),FEV20s=2, FEV20ns=2)) FEV.fit1 <- bugs(FEVdataBUGS, inits, parameters,"FEVWBModel.txt", working.directory = "H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=60000, n.thin=1, n.burnin=10000) print(FEV.fit1,digits=3) attach.bugs(FEV.fit1) ## Plot of Bayesian posterior means from a linear regression with interaction term and BCJ prior for the FEV data BayesSmoker <- FEV.fit\$mean\$meanFEVs BayesNonsmoker <- FEV.fit\$mean\$meanFEVns Age1 <- seq(10,19,1) plot(Age,FEV,xlab="Age", ylab="FEV", lwd=1.5) lines(Age1,BayesSmoker,type='l', lty=1, lwd=2,cex=1.2) lines(Age1,BayesNonsmoker,type='l', lty=2, lwd=2,cex=1.2) legend("bottomright", c("Smoker", "Nonsmoker"), lty=c(1,2), lwd=2) ## Compute CPO statistics and LPML. Could have put the appropriate lines in the code above. library(R2WinBUGS) mtilde0 <- c(2.8,3,4,3.3) cov0 <- diag(c(.04,.04,.04,.09)) prec0 <- solve(cov0) Xtilde <- matrix(c(1,11,0,0, 1,13,1,13, 1,16,0,0, 1,18,1,18),4,4,byrow=T) Xtildeinv <- solve(Xtilde) betastar <- c(t(Xtildeinv%*%mtilde0)) C0 <- Xtildeinv %*% cov0 %*% t(Xtildeinv) C0inv <- solve(C0) a <- 1.73 b <- 0.78 n <- length(FEV) r <- dim(Xtildeinv)[1] FEVdataBUGS <- list("n","r","betastar","C0inv","a","b","FEV","Age","Smoke") parameters <- c("beta","tau","CPOinv") inits <- list(list(tau=1,beta=c(0,0,0,0))) FEV.fit <- bugs(FEVdataBUGS, inits, parameters,"FEVWBModel1.txt", working.directory = "H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=12000, n.thin=1, n.burnin=2000) print(FEV.fit,digits=3) attach.bugs(FEV.fit) CPO <- 1/FEV.fit\$mean\$CPOinv LPML <- sum(log(CPO)) #### CODA for convergence diagnostics (not in book) mtilde0 <- c(2.8,3,4,3.3) cov0 <- diag(c(.04,.04,.04,.09)) prec0 <- solve(cov0) Xtilde <- matrix(c(1,11,0,0, 1,13,1,13, 1,16,0,0, 1,18,1,18),4,4,byrow=T) Xtildeinv <- solve(Xtilde) betastar <- c(t(Xtildeinv%*%mtilde0)) C0 <- Xtildeinv %*% cov0 %*% t(Xtildeinv) C0inv <- solve(C0) a <- 1.73 b <- 0.78 n <- length(FEV) r <- dim(Xtildeinv)[1] FEVdataBUGS <- list("n","r","betastar","C0inv","a","b","FEV","Age","Smoke") parameters <- c("beta","tau") inits <- list(list(tau=1, beta=c(0,0,0,0),FEV20s=2, FEV20ns=2), list(tau=.1, beta=c(10,10,10,10)), list(tau=5, beta=c(-10,-10,-10,-10)), list(tau=10, beta=c(10,-10,10,-10)), list(tau=15, beta=c(-10,10,-10,10)), list(tau=20, beta=c(10,-10,0,0)), list(tau=.5, beta=c(0,0,5,-10)), list(tau=1, beta=c(5,5,-10,0))) FEV.fit <- bugs(FEVdataBUGS, inits, parameters,"FEVWBModel.txt", codaPkg=TRUE, working.directory = "H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=8, n.iter=12000, n.thin=1, n.burnin=2000) attach.bugs(FEV.fit) setwd("H:/MyDocuments/BAYES/LinearRegression/CodeandData") library(coda) codamenu() ######## Section 9.6. Model diagnostics n <- length(FEV) LR <- lm(FEV~Age+Smoke+Age:Smoke) summary(LR) beta <- coef(LR) X <- cbind(rep(1,n),Age,Smoke,Age*Smoke) rstudsmoke <- rstudent(LR)[Smoke==1] rstudnonsmoke <- rstudent(LR)[Smoke==0] par(mfrow=c(2,2)) plot(fitted(LR),resid(LR),xlab="Fitted Values", ylab="Residuals",main="") lines(lowess(fitted(LR),resid(LR))) qqnorm(resid(LR),main="") qqline(resid(LR)) #hist(rstudent(LR),xlab="Residuals",main="") ### Figure 9.4 par(mfrow=c(2,2)) plot(fitted(LR)[Smoke==0],rstudnonsmoke,xlab="Fitted Values", ylab="Studentized Residuals",main="Nonsmokers \n Residual Plot") lines(lowess(fitted(LR)[Smoke==0],rstudnonsmoke)) qqnorm(rstudnonsmoke,main="Nonsmokers \n Normal Plot") qqline(rstudnonsmoke) plot(fitted(LR)[Smoke==1],rstudsmoke,xlab="Fitted Values", ylab="Studentized Residuals",main="Smokers \n Residual Plot") lines(lowess(fitted(LR)[Smoke==1],rstudsmoke)) qqnorm(rstudsmoke,main="Smokers \n Normal Plot") qqline(rstudsmoke) ## Plot of studentized residuals versus age plot(Age, rstudent(LR),xlab="Age", ylab="Residuals",main="") lines(lowess(Age,rstudent(LR))) ## Quadratic trend in age Age2 <- Age^2 LR1 <- lm(FEV~ Age + Age2 + Smoke + Age:Smoke + Age2:Smoke) summary(LR1) rstudsmoke1 <- rstudent(LR1)[Smoke==1] rstudnonsmoke1 <- rstudent(LR1)[Smoke==0] ### Figure 9.5 par(mfrow=c(1,2)) plot(fitted(LR1)[Smoke==0],rstudnonsmoke1,xlab="Fitted Values", ylab="Studentized Residuals",main="Nonsmokers \n Residual Plot") lines(lowess(fitted(LR1)[Smoke==0],rstudnonsmoke1)) #qqnorm(rstudnonsmoke1,main="Nonsmokers \n Normal Plot") #qqline(rstudnonsmoke1) plot(fitted(LR1)[Smoke==1],rstudsmoke1,xlab="Fitted Values", ylab="Studentized Residuals",main="Smokers \n Residual Plot") lines(lowess(fitted(LR1)[Smoke==1],rstudsmoke1)) #qqnorm(rstudsmoke1,main="Smokers \n Normal Plot") #qqline(rstudsmoke1) ## Exercise 9.19. Using the model for FEV with predictors A, A^2(1-S), S, SA. This is model 6 from section 9.7. A2NS <- Age*Age*(1-Smoke) AS <- Age*Smoke betahat <- c(-2.656, 0.733, -0.02, 5.046, -0.666) Y <- FEV n <- length(Y) p <- length(betahat) X <- matrix(c(rep(1,n),Age,A2NS,Smoke,AS),n,p) Yhat <- X %*% betahat H <- X %*% solve(t(X) %*% X) %*% t(X) rr <- Y - Yhat Hdiag <- diag(H) V <- diag(Hdiag) I <- diag(1, n, n) MSE <- c((t(Y-Yhat) %*% (Y-Yhat))/ (n-p)) r <- (solve(sqrt(MSE*(I-V))) %*% rr) D <- V %*% solve(I-V) %*% diag(c(r)) %*% r max(r) r==max(r) # case 132. max(r)=2.989. person is 12 yo and has FEV that is 342 highest out 345 max(D) D==max(D) # case 1. max(D)=0.71 #cbind(r,D) m6 <- lm(FEV~Age+A2NS+Smoke+AS) rstudent(m6) cooks.distance(m6) cbind(r,rstudent(m6)) # These are a little different because R uses least squares and I used the posterior mean with proper # reference prior cbind(D,cooks.distance(m6)) ## Have to divide D by p to get the same formula of Cook's D used by R # Is case 132 data fit well by the model? Answer. No R132 <- sqrt((n-p-1)*r[132]^2 / (n-p-r[132]^2)) alpha <- 0.05 abs(R132) > qt(1-alpha/2,n-p-1) # Calculate Ri for all i R <- sqrt((n-p-1)*r^2 / (n-p-r^2)) # Calculate the observed Rf and show that it is equal to Ri for case 132 m <- lm(FEV[-132]~Age[-132]+A2NS[-132]+Smoke[-132]+AS[-132]) b <- c(-2.47453, 0.70366, -0.01881, 4.86902, -0.63722) # Least squares estimate for beta: remove case 132 b <- c(-2.46, 0.7014, -0.01873, 4.854, -0.635) # Posterior mean for beta: remove case 132 MSE1 <- sum((FEV[-132]-X[-132,] %*% b)^2)/(n-1-p) Rf <- (FEV[132] - X[132,] %*% b) / sqrt(MSE1 * (1 + X[132,] %*% solve(t(X[-132,]) %*% X[-132,]) %*% X[132,])) R132 ## Close to equal to Rf Rf Posterior mean for beta: complete case analysis (-2.656, 0.733, -0.02, 5.046, -0.666) Posterior mean for beta: remove case 1 (highest Cook's D at 0.71; next highest is 0.38); a 19 yo NS with high FEV (ranked 341 out of 345). (-3.501,0.8727,-0.02558,5.896,-0.8063) ######## Section 9.7: Model Selection # Calculate LPML and DIC n <- length(FEV) parameters <- c("beta","tau","CPOinv") ## Model 1: (Smoke); Diffuse priors inits <- list(list(tau=1,beta=c(0,0))) FEVdataBUGS <- list("n","FEV","Smoke") FEV.fit1 <- bugs(FEVdataBUGS, inits, parameters,"Model1.txt", working.directory="H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=12000, n.thin=1, n.burnin=2000) CPO1 <- 1/FEV.fit1\$mean\$CPOinv LPML1 <- sum(log(CPO1)) print(FEV.fit1,digits=3) ## Model 2: (Age); Diffuse priors inits <- list(list(tau=1,beta=c(0,0))) FEVdataBUGS <- list("n","FEV","Age") FEV.fit2 <- bugs(FEVdataBUGS, inits, parameters,"Model2.txt", working.directory="H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=12000, n.thin=1, n.burnin=2000) CPO2 <- 1/FEV.fit2\$mean\$CPOinv LPML2 <- sum(log(CPO2)) print(FEV.fit2,digits=3) ## Model 3: (Age,Smoke); Diffuse priors inits <- list(list(tau=1,beta=c(0,0,0))) FEVdataBUGS <- list("n","FEV","Age","Smoke") FEV.fit3 <- bugs(FEVdataBUGS, inits, parameters,"Model3.txt", working.directory="H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=12000, n.thin=1, n.burnin=2000) CPO3 <- 1/FEV.fit3\$mean\$CPOinv LPML3 <- sum(log(CPO3)) print(FEV.fit3,digits=3) ## Model 4: (Age,Smoke,Age*Smoke); Diffuse priors inits <- list(list(tau=1,beta=c(0,0,0,0))) FEVdataBUGS <- list("n","FEV","Age","Smoke") FEV.fit4 <- bugs(FEVdataBUGS, inits, parameters,"Model4.txt", working.directory="H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=12000, n.thin=1, n.burnin=2000) CPO4 <- 1/FEV.fit4\$mean\$CPOinv LPML4 <- sum(log(CPO4)) print(FEV.fit4,digits=3) ## Model 5: (Age,Age2,Smoke,Age*Smoke,Age2*Smoke); Diffuse priors inits <- list(list(tau=1,beta=c(0,0,0,0,0,0))) FEVdataBUGS <- list("n","FEV","Age","Smoke") FEV.fit5 <- bugs(FEVdataBUGS, inits, parameters,"Model5.txt", working.directory="H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=12000, n.thin=1, n.burnin=2000) CPO5 <- 1/FEV.fit5\$mean\$CPOinv LPML5 <- sum(log(CPO5)) print(FEV.fit5,digits=3) ## Model 6: (Age,Smoke,Age*Smoke,Age2*(1-Smoke)); Diffuse priors inits <- list(list(tau=1,beta=c(0,0,0,0,0))) FEVdataBUGS <- list("n","FEV","Age","Smoke") FEV.fit6 <- bugs(FEVdataBUGS, inits, parameters,"Model6.txt", working.directory="H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData", n.chains=1, n.iter=12000, n.thin=1, n.burnin=2000) CPO6 <- 1/FEV.fit6\$mean\$CPOinv LPML6 <- sum(log(CPO6)) print(FEV.fit6,digits=3) ## Bayes factors PBF21 <- exp(LPML2-LPML1) PBF31 <- exp(LPML3-LPML1) PBF41 <- exp(LPML4-LPML1) PBF51 <- exp(LPML5-LPML1) PBF61 <- exp(LPML6-LPML1) PBF32 <- exp(LPML3-LPML2) PBF42 <- exp(LPML4-LPML2) PBF52 <- exp(LPML5-LPML2) PBF62 <- exp(LPML6-LPML2) PBF43 <- exp(LPML4-LPML3) PBF53 <- exp(LPML5-LPML3) PBF63 <- exp(LPML6-LPML3) PBF54 <- exp(LPML5-LPML4) PBF64 <- exp(LPML6-LPML4) PBF65 <- exp(LPML6-LPML5) c(PBF21,PBF31,PBF41,PBF51,PBF61,PBF32,PBF42,PBF52,PBF62,PBF43,PBF53,PBF63,PBF54,PBF64,PBF65) ## BIC statistics n <- length(FEV) m1 <- lm(FEV~Smoke) AIC(m1,k=2) AIC(m1,k=log(n)) #Gives BIC anova.lm(m1) sigma2 <- .57 tau <- 1/sigma2 r <- 2 BIC1 <- n*log(2*pi*sigma2) + (1/sigma2)*sum((FEV-(m1\$coef[1]+m1\$coef[2]*Smoke))^2) + (r+1)*log(n) BIC1a <- -n*log(tau)+n*log(2*pi) + tau*sum((FEV-(m1\$coef[1]+m1\$coef[2]*Smoke))^2) + (r+1)*log(n) BIC1 BIC1a m2 <- lm(FEV~Age) AIC(m2,k=log(n)) #Gives BIC m3 <- lm(FEV~Age+Smoke) AIC(m3,k=log(n)) m4 <- lm(FEV~Age+Smoke+Age*Smoke) AIC(m4,k=log(n)) anova.lm(m4) sigma2 <- .441 r <- 4 BIC4 <- n*log(2*pi*sigma2) + (1/sigma2)*sum((FEV-(m4\$coef[1]+m4\$coef[2]*Age+m4\$coef[3]*Smoke+m4\$coef[4]*Age*Smoke))^2) + (r+1)*log(n) BIC4 Age2 <- Age^2 AS <- Age*Smoke A2S <- Age*Age*Smoke A2NS <- Age*Age*(1-Smoke) m5 <- lm(FEV~Age+Age2+Smoke+AS+A2S) AIC(m5,k=log(n)) m6 <- lm(FEV~Age+A2NS+Smoke+AS) AIC(m6,k=log(n)) #### Section 9.8: Nonlinear regression # Dugong data x <- c(1.0, 1.5, 1.5, 1.5, 2.5, 4.0, 5.0, 5.0, 7.0, 8.0, 8.5, 9.0, 9.5, 9.5, 10.0, 12.0, 12.0, 13.0, 13.0, 14.5, 15.5, 15.5, 16.5, 17.0, 22.5, 29.0, 31.5) xstd <- (x-mean(x))/sd(x) y <- c(1.80, 1.85, 1.87, 1.77, 2.02, 2.27, 2.15, 2.26, 2.47, 2.19, 2.26, 2.40, 2.39, 2.41, 2.50, 2.32, 2.32, 2.43, 2.47, 2.56, 2.65, 2.47, 2.64, 2.56, 2.70, 2.72, 2.57) dugongdata <- data.frame(xstd,y) n <- 27 ## This function creates the vector mtilde (equation 3 in the Section 9.8) mt <- function(b1,b2,b3){ m <- matrix(c(b1-exp(b2),b1-exp(b2+b3),b1-exp(b2-b3)),3,1) return(m) } ## This function creates the matrix mtilde dot (equation 4 in Section 9.8) mtdot <- function(b2,b3){ mtd <- matrix(c(1,1,1,-exp(b2),-exp(b2+b3),-exp(b2-b3),0,-exp(b2+b3),exp(b2-b3)),3,3) return(mtd) } Xt <- matrix(c(1,1,1,0,1,-1),3,2) m0t <- c(2.4,2.6,2.0) delta <- 0.001 #b0old <- c(0,1,-1) #b0old <- c(2.65,-1.54,-1.08) # MLEs b0old <- c(10,10,-10) b0new <- b0old + solve(mtdot(b0old[2],b0old[3])) %*% (m0t-mt(b0old[1],b0old[2],b0old[3])) while(sum(abs(b0new-b0old))>delta){ b0old <- b0new b0new <- b0old + solve(mtdot(b0old[2],b0old[3])) %*% (m0t-mt(b0old[1],b0old[2],b0old[3])) } b0new # Answer is c(2.8000000, -0.9162907, -0.6931472) beta0 <- c(2.8000000, -0.9162907, -0.6931472) D <- diag(c(25,25,25)) cov0 <- solve(mtdot(beta0[2],beta0[3])) %*% D %*% t(solve(mtdot(beta0[2],beta0[3]))) prec0 <- solve(cov0) xstar <- seq(-1.3,2.6,0.1) K <- length(xstar) a <- 0.001 b <- 0.001 library(R2WinBUGS) DugongDataBUGS <- list("n","beta0","prec0","a","b","y","xstd","xstar","K") parameters <- c("beta","tau","thetastar") inits <- list(list(tau=1,beta=c(0,0,0))) Dugong.fit <- bugs(DugongDataBUGS, inits, parameters,"DugongWBcode.txt", working.directory = "H:\\MyDocuments\\BAYES\\LinearRegression\\CodeandData\\NonlinearRegression", n.chains=1, n.iter=60000, n.thin=1, n.burnin=10000) print(Dugong.fit,digits=3) attach.bugs(Dugong.fit) ## Figure 9.6. Plot estimated growth curve. Posterior mean and 95% pointwise band for theta Postm <- Dugong.fit\$mean\$thetastar Postl <- apply(thetastar,2,quantile,0.025) Postu <- apply(thetastar,2,quantile,0.975) plot(xstd,y,xlab="Standardized Age", ylab="Length", lwd=2) lines(xstar,Postm,type='l', lty=1, lwd=2,cex=1.2) lines(xstar,Postl,lty=2,lwd=2) lines(xstar,Postu,lty=2,lwd=2)