Ninf.test <- function(freq, K, alpha, delta){
#------------------------------------------------------------------------------------
# Examples of execution
#  
#    ldata <- c(14, 2, 1, 0, 
#                0, 0, 0, 0,
#                0, 0, 2, 0,
#                0, 0, 2, 0)
#    out <- Ninf.test(ldata, 3, 0.025, 0.1)
#
#> out
#[[1]]
#[[1]]$Score.Z.value
#         [,1]
#[1,] 1.753303
#
#[[1]]$Score.p.value
#           [,1]
#[1,] 0.03977503
#
#[[1]]$Wald.Z.value
#         [,1]
#[1,] 3.357821
#
#[[1]]$Wald.p.value
#             [,1]
#[1,] 0.0003927979
# (R version 2.13.1)
#
#    ldata <- c(9, 1, 0, 
#               3, 0, 0,
#               0, 1, 1)
#    out2 <- Ninf.test(ldata, 2, 0.025, 0.1)
#
#    Results of the out2 are omitted here.     
#------------------------------------------------------------------------------------
# PURPOSE
#     The non-inferiority score test and the non-inferiority Wald test for 
#     the difference in correlated proportions on the basis of the results from 
#     multiple raters are perfomred using this function.  
#
# ARUGUMENTS
#     freq  : a vector of elements in the (K+1) times (K+1) contingency table
#     K     : a value of the number of raters
#     alpha : a value of the level of significance at one-sided
#     delta : a value of non-inferiority margin
#
# OUTPUT
#     Score.Z.value : test statistic of the score test 
#     Score.p.value : p-value of the score test
#     Wald.Z.value  : test statistic of the Wald test
#     Wald.p.value  : p-value of the Wald test
#
#------------------------------------------------------------------------------------

n <- sum(freq)
new <- rep(K:0, each=(K+1))
std <- rep(K:0, times=(K+1))

dif <- new-std


#### Make each observed freqency
obsy <- NULL
obsx <- NULL

for (i in 1:K) {
  obsy <- append(obsy, assign(paste("ny", i, sep=""), t(freq)%*%ifelse(dif==(-(K-(K-i))), 1, 0)))
  obsx <- append(obsx, assign(paste("nx", i, sep=""), t(freq)%*%ifelse(dif==(K-(K-i)), 1, 0)))
  diag <- t(freq)%*%ifelse(dif==0, 1, 0)
}

opry <- (obsy)/n
oprx <- (obsx)/n

################################
pyv <- paste("py", K:1, sep="")
pyv

nyv <- paste("ny", K:1, sep="")
nyv

eq2 <- paste(paste(nyv, paste("log(", pyv,")", sep=""), sep="*"),collapse="+" )
eq2

pxv <- paste("px", (K-1):1, sep="")
pxv 

nxv <- paste("nx", (K-1):1, sep="")
nxv

eq3 <- paste(paste(nxv, paste("log(", pxv,")", sep=""), sep="*"),collapse="+" )
eq3

ak1 <- seq(1, 1/K, by = -1/K)
ak1
ak2 <- seq((K-1)/K, 1/K, by = -1/K)
ak2

tempA1 <- paste(paste(ak1, pyv ,sep="*"), collapse="+")
tempA1
tempA2 <- paste("(", paste(paste(ak2, pxv, sep="*"), collapse="+"), ")", sep="")
tempA2
capA <- paste(tempA1, tempA2, sep="-")
capA

nxK <- paste("nx", K, sep="")
nxK

eq1 <- paste(nxK, paste("log(", capA, "-delta)", sep=""), sep="*")
eq1

capB <- paste(pxv, collapse="+")
capB

capC <- paste(pyv, collapse="+")
capC

eq4n <- paste("(", paste("n", nxK, paste(nxv, collapse='-'), paste(nyv, collapse="-"), sep="-"), ")", sep="")
eq4p <- paste("log(1+delta", paste("(", capA, ")", sep=""), paste("(", capB, ")", sep=""), paste("(", capC, "))", sep=""), sep="-")
eq4 <- paste(eq4n, eq4p, sep="*")
eq4

ll <- paste(eq1, eq2, eq3, eq4, sep="+")
ll

llf <- evalq(parse(text=paste(ll, sep="")))

pv <- c(paste("py", 1:K, sep=""), paste("px", 1:(K-1), sep=""))
pv

drv <- NULL
for (i in 1:(2*K-1)) {
    drv <- append(drv, D(llf, pv[i]))
    }

#########################################
llh <- function(x){
    for (i in 1:K) {
    assign(paste("py", i, sep=""), x[i])
    }
    for (i in 1:(K-1)) {
    assign(paste("px", i, sep=""), x[i+K])
    }

 eval(parse(text=paste(ll, sep="")))
}


scr <- function(x){
    for (i in 1:K) {
    assign(paste("py", i, sep=""), x[i])
    }
    for (i in 1:(K-1)) {
    assign(paste("px", i, sep=""), x[i+K])
    }
    drv2 <- NULL
    for (i in 1:(2*K-1)) {
    drv2 <- append(drv2, eval(parse(text=drv[i])))
    }
return(drv2)
}

#### Estimating the MLEs using the quasi-Newton BFGS method with constraints 
fr <- function(x) {   #### The log-likelihood function
     return(llh(x))
}


grr <- function(x) {   #### The gradients of log-likelihood function
     return(scr(x))
}

#### Making a vector of initial values and matrices for constriction
inity <- NULL
initx <- NULL
for (i in 1:(K-1)) {
  inity <- append(inity, (opry[i]+oprx[i])/2+1e-12)
  initx <- append(initx, (opry[i]+oprx[i])/2+1e-12)
}
inityd <- (opry[K]+oprx[K])/2+delta*1.01+1e-12

initv <- c(inity, inityd, initx)
initv

um <- diag(2*K-1)
um

cy1 <- seq(1/K, K/K, by = 1/K)
cx1 <- seq(-1/K,(1-K)/K, by = -1/K)
cst1 <- c(cy1, cx1)
cst1

cy2 <- seq((-1-K)/K, (-K-K)/K, by = -1/K)
cx2 <- seq((1-K)/K,-1/K, by = 1/K)
cst2 <- c(cy2, cx2)
cst2

mui <- rbind(um, cst1, cst2)
mui

#### constrOptim function
solv <- constrOptim(initv, fr, grr, ui=mui, ci=c(numeric(2*K-1), delta, -1-delta), control=list(fnscale=-1), method='BFGS')
psolv <- solv[[1]]

#py1 <- psolv[1]
#py2  <- psolv[2]
#py3  <- psolv[3]
#px1  <- psolv[4]
#px2  <- psolv[5]

pymle <- psolv[1:K]
pymle
ppxmle <- psolv[(K+1):(2*K-1)]
ppxmle
wvec <- seq(1/K, K/K, by = 1/K)
wvec
wwvec <- wvec**2
wwvec
rxmle <- t(wvec)%*%pymle-t(wvec[1:(K-1)])%*%ppxmle-delta
rxmle
pxmle <- c(ppxmle, rxmle)
pxmle

#### The score test based on Z_ND (eq.(18))

d <- t(wvec)%*%oprx-t(wvec)%*%opry
v <- (t(wwvec)%*%pxmle+t(wwvec)%*%pymle-delta**2)/n
znd <- (d+delta)/sqrt(v)
pval_ZND <- (1-pnorm(znd))

#### The Wald test based on Zw (eq.(19))
s <- (t(wwvec)%*%oprx+t(wwvec)%*%opry-delta**2)/n
wald <- (d+delta)/sqrt(s)
pval_WALD <- (1-pnorm(wald))

list(Score.Z.value=znd, Score.p.value=pval_ZND, Wald.Z.value=wald, Wald.p.value=pval_WALD)

}