# generate data and fit the model
set.seed(5732)
x1 <- rnorm(100); x2 <- rnorm(100) 
y <- 1 + 3 * sin(2*pi*(x1-x2)-pi) + rnorm(x1) 
tpcubic.fit <- ssanova(y~x1*x2,method="m")

# evaluate the fit
range1 <- max(x1) - min(x1)
grid1 <- seq(min(x1)-.05*range1,max(x1)+.05*range1,length=25)
range2 <- max(x2) - min(x2)
grid2 <- seq(min(x2)-.05*range2,max(x2)+.05*range2,length=25)
new <- data.frame(x1=rep(grid1,25),x2=rep(grid2,rep(25,25)))
est <- predict(tpcubic.fit,newdata=new,se.fit=TRUE)
post.mean <- matrix(est$fit,25,25)
post.stdev <- matrix(est$se.fit,25,25)

# contour the posterior mean
filled.contour(grid1,grid2,post.mean,sub="REML Fit",
plot.axes={points(x1,x2);axis(1);axis(2)},asp=1)

# contour the posterior standard deviation
filled.contour(grid1,grid2,post.stdev,sub="Standard Deviation",
plot.axes={points(x1,x2);axis(1);axis(2)},asp=1)