#								  AGREE.COEFF3.RAW.R
#						 		    (June 23, 2017)
#Description: This script file contains a series of R functions for computing various agreement coefficients
#	      for multiple raters (2 or more) when the input data file is in the form of nxr matrix or data frame showing 
#             the actual ratings each rater (column) assigned to each subject (in row). That is n = number of subjects, and r = number of raters.
#             A typical table entry (i,g) represents the rating associated with subject i and rater g. 
#Author: Kilem L. Gwet, Ph.D. (gwet@agreestat.com)
#-----------------------------------------------------------------
#     EXAMPLES OF SIMPLE CALLS OF THE MAIN FUNCTIONS:
# > gwet.ac1.raw(YourRatings)       # to obtain gwet's AC1 coefficient
# > fleiss.kappa.raw(YourRatings)   # to obtain fleiss' unweighted generalized kappa coefficient
# > krippen.alpha.raw(YourRatings)  # to obtain krippendorff's unweighted alpha coefficient
# > conger.kappa.raw(YourRatings)   # to obtain conger's unweighted generalized kappa coefficient
# > bp.coeff.raw(YourRatings)       # to obtain Brennan-Prediger unweighted coefficient
#
#===========================================================================================
#gwet.ac1.raw: Gwet's AC1/Ac2 coefficient (Gwet(2008)) and its standard error for multiple raters when input 
#		   dataset is a nxr matrix of alphanumeric ratings from n subjects and r raters 
#-------------
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings
#from n subjects and r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Gwet, K. L. (2008). ``Computing inter-rater reliability and its variance in the presence of high
#		agreement." British Journal of Mathematical and Statistical Psychology, 61, 29-48.
#============================================================================================
gwet.ac1.raw <- function(ratings,weights="unweighted",conflev=0.95,N=Inf,print=TRUE){ 
  ratings.mat <- as.matrix(ratings)
  if (is.character(ratings.mat)){
    ratings.mat <- trim(toupper(ratings.mat))
    ratings.mat[ratings.mat==''] <- NA_character_
  }
  
  n <- nrow(ratings.mat) # number of subjects
  r <- ncol(ratings.mat) # number of raters
  f <- n/N # finite population correction 

  # creating a vector containing all categories used by the raters
 
  categ.init <- unique(na.omit(as.vector(ratings.mat)))
  categ <- sort(categ.init)
  q <- length(categ)

  # creating the weights matrix

  if (is.character(weights)){
     if (weights=="quadratic")
	      weights.mat<-quadratic.weights(categ)
     else if (weights=="ordinal")
	      weights.mat<-ordinal.weights(categ)
     else if (weights=="linear")
	      weights.mat<-linear.weights(categ)
     else if (weights=="radical")
	      weights.mat<-radical.weights(categ)
     else if (weights=="ratio")
	      weights.mat<-ratio.weights(categ)
     else if (weights=="circular")
	      weights.mat<-circular.weights(categ)
     else if (weights=="bipolar")
	      weights.mat<-bipolar.weights(categ)
     else weights.mat<-identity.weights(categ)
  }else weights.mat= as.matrix(weights)
  
  # creating the nxq agreement matrix representing the distribution of raters by subjects and category

  agree.mat <- matrix(0,nrow=n,ncol=q)
  for(k in 1:q){
    categ.is.k <- (ratings.mat==categ[k])
    agree.mat[,k] <- (replace(categ.is.k,is.na(categ.is.k),FALSE)) %*% rep(1,r)    
  }
  agree.mat.w <- t(weights.mat%*%t(agree.mat))
  # calculating gwet's ac1 coefficient

  ri.vec <- agree.mat%*%rep(1,q)
  sum.q <- (agree.mat*(agree.mat.w-1))%*%rep(1,q)
  n2more <- sum(ri.vec>=2)
  pa <- sum(sum.q[ri.vec>=2]/((ri.vec*(ri.vec-1))[ri.vec>=2]))/n2more

  pi.vec <- t(t(rep(1/n,n))%*%(agree.mat/(ri.vec%*%t(rep(1,q)))))
  pe <- sum(weights.mat) * sum(pi.vec*(1-pi.vec)) / (q*(q-1))
  gwet.ac1 <- (pa-pe)/(1-pe)

  # calculating variance, stderr & p-value of gwet's ac1 coefficient
  
  den.ivec <- ri.vec*(ri.vec-1)
  den.ivec <- den.ivec - (den.ivec==0) # this operation replaces each 0 value with -1 to make the next ratio calculation always possible.
  pa.ivec <- sum.q/den.ivec

  pe.r2 <- pe*(ri.vec>=2)
  ac1.ivec <- (n/n2more)*(pa.ivec-pe.r2)/(1-pe)
  pe.ivec <- (sum(weights.mat)/(q*(q-1))) * (agree.mat%*%(1-pi.vec))/ri.vec
  ac1.ivec.x <- ac1.ivec - 2*(1-gwet.ac1) * (pe.ivec-pe)/(1-pe)
  
  var.ac1 <- ((1-f)/(n*(n-1))) * sum((ac1.ivec.x - gwet.ac1)^2)
  stderr <- sqrt(var.ac1)# ac1's standard error
  p.value <- 2*(1-pt(abs(gwet.ac1/stderr),n-1))
  
  lcb <- gwet.ac1 - stderr*qt(1-(1-conflev)/2,n-1) # lower confidence bound
  ucb <- min(1,gwet.ac1 + stderr*qt(1-(1-conflev)/2,n-1)) # upper confidence bound
  if(print==TRUE) {
    if (weights=="unweighted") {
       cat("Gwet's AC1 Coefficient\n")
       cat('======================\n')	
       cat('Percent agreement:',pa,'Percent chance agreement:',pe,'\n')
       cat('AC1 coefficient:',gwet.ac1,'Standard error:',stderr,'\n')
       cat(conflev*100,'% Confidence Interval: (',lcb,',',ucb,')\n')
       cat('P-value: ',p.value,'\n')
    }
    else {
       cat("Gwet's AC2 Coefficient\n")
       cat('==========================\n')	
       cat('Percent agreement:',pa,'Percent chance agreement:',pe,'\n')
	     cat('AC2 coefficient:',gwet.ac1,'Standard error:',stderr,'\n')
       cat(conflev*100,'% Confidence Interval: (',lcb,',',ucb,')\n')
       cat('P-value: ',p.value,'\n')
	     cat('\n')
       if (!is.numeric(weights)) {
	        cat('Weights: ', weights,'\n')
          cat('---------------------\n') 	
	     }else{
	        cat('Weights: Custom Weights\n')
          cat('-----------------------\n')
       }
	     print(weights.mat)
    }
  }
  invisible(c(pa,pe,gwet.ac1,stderr,p.value))
}

