# Program 15
# Program to find IMSE-optimal design for  s > 2 support points for logistic link and MPL estimates
# The program requires input of the vector of parameters, betavec, the number of support points, s, and
# the vector of number of observations at the support points, nvec.  You need also to supply the values
# of qlo and qhi (the values of pi(x) that determine the values of x between which the integral is calculated;
# see Figure 4.8).
# Lastly, the number of subintervals used in Simpson's rule ("nsubintervals") must be defined.

#The program is currently set up for betavec = c(0,1), s = 4, nvec = c(2,2,2,2), qlo = 0.0001, qhi = 0.9999 
#and nsubintervals = 100

betavec <- c(0,1)
s <- 4
nvec <- c(2,2,2,2)

qlo <- 0.0001
qhi <- 0.9999
xlo <- (log(qlo/(1-qlo)) - betavec[1])/betavec[2] #equals -9.21024
xhi <- (log(qhi/(1-qhi)) - betavec[1])/betavec[2] #equals 9.21024

nsubintervals <- 100
nsteps <- 2*nsubintervals
h <- (xhi - xlo)/nsteps
xvalues <- seq(from=xlo,to=xhi,length=(nsteps+1))
temp <- 2 + 2*(((1:(nsteps-1)) %% 2) == 1)
multiplier <- c(1,temp,1)

library(brglm)

fx <- function(x)
{
 c(1,x)
}

mse <- function(x,betastarmat,outcomeprobs)
{#Calculates the MSE at a value x for support points in xvec
 pix <- 1/(1+exp(-sum(fx(x)*betavec)))
 pistarxmat <- 1/(1+exp(-betastarmat%*%fx(x)))
 mse <- sum(((pistarxmat-pix)^2)*outcomeprobs)
 mse
}

#Calculates all the possible vectors of responses at the s support points 
# (these are the rows of the matrix 'ymat'), and then calculates
# the corresponding values of the vector of estimates: betastar
p <- length(betavec)
nvecp1 <- nvec + 1
c1 <- c(1,nvecp1[1:(s-1)])
c2 <- cumprod(c1)
totalrows <- prod(nvecp1)
ymat <- matrix(0,totalrows,s)
probmat <- matrix(0,totalrows,s)
for(i in 1:s)
{
 ymat[,i] <- rep(0:nvec[i],times=c2[i],each=totalrows/(nvecp1[i]*c2[i]))
}
ymat

imse <- function(zvec)
{
 xvec <- 9.21024*cos(pi*zvec)
 betastarmat <- matrix(0,totalrows,p)
 for (i in 1:totalrows)
 {
  out <- brglm(cbind(ymat[i,],nvec - ymat[i,])~xvec,family=binomial("logit"))
  betastarmat[i,] <- out$coeff
 }
 betastarmat
 fxmat <- sapply(xvec,fx)
 etavec <- t(betavec)%*%fxmat
 probvec <- 1/(1+exp(-etavec))
 probmat <- matrix(0,totalrows,s)
 for(i in 1:s)
 {
  probmat[,i] <- rep(dbinom(0:nvec[i],nvec[i],probvec[i]),times=c2[i],each=totalrows/(nvecp1[i]*c2[i]))
 }
 probmat
 outcomeprobs <- apply(probmat,1,prod)
 zvalues <- sapply(xvalues,mse,betastarmat=betastarmat,outcomeprobs=outcomeprobs)
 imse <- sum(multiplier*zvalues)*h/3
 imse
}

##############################################################
#Searching for an IMSE-optimal design when n1=n2=n3=n4=2

nsims <- 20
minimse <- 1.0e08
for (i in 1:nsims)
{
 initial <- runif(s)
 out <- optim(initial,imse,NULL,method="Nelder-Mead")
 if (out$value < minimse) {minimse <- out$value
  design <- out$par}
}
transformeddesign <- 9.21024*cos(pi*design)
cat("Minimum value of IMSE is",minimse,"\n")
cat("design is", transformeddesign,"\n")