#Load the reshape package
require(reshape)

##################################################################################
####################### Multiple Intraclass Correlation ##########################
## Takes as arguments:
## Data   - the data set
## Yindex - The column index for the binary outcome
## Dindex - The column index for the first (lower) level clustering indices
## Hindex - The column index for the second (higher) level clustering indices
## alpha  - The significance level (defaulted at 0.05)
##
## Returns:
## Point estimates of the correlations at the first level B(A) and second level A
## Confidence intervals for each stage of correlation
##################################################################################
ICBA <- function(Data,Yindex, Dindex, Hindex,alpha=0.05){
  Y <- Data[,Yindex]
  D <- Data[,Dindex]
  H <- Data[,Hindex]
  N <- length(Y)
  p <- sum(Y) / N
  a <- length(unique(H))
  
  #Convert data format to counts
  bi <- rep(0,a)
  Ni <- rep(0,a)
  dat <- NULL #Empty storage space
  for(i in 1:a){
    Ysub <- Y[which(H==i)]
    Dsub <- as.vector(D[which(H==i)])
    temp <- as.data.frame(table(Ysub,Dsub))  
    temp.w <- reshape(temp,timevar="Ysub",idvar="Dsub",direction="wide")
    if(sum(Ysub)==0) temp.w <- cbind(temp.w,"Freq.1"=rep(0,nrow(temp.w)))
    #Add in the hospital index
    temp.w <- cbind((i*rep(1,nrow(temp.w))),temp.w,(temp.w[,2]+temp.w[,3]))
    dat <- rbind(dat,temp.w)
    
    #Obtain the bi's and Nis
    bi[i] <- nrow(temp.w)
    Ni[i] <- length(Ysub)
  }
  colnames(dat) <- c("Hosp","Doc","Freq0","Freq1","nij")
  
  B <- sum(bi)
  
  #Calculate S1 and initalize S2 and S3
  S1 <- sum(Y)
  S2 <- 0
  S3 <- 0
  
  #Calculate the sum of squared nij
  nij.sq <- dat[,5]**2 %*% rep(1,B)

  #Calculate sum of sum of squared nij/Ni
  #Also calculate vector of sums of nij(nij-1)
  nij.sq2 <- 0
  nnm1 <- NULL
  startsum <- 1
  for(i in 1:a){
    endsum <- sum(bi[seq(1:i)])
    nij.sq2 <- nij.sq2 + (sum(dat[(startsum:endsum),5]**2)/Ni[i])
    nnm1.temp <- sum(dat[(startsum:endsum),5]*(dat[(startsum:endsum),5]-1))
    nnm1 <- c(nnm1,nnm1.temp)
    
    #Calculate S2 and S3
    S2 <- S2 + (sum(dat[(startsum:endsum),4])**2)/Ni[i]
    S3 <- S3 + (sum((dat[(startsum:endsum),4]**2)/dat[(startsum:endsum),5]))
    
    startsum <- endsum + 1
  }
  
  #Sum of nij(nij-1)
  nnm1.sum <- sum(nnm1)
  
  #Sum of nij(nij-1) / Ni
  nnm1.sum2 <- sum(nnm1/Ni)
  
  #Calculate r1 and r2
  r1 <- (N-nij.sq2) / (B-a)
  r2 <- (nij.sq2 - (1/N)*nij.sq) / (a-1)
  r3 <- (N - (1/N)*sum(Ni**2))/(a-1)
  
  d1 <- (r3-r2) / r1
  d2 <- (r1*r3-r3-r1+r2) / r1
  
  #Calculate the correlation coefficients
  MSA <- (1/(a-1)) * (S2 - (1/N)*S1**2)
  MSB <- (1/(B-a)) * (S3 - S2)
  MSE <- (1/(N-B)) * (S1 - S3)
  
  #Calculate each of the sigma
  sigE <- MSE
  sigB <- (1/r1)*(MSB-MSE)
  sigA <- (1/r3)*(MSA - (r2/r1)*MSB - (1-r2/r1)*MSE)
  
  
  #Calculate rho(B(A))
  rhoB <- (sigB) / (sigA+sigB+sigE)
  if(is.na(rhoB)) rhoB = 0
  if(rhoB<0) rhoB = 0
  
  #Calculate rho(A)
  rhoA <- sigA / (sigA+sigB+sigE)
  if(is.na(rhoA)) rhoA=0
  if(rhoA<0) rhoA = 0
  
  #Calculate the variances
  rootfun <- function(x,y,corr){
    result <- sqrt(x*y*(1+corr*(x-1))*(1+corr*(y-1)))
    return(result)
  }
  
  #Calculate sum[sum[sqrt(nijnij'...)]]
  rootn <- 0
  #Calculate sum[1/Ni * sum[sqrt(nijnij'...)]]
  rootn2 <- 0
  startsum <- 1
  for(i in 1:a){
    endsum <- sum(bi[seq(1:i)])
    temp <- dat[(startsum:endsum),5]
    
    if(bi[i]>1){
      pairs <- combn((1:bi[i]), m=2, simplify=TRUE)
      root.out <- 0
      for(j in 1:(ncol(pairs))){
        root.out <- root.out + 2*rootfun(temp[pairs[1,j]],temp[pairs[2,j]],rhoB)
      }
      rootn <- rootn + root.out
      rootn2 <- rootn2 + root.out/Ni[i]
    }
      startsum <- endsum + 1
  }
  
  #Find the variance of S1
  VS1 <- N*p*(1-p) + rhoB*p*(1-p)*sum(dat[,5]*(dat[,5]-1)) + rhoA*p*(1-p)*rootn
  
  #Calculate the lambdas
  lambda1 <- -(1/(a-1))*N*p**2 + (d2/(N-B))*N*p + 
    (1/(a-1) - d1/(B-a))*(p**2) * (4*a + 4*p*(1-p)*rhoB*nnm1.sum2 + 4*p*rhoA*(1-p)*rootn2 - N) - 
    (d1/(B-a) - (d2/(N-B))) * ((1-p)*(B + N*rhoB - B*rhoB) + N*p**2)
  
  lambda2 <- -N*p/(N-B) - ((p**2)/(B-a)) * (4*a + 4*p*(1-p)*rhoB*nnm1.sum2 + 4*p*rhoA*(1-p)*rootn2 - N) +
    (1/(N-B) + 1/(B-a)) * (p*(1-p)*(B + N*rhoB - B*rhoB) + N*p**2)

  lambda3 <- -(N*p**2)/(a-1) + ((r2-r1)/(r1*(N-B)))*N*p + 
    (p**2)*(1/(a-1) + r2/(r1*(B-a))) * (4*a + 4*p*(1-p)*rhoB*nnm1.sum2 + 4*p*rhoA*(1-p)*rootn2 - N) -
    (r2/(r1*(B-a)) + (r2-r1)/(r1*(N-B))) * (p*(1-p)*(B + N*rhoB - B*rhoB) + N*p**2)
  
  #Calculate the variance of rhoB
  VrhoB <- (VS1/(lambda1**4)) * ((r3/r1)**2)  * (1 - 4*p + 4*p**2) *
    ((1/(N-B))**2) *(lambda1 - d2*lambda2)**2
  
  #Calculate the variance of rhoA
  VrhoA <- (VS1 / (lambda1**4)) * (1 - 4*p + 4*p**2) *
    (((r2-r1)/(r1*(N-B)))*lambda1 - (d2/(N-B))*lambda3)**2
  
  #Calculate the confidence intervals
  LUpperB <- log(rhoB)+qnorm(1-alpha)*sqrt(((1/rhoB)**2)*VrhoB / N)
  CLrhoB <- c(0,exp(LUpperB))
  
  LUpperA <- log(rhoA)+qnorm(1-alpha)*sqrt(((1/rhoA)**2)*VrhoA / N)
  CLrhoA <- c(0,exp(LUpperA))
  
  #Restrict the upper limits to be 1
  if(CLrhoB[2]>1) CLrhoB[2] <- 1
  if(CLrhoA[2]>1) CLrhoA[2] <- 1
  
  #Return the correlation estimates and the confidence intervals
  return(list("rhoA" = rhoA, "CLrhoA" = CLrhoA,  "rhoB" = rhoB , "CLrhoB" = CLrhoB))
}
#End Correlation function