#=====================================================================================
#fleiss.kappa.raw: This function computes Fleiss' generalized kappa coefficient (see Fleiss(1971)) and 
#		   its standard error for 3 raters or more when input dataset is a nxr matrix of alphanumeric 
#		   ratings from n subjects and r raters.
#-------------
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings
#from n subjects and r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Fleiss, J. L. (1981). Statistical Methods for Rates and Proportions. John Wiley & Sons.
#======================================================================================
fleiss.kappa.raw <- function(ratings,weights="unweighted",conflev=0.95,N=Inf,print=TRUE){ 
  ratings.mat <- as.matrix(ratings) 
  if (is.character(ratings.mat)){
    ratings.mat <- trim(toupper(ratings.mat))
    ratings.mat[ratings.mat==''] <- NA_character_
  }
  n <- nrow(ratings.mat) # number of subjects
  r <- ncol(ratings.mat) # number of raters
  f <- n/N # finite population correction 

  # creating a vector containing all categories used by the raters
 
  categ.init <- unique(na.omit(as.vector(ratings.mat)))
  categ <- sort(categ.init)
  q <- length(categ) # q = number of categories

  # creating the weights matrix

  if (is.character(weights)){
     if (weights=="quadratic")
	  weights.mat<-quadratic.weights(categ)
     else if (weights=="ordinal")
	  weights.mat<-ordinal.weights(categ)
     else if (weights=="linear")
	  weights.mat<-linear.weights(categ)
     else if (weights=="radical")
	  weights.mat<-radical.weights(categ)
     else if (weights=="ratio")
	  weights.mat<-ratio.weights(categ)
     else if (weights=="circular")
	  weights.mat<-circular.weights(categ)
     else if (weights=="bipolar")
	  weights.mat<-bipolar.weights(categ)
     else weights.mat<-identity.weights(categ)
  }else weights.mat= as.matrix(weights)
  
  # creating the nxq agreement matrix representing the distribution of raters by subjects and category

  agree.mat <- matrix(0,nrow=n,ncol=q)
  for(k in 1:q){
    categ.is.k <- (ratings.mat==categ[k])
    agree.mat[,k] <- (replace(categ.is.k,is.na(categ.is.k),FALSE)) %*% rep(1,r)
  }
  agree.mat.w <- t(weights.mat%*%t(agree.mat))

  # calculating fleiss's generalized kappa coefficient

  ri.vec <- agree.mat%*%rep(1,q)
  sum.q <- (agree.mat*(agree.mat.w-1))%*%rep(1,q)
  n2more <- sum(ri.vec>=2)
  pa <- sum(sum.q[ri.vec>=2]/((ri.vec*(ri.vec-1))[ri.vec>=2]))/n2more

  pi.vec <- t(t(rep(1/n,n))%*%(agree.mat/(ri.vec%*%t(rep(1,q)))))
  pe <- sum(weights.mat * (pi.vec%*%t(pi.vec)))
  fleiss.kappa <- (pa-pe)/(1-pe)

  # calculating variance, stderr & p-value of gwet's ac1 coefficient
  
  den.ivec <- ri.vec*(ri.vec-1)
  den.ivec <- den.ivec - (den.ivec==0) # this operation replaces each 0 value with -1 to make the next ratio calculation always possible.
  pa.ivec <- sum.q/den.ivec

  pe.r2 <- pe*(ri.vec>=2)
  kappa.ivec <- (n/n2more)*(pa.ivec-pe.r2)/(1-pe)
  pi.vec.wk. <- weights.mat%*%pi.vec
  pi.vec.w.k <- t(weights.mat)%*%pi.vec
  pi.vec.w <- (pi.vec.wk. + pi.vec.w.k)/2

  pe.ivec <- (agree.mat%*%pi.vec.w)/ri.vec
  kappa.ivec.x <- kappa.ivec - 2*(1-fleiss.kappa) * (pe.ivec-pe)/(1-pe)
  
  var.fleiss <- ((1-f)/(n*(n-1))) * sum((kappa.ivec.x - fleiss.kappa)^2)
  stderr <- sqrt(var.fleiss)# kappa's standard error
  p.value <- 2*(1-pt(abs(fleiss.kappa/stderr),n-1))
  
  lcb <- fleiss.kappa - stderr*qt(1-(1-conflev)/2,n-1) # lower confidence bound
  ucb <- min(1,fleiss.kappa + stderr*qt(1-(1-conflev)/2,n-1)) # upper confidence bound
  if(print==TRUE){
    cat("Fleiss' Kappa Coefficient\n")
    cat('==========================\n')	
    cat('Percent agreement:',pa,'Percent chance agreement:',pe,'\n')
    cat('Fleiss kappa coefficient:',fleiss.kappa,'Standard error:',stderr,'\n')
    cat(conflev*100,'% Confidence Interval: (',lcb,',',ucb,')\n')
    cat('P-value: ',p.value,'\n')
    cat('\n')
    if (weights != "unweighted") {
      if (!is.numeric(weights)) {
        cat('Weights: ', weights,'\n')
        cat('---------------------\n') 	
      }else{
        cat('Weights: Custom Weights\n')
        cat('-----------------------\n')
      }
      print(weights.mat)
    }
  }
  invisible(c(pa,pe,fleiss.kappa,stderr,p.value))
}
#=====================================================================================
#krippen.alpha.raw: This function computes Krippendorff's alpha coefficient (see Krippendorff(1970, 1980)) and 
#		   its standard error for 3 raters or more when input dataset is a nxr matrix of alphanumeric 
#		   ratings from n subjects and r raters.
#-------------
#The algorithm used to compute krippendorff's alpha is very different from anything that was published on this topic. Instead,
#it follows the equations presented by K. Gwet (2012)
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings
#from n subjects and r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Gwet, K. (2012). Handbook of Inter-Rater Reliability: The Definitive Guide to Measuring the Extent of Agreement Among 
#	Multiple Raters, 3rd Edition.  Advanced Analytics, LLC; 3rd edition (March 2, 2012)
#Krippendorff (1970). "Bivariate agreement coefficients for reliability of data." Sociological Methodology,2,139-150
#Krippendorff (1980). Content analysis: An introduction to its methodology (2nd ed.), New-bury Park, CA: Sage.
#======================================================================================
krippen.alpha.raw <- function(ratings,weights="unweighted",conflev=0.95,N=Inf,print=TRUE){ 
  ratings.mat <- as.matrix(ratings) 
  if (is.character(ratings.mat)){
    ratings.mat <- trim(toupper(ratings.mat))
    ratings.mat[ratings.mat==''] <- NA_character_
  }
  n <- nrow(ratings.mat) # number of subjects
  r <- ncol(ratings.mat) # number of raters
  f <- n/N # finite population correction 

  # creating a vector containing all categories used by the raters
 
  categ.init <- unique(na.omit(as.vector(ratings.mat)))
  categ <- sort(categ.init)
  q <- length(categ)

  # creating the weights matrix

  if (is.character(weights)){
     if (weights=="quadratic")
	  weights.mat<-quadratic.weights(categ)
     else if (weights=="ordinal")
	  weights.mat<-ordinal.weights(categ)
     else if (weights=="linear")
	  weights.mat<-linear.weights(categ)
     else if (weights=="radical")
	  weights.mat<-radical.weights(categ)
     else if (weights=="ratio")
	  weights.mat<-ratio.weights(categ)
     else if (weights=="circular")
	  weights.mat<-circular.weights(categ)
     else if (weights=="bipolar")
	  weights.mat<-bipolar.weights(categ)
     else weights.mat<-identity.weights(categ)
  }else weights.mat= as.matrix(weights)
  
  # creating the nxq agreement matrix representing the distribution of raters by subjects and category

  agree.mat <- matrix(0,nrow=n,ncol=q)
  for(k in 1:q){
    categ.is.k <- (ratings.mat==categ[k])
    agree.mat[,k] <- (replace(categ.is.k,is.na(categ.is.k),FALSE)) %*% rep(1,r)
  }
  agree.mat.w <- t(weights.mat%*%t(agree.mat))

  # calculating krippendorff's alpha coefficient

  ri.vec <- agree.mat%*%rep(1,q)
  agree.mat<-agree.mat[(ri.vec>=2),]
  agree.mat.w <- agree.mat.w[(ri.vec>=2),]
  ri.vec <- ri.vec[(ri.vec>=2)]
  ri.mean <- mean(ri.vec)
  n <- nrow(agree.mat)
  epsi <- 1/sum(ri.vec)
  sum.q <- (agree.mat*(agree.mat.w-1))%*%rep(1,q)
  paprime <-sum(sum.q/(ri.mean*(ri.vec-1)))/n
  pa <- (1-epsi)*paprime + epsi
  pi.vec <- t(t(rep(1/n,n))%*%(agree.mat/ri.mean))
  pe <- sum(weights.mat * (pi.vec%*%t(pi.vec)))
  krippen.alpha <- (pa-pe)/(1-pe)
  krippen.alpha.prime <- (paprime-pe)/(1-pe)

  # calculating variance, stderr & p-value of gwet's ac1 coefficient
  
  pa.ivec <- sum.q/(ri.mean*(ri.vec-1)) - pa*(ri.vec-ri.mean)/ri.mean
  krippen.ivec <- (pa.ivec-pe)/(1-pe)
  pi.vec.wk. <- weights.mat%*%pi.vec
  pi.vec.w.k <- t(weights.mat)%*%pi.vec
  pi.vec.w <- (pi.vec.wk. + pi.vec.w.k)/2

  pe.ivec <- (agree.mat%*%pi.vec.w)/ri.mean - pe * (ri.vec-ri.mean)/ri.mean
  krippen.ivec.x <- krippen.ivec - 2*(1-krippen.alpha.prime) * (pe.ivec-pe)/(1-pe)
  
  var.krippen <- ((1-f)/(n*(n-1))) * sum((krippen.ivec.x - krippen.alpha.prime)^2)
  stderr <- sqrt(var.krippen)# alpha's standard error
  p.value <- 2*(1-pt(abs(krippen.alpha/stderr),n-1))
  
  lcb <- krippen.alpha - stderr*qt(1-(1-conflev)/2,n-1) # lower confidence bound
  ucb <- min(1,krippen.alpha + stderr*qt(1-(1-conflev)/2,n-1)) # upper confidence bound
  if(print==TRUE){
    cat("Krippendorff's Alpha Coefficient\n")
    cat('==========================\n')	
    cat('Percent agreement:',pa,'Percent chance agreement:',pe,'\n')
    cat('Krippendorff alpha coefficient:',krippen.alpha,'Standard error:',stderr,'\n')
    cat(conflev*100,'% Confidence Interval: (',lcb,',',ucb,')\n')
    cat('P-value: ',p.value,'\n')
    cat('\n')
    if (weights != "unweighted") {
      if (!is.numeric(weights)) {
        cat('Weights: ', weights,'\n')
        cat('---------------------\n') 	
      }else{
        cat('Weights: Custom Weights\n')
        cat('-----------------------\n')
      }
      print(weights.mat)
    }
  }
  invisible(c(pa,pe,krippen.alpha,stderr,p.value))
}
#===========================================================================================
#conger.kappa.raw: Conger's kappa coefficient (see Conger(1980)) and its standard error for multiple raters when input 
#		   dataset is a nxr matrix of alphanumeric ratings from n subjects and r raters 
#-------------
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings
#from n subjects and r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Conger, A. J. (1980), ``Integration and Generalization of Kappas for Multiple Raters,"
#		Psychological Bulletin, 88, 322-328.
#======================================================================================
conger.kappa.raw <- function(ratings,weights="unweighted",conflev=0.95,N=Inf,print=TRUE){ 
  ratings.mat <- as.matrix(ratings) 
  if (is.character(ratings.mat)){
    ratings.mat <- trim(toupper(ratings.mat))
    ratings.mat[ratings.mat==''] <- NA_character_
  }
  n <- nrow(ratings.mat) # number of subjects
  r <- ncol(ratings.mat) # number of raters
  f <- n/N # finite population correction 

  # creating a vector containing all categories used by the raters
 
  categ.init <- unique(na.omit(as.vector(ratings.mat)))
  categ <- sort(categ.init)
  q <- length(categ)
  # creating the weights matrix

  if (is.character(weights)){
     if (weights=="quadratic")
	  weights.mat<-quadratic.weights(categ)
     else if (weights=="ordinal")
	  weights.mat<-ordinal.weights(categ)
     else if (weights=="linear")
	  weights.mat<-linear.weights(categ)
     else if (weights=="radical")
	  weights.mat<-radical.weights(categ)
     else if (weights=="ratio")
	  weights.mat<-ratio.weights(categ)
     else if (weights=="circular")
	  weights.mat<-circular.weights(categ)
     else if (weights=="bipolar")
	  weights.mat<-bipolar.weights(categ)
     else weights.mat<-identity.weights(categ)
  }else weights.mat= as.matrix(weights)
  
  # creating the nxq agreement matrix representing the distribution of raters by subjects and category

  agree.mat <- matrix(0,nrow=n,ncol=q)
  for(k in 1:q){
    categ.is.k <- (ratings.mat==categ[k])
    agree.mat[,k] <- (replace(categ.is.k,is.na(categ.is.k),FALSE)) %*% rep(1,r)
  }
  agree.mat.w <- t(weights.mat%*%t(agree.mat))
	
  # creating the rxq rater-category matrix representing the distribution of subjects by rater and category

  classif.mat <- matrix(0,nrow=r,ncol=q)
  for(k in 1:q){
  	with.mis <-(t(ratings.mat)==categ[k])
  	without.mis <- replace(with.mis,is.na(with.mis),FALSE)
  	classif.mat[,k] <- without.mis%*%rep(1,n)
  }
  # calculating conger's kappa coefficient

  ri.vec <- agree.mat%*%rep(1,q)
  sum.q <- (agree.mat*(agree.mat.w-1))%*%rep(1,q)
  n2more <- sum(ri.vec>=2)
  pa <- sum(sum.q[ri.vec>=2]/((ri.vec*(ri.vec-1))[ri.vec>=2]))/n2more
  ng.vec <- classif.mat%*%rep(1,q)
  pgk.mat <- classif.mat/(ng.vec%*%rep(1,q))
  p.mean.k <- (t(pgk.mat)%*%rep(1,r))/r 
  s2kl.mat <- (t(pgk.mat)%*%pgk.mat - r * p.mean.k%*%t(p.mean.k))/(r-1)
  pe <- sum(weights.mat * (p.mean.k%*%t(p.mean.k) -  s2kl.mat/r))
  conger.kappa <- (pa-pe)/(1-pe)

  # calculating variance, stderr & p-value of conger's kappa coefficient
  
  bkl.mat <- (weights.mat+t(weights.mat))/2
  pe.ivec1 <- r*(agree.mat%*%t(t(p.mean.k)%*%bkl.mat))
  pe.ivec2 = rep(0,n)
  
  lamda.ig.mat=matrix(0,n,r)
  if (is.numeric(ratings.mat)){
  	epsi.ig.mat <-1-is.na(ratings.mat)
  	epsi.ig.mat <- replace(epsi.ig.mat,is.na(epsi.ig.mat),FALSE)
  }else{
  	epsi.ig.mat <- 1-(ratings.mat=="")
  	epsi.ig.mat <- replace(epsi.ig.mat,is.na(epsi.ig.mat),FALSE)
  }
  for(k in 1:q){
  	lamda.ig.kmat=matrix(0,n,r)
  	for(l in 1:q){
  		delta.ig.mat <- (ratings.mat==categ[l])
  		delta.ig.mat <- replace(delta.ig.mat,is.na(delta.ig.mat),FALSE)
  		lamda.ig.kmat <- lamda.ig.kmat + weights.mat[k,l] * (delta.ig.mat - (epsi.ig.mat - rep(1,n)%*%t(ng.vec/n)) * (rep(1,n)%*%t(pgk.mat[,l])))
  	}
  	lamda.ig.kmat = lamda.ig.kmat*(rep(1,n)%*%t(n/ng.vec))
  	lamda.ig.mat = lamda.ig.mat+ lamda.ig.kmat*(r*mean(pgk.mat[,k]) - rep(1,n)%*%t(pgk.mat[,k]))
  }
  pe.ivec <- (lamda.ig.mat%*%rep(1,r)) / (r*(r-1))
  den.ivec <- ri.vec*(ri.vec-1)
  den.ivec <- den.ivec - (den.ivec==0) # this operation replaces each 0 value with -1 to make the next ratio calculation always possible.
  pa.ivec <- sum.q/den.ivec
  pe.r2 <- pe*(ri.vec>=2)
  conger.ivec <- (n/n2more)*(pa.ivec-pe.r2)/(1-pe) 
  conger.ivec.x <- conger.ivec - 2*(1-conger.kappa) * (pe.ivec-pe)/(1-pe)
  
  var.conger <- ((1-f)/(n*(n-1))) * sum((conger.ivec.x - conger.kappa)^2)
  stderr <- sqrt(var.conger)# conger's kappa standard error
  p.value <- 2*(1-pt(abs(conger.kappa/stderr),n-1))
  
  lcb <- conger.kappa - stderr*qt(1-(1-conflev)/2,n-1) # lower confidence bound
  ucb <- min(1,conger.kappa + stderr*qt(1-(1-conflev)/2,n-1)) # upper confidence bound
  if(print==TRUE) {
    cat("Conger's Kappa Coefficient\n")
    cat('==========================\n')	
    cat('Percent agreement: ',pa,'Percent chance agreement: ',pe,'\n')
    cat("Conger's kappa coefficient: ",conger.kappa,'Standard error:',stderr,'\n')
    cat(conflev*100,'% Confidence Interval: (',lcb,',',ucb,')\n')
    cat('P-value: ',p.value,'\n')
    cat('\n')
    if (weights != "unweighted") {
      if (!is.numeric(weights)) {
        cat('Weights: ', weights,'\n')
        cat('---------------------\n') 	
      }else{
        cat('Weights: Custom Weights\n')
        cat('-----------------------\n')
      }
      print(weights.mat)
    }
  }
  invisible(c(pa,pe,conger.kappa,stderr,p.value))
}

