Ninf.ssize <- function(freq, K, alpha, beta, delta){
#------------------------------------------------------------------------------------
# Examples of execution
#  
#    ldata <- c(0.05, 0.02, 0.05,  ### (qN3,  qN2,  qN1,
#               0.05, 0.02, 0.05)  ###  qS3,  qS2,  qS1)                            
#    out <- Ninf.ssize(ldata, 3, 0.025, 0.2, 0.1)
#
#> out
#[[1]]
#[[1]]$Sample.size
#[1] 119.2904
# (R version 2.13.1)
#
#    ldata <- c(0.05, 0.05,  ### (qN2,  qN1,
#               0.05, 0.05)  ###  qS2,  qS1)                            
#    out2 <- Ninf.ssize(ldata, 2, 0.025, 0.2, 0.1)
#
#    Result of the out2 is omitted here. 
#------------------------------------------------------------------------------------
# PURPOSE
#     The approximate sample size required for a one-sided non-inferiority score test 
#     for the difference in correlated proportions on the basis of the expected values 
#     from multiple raters are estimated using this function.  
#
# ARUGUMENTS
#     freq  : a vector of expected proportions (qNk and qSk)
#     K     : a value of the number of raters
#     alpha : a value of the level of significance at one-sided
#     beta  : a value of the 1-power
#     delta : a value of non-inferiority margin
#
# OUTPUT
#     Sample.size : estimated sample size
#
# NOTE
#   This function is available when the number of rater is two or three.
#   The generalized function adapted for any number of raters will be available near 
#   future. 
#
#------------------------------------------------------------------------------------

#### K=2
Ninf.ssize2 <- function(freq, K, alpha, beta, delta){

ipn2 <- freq[1]
ipn1 <- freq[2]
ips2 <- freq[3]
ips1 <- freq[4]

ipdiag <- 1-ipn2-ipn1-ips2-ips1

#### Difference of probabilities
dif <- ((ipn2+ipn1/2)-(ips2+ips1/2))

#### Estimating the MLEs based on null hypothesis by the quasi-Newton BFGS method with constraints 
psolv <- NULL

fr <- function(x) {   #### The log-likelihood function
    py2  <- x[1]
    px1  <- x[2]
    py1  <- x[3]
    ipn2*log(py2+0.5*py1-0.5*px1-delta)+ips2*log(py2)+ipn1*log(px1)+ips1*log(py1)+(1-ipn2-ips2-ipn1-ips1)*log(1-2*py2-0.5*px1-1.5*py1+delta)
}

grr <- function(x) {  #### The gradients of log-likelihood function
    py2  <- x[1]
    px1  <- x[2]
    py1  <- x[3]
    c((ipn2/(py2+0.5*py1-0.5*px1-delta)+ips2/py2-2*(1-ipn2-ips2-ipn1-ips1)/(1-2*py2-0.5*px1-1.5*py1+delta)),
       ((-1*ipn2)/2/(py2+0.5*py1-0.5*px1-delta)+ipn1/px1-(1-ipn2-ips2-ipn1-ips1)/2/(1-2*py2-0.5*px1-1.5*py1+delta)),
       (ipn2/2/(py2+0.5*py1-0.5*px1-delta)+ips1/py1-3*(1-ipn2-ips2-ipn1-ips1)/2/(1-2*py2-0.5*px1-1.5*py1+delta)))
}

#### constrOptim function
solv <- constrOptim(c((ipn2+ips2)/2+delta, (ipn1+ips1)/2, (ipn1+ips1)/2), fr, grr, ui=rbind(c(1, 0, 0),c(0, 1, 0),c(0, 0, 1),c(1, -0.5, 0.5), c(-2, -0.5, -1.5)), ci=c(0, 0, 0, delta, -1-delta), control=list(fnscale=-1), method='BFGS')
psolv <- solv[[1]]

py2  <- psolv[1]
px1  <- psolv[2]
py1  <- psolv[3]

####Calculation of the sample size (eq.(23))
Rval <- sqrt(1/4*(8*py2+3*py1-px1-4*delta-4*delta**2))
Sval  <- sqrt(1/4*(8*ips2+3*ips1-ipn1+4*dif-4*dif**2))

zalpha <- qnorm(1-alpha)
zbeta <- qnorm(1-beta)

ssize <- ((Rval*zalpha+Sval*zbeta)/(dif+delta))**2
ssize
list(Sample.size=ssize)
}

#### K=3
Ninf.ssize3 <- function(freq, K, alpha, beta, delta){
ipn3 <- freq[1]
ipn2 <- freq[2]
ipn1 <- freq[3]
ips3 <- freq[4]
ips2 <- freq[5]
ips1 <- freq[6]

ipdiag <- 1-ipn3-ipn2-ipn1-ips3-ips2-ips1

#### Difference of probabilities
dif <- ((ipn3+ipn2*2/3+ipn1/3)-(ips3+ips2*2/3+ips1/3))

#### Estimating the MLEs based on null hypothesis by the quasi-Newton BFGS method with constraints 
psolv <- NULL

fr <- function(x) {   #### The log-likelihood function
    py3  <- x[1]
    px2  <- x[2]
    py2  <- x[3]
    px1  <- x[4]
    py1  <- x[5]
   ipn3*log(py3+2/3*py2+1/3*py1-2/3*px2-1/3*px1-delta)+ips3*log(py3)+ipn2*log(px2)+ips2*log(py2)+ipn1*log(px1)+ips1*log(py1)+(1-ipn3-ips3-ipn2-ips2-ipn1-ips1)*log(1-2*py3-5/3*py2-4/3*py1-1/3*px2-2/3*px1+delta)
}

grr <- function(x) {  #### The gradients of log-likelihood function
    py3  <- x[1]
    px2  <- x[2]
    py2  <- x[3]
    px1  <- x[4]
    py1  <- x[5]
    c(
       (ipn3/(py3+2/3*py2+1/3*py1-2/3*px2-1/3*px1-delta)+ips3/py3-2*(1-ipn3-ips3-ipn2-ips2-ipn1-ips1)/(1-2*py3-5/3*py2-4/3*py1-1/3*px2-2/3*px1+delta)),  # gradient of py3
       (-2*ipn3/(3*(py3+2/3*py2+1/3*py1-2/3*px2-1/3*px1-delta))+ipn2/px2-(1-ipn3-ips3-ipn2-ips2-ipn1-ips1)/(3*(1-2*py3-5/3*py2-4/3*py1-1/3*px2-2/3*px1+delta))),  # gradient of px2
       (2*ipn3/(3*(py3+2/3*py2+1/3*py1-2/3*px2-1/3*px1-delta))+ips2/py2-5*(1-ipn3-ips3-ipn2-ips2-ipn1-ips1)/(3*(1-2*py3-5/3*py2-4/3*py1-1/3*px2-2/3*px1+delta))),  # gradient of py2
       (-ipn3/(3*(py3+2/3*py2+1/3*py1-2/3*px2-1/3*px1-delta))+ipn1/px1-2*(1-ipn3-ips3-ipn2-ips2-ipn1-ips1)/(3*(1-2*py3-5/3*py2-4/3*py1-1/3*px2-2/3*px1+delta))),  # gradient of px1
       (ipn3/(3*(py3+2/3*py2+1/3*py1-2/3*px2-1/3*px1-delta))+ips1/py1-4*(1-ipn3-ips3-ipn2-ips2-ipn1-ips1)/(3*(1-2*py3-5/3*py2-4/3*py1-1/3*px2-2/3*px1+delta)))   # gradient of py1
      )
}


#### constrOptim function
solv <- constrOptim(c((ipn3+ips3)/2+delta, (ipn2+ips2)/2, (ipn2+ips2)/2, (ipn1+ips1)/2, (ipn1+ips1)/2), fr, grr, 
ui=rbind(c(1, 0, 0, 0, 0), c(0, 1, 0, 0, 0), c(0, 0, 1, 0, 0), c(0, 0, 0, 1, 0), c(0, 0, 0, 0, 1),  
c(1, -2/3, 2/3, -1/3, 1/3), c(-2, -1/3, -5/3, -2/3, -4/3)), 
ci=c(0, 0, 0, 0, 0, delta, -1-delta), control=list(fnscale=-1), method='BFGS')
psolv <- solv[[1]]

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

####Calculation of the sample size (eq.(23))
Rval <- sqrt((18*py3+10*py2+4*py1-2*px1-2*px2-9*delta-9*delta**2)/9)
Sval  <- sqrt((18*ips3+10*ips2+4*ips1-2*ipn1-2*ipn2+9*dif-9*dif**2)/9)

zalpha <- qnorm(1-alpha)
zbeta <- qnorm(1-beta)

ssize <- ((Rval*zalpha+Sval*zbeta)/(dif+delta))**2
ssize
list(Sample.size=ssize)
}

#### Coding for the selection of the functions
if (K==2){
res <- Ninf.ssize2(freq, K, alpha, beta, delta)
}else{
res <- Ninf.ssize3(freq, K, alpha, beta, delta)
}
list(res)
}