#################################################################################
# Bayesian model selection and statistical modeling 
# Chapman & Hall/CRC Taylor and Francis Group 
#
# Chapter7: Section7.4.2
#
# P-spline logistic regression modeling with MBIC
#
# Author : Tomohiro Ando
#
################################################################################
library(splines)
library(rpart)

# B-Splines design matrix

bspline <- function(x,xleft,xright,ndx,bdeg){
	dx <- (xright - xleft)/ndx
	knots <- seq(xleft - bdeg * dx , xright + bdeg * dx , by = dx)
	B <- spline.des(knots , x , bdeg +1 , 0 * x)$design
B
}


# Data

x <- kyphosis$Age
y <- as.numeric(kyphosis$Kyphosis)-1
n <- length(y)

#Settings

m <- 13

#lambda <- 10
#lambda <- 0.1
#lambda <- 0.001
lambda <- 0.00001

xleft  <- min(x)-0.0001
xright <- max(x)+0.0001
ndx    <- p-3
bdeg   <- 3

B <- bspline(x,xleft,xright,ndx,bdeg)

Dk <- diff(diag(1,m+3),differences=2)
K  <- t(Dk)%*%Dk

Beta <- runif(m+3,-0.0001,0.001)
for(k in 1:100){
Beta.old <- Beta
O <- B%*%Beta
Pi <- as.vector(exp(O)/(exp(O)+1))
W <- diag(Pi*(1-Pi))
Eta <- (y-Pi)/(Pi*(1-Pi))+O
Beta <- solve(t(B)%*%W%*%B+n*lambda*K)%*%t(B)%*%W%*%Eta
if(sum(abs(Beta-Beta.old))<=10^(-8)){break}
}

O <- B%*%Beta
Pi <- exp(O)/(exp(O)+1)
H <- B%*%solve(t(B)%*%W%*%B+n*lambda*K)%*%t(B)%*%W
MBIC <- 2*sum(log(exp(O)+1))-2*sum(y*O)+log(n)*sum(diag(H))
print(MBIC)


#Plot

par(cex.lab=1.2)
par(cex.axis=1.2)
plot(x,y,xlab="Age (month)",ylab="Kyphosis")
a <- cbind(x,Pi); a <- a[order(a[,1]),1:2]
lines(a[,1],a[,2],lwd=3)


#Plot

pairs(kyphosis)