#===========================================================================================
#bp.coeff.raw: Brennan-Prediger coefficient (see Brennan & Prediger(1981)) and its standard error for multiple raters when input 
#		   dataset is a nxr matrix of alphanumeric ratings from n subjects and r raters 
#-------------
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings
#from n subjects and r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Brennan, R.L., and Prediger, D. J. (1981). ``Coefficient Kappa: some uses, misuses, and alternatives."
#           Educational and Psychological Measurement, 41, 687-699.
#======================================================================================
bp.coeff.raw <- function(ratings,weights="unweighted",conflev=0.95,N=Inf,print=TRUE){ 
  ratings.mat <- as.matrix(ratings) 
  if (is.character(ratings.mat)){
    ratings.mat <- trim(toupper(ratings.mat))
    ratings.mat[ratings.mat==''] <- NA_character_
  }
  n <- nrow(ratings.mat) # number of subjects
  r <- ncol(ratings.mat) # number of raters
  f <- n/N # finite population correction 

  # creating a vector containing all categories used by the raters
 
  categ.init <- unique(na.omit(as.vector(ratings.mat)))
  categ <- sort(categ.init)
  q <- length(categ)

  # creating the weights matrix

  if (is.character(weights)){
     if (weights=="quadratic")
	  weights.mat<-quadratic.weights(categ)
     else if (weights=="ordinal")
	  weights.mat<-ordinal.weights(categ)
     else if (weights=="linear")
	  weights.mat<-linear.weights(categ)
     else if (weights=="radical")
	  weights.mat<-radical.weights(categ)
     else if (weights=="ratio")
	  weights.mat<-ratio.weights(categ)
     else if (weights=="circular")
	  weights.mat<-circular.weights(categ)
     else if (weights=="bipolar")
	  weights.mat<-bipolar.weights(categ)
     else weights.mat<-identity.weights(categ)
  }else weights.mat= as.matrix(weights)
  
  # creating the nxq agreement matrix representing the distribution of raters by subjects and category

  agree.mat <- matrix(0,nrow=n,ncol=q)
  for(k in 1:q){
    categ.is.k <- (ratings.mat==categ[k])
    agree.mat[,k] <- (replace(categ.is.k,is.na(categ.is.k),FALSE)) %*% rep(1,r)
  }
  agree.mat.w <- t(weights.mat%*%t(agree.mat))

  # calculating gwet's ac1 coefficient

  ri.vec <- agree.mat%*%rep(1,q)
  sum.q <- (agree.mat*(agree.mat.w-1))%*%rep(1,q)
  n2more <- sum(ri.vec>=2)
  pa <- sum(sum.q[ri.vec>=2]/((ri.vec*(ri.vec-1))[ri.vec>=2]))/n2more

  pi.vec <- t(t(rep(1/n,n))%*%(agree.mat/(ri.vec%*%t(rep(1,q)))))
  pe <- sum(weights.mat) / (q^2)
  bp.coeff <- (pa-pe)/(1-pe)

  # calculating variance, stderr & p-value of gwet's ac1 coefficient
  
  den.ivec <- ri.vec*(ri.vec-1)
  den.ivec <- den.ivec - (den.ivec==0) # this operation replaces each 0 value with -1 to make the next ratio calculation always possible.
  pa.ivec <- sum.q/den.ivec

  pe.r2 <- pe*(ri.vec>=2)
  bp.ivec <- (n/n2more)*(pa.ivec-pe.r2)/(1-pe)
  var.bp <- ((1-f)/(n*(n-1))) * sum((bp.ivec - bp.coeff)^2)
  stderr <- sqrt(var.bp)# BP's standard error
  p.value <- 2*(1-pt(abs(bp.coeff/stderr),n-1))
  
  lcb <- bp.coeff - stderr*qt(1-(1-conflev)/2,n-1) # lower confidence bound
  ucb <- min(1,bp.coeff + stderr*qt(1-(1-conflev)/2,n-1)) # upper confidence bound
  if(print==TRUE) {
    cat("Brennan-Prediger Coefficient\n")
    cat('============================\n')	
    cat('Percent agreement:',pa,'Percent chance agreement:',pe,'\n')
    cat('B-P coefficient:',bp.coeff,'Standard error:',stderr,'\n')
    cat(conflev*100,'% Confidence Interval: (',lcb,',',ucb,')\n')
    cat('P-value: ',p.value,'\n')
    cat('\n')
    if (weights != "unweighted") {
      if (!is.numeric(weights)) {
        cat('Weights: ', weights,'\n')
        cat('---------------------\n') 	
      }else{
        cat('Weights: Custom Weights\n')
        cat('-----------------------\n')
      }
      print(weights.mat)
    }
  }
  invisible(c(pa,pe,bp.coeff,stderr,p.value))
}

