###############################################################################################################
# --- Power Calculation for Case-Cohort Studies with Nonrare Events
# --- Function: Power_CaseCohort_Nonrare <-function(alpha, n, pD, ratio, HR, pC, p0, q0)
#               alpha------Test size 
#               n    ------Full cohort size 
#               pD   ------Observed failure rate
#               ratio------Ratio of sizes of exposure group to control in the population
#               HR   ------Hazard ratio of exposure group to control to be detected 
#               pC   ------Premature drop-out rate (Input "NA" in case unspecified)
#               p0   ------Subcohort sampling probability (Input "NA" in case unspecified)
#               q0   ------Failure sampling probability (Input "NA" in case unspecified)
# --- Reference:  Cai, J. and Zeng, D. (2007) Power Calculation for Case-Cohort Studies with Nonrare Events
#                 BIOMETRICS 10, 1111-1118
###############################################################################################################

Power_CaseCohort_Nonrare <-function(alpha, n, pD, ratio, HR, pC, p0, q0)
{
 D=round(n*pD, digits = 0) # D : Total number of failure
 beta<-1                   # Assume same censoring patterns in two groups
 gamma<-ratio/(1+ratio)    # gamma = Proportion of the first group (exposure group)
 theta<-log(HR)            # HR = Hazard ratio of exposure group to control to be detected
 Z_a<-qnorm(alpha)

 if (is.na(pC))            #(II) When pC is not specified by input
 {
  pC_value<-seq(0,1,0.1)    # Produce powers by pC (0.1~1 by 0.1) attention!
  nC=n*pC_value             # nC : Total number of dropouts   pC_value='Premature drop-out rate      (pC)'

  if (is.na(p0))
  {
   p<-seq(0.1,1,0.1)          #p='Subcohort sampling probability (p)'

   if (is.na(q0))              # (4) Produce powers by p (0.1~1 by 0.1) and q (0.1~1 by 0.1)
   {
    q<-seq(0.1,1,0.1)           #q='Failure sampling probability (q)'
   }
   else                     # (3) Input q and produce powers by p (0.1~1 by 0.1)
   {
    q<-q0
   }
   
  }
  else
  {
   p<-p0

   if (is.na(q0))              # (2) Input p and produce powers by q (0.1~1 by 0.1)
   {
    q<-seq(0.1,1,0.1)
   }
   else
   {
    q<-q0                     # (1) Input p and q and produce a power
   }
   
  }
  
 }
 else       ##(I) When pC is specified by input
                  #(I) When pC is specified by input
 {
  pC_value=pC
  nC=n*pC_value             #  nC : Total number of dropouts

  if (is.na(p0))
  {
   p<-seq(0.1,1,0.1)

   if (is.na(q0))              # (4) Produce powers by p (0.1~1 by 0.1) and q (0.1~1 by 0.1)
   {
    q<-seq(0.1,1,0.1)
   }
   else               # (3) Input q and produce powers by p (0.1~1 by 0.1)
   {
    q=q0
   }

  }
  else
  {
   p=p0

   if (is.na(q0))              # (2) Input p and produce powers by q (0.1~1 by 0.1)
   {
    q<-seq(0.1,1,0.1)
   }
   else
   {
    q=q0                      # (1) Input p and q and produce a power
   }

  }
 
 }

 power_data<-expand.grid(pC_value,q,p)
 power_data<-data.frame(beta,gamma,theta,Z_a,nC,power_data)
 colnames(power_data)<-c("beta","gamma","theta","Z_a","nC","pC_value","q","p")
 err_ind=0
 
 if (is.na(pC))
 {
  power_data_new<-subset(power_data,nC<(n-D))
  m_n<-matrix(n,nrow(power_data_new),1)
  m_D<-matrix(D,nrow(power_data_new),1)
  nk_L<-m_n-power_data_new[5]-m_D
  NoErrSet<-data.frame(power_data_new,nk_L)
  colnames(NoErrSet)<-c("beta","gamma","theta","Z_a","nC","pC_value","q","p","nk_L")
 }
 else
 {
  nk_L = n-nC-D
  if ((nk_L<=0) & (!is.na(pC))) 
  {
   err_nk<-"Lower bound of size of risk sets is not positive for the given" 
   err_ind<-1
   ErrSet<-data.frame(power_data,nk_L,err_nk,err_ind)
  }
  if (nk_L > 0)
  {
   NoErrSet<-data.frame(power_data,nk_L)
   }
 }



#/**** (1) With error ****/
if (err_ind == 1)
{
print("ERROR",quote=F)
print("Full cohort size (n)")
print(n)
print("Failure proportion (pD)")
print(pD)
print("Total # of failure (D)")
print(D)
print("Ratio of exposure to control (ratio)")
print(ratio)
print("Hazard Ratio (HR)")
print(HR)
print("ERROR MESSAGE")
print(err_nk)
print("Premature drop-out rate (pC)")
print(pC_value)
print("The number of premature drop-out(nC)")
print(nC)
ErrSet
}else #/**** (2) Without error ****/
{
n_power=nrow(NoErrSet)
k=seq(0,D,1)
id<-seq(1,n_power,1)
for (i in 1:n_power)
{
nk_L<-NoErrSet[i,9]
nk_inv<-matrix(1/nk_L,D+1,D)
sum1<-matrix(0,D+1,D)
for (j in 1:D)
{
nk_inv[(j+1):(D+1),j]<-1/(n-j+1)
if (j==1)
{
sum1[,j]<-nk_inv[,j]
}else
{
sum1[,j] = sum1[,(j-1)]+nk_inv[,j]
}
}
sum1_sq<-(1-sum1)*(1-sum1)
J1 = matrix(1,D,1)
if (i==1)
{
sum2<-sum1_sq%*%J1
}
else
{
sum2<-c(sum2,sum1_sq%*%J1)
}
}
id_k_sum2<-expand.grid(k,id)
id_k_sum2<-data.frame(id_k_sum2,sum2)
sum_all<-id_k_sum2
sum_all1<-sum_all
sum_all[,1]<-sum_all1[,2]
sum_all[,2]<-sum_all1[,1]
colnames(sum_all)<-c("id","k","sum2")
}

if (err_ind == 0)
{

a<-dim(NoErrSet)
id<-seq(1,a[[1]],1)
NoErrSet<-data.frame(NoErrSet,id)

Chi_all<-merge(NoErrSet,sum_all,by.NoErrSet="id",by.sum_all="id")
for (i in 1:nrow(NoErrSet))
{
temp<-pD/NoErrSet[i,8]-( (1-NoErrSet[i,8])/NoErrSet[i,8] - (1-NoErrSet[i,8])*(1-NoErrSet[i,7])/NoErrSet[i,7])/n*sum2[((i-1)*(D+1)+1):((D+1)*i)]
for (j in 1:(D+1))
{
if (temp[j]<=0)
{
temp[j]<-0
}
temp_min<-min(temp)
temp_max<-max(temp)
}
if (i==1)
{
Chi<-temp
Chi_L<-temp_min
Chi_U<-temp_max
}
else
{
Chi<-c(Chi,temp)
Chi_L<-c(Chi_L,temp_min)
Chi_U<-c(Chi_U,temp_max)
}
}

Chi_all<-data.frame(Chi_all,Chi)

Chi_LU<-data.frame(NoErrSet,Chi_L,Chi_U)

numer<-sqrt(n*beta*gamma*(1-gamma))*theta*pD
for (i in 1:nrow(Chi_LU))
{
if (i==1)
{
if (Chi_LU[i,12]==0)
{
Z_L<-Z_a + numer/sqrt(Chi_LU[i,12]+1E-20)
}
else
{
Z_L<-Z_a + numer/sqrt(Chi_LU[i,12])
}
if (Chi_LU[i,11]==0)
{
Z_U<-Z_a + numer/sqrt(Chi_LU[i,11]+1E-20)
}
else
{
Z_U<-Z_a + numer/sqrt(Chi_LU[i,11])
}
}else
{
if (Chi_LU[i,12]==0)
{
Z_L<-c(Z_L,Z_a + numer/sqrt(Chi_LU[i,12]+1E-20))
}
else
{
Z_L<-c(Z_L,Z_a + numer/sqrt(Chi_LU[i,12]))
}
if (Chi_LU[i,11]==0)
{
Z_U<-c(Z_U,Z_a + numer/sqrt(Chi_LU[i,11]+1E-20))
}
else
{
Z_U<-c(Z_U,Z_a + numer/sqrt(Chi_LU[i,11]))
}
}
}
Power_L=pnorm(Z_L)
Power_U=pnorm(Z_U)

for (i in 1:nrow(Chi_LU))
{
p_v<-Chi_LU[i,8]
if (i==1)
{
den_Rare<-p_v+(1-p_v)*pD
if (den_Rare==0)
{
Z_Rare<-Z_a + sqrt(n*p_v)*theta*sqrt( beta*gamma*(1-gamma)*pD/(p_v+(1-p_v)*pD+1E-20)) 
}
else
{
Z_Rare<-Z_a + sqrt(n*p_v)*theta*sqrt( beta*gamma*(1-gamma)*pD/(p_v+(1-p_v)*pD) )
}
}else
{
den_Rare<-c(den_Rare,p_v+(1-p_v)*pD)
if (den_Rare[i]==0)
{
Z_Rare<-c(Z_Rare,Z_a + sqrt(n*p_v)*theta*sqrt( beta*gamma*(1-gamma)*pD/(p_v+(1-p_v)*pD+1E-20))) 
}
else
{
Z_Rare<-c(Z_Rare,Z_a + sqrt(n*p_v)*theta*sqrt( beta*gamma*(1-gamma)*pD/(p_v+(1-p_v)*pD)))
}
}
}
Power_Rare=pnorm(Z_Rare)

power_LU<-data.frame(Chi_LU,numer,Z_L,Z_U,Power_L,Power_U,den_Rare,Z_Rare,Power_Rare) #label Power_L='Lower Bound Power' Power_U='Upper Bound Power' Power_Rare='Rare Disease Power'


print("<Power for Case-Cohort design with Nonrare Events>")
print(alpha)
print("Full cohort size (n)")
print(n)
print("Failure proportion (pD)")
print(pD)
print("Total # of failure (D)")
print(D)
print("Ratio of exposure to control (ratio)")
print(ratio)
print("Hazard Ratio (HR)")
print(HR)
power_LU_final<-data.frame(power_LU[,8],power_LU[,7],power_LU[,6],power_LU[,16],power_LU[,17],power_LU[,20])
colnames(power_LU_final)<-c("Subcohort Sampling Probability (p)","Failure Sampling Probability (q)","Premature drop-out rate (pC)","Lower Bound Power","Upper Bound Power","Rare Disease Power")
power_LU_final

}
}