barsVSgss.R

library(BARS)

library(gss)

x <- c(1:110)

mu <- 20*abs(exp(-x/100)*cos(x/5))

## replicate example of "bars"

set.seed(20061001,"Mersenne-Twister")

y <- rnorm(110,mu,1)

system.time(out <- bars(x,y,prior="uniform",priorparam=c(1,15),fits=TRUE))

system.time(gssfit1 <- ssanova0(y~x))

pred.gss <- predict(gssfit1,newdata=data.frame(x=x),se=TRUE)

upper <- apply(out$sampfits,2,quantile,0.975)

lower <- apply(out$sampfits,2,quantile,0.025)

upr <- pred.gss$fit+1.96*pred.gss$se.fit

lwr <- pred.gss$fit-1.96*pred.gss$se.fit

bars.out <- sum(mu < lower | upper < mu)

gss.out <- sum(mu < lwr | upr < mu)

layout(matrix(1:2,nrow=1))

m <- min(c(lower,lwr))

M <- max(c(upper,upr))

plot(x,out$postmeans,

type="n",xlab="x",ylab="y",

main=paste("BARS result (",bars.out,"out of",length(x),")"),

ylim=c(m,M))

polygon(c(x,rev(x)),c(lower,rev(upper)),col="grey30",lty=0)

lines(x,mu,col=2,lwd=2)

plot(x,pred.gss$fit,type="n",

xlab="x",ylab="y",

main=paste("GSS result (",gss.out,"out of",length(x),")"),

ylim=c(m,M))

polygon(c(x,rev(x)),c(lwr,rev(upr)),col="grey30",lty=0)

lines(x,mu,col=2,lwd=2)

## replicate example of "barsP"

set.seed(20061001,"Mersenne-Twister")

y <- rpois(110,mu)

system.time(out <- barsP(x,y,prior="uniform",priorparam=c(1,15),fits=TRUE))

system.time(gssfit1 <- gssanova0(y~x,family="poisson"))

pred.gss <- predict(gssfit1,newdata=data.frame(x=x),se=TRUE)

upper <- apply(out$sampfits,2,quantile,0.975)

lower <- apply(out$sampfits,2,quantile,0.025)

upr <- exp(pred.gss$fit+1.96*pred.gss$se.fit)

lwr <- exp(pred.gss$fit-1.96*pred.gss$se.fit)

bars.out <- sum(mu < lower | upper < mu)

gss.out <- sum(mu < lwr | upr < mu)

X11();layout(matrix(1:2,nrow=1))

m <- min(c(lower,lwr))

M <- max(c(upper,upr))

plot(x,out$postmeans,

type="n",xlab="x",ylab="y",

main=paste("BARS result (",bars.out,"out of",length(x),")"),

ylim=c(m,M))

polygon(c(x,rev(x)),c(lower,rev(upper)),col="grey30",lty=0)

lines(x,mu,col=2,lwd=2)

plot(x,pred.gss$fit,type="n",

xlab="x",ylab="y",

main=paste("GSS result (",gss.out,"out of",length(x),")"),

ylim=c(m,M))

polygon(c(x,rev(x)),c(lwr,rev(upr)),col="grey30",lty=0)

lines(x,mu,col=2,lwd=2)