# load data
data(wesdr)

# fit linear spline with full interactions and obtain diagnostics
wesdr.fit0 <- gssanova(ret~dur*bmi*gly,data=wesdr,
                       family="binomial",type="linear")
sum.fit0 <- summary(wesdr.fit0,diag=TRUE)
round(sum.fit0$pi,2)
round(sum.fit0$cos,2)

# fit additive cubic spline and obtain diagnostics
wesdr.fit <- gssanova(ret~dur+bmi+gly,data=wesdr,family="binomial")
sum.fit <- summary(wesdr.fit,diag=TRUE)
round(sum.fit$kappa,2)
round(sum.fit$pi,2)
round(sum.fit$cos,2)

# evaluate the fit and plot
est.dur <- predict(wesdr.fit,wesdr,se=TRUE,inc="dur")
est.gly <- predict(wesdr.fit,wesdr,se=TRUE,inc="gly")
est.bmi <- predict(wesdr.fit,wesdr,se=TRUE,inc="bmi")

# plot the fit
par(mar=c(5,5,2,2)+.1,mfrow=c(1,3),mex=.6,pty="s")

ind <- order(wesdr$dur)
plot(wesdr$dur[ind],est.dur$fit[ind],type="l",ylim=c(-4,4),cex.lab=1.2,cex.axis=1,
xlab="duration",ylab="")
lines(wesdr$dur[ind],(est.dur$fit+1.96*est.dur$se)[ind],lty=2)
lines(wesdr$dur[ind],(est.dur$fit-1.96*est.dur$se)[ind],lty=2)
rug(wesdr$dur,ticksize=.05,lwd=.005)

ind <- order(wesdr$gly)
plot(wesdr$gly[ind],est.gly$fit[ind],type="l",ylim=c(-4,4),cex.lab=1.2,cex.axis=1,
xlab="gly. hemoglobin",ylab="")
lines(wesdr$gly[ind],(est.gly$fit+1.96*est.gly$se)[ind],lty=2)
lines(wesdr$gly[ind],(est.gly$fit-1.96*est.gly$se)[ind],lty=2)
rug(wesdr$gly,ticksize=.05,lwd=.005)

ind <- order(wesdr$bmi)
plot(wesdr$bmi[ind],est.bmi$fit[ind],type="l",ylim=c(-4,4),cex.lab=1.2,cex.axis=1,
xlab="body mass index",ylab="")
lines(wesdr$bmi[ind],(est.bmi$fit+1.96*est.bmi$se)[ind],lty=2)
lines(wesdr$bmi[ind],(est.bmi$fit-1.96*est.bmi$se)[ind],lty=2)
rug(wesdr$bmi,ticksize=.05,lwd=.005)