R/cor.cohen.kappa.onesample.1971.fleiss.R

In burrm/lolcat: Miscellaneous Useful Functions

#
# Fleiss. Measuring Nominal Scale Agreement Among Many Raters. 
# Psychological Bulletin. 1971 76-5. pg 378-382 
#

#' Cohen's Kappa   
#' 
#' Calculate Fleiss' extension of Cohen's Kappa based on Fleiss' 1971 paper.
#'
#' @param subject A vector identifying subjects.
#' @param rater A vector identifying raters. Not used in this calculation.
#' @param rating A vector identifying ratings.
#' @param weight A vector identifying weights.
#' @param alternative The alternative hypothesis to use for the test computation.
#' @param conf.level The confidence level for this test, between 0 and 1.
#'
#' @return The results of the statistical test.

cor.cohen.kappa.onesample.1971.fleiss <- function(
    subject #
   ,rater  #not used, standardized data input for kappa and Light's G
   ,rating
   ,weight = rep(1, length(rating))
   ,alternative = c("two.sided", "greater", "less")
   ,conf.level = .95
) {
  validate.htest.alternative(alternative = alternative)
  
  if (any(weight < 0)) {
    weight[which(weight < 0)] <- 0
  }

  if (!is.factor(subject)) {
    subject <- factor(subject)
  }

  #if (!is.factor(rater)) {
  #  rater <- factor(rater)
  #}

  if (!is.factor(rating)) {
    rating <- factor(rating)
  }

  # N <- number of subjects
  N <- length(levels(subject))
  total_assignments <- sum(weight)

  #n_i  <- number of ratings for a subject
  #n_j  <- number of ratings assigned to a category
  #n_ij <- number of raters who assigned the ith subject to the jth category

  n_i <- list()
  n_j <- list()
  n_ij <- list()

  unique_subjects <- levels(subject)
  unique_ratings <- levels(rating)

  for (i in unique_subjects) {
    n_ij[[i]] <- list()

    for (j in unique_ratings) {
      n_ij[[i]][[j]] <- 0
    }
  }

  for (i in unique_subjects) {
    n_i[[i]] <- 0
  }

  for (j in unique_ratings) {
      n_j[[j]] <- 0
  }

  for (idx in 1:length(weight)) {
    this_subject <- unique_subjects[subject[idx]]
    this_rating  <- unique_ratings[rating[idx]]

    n_ij[[this_subject]][[this_rating]] <- n_ij[[this_subject]][[this_rating]] + weight[idx]
    n_i[[this_subject]] <- n_i[[this_subject]] + weight[idx]
    n_j[[this_rating]]  <- n_j[[this_rating]]  + weight[idx] 
  }

  #p_j <- proportion of all assignments made to the jth category
  #p_j <- (1/(Nn))*sum[over subjects, i](n_ij), use n_j and total_assignments instead 

  p_j <- list()
  for (j in unique_ratings) {
      p_j[[j]] <- n_j[[j]] / total_assignments
  }

  #P_i <- extent of agrement among n raters, proportion of agreeing pairs out of (n(n-1 possible))

  P_i <- list()
  P_bar <- 0

  for (i in unique_subjects) {
    P_i[[i]] <- 0

    for (j in unique_ratings) {
      P_i[[i]] <- P_i[[i]] + (n_ij[[i]][[j]] * (n_ij[[i]][[j]]-1))
    }

    P_i[[i]] <- P_i[[i]] / (n_i[[i]] * (n_i[[i]] - 1))

    P_bar <- P_bar + P_i[[i]]
  }

  P_bar <- P_bar / N

  #P_bar_c <- mean proportion agreement expected by chance
  #sum_p_j_cubed <- used for variance
  P_bar_c <- 0
  sum_p_j_cubed <- 0

  for (j in unique_ratings) {
     P_bar_c <- P_bar_c + (p_j[[j]])^2
     sum_p_j_cubed <- sum_p_j_cubed + (p_j[[j]])^3
  }

  kappa <- (P_bar - P_bar_c) / (1 - P_bar_c)
  #n for variance calculation, all subjects should have equal number of raters, but...
  n_var <- mean(rmnames(unlist(n_i)))
  se.kappa <- sqrt((2/(N*n_var*(n_var-1)) * 
                   ((P_bar_c - ((2*n_var - 3) * (P_bar_c^2)) + ((2*(n_var - 2)*sum_p_j_cubed)))/
                   ((1-P_bar_c)^2))))

  

  
  #SE Kappa and CI

  z <- kappa/se.kappa
 
  cv      <- qnorm(conf.level+(1-conf.level)/2)
  ci.add  <- cv*se.kappa
  
  ci.lower <- max(-1,kappa-ci.add)
  ci.upper <- min(1,kappa+ci.add)
  
  p.value <- if (alternative[1] == "two.sided") {
    tmp<-pnorm(z)
    min(tmp,1-tmp)*2
  } else if (alternative[1] == "greater") {
    pnorm(z,lower.tail = FALSE)
  } else if (alternative[1] == "less") {
    pnorm(z,lower.tail = TRUE)
  } else {
    NA
  }
  
  
  

  retval<-list(data.name   = "subjects and ratings",
               statistic   = z, 
               estimate    = c(kappa = kappa 
                               ,se.kappa = se.kappa
               ),
               parameter   = 0,
               p.value     = p.value,
               null.value  = 0,
               alternative = alternative[1],
               method      = "Cohen's Kappa (Fleiss 1971 Extension)",
               conf.int    = c(ci.lower, ci.upper)
  )
  
  #names(retval$estimate) <- c("sample mean")
  names(retval$statistic) <- "z"
  names(retval$null.value) <- "kappa"
  names(retval$parameter) <- "null hypothesis kappa"
  attr(retval$conf.int, "conf.level")  <- conf.level
  
  class(retval)<-"htest"
  retval
}