# # # R CODE FOR REPRODUCING THE FIGURES AND ANALYSES IN JW'S # "BAYESIAN AND FREQUENTIST REGRESSION ANALYSIS" CHAPTER 11 # CODE WRITTEN BY JON WAKEFIELD, UNLESS OTHERWISE STATED # #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # # # Polynomials for the LIDAR data # library(SemiPar) data(lidar) attach(lidar) xvals <- seq(min(range),max(range),.1) lidar2 <- lm(logratio~poly(range,deg=2,raw=T)) lidar3 <- lm(logratio~poly(range,deg=3,raw=T)) lidar4 <- lm(logratio~poly(range,deg=4,raw=T)) lidar8 <- lm(logratio~poly(range,deg=8,raw=T)) # # Fig 11.1(a) # pdf("lidardeg2fig.pdf",height=4.5,width=4.5) plot(logratio~range,xlab="Range (m)",ylab="Log Ratio",col="grey") lines(xvals,lidar2$coeff[1]+lidar2$coeff[2]*xvals+lidar2$coeff[3]*xvals*xvals) dev.off() # # Fig 11.1(b) # pdf("lidardeg3fig.pdf",height=4.5,width=4.5) plot(logratio~range,xlab="Range (m)",ylab="Log Ratio",col="grey") lines(xvals,lidar3$coeff[1]+lidar3$coeff[2]*xvals+lidar3$coeff[3]*xvals*xvals+ lidar3$coeff[4]*xvals*xvals*xvals,lty=1) dev.off() # # Fig 11.1(c) # pdf("lidardeg4fig.pdf",height=4.5,width=4.5) plot(logratio~range,xlab="Range (m)",ylab="Log Ratio",col="grey") lines(xvals,lidar4$coeff[1]+lidar4$coeff[2]*xvals+lidar4$coeff[3]*xvals*xvals+ lidar4$coeff[4]*xvals*xvals*xvals+lidar4$coeff[5]*xvals*xvals*xvals*xvals, lty=1) dev.off() # # Fig 11.1(d) # pdf("lidardeg8fig.pdf",height=4.5,width=4.5) plot(logratio~range,xlab="Range (m)",ylab="Log Ratio",col="grey") lines(xvals,lidar8$coeff[1]+lidar8$coeff[2]*xvals+lidar8$coeff[3]*xvals^2+ lidar8$coeff[4]*xvals^3+lidar8$coeff[5]*xvals^4+lidar8$coeff[6]*xvals^5+ lidar8$coeff[7]*xvals^6+lidar8$coeff[8]*xvals^7+lidar8$coeff[9]*xvals^8, lty=1) dev.off() # # Illustration of splines with 2 knots for the lidar data # y <- logratio x <- range xi1 <- min(x) + (range(x)[2]-range(x)[1])/3 xi2 <- min(x) + 2*(range(x)[2]-range(x)[1])/3 # Horizontal lines xtruncflat1 <- ifelse(x>xi1,1,0) xtruncflat2 <- ifelse(x>xi2,1,0) pieceflat <- lm(y~1+xtruncflat1+xtruncflat2) # # Fig 11.2(a) # pdf("lidarpiece0.pdf",height=4.5,width=4.5) plot(x,y,ylim=c(min(y),max(y)),ylab="y",type="n") points(x,y,col="grey") segments(x0=min(x),x1=xi1,y0=pieceflat$coeff[1],y1=pieceflat$coeff[1]) segments(x0=xi1,x1=xi2,y0=pieceflat$coeff[1]+pieceflat$coeff[2], y1=pieceflat$coeff[1]+pieceflat$coeff[2]) segments(x0=xi2,x1=max(x),y0=pieceflat$coeff[1]+pieceflat$coeff[2]+ pieceflat$coeff[3],y1=pieceflat$coeff[1]+pieceflat$coeff[2]+ pieceflat$coeff[3]) abline(v=xi1,lty=2) abline(v=xi2,lty=2) dev.off() # xtrunclin1 <- ifelse(x>xi1,x-xi1,0) xtrunclin2 <- ifelse(x>xi2,x-xi2,0) piecelin <- lm(y~1+x+xtrunclin1+xtrunclin2) # # Fig 11.2(b) # pdf("lidarpiece1.pdf",height=4.5,width=4.5) plot(x,predict(piecelin),type="l",ylim=c(min(y),max(y)),ylab="y") points(x,y,col="grey") abline(v=xi1,lty=2) abline(v=xi2,lty=2) dev.off() # Now for quadratic x2 <- x*x xtruncquad1 <- ifelse(x>xi1,(x-xi1)^2,0) xtruncquad2 <- ifelse(x>xi2,(x-xi2)^2,0) piecequad <- lm(y~1+x+x2+xtruncquad1+xtruncquad2) # # Fig 11.2(c) # pdf("lidarpiece2.pdf",height=4.5,width=4.5) plot(x,predict(piecequad),type="l",ylim=c(min(y),max(y)),ylab="y") points(x,y,col="grey") abline(v=xi1,lty=2) abline(v=xi2,lty=2) dev.off() # Now for cubic x3 <- x2*x xtrunccub1 <- ifelse(x>xi1,(x-xi1)^3,0) xtrunccub2 <- ifelse(x>xi2,(x-xi2)^3,0) piececub <- lm(y~1+x+x2+x3+xtrunccub1+xtrunccub2) # # Fig 11.2(d) # pdf("lidarpiece3.pdf",height=4.5,width=4.5) plot(x,predict(piececub),type="l",ylim=c(min(y),max(y)),ylab="y") points(x,y,col="grey") abline(v=xi1,lty=2) abline(v=xi2,lty=2) dev.off() # # Piecewise splines # xi1 <- 1/3 xi2 <- 2/3 x <- y <- seq(0,1,.01) # # Fig 11.3 # pdf("piecewiselin.pdf",height=4.5,width=4.5) par(mfrow=c(1,1),pty="s") plot(x,y,type="l",lty=1,ylab="f(x)") lines(c(0,1),c(1,1),lty=1) lines(c(xi1,1),c(0,2*xi1),lty=2) lines(c(xi2,1),c(0,xi1),lty=2) text(.31,.07,expression(xi[1])) lines(x=c(.33,.33),y=c(-.01,.01)) text(.64,.07,expression(xi[2])) lines(x=c(.66,.66),y=c(-.01,.01)) dev.off() # # Basis functions for piecewise cubic # # Fig 11.4(a) # pdf("piecewisecuball1.pdf",height=4.5,width=4.5) par(mfrow=c(1,1),pty="s") plot(x,y,type="n",lty=1,ylab="f(x)") lines(c(0,1),c(1,1),lty=1) lines(x,y,lty=1) lines(seq(0,1,.01),seq(0,1,.01)^2,lty=1) lines(seq(0,1,.01),seq(0,1,.01)^3,lty=1) dev.off() beta3 <- 3.3 beta4 <- 12 x3 <- ifelse(x>xi1,(x-xi1)^3,0) x4 <- ifelse(x>xi2,(x-xi2)^3,0) # # Fig 11.4(b) # pdf("piecewisecuball2.pdf",height=4.5,width=4.5) par(mfrow=c(1,1),pty="s") plot(x,y,type="n",lty=1,ylab="f(x)") lines(x,x3*beta3,lty=2) text(.31,.07,expression(xi[1])) lines(x=c(.33,.33),y=c(-.01,.01)) text(.69,.07,expression(xi[2])) lines(x=c(.66,.66),y=c(-.01,.01)) lines(x,x4*beta4,lty=2) text(.31,.07,expression(xi[1])) lines(x=c(.33,.33),y=c(-.01,.01)) text(.69,.07,expression(xi[2])) lines(x=c(.66,.66),y=c(-.01,.01)) dev.off() # # Lidar with natural smoothing spline and CV # library(pspline) library(SemiPar) data(lidar) attach(lidar) dfval <- seq(4,20,.1) ocvval <- gcvval <- rocvval <- rgcvval <- NULL for (i in 1:length(dfval)){ origmod <- smooth.Pspline(x=range,y=logratio,method=2,df=dfval[i]) ocvval[i] <- origmod$cv gcvval[i] <- origmod$gcv } # # Fig 11.5 # pdf("lidarnatsmoothcvplot.pdf",height=4.5,width=4.5) plot(ocvval~dfval,type="l",lty=1,ylab="Cross-Validation Score", xlab="Effective Degrees of Freedom") lines(gcvval~dfval,lty=2) legend("topright",legend=c("Ordinary Cross-Validation", "Generalized Cross-Validation"),lty=1:2,bty="n") dev.off() dfminocvorig <- dfval[which.min(ocvval)] dfmingcvorig <- dfval[which.min(gcvval)] cat("Min DF for original and OCV = ",dfminocvorig,"\n") # 9.3 cat("Min DF for original and GCV = ",dfmingcvorig,"\n") # 9.4 # # Lidar with smoothing splines and mixed model # library(SemiPar) data(lidar) attach(lidar) knots.pos <- seq(400,700,length=20) fit <- spm(lidar$logratio~f(lidar$range,basis="trunc.poly", degree=3,knots=knots.pos)) summary(fit) par(mfrow=c(1,1)) # # Fig 11.6 # pdf("lidarmixed.pdf",height=6,width=7) plot(logratio~range,col="grey",xlab="Range (m)",ylab="Log Ratio") lines(sm.spline(x=range,y=logratio,df=dfmingcvorig),lty=1) lines(x=lidar$range,y=fit$fit$fitted,lty=2) legend("bottomleft",legend=c("GCV natural cubic spline","Mixed model cubic spline"),lty=1:2,bty="n") dev.off() # # B-splines figure # library(splines) x <- seq(0,1,.01) K <- 9 # K=no of knots, so basis is of dimension K+4 (inc intercept). knots <- seq(1,K)/(K+1) y <- bs(x,knots=knots,intercept=T) nbasis <- K+4 # # Fig 11.7 # pdf("Bsplinefig.pdf",h=6,w=4.5) plot(y[,1]~x,type="l",xlab="x",ylab="B-spline",ylim=c(0,max(y))) for (j in 2:nbasis){ lines(y[,j]~x,lty=j) } points(x=knots,y=rep(0,K)) dev.off() # # LIDAR example at the end of Section 11.2.7 # # # Tatiana Maravina's code for carrying out confidence band calculation. # # First get function that computes knots locations # source("http://matt-wand.utsacademics.info/webspr/default.knots.r") # # Fit using mixed models representations # for later comparison with my 'by hand' results # library(SemiPar) data(lidar) attach(lidar) fit=spm(logratio~f(range, basis="trunc.poly", degree=3)) plot(fit) points(range, logratio, pch=20) # # 'By hand' fitting of penalized splines and computation of CIs # Define some functions first # Tprimenorm <- function (x, knots, W, p) { # computes ||T'(x)|| - function to be integrated to get kappa0 # x can be a vector K <- length(knots) m <- length(x) S <- Sprime <- matrix(nrow=m, ncol=1+p+K) for (i in 1:m) { S[i,] <- c(x[i]^(0:p), ((x[i]-knots)^p)*(x[i]>knots)) Sprime[i,] <- c(0, 1:p*(x[i]^(0:(p-1))), (3*(x[i]-knots)^2)*(x[i]>knots)) } S <- S%*%W # m-by-n Sprime <- Sprime%*%W # m-by-n Snorm2 <- apply(S^2, 1, sum) numer <- sqrt(Snorm2*apply(Sprime^2, 1, sum)- apply((S*Sprime)^2,1,sum)) # m-by-1 return(numer/Snorm2) } # f <- function (x, kappa0, alpha) { # solution of f(x)=0 gives c 2*(1-pnorm(x))+kappa0*exp(-x^2/2)/pi-alpha } # ## data # x <- lidar$range y <- lidar$logratio ## Standardize x # Otherwise, I kept running into "system is computationally singular" problem # when inverting (t(X)%*%X + lambda*D) mu <- mean(x) sig <- sd(x) xs <- (x-mu)/sig # n <- length(y) K <- 35 # number of knots knots <- default.knots(xs, K) # knots p <- 3 # order (3 for cubic) X <- (matrix(nrow=n, ncol=K, data=xs)-matrix(nrow=n, ncol=K, data=knots, byrow=TRUE)) X <- (X^p)*(X>0) # use truncated power basis of order p Xmat <- cbind(model.matrix(~xs+I(xs^2)+I(xs^3)), X) # EDIT TO MATCH ORDER! D <- diag(c(rep(0,p+1), rep(1,K))) # ## Find lambda that minimizes GCV # lambda <- seq(0,1, by=0.01) N <- length(lambda) GCV <- df <- numeric(N) for (i in 1:N) { # following the Semiparametric Regression book and # http://www.uow.edu.au/~mwand/SPmanu.pdf # the multiplier is lambda^(2*p) W <- solve(t(Xmat)%*%Xmat+(lambda[i]^(2*p))*D)%*%t(Xmat) beta.hat <- W%*%y y.hat <- Xmat%*%beta.hat S <- Xmat%*%W Sii <- diag(S) df[i] <- sum(Sii) GCV[i] <- sum((y-y.hat)^2)*n/(n-df[i])^2 } plot(lambda, GCV, type="l", xlab=expression(lambda)) lambdaopt <- lambda[which.min(GCV)] # optimal lambda is 0.41 abline(v=lambdaopt, lty=2, col="darkgrey") # ## Fit model with optimal value of lambda # W <- solve(t(Xmat)%*%Xmat+(lambdaopt^(2*p))*D)%*%t(Xmat) beta.hat <- W%*%y y.hat <- Xmat%*%beta.hat ## Find kappa0 and c a <- min(xs) b <- max(xs) kappa0 <- integrate(Tprimenorm, a, b, knots=knots, W=W, p=p)$value kappa0 # 15.44 cc <- uniroot(f, c(0, 10), kappa0=19, alpha=0.05)$root # c=3.11 as in text ## Variance estimation # Clearly, the data are heteroscedastic resids <- y-y.hat z <- log(resids^2) plot(x, z, pch=20) # looks linear so I will simply use lm() instead of nonparametric methods # q=spm(z~f(x))$fit$fitted q <- lm(z~x) # evaluate approximate standard errors at the observed range values sigma.e <- sqrt(exp(predict(q))) S <- Xmat%*%W Snorm <- sqrt(apply(S^2, 1, sum)) se <- sigma.e*Snorm # what if assumed constant variance? v <- sum(diag(S)) # trace(S) vv <- sum(diag(t(S)%*%S)) # trace(t(S)%*%S) sigma.e.homosc <- sqrt(sum(resids^2)/(n-2*v-vv)) sigma.e.homosc # compare with fit$fit$sigma 0.083 vs 0.080 se.homosc <- sigma.e.homosc*Snorm # compare with predict(fit, newdata=lidar, se=TRUE)$se # # Confidence bands # # Fig 11.8(a) # pdf("lidarsim2.pdf",h=6,w=6) par(mfrow=c(1,1)) plot(x, y, pch=20, xlab="Range (m)", ylab="Log Ratio",col="grey") lines(x, y.hat, lwd=2) lines(x, y.hat-cc*se, lty=2) lines(x, y.hat+cc*se, lty=2) lines(x, y.hat-1.96*se, lty=3) lines(x, y.hat+1.96*se, lty=3) legend("bottomleft", legend=c("Pointwise", "Simultaneous"),lty=3:2, ,bty="n") dev.off() # # Fig 11.8(b) # pdf("lidarsim1.pdf",h=6,w=6) plot(x, y, pch=20, xlab="Range (m)", ylab="Log Ratio", col="grey") lines(x, y.hat, lwd=2) lines(x, y.hat-cc*se.homosc, lty=2) lines(x, y.hat+cc*se.homosc, lty=2) lines(x, y.hat-1.96*se.homosc, lty=3) lines(x, y.hat+1.96*se.homosc, lty=3) legend("bottomleft", legend=c("Pointwise", "Simultaneous"), lty=3:2, bty="n") dev.off() # # Mixed model LIDAR spline example at the end of Section 11.2.7 # knots.pos <- seq(400,700,length=20) fit <- spm(lidar$logratio~f(lidar$range,basis="trunc.poly", degree=3,knots=knots.pos)) xseq <- seq(400,700,1) cubfit <- fit$fit$coef$fixed[1] + fit$fit$coef$fixed[2]*xseq + fit$fit$coef$fixed[3]*xseq^2 +fit$fit$coef$fixed[4]*xseq^3 xbas <- matrix(0,nrow=20,ncol=301) # totfit <- cubfit+randomfit randomfit <- 0 rand <- matrix(0,nrow=20,ncol=301) figname <- c("lidars1.pdf","lidars2.pdf","lidars3.pdf","lidars4.pdf","lidars5.pdf","lidars6.pdf","lidars7.pdf","lidars8.pdf","lidars9.pdf","lidars10.pdf","lidars11.pdf","lidars12.pdf","lidars13.pdf","lidars14.pdf","lidars15.pdf","lidars16.pdf","lidars17.pdf","lidars18.pdf","lidars19.pdf","lidars20.pdf") caps <- c("1st","2nd","3rd","4th","5th","6th","7th","8th","9th","10th","11th","12th","13th","14th","15th","16th","17th","18th","19th","20th") for (i in 1:20){ xbas[i,] <- ifelse(xseq > knots.pos[i],(xseq-knots.pos[i])^3,0) rand[i,] <- xbas[i,]*fit$fit$coef$random[i] # Fig 11.9 -- create 20 different files pdf(figname[i],height=5,width=5) par(mfrow=c(1,1),pty="s") plot(logratio~range,col="grey",xlab="",ylab="",ylim=c(-1.5,1.5), main=caps[i],cex.main=2) lines(xseq,cubfit,lty=1) lines(xseq,rand[i,],lty=2) randomfit <- randomfit + xbas[i,]*fit$fit$coef$random[i] abline(v=knots.pos[i],lty=3) # totfit <- cubfit+randomfit # lines(xseq,totfit,lty=3) dev.off() } # #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # Putting prior on the effective degrees of freedom in Section 11.2.9 # and then analyzing the SBMD data -- Code written by Youyi Fong #-------------------------------------------------------------------- # Code to evaluate induced prior distribution on the degrees of freedom # DFfun <- function(m,n,Cmat,sigmae2,sigmab2){ Lambda <- diag(c(0,rep(sigmae2/sigmab2,m))) Totdf <- tr( solve(crossprod(Cmat) + Lambda) %*% crossprod(Cmat) ) } # create.matrix=function (n, m){ out = matrix(0,n*m,m) for (i in 1:m) out[(i-1)*n+1:n,i]=1 out } # # Simple one-way ANOVA, example at end of Section 11.2.9 # m <- 10 # Number of groups n <- 5 # Number of observations per group # Xmat <- matrix(1,nrow=m*n,ncol=1) Zmat <- create.matrix(n,m) tr <- function(M) sum(diag(M)) Cmat <- cbind(Xmat,Zmat) sigmae2 <- 1 nsim <- 10000 sigmab2 <- 1/rgamma(nsim,.5,.005) par(mfrow=c(1,1)) cuto <- 2.5 # # Fig 11.10(a) # pdf("DFPrior1.pdf",w=4,h=4) hist(sqrt(sigmab2[sqrt(sigmab2)