############################################ # Spline basis by hand & using splines package ############################################ p=function(x){ # mother B-spline le=length(x) f=rep(0,le) for(i in 1:le){ if(x[i]>=0){f[i]=0.5*x[i]^2} if(x[i]>=1){f[i]=0.75-(x[i]-1.5)^2} if(x[i]>=2){f[i]=0.5*(3-x[i])^2} if(x[i]>3){f[i]=0} } f } b=function(x,m,d,j){p((x-m)/d+3-j)} xl=0; xr=10; J=5; d=(xr-xl)/(J-2) x=seq(xl-2*d,xr+2*d,0.1) plot(x,b(x,xl,d,1),type="l") lines(x,b(x,xl,d,2)) lines(x,b(x,xl,d,3)) lines(x,b(x,xl,d,4)) lines(x,b(x,xl,d,5)) library(splines) # included w/ R installation xm=bs(x,knots=seq(xl-d,xr+d,d),degree=2,Boundary.knots=c(xl-2*d,xr+2*d))[,2:(J+1)] points(x,xm[,1]) points(x,xm[,2]) points(x,xm[,3]) points(x,xm[,4]) points(x,xm[,5]) ############################################ # Ethanol data, NOx vs. E only # Unpenalized vs. penalized likelihood # for different lambda ############################################ library(lattice) data(ethanol) attach(ethanol) plot(ethanol) a=min(E); b=max(E) x=seq(a,b,length=100) n=length(E) fit=loess(NOx~E) pred=predict(fit,x,se=TRUE) plot(x,pred$fit,type="l",xlab="equivalence ratio",ylab="NOx",main="Lowess Fit") lines(x,pred$fit-1.96*pred$se.fit,lty=3) lines(x,pred$fit+1.96*pred$se.fit,lty=3) points(E,NOx) library(splines) # included w/ R installation J=20 d=(b-a)/(J-2) X=bs(E,knots=seq(a-d,b+d,d),degree=2, Boundary.knots=c(a-2*d,b+2*d))[,2:(J+1)] Xp=bs(x,knots=seq(a-d,b+d,d),degree=2, Boundary.knots=c(a-2*d,b+2*d))[,2:(J+1)] fit=lm(NOx~X-1) # fit model, create 95% CI’s lo=x; hi=x for(i in 1:100){ lo[i]=Xp[i,]%*%fit$coef-qt(0.975,n-20)*Xp[i,]%*%summary(fit)$cov%*%Xp[i,] hi[i]=Xp[i,]%*%fit$coef+qt(0.975,n-20)*Xp[i,]%*%summary(fit)$cov%*%Xp[i,] } # fit B-spline model, display fitted regression line and 95% pointwise CI’s par(mfrow=c(1,1)) plot(E,NOx,main="Unpenalized B-spline, J=20",ylim=c(-1,5)) lines(x,Xp%*%fit$coef) lines(x,lo,lty=3) lines(x,hi,lty=3) ######################## # Penalized B-spline fit ######################## D=matrix(0,nrow=J-2,ncol=J) for(i in 1:(J-2)){D[i,i]=1; D[i,i+1]=-2; D[i,i+2]=1} par(mfrow=c(1,1)) lambda=1000 # neg. penalized log-likelihood function penlike=function(xi){ sum((NOx-X%*%xi)^2)+lambda*sum((d*D%*%xi)^2) } f=optim(fit$coef,penlike,hessian=T) cov=solve(f$hessian) for(i in 1:100){ lo[i]=Xp[i,]%*%f$par-qnorm(0.975)*Xp[i,]%*%cov%*%Xp[i,] hi[i]=Xp[i,]%*%f$par+qnorm(0.975)*Xp[i,]%*%cov%*%Xp[i,] } # display fitted penalized regression line and 95% pointwise CI’s plot(E,NOx,main=bquote("J=20,"~lambda==.(lambda)),ylim=c(-1,5)) lines(x,Xp%*%f$par) lines(x,lo,lty=3) lines(x,hi,lty=3) ############################################ # Ache hunting via JAGS, additive ############################################ library(R2jags) library(splines) # included w/ R installation d=read.table("http://people.stat.sc.edu/hansont/stat740/ache.txt",header=T) attach(d) # makes 'age', 'days', and 'monkeys' available to R n=length(age) id=factor(1:47) J=20 al=min(age); ar=max(age); d=(ar-al)/(J-2) xbs=bs(age,knots=seq(al-d,ar+d,d),degree=2, Boundary.knots=c(al-2*d,ar+2*d))[,2:(J+1)] D=matrix(0,nrow=J-2,ncol=J) for(i in 1:(J-2)){D[i,i]=1; D[i,i+1]=-2; D[i,i+2]=1} M=t(D)%*%solve(D%*%t(D)) # starting values for JAGS start=list(list("b"=c(0,0),"tau"=1,"u"=rep(0,n),"e"=rep(0,J-2))) model=function(){ # model for JAGS for(i in 1:n){ monkeys[i]~dpois(mu[i]) mu[i]<-days[i]*exp(inprod(xbs[i,],xi[])+u[i]) u[i]~dnorm(0,tau) } xi<-b[1]+b[2]*(1:J-J/2)+M%*%e for(j in 1:(J-2)){e[j]~dnorm(0,tau)} for(i in 1:2){b[i]~dnorm(0,0.0001)} tau~dgamma(0.0001,0.0001) } kills=list("monkeys","xbs","days","n","M","J") # data for JAGS param=c("b","xi","tau") f=jags(data=kills,inits=start,parameters.to.save=param, n.chains=1,n.iter=52000,n.burnin=2000,model.file=model, n.thin=5) print(f) traceplot(f) ############################################ # Ache hunting via JAGS, additive # kill rate vs. age for median-skilled hunter ############################################ grid=seq(al,ar,length=100) # age grid xrate=bs(grid,knots=seq(al-d,ar+d,d),degree=2, Boundary.knots=c(al-2*d,ar+2*d))[,2:(J+1)] model=function(){ # model for JAGS for(i in 1:n){ monkeys[i]~dpois(mu[i]) mu[i]<-days[i]*exp(inprod(xbs[i,],xi[])+u[i]) u[i]~dnorm(0,tau) } # this next line is new for(i in 1:100){rate[i]<-inprod(xrate[i,],xi[])} xi<-b[1]+b[2]*(1:J-J/2)+M%*%e for(j in 1:(J-2)){e[j]~dnorm(0,tau)} for(i in 1:2){b[i]~dnorm(0,0.0001)} tau~dgamma(0.0001,0.0001) } kills=list("monkeys","xbs","days","n","M","J","xrate") # data for JAGS param=c("rate") f=jags(data=kills,inits=start,parameters.to.save=param, n.chains=1,n.iter=52000,n.burnin=2000,model.file=model, n.thin=5) f2=print(f) # get posterior summary statistics med=f2$summary[2:101,5] low=f2$summary[2:101,3] upp=f2$summary[2:101,7] plot(grid,med,type="l",lty=1,ylim=c(min(low),max(upp)), xlab="Age (years)",ylab="log kill rate",main="Monkey Hunting") polygon(c(rev(grid),grid),c(rev(low),upp),col='grey80',border=NA) lines(grid,low,lty=2) lines(grid,upp,lty=2) lines(grid,med,lty=1) ############################################ # Ache hunting via BayesX, additive # sx(age) is penalized canonical B-spline ############################################ library(R2BayesX) f=bayesx(monkeys~sx(age)+sx(id,bs="re"),family="poisson",offset=log(days),method="REML") summary(f) plot(f) plot(f,term="sx(age)") bayesx.term.options(bs="ps",method="REML") # BayesX also can do variable selection choosing linear and/or smooth terms, # spatial terms, etc. Also method="MCMC" works but is buggy still. ############################################ # Ethanol via BayesX, additive ############################################ f=bayesx(NOx~sx(E)+sx(C),method="REML") summary(f) plot(f) plot(f,term="sx(E)") plot(f,term="sx(C)") f=bayesx(NOx~sx(E)+sx(C),method="MCMC") summary(f) plot(f) plot(f,term="sx(E)") plot(f,term="sx(C)") ############################################ # Orings via BayesX, additive ############################################ library(DPpackage) data(orings) attach(orings) f=bayesx(ThermalDistress~sx(Temperature)+sx(Pressure),method="REML") summary(f) plot(f) plot(f,term="sx(Temperature)") plot(f,term="sx(Pressure)") f=bayesx(ThermalDistress~sx(Temperature)+sx(Pressure),method="MCMC") summary(f) plot(f) plot(f,term="sx(Temperature)") plot(f,term="sx(Pressure)")