#
#-----  Additional functions needed to run the main functions. If the main functions must be included in another R script, then
# the user will need to add these additional functions to the new script file.
#
# ==============================================================
# trim(x): This is an r function for trimming leading and trealing blanks
# ==============================================================
trim <- function( x ) {
  gsub("(^[[:space:]]+|[[:space:]]+$)", "", x) 
}
# ==============================================================
# The following functions generate various weight matrices used 
# in the weighted or unweighted analyses.
# ==============================================================
identity.weights<-function(categ){
  weights<-diag(length(categ))
  return (weights)
}

quadratic.weights<-function(categ){
  q<-length(categ)
  weights <- diag(q)
  if (is.numeric(categ)) { 
    categ.vec <- sort(categ)
  }
  else {
    categ.vec<-1:length(categ)
  }
  xmin<-min(categ.vec)
  xmax<-max(categ.vec)
  for(k in 1:q){
    for(l in 1:q){
      weights[k,l] <- 1-(categ.vec[k]-categ.vec[l])^2/(xmax-xmin)^2 
    }
  }
  return (weights)
}
linear.weights<-function(categ){
  q<-length(categ)
  weights <- diag(q)
  if (is.numeric(categ)) { 
    categ.vec <- sort(categ)
  }
  else {
    categ.vec<-1:length(categ)
  }
  xmin<-min(categ.vec)
  xmax<-max(categ.vec)
  for(k in 1:q){
    for(l in 1:q){
      weights[k,l] <- 1-abs(categ.vec[k]-categ.vec[l])/abs(xmax-xmin)
    }
  }
  return (weights)
}
#--------------------------------
radical.weights<-function(categ){
  q<-length(categ)
  weights <- diag(q)
  if (is.numeric(categ)) { 
    categ.vec <- sort(categ)
  }
  else {
    categ.vec<-1:length(categ)
  }
  xmin<-min(categ.vec)
  xmax<-max(categ.vec)
  for(k in 1:q){
    for(l in 1:q){
      weights[k,l] <- 1-sqrt(abs(categ.vec[k]-categ.vec[l]))/sqrt(abs(xmax-xmin))
    }
  }
  return (weights)
}
#--------------------------------
ratio.weights<-function(categ){
  q<-length(categ)
  weights <- diag(q)
  if (is.numeric(categ)) { 
    categ.vec <- sort(categ)
  }
  else {
    categ.vec<-1:length(categ)
  }
  xmin<-min(categ.vec)
  xmax<-max(categ.vec)
  for(k in 1:q){
    for(l in 1:q){
      weights[k,l] <- 1-((categ.vec[k]-categ.vec[l])/(categ.vec[k]+categ.vec[l]))^2 / ((xmax-xmin)/(xmax+xmin))^2
    }
  }
  return (weights)
}

