### Code for Section 4.1: calculation of unbiased and bias-adjusted estimators
### for the MUSEC trial reported by Bauer et al. and Zajicek et al.

################################################################################

library(gsDesign)
library(mvtnorm)
library(fda.usc)

### Interim data

n1_1 = 101   # Number of subjects randomized to CE arm
n0_1 = 97    # Number of subjects randomized to placebo

s1_1 = 27    # Number of subjects with relief from muscle stiffness in CE arm
s0_1 = 12    # Number of subjects with relief from muscle stiffness in placebo

ptilde_1 = (s0_1 + s1_1)/(n0_1 + n1_1)
Ihat_1 = 1/(ptilde_1*(1-ptilde_1)*(1/n0_1 + 1/n1_1))
Z1 = (s1_1/n1_1 - s0_1/n0_1)*sqrt(Ihat_1)   # Wald test statistic

### Final data

n1_2 = 143   # Number of subjects randomized to CE arm
n0_2 = 134   # Number of subjects randomized to placebo

s1_2 = 42    # Number of subjects with relief from muscle stiffness in CE arm
s0_2 = 21    # Number of subjects with relief from muscle stiffness in placebo

ptilde_2 = (s0_2 + s1_2)/(n0_2 + n1_2)
Ihat_2 = 1/(ptilde_2*(1-ptilde_2)*(1/n0_2 + 1/n1_2))
Z2 = (s1_2/n1_2 - s0_2/n0_2)*sqrt(Ihat_2)   # Wald test statistic

### O'Brien-Fleming stopping boundaries

original = gsDesign(k=2, test.type = 1, beta = 0.2, sfu='OF', n.fix = 400)


################################################################################
### Calculate estimators

# Naive (MLE)
MLE = s1_2/n1_2 - s0_2/n0_2
MLE_se = sqrt(1/Ihat_2)

# First stage MLE
MLE1 = s1_1/n1_1 - s0_1/n0_1
MLE1_se = sqrt(1/Ihat_1)

# Second stage MLE
MLE2 = (s1_2 - s1_1)/(n1_2 - n1_1) - (s0_2 - s0_1)/(n0_2 - n0_1)

ptilde2_star = ((s1_2-s1_1)+(s0_2-s0_1))/(n1_2+n0_2-n1_1-n0_1)
MLE2_se = sqrt(ptilde2_star*(1-ptilde2_star)*(1/(n1_2-n1_1)+1/(n0_2-n0_1)))

# MUE (unconditional)

cov_matrix = matrix(c(1, sqrt(Ihat_1/Ihat_2), sqrt(Ihat_1/Ihat_2), 1),
                    nrow = 2, byrow = TRUE)

e1 = original$upper$bound[1]
f1 = -Inf

pval_fn = function(delta) {
  
  mu = delta*sqrt(c(Ihat_1,Ihat_2))
  
  pnorm(delta*sqrt(Ihat_1) - e1) + 
    pmvnorm(lower=c(f1, Z2), upper=c(e1, Inf), mean = mu, sigma = cov_matrix)[1]
}

MUEobj = function(p){
  L = (pval_fn(delta=p) - 0.5)^2
}

MUE = optimize(MUEobj,interval = c(-1,1))$minimum


# Bias-corrected MLE (unconditional)

Ihat_2star = Ihat_2 - Ihat_1

bias = function(delta) {
  
  (Ihat_2star/(Ihat_2*sqrt(Ihat_1)))*(dnorm(e1 - delta*sqrt(Ihat_1)) - 
                                        dnorm(f1 - delta*sqrt(Ihat_2)))
}

ubcMLE_obj = function(delta){
  L = (delta + bias(delta) - MLE)^2
}

ubcMLE = optimize(ubcMLE_obj, interval = c(-1,1))$minimum


# UMVUE (unconditional)

UMVUE = MLE - sqrt(Ihat_2star/(Ihat_1*Ihat_2))*
  (dnorm((e1-sqrt(Ihat_1/Ihat_2)*Z2)/sqrt(Ihat_2star/Ihat_2)) - 
     dnorm((f1-sqrt(Ihat_1/Ihat_2)*Z2)/sqrt(Ihat_2star/Ihat_2)))/
  (pnorm((e1-sqrt(Ihat_1/Ihat_2)*Z2)/sqrt(Ihat_2star/Ihat_2)) - 
     pnorm((f1-sqrt(Ihat_1/Ihat_2)*Z2)/sqrt(Ihat_2star/Ihat_2)))


# Bias-corrected MLE (conditional)

bias_c = function(delta) {
  
  (1/sqrt(Ihat_1))*(dnorm(delta*sqrt(Ihat_1)-f1)-dnorm(delta*sqrt(Ihat_1)-e1))/
    (pnorm(e1 - delta*sqrt(Ihat_1)) - pnorm(f1 - delta*sqrt(Ihat_1)))
}


