R/wquad.conc.R

In raters: A Modification of Fleiss' Kappa in Case of Nominal and Ordinal Variables

wquad.conc <-
  function(db, test = "Default", B = 1000, alpha = 0.05) {
    C <- ncol(db)
    R <- nrow(db)
    row.sum <- rowSums(db)
    db2 <- db^2
    row.sum2 <- rowSums(db2)

    w <- matrix(, nrow = C, ncol = C)
    for (j in 1:C) {
      for (k in 1:C) {
        w[j, k] <- 1 - ((abs(j - k))^2 / (C - 1)^2)
      }
    }

    ww <- rep(0, R)

    for (i in 1:R) {
      for (j in 1:(C - 1)) {
        l <- j + 1
        for (k in l:C) {
          ww[i] <- ww[i] + (db[i, j] * db[i, k] * w[j, k])
        }
      }
    }

    xi <- c()
    pi <- c()
    for (h in 1:R) {
      xi[h] <- (0.5 * row.sum2[h]) + ww[h] - (0.5 * row.sum[h])
      pi[h] <- (2 * xi[h]) / (row.sum[h] * (row.sum[h] - 1))
    }

    p.avg <- mean(pi)
    pe <- (1 / (C^2)) * sum(w)

    s.star <- (p.avg - pe) / (1 - pe)


    obj.min <- c()
    for (i in 1:R) {
      obj.min[i] <- (row.sum[i] - 2) / (row.sum[i] - 1)
    }
    sum.obj.min <- sum(obj.min) * ((3 * C) / (2 * R))

    min.s <- (sum.obj.min - 2 * C + 1) / (C + 1)


    s.star
    pi.boot.w <- list()
    p.boot.w <- c()
    s.boot.w <- c()
    for (i in 1:B) {
      pi.boot.w[[i]] <- sample(pi, size = nrow(db), replace = TRUE)
      p.boot.w[i] <- mean(pi.boot.w[[i]])
      s.boot.w[i] <- (p.boot.w[i] - pe) / (1 - pe)
    }
    s.boot.ci.w <- quantile(s.boot.w, probs = c(alpha / 2, 1 - alpha / 2))


    Default <- function(db) {
      s.res <- c(s.star, min.s, s.boot.ci.w)
      names(s.res) <- c("S*", "min", "LCL", "UCL")
      s.res
    }

    MC <- function(db) {
      matrix.mc <- list()
      matrix.mc2 <- list()
      R <- nrow(db)
      C <- ncol(db)
      w.sum <- list()
      rs.mc <- list()
      rs2.mc <- list()
      xi.mc <- list()
      pi.mc <- list()

      for (h in 1:B) {
        matrix.mc[[h]] <- matrix(,
          nrow = R, ncol = C,
          byrow = T
        )
        matrix.mc2[[h]] <- matrix(, nrow = R, ncol = C, byrow = T)
        w.sum[[h]] <- rep(0, R)
        xi.mc[[h]] <- rep(0, R)
        pi.mc[[h]] <- rep(0, R)
      }

      for (k in 1:B) {
        for (j in 1:R) {
          matrix.mc[[k]][j, ] <- t(rmultinom(1, size = rowSums(db)[j], prob = rep(1 / C, C)))
          matrix.mc2[[k]] <- (matrix.mc[[k]])^2
          rs.mc[[k]] <- rowSums(matrix.mc[[k]])
          rs2.mc[[k]] <- rowSums(matrix.mc2[[k]])
        }
      }


      for (g in 1:B) {
        for (i in 1:R) {
          for (j in 1:(C - 1)) {
            l <- j + 1
            for (k in l:C) {
              w.sum[[g]][i] <- w.sum[[g]][i] + (matrix.mc[[g]][i, j] * matrix.mc[[g]][i, k] * w[j, k])
            }
          }
        }
      }

      p.avg.mc <- c()
      pe.mc <- c()
      s.star.mc <- c()
      for (g in 1:B) {
        for (h in 1:R) {
          xi.mc[[g]][h] <- (0.5 * rs2.mc[[g]][h]) + w.sum[[g]][h] - (0.5 * rs.mc[[g]][h])
          pi.mc[[g]][h] <- (2 * xi.mc[[g]][h]) / (rs.mc[[g]][h] * (rs.mc[[g]][h] - 1))
          p.avg.mc[g] <- mean(pi.mc[[g]])
          pe.mc[g] <- (1 / (C^2)) * sum(w)
          s.star.mc[g] <- (p.avg.mc[g] - pe.mc[g]) / (1 - pe.mc[g])
        }
      }


      crit.s.mc <- quantile(s.star.mc, 0.95)
      binary <- c()
      for (i in 1:length(s.star.mc)) {
        binary[i] <- (s.star.mc[i] >= s.star)
      }
      pvalue <- sum(binary) / B
      s.res <- c(s.star, min.s, s.boot.ci.w, pvalue)
      names(s.res) <- c("S*", "min", "LCL", "UCL", "pvalue")
      s.res
    }
    switch(test,
      Default = Default(db),
      MC = MC(db)
    )
  }