#--------------------------------
circular.weights<-function(categ){
  q<-length(categ)
  weights <- diag(q)
  if (is.numeric(categ)) { 
    categ.vec <- sort(categ)
  }
  else {
    categ.vec<-1:length(categ)
  }
  xmin<-min(categ.vec)
  xmax<-max(categ.vec)
  U = xmax-xmin+1
  for(k in 1:q){
    for(l in 1:q){
      weights[k,l] <- (sin(pi*(categ.vec[k]-categ.vec[l])/U))^2
    }
  }
  weights <- 1-weights/max(weights)
  return (weights)
}

#--------------------------------
bipolar.weights<-function(categ){
  q<-length(categ)
  weights <- diag(q)
  if (is.numeric(categ)) { 
    categ.vec <- sort(categ)
  }
  else {
    categ.vec<-1:length(categ)
  }
  xmin<-min(categ.vec)
  xmax<-max(categ.vec)
  for(k in 1:q){
    for(l in 1:q){
      if (k!=l)
        weights[k,l] <- (categ.vec[k]-categ.vec[l])^2 / (((categ.vec[k]+categ.vec[l])-2*xmin)*(2*xmax-(categ.vec[k]+categ.vec[l])))
      else weights[k,l] <- 0
    }
  }
  weights <- 1-weights/max(weights)
  return (weights)
}
#--------------------------------
ordinal.weights<-function(categ){
  q<-length(categ)
  weights <- diag(q)
  categ.vec<-1:length(categ)
  for(k in 1:q){
    for(l in 1:q){
      nkl <- max(k,l)-min(k,l)+1
      weights[k,l] <- nkl * (nkl-1)/2
    }
  }
  weights <- 1-weights/max(weights)
  return (weights)
}