#
# 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
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.