cbcMLE_obj = function(delta){
  L = (delta + bias_c(delta)*Ihat_1/Ihat_2 - MLE)^2
}

cbcMLE = optimize(cbcMLE_obj, interval = c(-1,1))$minimum


# UMVCUE (conditional)

Ihat_2star = Ihat_2 - Ihat_1

w1 = (Ihat_2star)^(-1)/sqrt(Ihat_1^(-1) + (Ihat_2star)^(-1))
w2 = sqrt(Ihat_1^(-1) + (Ihat_2star)^(-1))/(Ihat_1^(-1))

UMVCUE = MLE - w1*(dnorm(w2*(MLE - f1/sqrt(Ihat_1))) - 
                     dnorm(w2*(MLE - e1/sqrt(Ihat_1))))/
  (pnorm(w2*(MLE - f1/sqrt(Ihat_1))) - 
     pnorm(w2*(MLE - e1/sqrt(Ihat_1))))


# MUE (conditional)

# cMUE_fn = function(delta) {
#   
#   mu = delta*sqrt(c(Ihat_1,Ihat_2))
#   
#   log_S2 = log(pnorm(e1, mu[1]) - pnorm(f1, mu[1]))
#   
#   z_int   = seq(Z2, 10, length.out = 10000)
#   
#   dfz = -0.5*(z_int - delta*sqrt(Ihat_2))^2 +
#     log((pnorm(e1, z_int*sqrt(Ihat_1/Ihat_2), sqrt((Ihat_2 - Ihat_1)/Ihat_2)) -
#            pnorm(f1, z_int*sqrt(Ihat_1/Ihat_2), sqrt((Ihat_2 - Ihat_1)/Ihat_2)))) - 
#     log(sqrt(2*pi))
#   
#   cdfz    = exp(dfz - log_S2)
#   int     = int.simpson2(z_int, cdfz)
# }
# 
# cMUE_obj = function(p){
#   L = (cMUE_fn(delta=p) - 0.5)^2
# }
# 
# cMUE = optimize(cMUE_obj,interval = c(-1,1))$minimum

f2                  <- function(z, theta, I1, I2, l, u, log = FALSE) {
  if (!log) {
    exp(-0.5*(z - theta*sqrt(I2))^2)*
      (pnorm(u, z*sqrt(I1/I2), sqrt((I2 - I1)/I2)) -
         pnorm(l, z*sqrt(I1/I2), sqrt((I2 - I1)/I2)))/sqrt(2*pi)
  } else {
    if (is.infinite(u)) {
      -0.5*(z - theta*sqrt(I2))^2 +
        pnorm((z*sqrt(I1/I2) - l)/sqrt((I2 - I1)/I2), log.p = TRUE) -
        log(sqrt(2*pi))
    } else {
      -0.5*(z - theta*sqrt(I2))^2 +
        log((pnorm(u, z*sqrt(I1/I2), sqrt((I2 - I1)/I2)) -
               pnorm(l, z*sqrt(I1/I2), sqrt((I2 - I1)/I2)))) - log(sqrt(2*pi))
    }
  }
}

cmue_optim          <- function(theta, k, z, I1, I2, l, u) {
  if (k == 1) {
    S1       <- pnorm(l, theta*sqrt(I1), 1) + pnorm(u, theta*sqrt(I1), 1,
                                                    lower.tail = FALSE)
    if (z >= u) {
      num    <- pnorm(l, theta*sqrt(I1)) + pnorm(z, theta*sqrt(I1)) -
        pnorm(u, theta*sqrt(I1))
    } else {
      num    <- pnorm(z, theta*sqrt(I1))
    }
    int      <- num/S1
  } else {
    if (is.infinite(u)) {
      log_S2 <- pnorm(theta*sqrt(I1) - l, log.p = TRUE)
    } else {
      log_S2 <- log(pnorm(u, theta*sqrt(I1)) - pnorm(l, theta*sqrt(I1)))
    }
    z_int    <- seq(z, qnorm(1 - 0.0001, z), length.out = 1000)
    cf2z     <- exp(f2(z_int, theta, I1, I2, l, u, log = TRUE) - log_S2)
    int      <- int.simpson2(z_int, cf2z)
  }
  (int - 0.5)^2
}


cMUE = optim(MLE2, cmue_optim, k = 2, z = Z2,
             I1 = Ihat_1, I2 = Ihat_2, l = f1, u = e1,
             method = "Nelder-Mead")$par

### Print estimates for Table 3

results = matrix(c(MLE, MLE1, MLE2, MUE, ubcMLE, UMVUE, cMUE, cbcMLE, UMVCUE))

rownames(results) = c('MLE', 'MLE1', 'MLE2', 'MUE', 'Bias-corrected MLE',
                      'UMVUE', 'MUE (conditional)',
                      'Bias-corrected MLE (conditional)', 'UMVCUE')

print(results)

relative_diff = results/MLE
print(relative_diff)