## rhoR coverage - 100.00%

1 |  |
    ###
2 |  |
    #' @include createSimulatedCodeSet.R
3 |  |
    #' @title Get Test Set
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function gets a \emph{testSet} from a larger \code{\link{codeSet}} given certain sampling parameters.
7 |  |
    #'
8 |  |
    #' @details
9 |  |
    #' A \emph{testSet} is a \code{\link{codeSet}} that is a subset of a larger \code{\link{codeSet}} with a given set of properties.  A \emph{testSet} is constructed by sampling (without replacement) P rows from rows in the larger \code{\link{codeSet}} where the first rater's code was 1, and then appending an additional sample (without replacement) of R rows taken at random from the larger \code{\link{codeSet}} excluding rows included in the first P rows sampled. P is computed as the minbaserate * length of the \emph{testset}. R is computed as testSetLength - P. The result of this sampling procedure is to create a sample with a minimum baserate regardless of the baserate of the larger \code{\link{codeSet}}.If \emph{testSetBaserateInflation} is set to zero, the function selects rows at random.
10 |  |
    #'
11 |  |
    #'
12 |  |
    #'
13 |  |
    #' @export
14 |  |
    #'
15 |  |
    #' @param set The \code{\link{codeSet}} from which the \emph{testSet} is taken
16 |  |
    #' @param testSetLength The length of the \emph{testSet} to be taken
17 |  |
    #' @param testSetBaserateInflation The minimum guaranteed \code{\link{baserate}} of the \emph{testSet}. Default to 0
18 |  |
    #'
19 |  |
    #' @keywords testSet
20 |  |
    #'
21 |  |
    #' @return A \code{\link{codeSet}} with the properties specified
22 |  |
    #'
23 |  |
    ###
24 |  |
25 |  |
    getTestSet = function(set, testSetLength, testSetBaserateInflation = 0){
26 | 1x |
      if(length(unlist(set)) < testSetLength){stop("testSetLength must be less than the length of set")}
27 | 1x |
      if(testSetLength %% 1 | testSetLength < 1) {stop("testSetLength value must be a positive integer.")}
28 | 1x |
      if(testSetBaserateInflation < 0){stop("testSetBaserateInfation must be positive")}
29 | 1x |
      if(testSetBaserateInflation >= 1){stop("testSetBaserateInflation must be below 1")}
30 | 5x |
      if(length(which(set[,1] == 1 & set[,2] == 1)) == nrow(set)){
31 | 1x |
        return(set[1:testSetLength,])
32 | 4x |
      }else if(length(which(set[,1] == 0 & set[,2] == 0)) == nrow(set)){
33 | 2x |
        if(testSetBaserateInflation == 0){
34 | 1x |
          return(set[1:testSetLength,])
35 |  |
        }else{
36 | 1x |
          stop("Not enough positives in first rater to inflate to this level")
37 |  |
        }
38 |  |
      }
39 | 2x |
      return(getHandSet(set = set, handSetLength = testSetLength, handSetBaserate = testSetBaserateInflation, returnSet = T))
40 |  |
    }

1 |  |
    ###
2 |  |
    #' @include kappaCT.R
3 |  |
    #' @title Rho Min
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function calculates the minimum testSetLength where it is possible to get a rho less than alpha for the given parameters of rho.
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #' @param baserate A \code{\link{baserate}}
10 |  |
    #' @param alpha The threshold of significance for rho (similar to an alpha level for a p value), defaulted to 0.05
11 |  |
    #' @param inc An integer indicating by how much the testSetLength should increase each iteration
12 |  |
    #' @param printInc A boolean indicating whether to print out each increment value with it's corresponding significance for rho
13 |  |
    #' @param ... Any additional parameters passed into \code{\link{rho}}
14 |  |
    #' @keywords rhoMin
15 |  |
    #'
16 |  |
    #' @examples
17 |  |
    #' #Add testSetBaserateInflation as an additional parameter
18 |  |
    #' rhoMin(0.2, testSetBaserateInflation = 0.33)
19 |  |
    #'
20 |  |
    #' #Add testSetBaserateInflation as well as changing inc and selecting printInc
21 |  |
    #' rhoMin(0.2, inc = 5, printInc = TRUE, testSetBaserateInflation = 0.33)
22 |  |
    #'
23 |  |
    #' @return
24 |  |
    #' The minimum length of testSet, to the nearest multiple of inc, greater than the minimum length, that would give a value where rho less than alpha becomes mathematically possible.
25 |  |
    ###
26 |  |
    rhoMin = function(baserate, alpha = 0.05, inc = 10, printInc = FALSE, ...){
27 | 4x |
      rhoVal = 1
28 | 4x |
      testSetLength = 0
29 | 1x |
      if(inc%%1 != 0 | inc <= 0) stop("Inc value must be a positive integer.")
30 | 1x |
      if(alpha > 1 | alpha < 0) stop("Alpha must be between 0 and 1.")
31 |  |
32 |  |
33 | 2x |
      while(rhoVal > alpha){
34 | 6x |
        testSetLength = testSetLength + inc
35 | 6x |
        rhoVal = rho(1, baserate, testSetLength, ...)
36 | 6x |
        if(printInc) {
37 | 3x |
          cat(testSetLength, rhoVal,"\n")
38 |  |
        }
39 |  |
      }
40 | 2x |
      return(testSetLength)
41 |  |
42 |  |
    }

1 |  |
    ###
2 |  |
    #' @include rhoK.R
3 |  |
    #' @title Rho (set)
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function calculates rho and kappa for a given \code{\link[=getTestSet]{testSet}}, and returns a list containing both values. Called by \code{\link{rho}}.
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #'
10 |  |
    #' @template rho-params
11 |  |
    #'
12 |  |
    #' @return A list of the format:\describe{
13 |  |
    #' \item{rho}{The rho of the \code{\link{codeSet}}}
14 |  |
    #' \item{kappa}{The Cohen's Kappa of the \code{\link{codeSet}}}
15 |  |
    #' }
16 |  |
    ###
17 |  |
    rhoSet = function(
18 |  |
      x,
19 |  |
      OcSBaserate = NULL,
20 |  |
      testSetBaserateInflation = 0,
21 |  |
      OcSLength = 10000,
22 |  |
      replicates = 800,
23 |  |
      ScSKappaThreshold = .65,
24 |  |
      ScSKappaMin = .40,
25 |  |
      ScSPrecisionMin = 0.6,
26 |  |
      ScSPrecisionMax = 1
27 |  |
    ) {
28 | 4x |
      kappa = kappaSet(x);
29 |  |
30 | 4x |
      row.count = nrow(x);
31 | 4x |
      if(is.null(OcSBaserate)) {
32 | 3x |
        OcSBaserate = length(which(x[,1] == 1)) / row.count;
33 |  |
      }
34 |  |
35 | 4x |
      return(list(
36 | 4x |
        rho = rhoK(
37 | 4x |
          kappa,
38 | 4x |
          OcSBaserate,
39 | 4x |
          row.count,
40 | 4x |
          testSetBaserateInflation,
41 | 4x |
          OcSLength,
42 | 4x |
          replicates,
43 | 4x |
          ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax
44 |  |
        ),
45 | 4x |
        kappa = kappa
46 |  |
      ));
47 |  |
    }

1 |  |
    ###
2 |  |
    #' @include rhoCT.R
3 |  |
    #' @title Calculate kappa
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function calculates Cohen's kappa on a \code{\link{contingencyTable}} or a \code{\link{codeSet}}
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #' @param data A \code{\link{contingencyTable}} or a \code{\link{codeSet}}
10 |  |
    #'
11 |  |
    #' @keywords kappa
12 |  |
    #'
13 |  |
    #' @seealso \code{\link{kappaSet}} and \code{\link{kappaCT}}
14 |  |
    #'
15 |  |
    #' @examples
16 |  |
    #' #Given a code set
17 |  |
    #' kappa(data = codeSet)
18 |  |
    #'
19 |  |
    #' #Given a contingency Table
20 |  |
    #' kappa(data = contingencyTable)
21 |  |
    #'
22 |  |
    #' @return
23 |  |
    #' The kappa of the \code{\link{contingencyTable}} or \code{\link{codeSet}}
24 |  |
    ###
25 |  |
    kappa <- function(data) {
26 | 61x |
      if (all(dim(data) == 2)) {
27 | 5x |
        return(kappaCT(data))
28 |  |
      } else {
29 | 56x |
        return(kappaSet(data))
30 |  |
      }
31 |  |
    }

1 |  |
    ###
2 |  |
    #' @title Get Handset
3 |  |
    #' @description This function is to get a handset of a set and calculate the kappa
4 |  |
    #' @param set This is the set to take a handset of
5 |  |
    #' @param handSetLength This is the length of the handset to take
6 |  |
    #' @param handSetBaserate This is the minimum baserate to inflate the handset to
7 |  |
    #' @param returnSet If TRUE, then return the handSet if FALSE, return the kappa of the handSet
8 |  |
    #' @keywords handset, hand, inflation
9 |  |
    #' @return The function returns the handSet if returnSet is TRUE or the kappa of the handSet if not
10 |  |
    #' @export
11 |  |
    ###
12 |  |
    getHandSet = function(set, handSetLength, handSetBaserate, returnSet=FALSE) {
13 |  |
      #positives is the minimum number of positive pairs in the first rater
14 | 28305x |
      positives = ceiling(handSetLength * handSetBaserate);
15 | 28305x |
      posInd = which(set[,1] == 1);
16 |  |
17 | 1x |
      if(positives > length(posInd)){stop("Not enough positives in first rater to inflate to this level")}
18 |  |
      #if there is an inflated baserate positives will be > 0, so if handSetBaserate > 0, positives > 0
19 | 28304x |
      if (positives > 0) {
20 |  |
        #goes through and picks a set that fits the number of positives first then randomly samples from the rest
21 |  |
        #if the set created is not a valid set, then it reruns the process until a valid set is found
22 | 11502x |
        positiveIndices = posInd[sample.int(length(posInd),size=positives,replace=FALSE)];
23 | 11502x |
        others = set[!(1:nrow(set) %in% positiveIndices),];
24 | 11502x |
        otherIndices = sample.int(nrow(others),size=(handSetLength - positives),replace=FALSE);
25 | 11502x |
        this.set = rbind(set[positiveIndices,], others[otherIndices,]);
26 | 16802x |
      } else if (positives == 0){
27 |  |
        #if there is no restriction on positives, then the set is generated by randomly sampling from the entire set
28 |  |
        #this.set=matrix(0,2,2);
29 |  |
        #while ((sum(this.set[,2]) == 0 | sum(this.set[,2]) == nrow(this.set))) {
30 | 16802x |
          theseIndices = sample.int(nrow(set),size=handSetLength,replace=FALSE);
31 | 16802x |
          this.set = set[theseIndices,];
32 |  |
        #}
33 |  |
      }
34 |  |
35 |  |
      #if return set is true return set rather than kappa
36 | 28304x |
      if(returnSet){
37 | 3x |
        return(this.set);
38 |  |
      }
39 |  |
40 |  |
      #otherwise return the kappa of the set
41 | 28301x |
      handSetKappa = calcKappa(this.set);
42 | 28301x |
      return(handSetKappa);
43 |  |
    }
44 |  |
45 |  |
    ###
46 |  |
    #' @title Get Handset
47 |  |
    #' @description This function is to get a handset of a set and calculate the kappa
48 |  |
    #'
49 |  |
    #' @param full.ct This is the set to take a handset of
50 |  |
    #' @param handSetLength This is the length of the handset to take
51 |  |
    #' @param handSetBaserate This is the minimum baserate to inflate the handset to
52 |  |
    #' @param as_kappa If FALSE then return the handSet, if TRUE (default) return the kappa of the handSet
53 |  |
    #'
54 |  |
    #' @return The function returns the handSet if returnSet is TRUE or the kappa of the handSet if not
55 |  |
    #' @export
56 |  |
    ###
57 |  |
    getHandCT <- function(full.ct, handSetLength, handSetBaserate, as_kappa = TRUE) {
58 | 3x |
      positives <- ceiling(handSetLength * handSetBaserate)
59 |  |
60 | 3x |
      if (positives > sum(full.ct[1,])) {
61 | 1x |
        stop("Not enough positives in first rater to inflate to this level")
62 |  |
      }
63 |  |
64 | 2x |
      positiveIndices <- sample_ct(full.ct[1,], positives)
65 |  |
66 | 2x |
      other_ct <- full.ct
67 | 2x |
      other_ct[1, ] <- other_ct[1, ] - positiveIndices
68 |  |
69 | 2x |
      other_ct_vec <- as.vector(other_ct)
70 | 2x |
      this_ct <- sample_ct(other_ct_vec,  handSetLength - positives)
71 | 2x |
      this_ct[1, ] <- this_ct[1, ] + positiveIndices
72 |  |
73 | 2x |
      if (!as_kappa) {
74 | 1x |
        return(this_ct)
75 |  |
      }
76 |  |
77 | 1x |
      return(kappa_ct(this_ct))
78 |  |
    }
79 |  |
80 |  |
    sample_ct <- function(x, size) {
81 | 4x |
      x_vec <- as.vector(x)
82 | 4x |
      x_sampled <- matrix(rep(0, length(x)), ncol = length(x_vec) / 2)
83 | 4x |
      for(i in seq(size)) {
84 | 210x |
        s <- sample(seq(x_vec), 1, replace = F, prob = x_vec / sum(x_vec))
85 | 210x |
        x_sampled[[s]] = x_sampled[[s]] + 1
86 | 210x |
        x_vec[s] = x_vec[s] - 1
87 |  |
      }
88 |  |
89 | 4x |
      return(x_sampled)
90 |  |
    }

1 |  |
    ###
2 |  |
    #' @include rhoSet.R
3 |  |
    #' @title Rho (contingency Table)
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function calculates rho and kappa for a given \code{\link{contingencyTable}}, and returns a list containing both values. Called by \code{\link{rho}}.
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #'
10 |  |
    #' @template rho-params
11 |  |
    #'
12 |  |
    #' @return A list of the format:\describe{
13 |  |
    #' \item{rho}{The rho of the \code{\link{contingencyTable}}}
14 |  |
    #' \item{kappa}{The Cohen's Kappa of the \code{\link{contingencyTable}}}
15 |  |
    #' }
16 |  |
    ###
17 |  |
    rhoCT = function(
18 |  |
      x,
19 |  |
      OcSBaserate = NULL,
20 |  |
      testSetBaserateInflation = 0,
21 |  |
      OcSLength = 10000,
22 |  |
      replicates = 800,
23 |  |
      ScSKappaThreshold = .65,
24 |  |
      ScSKappaMin = .40,
25 |  |
      ScSPrecisionMin = 0.6,
26 |  |
      ScSPrecisionMax = 1
27 |  |
    ){
28 | 2x |
      if (any(x < 0)) {
29 | 1x |
        stop("Values in Contingency Table must be positive")
30 |  |
      }
31 | 1x |
      kappa = kappa_ct(x);
32 |  |
33 | 1x |
      if (is.null(OcSBaserate)) {
34 | 1x |
        OcSBaserate = (sum(x[1,])) / (sum(x))
35 |  |
      }
36 |  |
37 | 1x |
      rho_result <- rhoK(
38 | 1x |
        kappa,
39 | 1x |
        OcSBaserate,
40 | 1x |
        sum(x),
41 | 1x |
        testSetBaserateInflation,
42 | 1x |
        OcSLength,
43 | 1x |
        replicates,
44 | 1x |
        ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax
45 |  |
      )
46 |  |
47 | 1x |
      return(list(
48 | 1x |
        rho = rho_result,
49 | 1x |
        kappa = kappa
50 |  |
      ))
51 |  |
    }

1 |  |
    ###
2 |  |
    #' @include codeSet.R
3 |  |
    #' @title Create Simulated codeSet
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' Creates a simulated \code{\link{codeSet}} with the given parameters
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #'
10 |  |
    #' @param length the length of the simulated \code{\link{codeSet}} to be created
11 |  |
    #' @param baserate the \code{\link{baserate}} of the simulated \code{\link{codeSet}}
12 |  |
    #' @param kappaMin the minimum kappa of the simulated \code{\link{codeSet}}
13 |  |
    #' @param kappaMax the maximum kappa of the simulated \code{\link{codeSet}}
14 |  |
    #' @param precisionMin the minimum precision of the simulated \code{\link{codeSet}}
15 |  |
    #' @param precisionMax the maximum precision of the simulated \code{\link{codeSet}}
16 |  |
    #' @param ... Parameters passed to createRandomSet (e.g. type = "set" or type = "ct")
17 |  |
    #' @param tries the maximum number of tries to generate a valid set, smaller set lengths may require an increased number of tries
18 |  |
    #'
19 |  |
    #' @keywords codeSet, create
20 |  |
    #'
21 |  |
    #' @details
22 |  |
    #' \code{\link{codeSet}}s are generated by first picking a random kappa within its range and a random precision within its range.  If the random kappa, random precision, and baserate are not mathematically possible, then the precision is resampled from a range of mathematically possible values within its range.  A unique simulated \code{\link{codeSet}} is then constructed given these parameters.
23 |  |
    #'
24 |  |
    #' @return A \code{\link{codeSet}} that fulfills the given parameters
25 |  |
    #'
26 |  |
    ###
27 |  |
    createSimulatedCodeSet = function(length, baserate, kappaMin, kappaMax, precisionMin, precisionMax, ..., tries = 50){
28 | 1x |
      if(length <= 0){stop("length must be positive")}
29 | 1x |
      if(baserate < 0){stop("baserate must be positive")}
30 | 2x |
      if(kappaMin < 0 | kappaMin > 1) stop("kappaMin must be between 0 and 1.")
31 | 2x |
      if(kappaMax < 0 | kappaMax > 1) stop("kappaMax must be between 0 and 1.")
32 | 1x |
      if(kappaMin > kappaMax){stop("kappaMin must be less than kappaMax")}
33 | 2x |
      if(precisionMin < 0 | precisionMin > 1) stop("precisionMin must be between 0 and 1.")
34 | 2x |
      if(precisionMax < 0 | precisionMax > 1) stop("precisionMax must be between 0 and 1.")
35 | 1x |
      if(precisionMin > precisionMax){stop("precisionMin must be less than precisionMax")}
36 | 1x |
      if(tries < 1){stop("tries must be greater than 1")}
37 |  |
38 | 2x |
      setFound = NULL
39 | 2x |
      validSet = F
40 | 2x |
      triesLeft = tries
41 | 2x |
      Ks = c()
42 | 2x |
      while(!validSet && triesLeft > 0) {
43 | 51x |
        set = createRandomSet(setLength = length, baserate = baserate, kappaMin = kappaMin, kappaMax = kappaMax, minPrecision = precisionMin, maxPrecision = precisionMax)
44 | 51x |
        k = kappa(set)
45 |  |
46 | 51x |
        if ((k > kappaMin && k < kappaMax) || any_equal(k, c(kappaMin, kappaMax))) {
47 | 1x |
          setFound = set
48 | 1x |
          validSet = T
49 |  |
        }
50 |  |
51 | 51x |
        Ks = c(Ks, k)
52 | 51x |
        triesLeft = triesLeft - 1
53 |  |
      }
54 |  |
55 | 2x |
      if (is.null(setFound)) {
56 | 1x |
        Ks_under_min <- Ks[which(Ks < kappaMin)]
57 | 1x |
        Kdist.to.min <- kappaMin - Ks_under_min
58 | 1x |
        Kdist.to.min.low <- ifelse(length(Ks_under_min) > 0, min(Kdist.to.min), NA)
59 |  |
60 | 1x |
        Ks_over_min  <- Ks[which(Ks > kappaMax)]
61 | 1x |
        Kdist.to.max <- Ks_over_min - kappaMax
62 | 1x |
        Kdist.to.max.low <- ifelse(length(Ks_over_min) > 0, min(Kdist.to.max), NA)
63 |  |
64 | 1x |
        Kdist.all_mins <- c(Kdist.to.min, Kdist.to.max)
65 |  |
66 | 1x |
        both_mins <- c(Kdist.to.min.low, Kdist.to.max.low)
67 | 1x |
        both_mins <- both_mins[!is.na(both_mins)]
68 | 1x |
        closest <- c(Ks_under_min, Ks_over_min) [which(Kdist.all_mins == min(both_mins))][1]
69 |  |
70 | 1x |
        stop(paste0(
71 | 1x |
          "Unable to create a set of the provided length (", length, ") ",
72 | 1x |
            "in the provided kappa range (", kappaMin, "\U2014", kappaMax, "). ",
73 | 1x |
            "The closest kappa obtained after ", tries, " tries was: ",
74 | 1x |
            round(x = closest, digits = 5), "\n",
75 | 1x |
            "\tTry increasing the set length or the range of valid kappas"
76 |  |
        ))
77 |  |
      }
78 |  |
79 | 1x |
      return(setFound)
80 |  |
81 |  |
      # return(createRandomSet(setLength = length, baserate = baserate, kappaMin = kappaMin, kappaMax = kappaMax, minPrecision = precisionMin, maxPrecision = precisionMax, ...))
82 |  |
    }

1 |  |
    ###
2 |  |
    # Generate Precision Combo
3 |  |
    #
4 |  |
    # This function generates a precision that fits within the range and will create a valid set
5 |  |
    # @param kappa This is the kappa of the desired pair
6 |  |
    # @param Pmin This is the minimum precision desired
7 |  |
    # @param Pmax This is the maximum precision desired
8 |  |
    # @param BR This is the baserate of the desired pair
9 |  |
    # @keywords generate, combo, kappa, precision, baserate
10 |  |
    # @return This returns a precision value that will generate a valid set with the other parameters
11 |  |
    ####
12 |  |
13 |  |
    genPcombo = function(kappa, Pmin, Pmax, BR){
14 |  |
      #gets a random precision in the precision range
15 | 821x |
      currPrecision = stats::runif(1, Pmin, Pmax);
16 |  |
      #checks is the combination of precision, kappa, and baserate is valid.  If it is vaid, return the precision, if not, recall the function until a valid precision is found
17 | 821x |
      if(!checkBRPKcombo(BR, currPrecision, kappa)){
18 | 819x |
        genPcombo(kappa, Pmin, Pmax, BR);
19 |  |
      }else{
20 | 1x |
        return(currPrecision);
21 |  |
      }
22 |  |
    }

1 |  |
    ###
2 |  |
    # Precision Recall Set
3 |  |
    #
4 |  |
    # This creates a set with a given recall and precision
5 |  |
    # @param precision This is the precision of the given set
6 |  |
    # @param recall This is the recall of the given set
7 |  |
    # @param length This is the length of the set
8 |  |
    # @param baserate This is the baserate of the set
9 |  |
    # @keywords prset, precision, recall, generation, set
10 |  |
    # @return returns a set with the given parameters
11 |  |
    ###
12 |  |
    prset = function(precision, recall, length, baserate) {
13 |  |
      #generates a set given the different 4 quadrants
14 | 28354x |
      gold1s = max(round(baserate * length , 0), 1); # TP + FN
15 | 28354x |
      gold0s = length - gold1s; # FP + FN
16 | 28354x |
      TP = max(round(gold1s * recall , 0) , 1); # TP
17 | 28354x |
      FP = min(gold0s, max(round(TP * (1 - precision) / precision , 0), 1)); # FP
18 | 28354x |
      gold = c(rep(1, gold1s), rep(0, gold0s));
19 | 28354x |
      silver = c(rep(1, TP),rep(0, gold1s - TP), rep(1, FP), rep(0, gold0s - FP));
20 |  |
21 | 28354x |
      return(cbind(gold, silver));
22 |  |
    }
23 |  |
24 |  |
    pr.ct = function(precision, recall, length, baserate) {
25 |  |
      #generates a set given the different 4 quadrants
26 | 1x |
      gold1s = max(round(baserate * length , 0), 1); # TP + FN
27 | 1x |
      gold0s = length - gold1s; # FP + FN
28 | 1x |
      TP = max(round(gold1s * recall , 0) , 1); # TP
29 | 1x |
      FP = min(gold0s, max(round(TP * (1 - precision) / precision , 0), 1)); # FP
30 |  |
31 | 1x |
      ct = matrix(0, ncol = 2, nrow = 2)
32 | 1x |
      ct[1,1] = TP
33 | 1x |
      ct[1,2] = gold1s - TP
34 | 1x |
      ct[2,1] = FP
35 | 1x |
      ct[2,2] = gold0s - FP
36 | 1x |
      return(ct);
37 |  |
    }

1 |  |
    ###
2 |  |
    #' @include rhoR.R
3 |  |
    #'
4 |  |
    #' @title Rho
5 |  |
    #'
6 |  |
    #' @description
7 |  |
    #' This function calculates rho for a \code{\link[=getTestSet]{testSet}}, \code{\link{contingencyTable}}, or an observed kappa value with associated set parameters (testSetLength and OcSBaserate).
8 |  |
    #'
9 |  |
    #' @details
10 |  |
    #' Rho is a Monte Carlo rejective method of interrater reliability statistics, implemented here for Cohen's Kappa. Rho constructs a collection of data sets in which kappa is below a specified threshold, and computes the empirical distribution on kappa based on the specified sampling procedure. Rho returns the percent of the empirical distribution greater than or equal to an observed kappa. As a result, Rho quantifies the type 1 error in generalizing from an observed test set to a true value of agreement between two raters.
11 |  |
    #'
12 |  |
    #' Rho starts with an observed kappa value, calculated on a subset of a \code{\link{codeSet}}, known as an observed \code{\link[=getTestSet]{testSet}}, and a \emph{kappa threshold} which indicates what is considered significant agreement between raters.
13 |  |
    #'
14 |  |
    #' It then generates a collection of fully-coded, simulated
15 |  |
    #' \code{\link[=getTestSet]{codeSets}} (ScS), further described in
16 |  |
    #' \code{\link{createSimulatedCodeSet}}, all of which have a kappa value below the kappa
17 |  |
    #' threshold and similar properties as the original \code{\link{codeSet}}.
18 |  |
    #'
19 |  |
    #' Then, kappa is calculated on a \code{\link[=getTestSet]{testSet}} sampled from each of the ScSs in the
20 |  |
    #' collection to create a null hypothesis distribution. These \code{\link[=getTestSet]{testSets}} mirror the observed \code{\link[=getTestSet]{testSets}} in their size and sampling method.  How these \code{\link[=getTestSet]{testSets}} are sampled is futher described in \code{\link{getTestSet}}.
21 |  |
    #'
22 |  |
    #' The null hypothesis is that the observed \code{\link[=getTestSet]{testSet}}, was sampled from a data set, which, if both raters were to code in its entirety, would result in a level of agreement below the kappa threshold.
23 |  |
    #'
24 |  |
    #' For example, using an alpha level of 0.05, if the observed kappa is greater than 95 percent of the kappas in the null hypothesis distribution, the null hypothesis is rejected. Then one can conclude that the two raters would have acceptable agreement had they coded the entire data set.
25 |  |
    #'
26 |  |
    #' @export
27 |  |
    #'
28 |  |
    #' @template rho-params
29 |  |
    #' @param testSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
30 |  |
    #'
31 |  |
    #' @seealso \code{\link{rhoK}}
32 |  |
    #'
33 |  |
    #' @examples
34 |  |
    #' # Given an observed kappa value
35 |  |
    #' rho(x = 0.88, OcSBaserate = 0.2, testSetLength = 80)
36 |  |
    #'
37 |  |
    #' # Given a test Set
38 |  |
    #' rho(x = codeSet)
39 |  |
    #'
40 |  |
    #' # Given a contingency Table
41 |  |
    #' rho(x = contingencyTable)
42 |  |
    #'
43 |  |
    #' @return rho and kappa for the given data and parameters (unless kappa is given)
44 |  |
    ###
45 |  |
46 |  |
    rho = function(
47 |  |
      x,
48 |  |
      OcSBaserate = NULL,
49 |  |
      testSetLength = NULL,
50 |  |
      testSetBaserateInflation = 0,
51 |  |
      OcSLength = 10000,
52 |  |
      replicates = 800,
53 |  |
      ScSKappaThreshold = .65,
54 |  |
      ScSKappaMin = .40,
55 |  |
      ScSPrecisionMin = 0.6,
56 |  |
      ScSPrecisionMax = 1
57 |  |
    ) {
58 |  |
59 |  |
      ### Calculate rho given a kappa value
60 | 10x |
      if (class(x) == "numeric") {
61 | 1x |
        if (is.null(OcSBaserate)) { stop("Must give a baserate when a kappa value is given") }
62 | 1x |
        if (is.null(testSetLength)) { stop("Must give a testSetLength when a kappa value is given") }
63 |  |
64 | 6x |
        result <- rhoK(x, OcSBaserate,
65 | 6x |
                   testSetLength,
66 | 6x |
                   testSetBaserateInflation,
67 | 6x |
                   OcSLength,
68 | 6x |
                   replicates,
69 | 6x |
                   ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
70 |  |
      }
71 |  |
72 |  |
      ### Calculate rho given a code set
73 | 2x |
      else if (nrow(x) != 2) {
74 | 1x |
        result <- rhoSet(x, OcSBaserate,
75 | 1x |
                      testSetBaserateInflation,
76 | 1x |
                      OcSLength,
77 | 1x |
                      replicates,
78 | 1x |
                      ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
79 |  |
      }
80 |  |
81 |  |
      ### Calculate rho given a contingency table
82 |  |
      else {
83 | 1x |
        result <- rhoCT(x, OcSBaserate,
84 | 1x |
                     testSetBaserateInflation,
85 | 1x |
                     OcSLength,
86 | 1x |
                     replicates,
87 | 1x |
                     ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
88 |  |
      }
89 |  |
90 | 8x |
      return(result)
91 |  |
    }

1 |  |
    ###
2 |  |
    #' @include rhoMin.R
3 |  |
    #' @title Calculate Baserate
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function calculates the baserate of the first rater, second rater, and the average of both the raters.
7 |  |
    #'
8 |  |
    #' @details
9 |  |
    #' A baserate is the percentage, as a decimal, that a positive code appears in data (either a \code{\link{codeSet}} or \code{\link{contingencyTable}}) for a given rater.  It is assumed that the first rater is more experienced and thus provides a better estimation of the actual baserate for a given code, so the first rater's baserate is often used as if it is the actual baserate.  If the raters are assumed to have the same experience level, the average baserate may give a better estimation.  If the second rater is more experienced, the second rater's baserate may give a better estimation.  Functions assume that the first rater is the more experienced rater and thus uses the first rater's baserate as the overall baserate estimation.
10 |  |
    #'
11 |  |
    #' @export
12 |  |
    #'
13 |  |
    #' @param data The \code{\link[=getTestSet]{testSet}} or \code{\link{contingencyTable}} for which the baserate is calculated
14 |  |
    #'
15 |  |
    #' @keywords baserate
16 |  |
    #'
17 |  |
    #' @seealso \code{\link{baserateSet}} and \code{\link{baserateCT}}
18 |  |
    #'
19 |  |
    #' @examples
20 |  |
    #' #Given a code set
21 |  |
    #' baserate(data = codeSet)
22 |  |
    #'
23 |  |
    #' #Given a contingency Table
24 |  |
    #' baserate(data = contingencyTable)
25 |  |
    #'
26 |  |
    #' @return A list of the format:\describe{
27 |  |
    #' \item{firstBaserate}{The percentage of the data for which a positive code, or a 1, appears in the first rater}
28 |  |
    #' \item{secondBaserate}{The percentage of the data for which a positive code, or a 1, appears in the second rater}
29 |  |
    #' \item{averageBaserate}{The average of the firstBaserate and secondBaserate.}
30 |  |
    #' }
31 |  |
    #'
32 |  |
    ###
33 |  |
    baserate = function(data){
34 | 7x |
      if(nrow(data) != 2){
35 | 4x |
        return(baserateSet(data))
36 |  |
      }else{
37 | 3x |
        return(baserateCT(data))
38 |  |
      }
39 |  |
    }

1 |  |
    ###
2 |  |
    #' @include rho.R
3 |  |
    #'
4 |  |
    #' @title Rho (kappa)
5 |  |
    #'
6 |  |
    #' @description
7 |  |
    #' This function calculates rho for an observed kappa value with associated set parameters
8 |  |
    #' (testSetLength and OcSBaserate). Called by \code{\link{rho}}. A p-value is returned and if
9 |  |
    #' this value is less than 0.05, it is said that the handset does generalize to the entire set
10 |  |
    #'
11 |  |
    #' @export
12 |  |
    #'
13 |  |
    #' @template rho-params
14 |  |
    #' @param testSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
15 |  |
    #' @param method set to "c" to calculate using the C++ implmentation. Defaults to "standard"
16 |  |
    #'
17 |  |
    #' @return rho for the given parameters
18 |  |
    ###
19 |  |
    rhoK = function(
20 |  |
      x,
21 |  |
      OcSBaserate,
22 |  |
      testSetLength,
23 |  |
      testSetBaserateInflation = 0,
24 |  |
      OcSLength = 10000,
25 |  |
      replicates = 800,
26 |  |
      ScSKappaThreshold = .65,
27 |  |
      ScSKappaMin = .40,
28 |  |
      ScSPrecisionMin = 0.6,
29 |  |
      ScSPrecisionMax = 1,
30 |  |
      method = "standard"
31 |  |
    ) {
32 |  |
      if (
33 | 45x |
        (x > 1 | x < 0) ||
34 | 45x |
        (OcSBaserate > 1 | OcSBaserate < 0)
35 | 3x |
      ) stop("x and OcSBaserate must be between 0 and 1.")
36 |  |
37 | 42x |
      if (OcSBaserate < 0 || ScSKappaThreshold < 0 || ScSKappaMin < 0)
38 | 2x |
        stop("OcSBaserate, ScSKappaThreshold, and ScSKappaMin must be positive.")
39 |  |
40 |  |
      if (
41 | 40x |
        testSetLength %% 1 != 0 ||
42 | 40x |
        OcSLength %% 1 != 0 ||
43 | 40x |
        replicates %% 1 != 0
44 | 3x |
      ) stop("testSetLength, OcSLength, and replicates values must be a positive integer.")
45 |  |
46 | 37x |
      if (ScSKappaThreshold < ScSKappaMin)
47 | 1x |
        stop("ScSKappaThreshold must be greater than ScSKappaMin.")
48 |  |
49 |  |
      if (
50 | 36x |
        (ScSPrecisionMin < 0 | ScSPrecisionMin > 1) ||
51 | 36x |
        (ScSPrecisionMax < 0 | ScSPrecisionMax > 1)
52 | 4x |
      ) stop("ScSPrecisionMin and ScSPrecisionMax must be between 0 and 1.")
53 |  |
54 | 32x |
      if (ScSPrecisionMax < ScSPrecisionMin)
55 | 1x |
        stop("ScSPrecisionMax must be greater than ScSPrecisionMin.")
56 |  |
57 | 31x |
      result <- NULL
58 | 31x |
      if (method == "standard") {
59 | 21x |
        result <- calcRho(  x, OcSBaserate, testSetLength, testSetBaserateInflation, OcSLength, replicates, ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
60 |  |
      }
61 |  |
      else {
62 | 10x |
        result <- calcRho_c(x, OcSBaserate, testSetLength, testSetBaserateInflation, OcSLength, replicates, ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax)
63 |  |
      }
64 | 31x |
      return(result)
65 |  |
    }

1 |  |
    contingencyToSet = function(TP, FP, FN, TN){
2 | 12x |
      setLength = TP + FP + FN + TN;
3 |  |
4 | 12x |
      gold1s = TP + FN;
5 | 12x |
      gold0s = setLength - gold1s;
6 |  |
7 | 12x |
      gold = c(rep(1, gold1s), rep(0, gold0s));
8 | 12x |
      silver = c(rep(1, TP),rep(0, gold1s - TP), rep(1, FP), rep(0, gold0s - FP));
9 |  |
10 | 12x |
      return(cbind(gold, silver));
11 |  |
    };

1 |  |
    genPKcombo = function(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate){
2 |  |
      # currKappa = sample(kappaDistribution, 1, replace = FALSE, prob = kappaProbability)
3 |  |
      # currPrecision = sample(precisionDistribution, 1, replace = FALSE, prob = precisionProbability)
4 | 36803x |
      currKappa = kappaDistribution[sample.int(length(kappaDistribution), 1, replace = FALSE, prob = kappaProbability)]
5 | 36803x |
      currPrecision = precisionDistribution[sample.int(length(precisionDistribution), 1, replace = FALSE, prob = precisionProbability)]
6 |  |
7 | 36803x |
      if(!checkBRPKcombo(baserate, currPrecision, currKappa)){
8 | 1850x |
        precisionMin = (2*baserate*currKappa - 2*baserate - currKappa)/(currKappa - 2)
9 | 1850x |
        indicies = which(precisionDistribution > precisionMin)
10 |  |
11 | 1850x |
        if (length(indicies) == 0) {
12 | 1838x |
          return(genPKcombo(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate))
13 |  |
        }
14 |  |
15 | 12x |
        precisionDistribution = precisionDistribution[indicies]
16 | 12x |
        precisionProbability = precisionProbability[indicies]
17 |  |
        # currPrecision = sample(precisionDistribution, 1, replace = FALSE, prob = precisionProbability)
18 | 12x |
        currPrecision = precisionDistribution[sample.int(length(precisionDistribution), 1, replace = FALSE, prob = precisionProbability)]
19 | 12x |
        return(list(precision = currPrecision, kappa = currKappa))
20 |  |
      } else {
21 | 34953x |
        return(list(precision = currPrecision, kappa = currKappa))
22 |  |
      }
23 |  |
    }

1 |  |
    #' Generate a Handset
2 |  |
    #'
3 |  |
    #' @description Generate a vector representing indices of set, using the handSetBaserate
4 |  |
    #' to determine the minimum number of indices that are positive
5 |  |
    #'
6 |  |
    #' @param set matrix of two columns
7 |  |
    #' @param handSetLength number of indices to find
8 |  |
    #' @param handSetBaserate number between 0 and 1 to use as a minimum number of positive indices
9 |  |
    #'
10 |  |
    #' @return vector of indices from set
11 |  |
    #' @export
12 |  |
    getHandSetIndices = function(set, handSetLength = 20, handSetBaserate = 0.2) {
13 | 4x |
      positives = ceiling(handSetLength * handSetBaserate);
14 | 4x |
      posInd = which(set[, 1] == 1);
15 |  |
16 | 4x |
      if (positives > length(posInd)) {
17 | 1x |
        stop("Not enough positives in first rater to inflate to this level")
18 |  |
      }
19 |  |
20 | 3x |
      this.set = NULL
21 | 3x |
      if (positives > 0) {
22 | 2x |
        positiveIndices = posInd[sample.int(length(posInd), size=positives, replace=FALSE)]
23 | 2x |
        others = which(!(1:nrow(set) %in% positiveIndices))
24 | 2x |
        otherIndices = sample.int(length(others), size = (handSetLength - positives), replace = FALSE)
25 | 2x |
        this.set = c(positiveIndices, others[otherIndices])
26 |  |
      }
27 | 1x |
      else if (positives == 0) {
28 | 1x |
        theseIndices = sample.int(length(set),size=handSetLength,replace=FALSE)
29 | 1x |
        this.set = theseIndices
30 |  |
      }
31 |  |
32 | 3x |
      return(sample(this.set))
33 |  |
    }

1 |  |
    ###
2 |  |
    #' @include kappa.R
3 |  |
    #' @title Calculate kappa (Set)
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function calculates Cohen's kappa for a given \code{\link{codeSet}}. Called by \code{\link{kappa}}.
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #' @param set A \code{\link{codeSet}}
10 |  |
    #'
11 |  |
    #' @keywords kapppaSet
12 |  |
    #'
13 |  |
    #' @seealso \code{\link{kappa}} and \code{\link{kappaCT}}
14 |  |
    #'
15 |  |
    #' @return
16 |  |
    #' The kappa of the \code{\link{codeSet}}
17 |  |
    #'
18 |  |
    ###
19 |  |
    kappaSet = function(set){
20 | 62x |
      if(any(set != 1 & set != 0)) {
21 | 1x |
        stop("Each set must only consist of zeros and ones.")
22 |  |
      }
23 |  |
24 | 61x |
      return(calcKappa(set,isSet = TRUE,kappaThreshold = NULL));
25 |  |
    }
26 |  |

1 |  |
    ###
2 |  |
    # Get Boot P Value
3 |  |
    #
4 |  |
    # This function gets the p value of the result in the distribution
5 |  |
    # @param distribution The distribution that the result falls into
6 |  |
    # @param result The value that you are evaluating where in the distribution it falls
7 |  |
    # @keywords distribution, rho, p, boot
8 |  |
    # @return This returns the p value from the result into the distribution
9 |  |
    ###
10 |  |
    getBootPvalue = function(distribution, result) {
11 |  |
      #returns the percentage of the time that the distribution was greater or equal to the observed kappa
12 | 94x |
      N = length(distribution);
13 | 94x |
      if(result < mean(distribution, na.rm = T)) {
14 |  |
        #if the result is less than the mean of the distribution, than the p value is 1
15 | 1x |
        return(1);
16 |  |
      } else {
17 |  |
        #return the number of times that the distribution is greater than the result as a percentage of the total number of items in the distribution
18 | 93x |
        return(length(which(distribution >= result)) / N);
19 |  |
      }
20 |  |
    }

1 |  |
    ###
2 |  |
    #' @include baserateSet.R
3 |  |
    #' @title Calculate Baserate (CT)
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function calculates the baserate of the first rater, second rater, and the average of both the raters. Called by \code{\link{baserate}}.
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #'
10 |  |
    #' @param CT The \code{\link{contingencyTable}} for which the baserate is calculated
11 |  |
    #'
12 |  |
    #' @keywords baserateCT
13 |  |
    #'
14 |  |
    #' @seealso \code{\link{baserate}} and \code{\link{baserateSet}}
15 |  |
    #'
16 |  |
    #' @return A list of the format:\describe{
17 |  |
    #' \item{firstBaserate}{The percentage of the data for which a positive code, or a 1, appears in the first rater}
18 |  |
    #' \item{secondBaserate}{The percentage of the data for which a positive code, or a 1, appears in the second rater}
19 |  |
    #' \item{averageBaserate}{The average of the firstBaserate and secondBaserate.}
20 |  |
    #' }
21 |  |
    #'
22 |  |
    ###
23 |  |
24 |  |
    baserateCT = function(CT){
25 | 1x |
      if(any(CT < 0)){stop("Values in Contingency Table must be positive")}
26 | 3x |
      firstBaserate = sum(CT[1,])/sum(CT)
27 | 3x |
      secondBaserate = sum(CT[,1])/sum(CT)
28 | 3x |
      averageBaserate = mean(c(firstBaserate,secondBaserate))
29 | 3x |
      return(list(firstBaserate = firstBaserate, secondBaserate = secondBaserate, averageBaserate = averageBaserate))
30 |  |
    }

1 |  |
    ###
2 |  |
    #' @include kappaSet.R
3 |  |
    #' @title Calculate kappa (contingency Table)
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function calculates Cohen's kappa on a \code{\link{contingencyTable}}. Called by \code{\link{kappa}}.
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #' @param ct A \code{\link{contingencyTable}}
10 |  |
    #'
11 |  |
    #' @keywords kappaCT
12 |  |
    #'
13 |  |
    #' @seealso \code{\link{kappa}} and \code{\link{kappaSet}}
14 |  |
    #'
15 |  |
    #' @return
16 |  |
    #' The kappa of the \code{\link{contingencyTable}}
17 |  |
    ###
18 |  |
    kappaCT = function(ct){
19 | 1x |
      if(any(ct < 0)){stop("Values in Contingency Table must be positive")}
20 | 1x |
      if(ncol(ct) != 2 || nrow(ct) != 2) stop("Incorrect number of dimensions: Contingency table must be 2x2.")
21 | 7x |
      for(i in 1:nrow(ct)) {
22 | 13x |
        for(j in 1:ncol(ct)) {
23 | 1x |
          if (ct[i,j]%%1 != 0) stop("Contingency table values must be positive integers.")
24 |  |
        }
25 |  |
      }
26 |  |
27 | 6x |
      return(calcKappa(ct,isSet = FALSE,kappaThreshold = NULL));
28 |  |
    }

1 |  |
    generateKPs = function(
2 |  |
      numNeeded, baserate,
3 |  |
      kappaMin, kappaMax,
4 |  |
      precisionMin, precisionMax,
5 |  |
      distributionType = 'FLAT',
6 |  |
      distributionLength = 100000
7 |  |
    ) {
8 | 89x |
      kappaDistribution = seq(kappaMin, kappaMax, length = distributionLength)
9 |  |
10 | 89x |
      if (distributionType == 'FLAT') {
11 |  |
        #probability for a flat distribution
12 | 88x |
        kappaProbability = NULL
13 |  |
      }
14 | 1x |
      else if (distributionType == 'BELL') {
15 |  |
        #probability for a bell curve
16 | 1x |
        kappaProbability = stats::dnorm(kappaDistribution, mean = 0.65, sd = 0.1)
17 |  |
      }
18 |  |
19 | 89x |
      precisionDistribution = seq(precisionMin, precisionMax, length = 10000)
20 | 89x |
      precisionProbability = NULL
21 |  |
22 | 89x |
      KPs = tryCatch({
23 | 89x |
        KPs = list()
24 | 89x |
        for (i in 1:numNeeded) {
25 | 34965x |
          KPs[[length(KPs) + 1]] = genPKcombo(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate)
26 |  |
        }
27 | 88x |
        KPs
28 | 89x |
      },error=function(e){
29 | 1x |
        NULL
30 |  |
      })
31 | 1x |
      if(is.null(KPs)){stop("Could not generate a valid set given ranges of kappa and precision")}
32 |  |
33 | 88x |
      return (KPs);
34 |  |
    }

1 |  |
    ###
2 |  |
    # This function calculates rho given the parameters of the initial set and the kappa of the observed handset and is
3 |  |
    # called by rhoK()
4 |  |
    ###
5 |  |
    calcRho = function(
6 |  |
      x,
7 |  |
      OcSBaserate,
8 |  |
      testSetLength,
9 |  |
      testSetBaserateInflation = 0,
10 |  |
      OcSLength = 10000,
11 |  |
      replicates = 800,
12 |  |
      ScSKappaThreshold = 0.65,
13 |  |
      ScSKappaMin = 0.40,
14 |  |
      ScSPrecisionMin = 0.6,
15 |  |
      ScSPrecisionMax = 1,
16 |  |
      KPs = NULL
17 |  |
    ) {
18 |  |
19 | 91x |
      if(is.null(KPs)) {
20 | 31x |
        KPs = generateKPs(replicates, OcSBaserate, ScSKappaMin, ScSKappaThreshold, ScSPrecisionMin, ScSPrecisionMax)
21 |  |
      }
22 | 91x |
      if(length(KPs) < replicates) {
23 | 10x |
        replicates = length(KPs)
24 |  |
      }
25 |  |
26 |  |
      #creates distribution from random kappa sets
27 | 91x |
      savedKappas = rep(0,replicates);
28 |  |
29 |  |
      # cat("Runs: ", runs, "\n")
30 |  |
      # cat("KPs length:, ", length(KPs), "\n")
31 | 91x |
      for (i in 1:replicates) {
32 | 28300x |
        KP = KPs[i]
33 | 28300x |
        kappa = KP[[1]]$kappa
34 | 28300x |
        precision = KP[[1]]$precision
35 | 28300x |
        recall = getR(kappa, OcSBaserate, precision)
36 | 28300x |
        fullSet = prset(precision = precision, recall = recall, length = OcSLength, baserate = OcSBaserate);
37 |  |
38 |  |
        #computes kappa of handset to add to distribution
39 | 28300x |
        savedKappas[i] = getHandSet(set = fullSet, handSetLength = testSetLength, handSetBaserate = testSetBaserateInflation);
40 |  |
      }
41 |  |
42 |  |
      # compare observedKappa to distribution
43 | 91x |
      return(getBootPvalue(savedKappas, x));
44 |  |
    }

1 |  |
    ###
2 |  |
    # Create A Random Set
3 |  |
    #
4 |  |
    # This creates a set with given parameters
5 |  |
    # @param setLength The length of the set to be created
6 |  |
    # @param baserate The baserate of the set to be created
7 |  |
    # @param kappa The kappa of the set to be created
8 |  |
    # @param minPrecision The minimum precision accepted in this set. Defaults to 0
9 |  |
    # @param maxPrecision The maximum precision accepted in this set. Defaults to 1
10 |  |
    # @keywords random, set, generation, kappa, ranges
11 |  |
    # @return This returns a set with the given parameters
12 |  |
    ###
13 |  |
    createRandomSet = function(setLength, baserate, kappaMin, kappaMax, minPrecision = 0, maxPrecision = 1, type = "set"){
14 |  |
      #generates a valid precision with the given kappa and baserate
15 | 54x |
      KP = generateKPs(1, baserate, kappaMin, kappaMax, minPrecision, maxPrecision)[[1]]
16 | 54x |
      kappa = KP$kappa
17 | 54x |
      precision = KP$precision
18 | 54x |
      recall = getR(kappa, baserate, precision)
19 |  |
20 | 54x |
      if(type == "ct") {
21 | 1x |
        fullSet = pr.ct(precision = precision, recall = recall, length = setLength, baserate = baserate);
22 |  |
      }
23 |  |
      else {
24 | 53x |
        fullSet = prset(precision = precision, recall = recall, length = setLength, baserate = baserate);
25 |  |
      }
26 |  |
27 | 54x |
      return(fullSet);
28 |  |
    }

1 |  |
    ###
2 |  |
    #' @include rhoMin.R
3 |  |
    #' @title Calculate Baserate (Set)
4 |  |
    #'
5 |  |
    #' @description
6 |  |
    #' This function will calculate the baserate of the first rater, second rater, and the average of both the raters. Called by \code{\link{baserate}}.
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #'
10 |  |
    #' @param set The \code{\link{codeSet}} for which the baserate is calculated
11 |  |
    #'
12 |  |
    #' @keywords baserateSet
13 |  |
    #'
14 |  |
    #' @seealso \code{\link{baserate}} and \code{\link{baserateCT}}
15 |  |
    #'
16 |  |
    #' @return A list of the format:\describe{
17 |  |
    #' \item{firstBaserate}{The percentage that a positive code, or a 1, appears in the first rater}
18 |  |
    #' \item{secondBaserate}{The percentage that a positive code, or a 1, appears in the second rater}
19 |  |
    #' \item{averageBaserate}{The average percentage that a positive code, or a 1, appears in either of the two raters}
20 |  |
    #' }
21 |  |
    #'
22 |  |
    ###
23 |  |
    baserateSet = function(set){
24 | 4x |
      firstBaserate = length(which(set[,1] == 1))/nrow(set)
25 | 4x |
      secondBaserate = length(which(set[,2] == 1))/nrow(set)
26 | 4x |
      averageBaserate = mean(c(firstBaserate,secondBaserate))
27 | 4x |
      return(list(firstBaserate = firstBaserate, secondBaserate = secondBaserate, averageBaserate = averageBaserate))
28 |  |
    }

1 |  |
    ###
2 |  |
    #' @title Rho using a file
3 |  |
    #'
4 |  |
    #' @description
5 |  |
    #' This function calculates rho and kappa for a given \code{\link[=getTestSet]{testSet}} as defined by the file
6 |  |
    #' and columns (col1, col2), and returns a list containing both values. Called by \code{\link{rho}}.
7 |  |
    #'
8 |  |
    #' @export
9 |  |
    #'
10 |  |
    #' @param col1 The first column from file
11 |  |
    #' @param col2 The second column from file
12 |  |
    #' @template rho-params
13 |  |
    #'
14 |  |
    #' @return A list of the format:\describe{
15 |  |
    #' \item{rho}{The rho of the \code{\link{codeSet}}}
16 |  |
    #' \item{kappa}{The Cohen's Kappa of the \code{\link{codeSet}}}
17 |  |
    #' }
18 |  |
    ###
19 |  |
    rho.file <- function(
20 |  |
      x, col1, col2,
21 |  |
      OcSBaserate = NULL,
22 |  |
      testSetBaserateInflation = 0,
23 |  |
      OcSLength = 10000,
24 |  |
      replicates = 800,
25 |  |
      ScSKappaThreshold = .65,
26 |  |
      ScSKappaMin = .40,
27 |  |
      ScSPrecisionMin = 0.6,
28 |  |
      ScSPrecisionMax = 1
29 |  |
    ) {
30 | 3x |
      setFile = utils::read.csv(x)
31 | 3x |
      set = as.matrix(setFile[,c(col1, col2)])
32 |  |
33 | 3x |
      if (!any(is.na(set))) {
34 | 2x |
        rhoSet(
35 | 2x |
          set,
36 | 2x |
          OcSBaserate,
37 | 2x |
          testSetBaserateInflation,
38 | 2x |
          OcSLength, replicates,
39 | 2x |
          ScSKappaThreshold, ScSKappaMin,
40 | 2x |
          ScSPrecisionMin, ScSPrecisionMax
41 |  |
        )
42 |  |
      }
43 |  |
      else {
44 | 1x |
        stop("The columns provided for col1 and col2 must not contain NA")
45 |  |
      }
46 |  |
    }

1 |  |
    any_equal <- function(target, current, tolerance = sqrt(.Machine$double.eps)) {
2 | 54x |
      any(sapply(target, function(targ) abs(current - targ) < tolerance ))
3 |  |
    }
4 |  |
    all_equal <- function(target, current, tolerance = sqrt(.Machine$double.eps)) {
5 | 5x |
      all(sapply(target, function(targ) abs(current - targ) < tolerance ))
6 |  |
    }
7 |  |
8 |  |
9 |  |
    #' Convert a codeset to a contingency table
10 |  |
    #'
11 |  |
    #' @param x codeset
12 |  |
    #'
13 |  |
    #' @return contingency table as a 2x2 matrix
14 |  |
    #'
15 |  |
    #' @export
16 |  |
    as.contingency.table <- function(x) {
17 | 1x |
      tp = length(which(x[,1] == 1 & x[,2] == 1))
18 | 1x |
      tn = length(which(x[,1] == 0 & x[,2] == 0))
19 | 1x |
      fp = length(which(x[,1] == 0 & x[,2] == 1))
20 | 1x |
      fn = length(which(x[,1] == 1 & x[,2] == 0))
21 |  |
22 | 1x |
      return(matrix(c(tp,fp,fn,tn), ncol = 2, nrow = 2))
23 |  |
    }
24 |  |
25 |  |
    #' Convert codeset to contingency table
26 |  |
    #'
27 |  |
    #' @param x matrix contingency table (2x2)
28 |  |
    #'
29 |  |
    #' @return 2-column matrix representation of the contingency table
30 |  |
    #'
31 |  |
    #' @export
32 |  |
    as.codeset <- function(x) {
33 | 1x |
      return(matrix(c(rep(1, x[1,1]*2), rep(1:0, x[1,2]), rep(0:1, x[2,1]), rep(0, x[2,2]*2)), byrow=TRUE, ncol=2))
34 |  |
    }

1 |  |
    ###
2 |  |
    # Calculate Kappa
3 |  |
    #
4 |  |
    # This function calculates kappa given a set
5 |  |
    # @param set This is the set to calculate kappa from
6 |  |
    # @param isSet TRUE if set FALSE if contigency table
7 |  |
    # @param kappaThreshold if null then return kappa otherwise return list of kappa and above or not
8 |  |
    # @keywords kappa, calculate
9 |  |
    # @return This returns the kappa of the set given
10 |  |
    ###
11 |  |
    calcKappa <- function(set, isSet = TRUE, kappaThreshold = NULL) {
12 | 28374x |
      if (!isSet) {
13 | 12x |
        set <- contingencyToSet(set[1, 1], set[2, 1], set[1, 2], set[2, 2])
14 |  |
      }
15 |  |
16 |  |
      if (
17 | 28374x |
        length(which(set[, 1] == set[, 2] & set[, 1] == 1)) == nrow(set) |
18 | 28374x |
        length(which(set[, 1] == set[, 2] & set[, 1] == 0)) == nrow(set)
19 |  |
      ) {
20 | 939x |
        return(1)
21 |  |
      }
22 |  |
23 |  |
      #calculates kappa by calculating the agreement and the baserates and then creating the adjacenty matrix
24 | 27435x |
      agreement <- length(which(set[, 1] == set[, 2])) / nrow(set);
25 | 27435x |
      baserate2 <- length(which(1 == set[, 2])) / nrow(set);
26 | 27435x |
      baserate1 <- length(which(1 == set[, 1])) / nrow(set);
27 | 27435x |
      randomAgreement <- baserate1 * baserate2 + (1 - baserate1) * (1 - baserate2);
28 | 27435x |
      kappa <- (agreement - randomAgreement) / (1 - randomAgreement);
29 |  |
30 | 27435x |
      if (is.null(kappaThreshold)) {
31 | 27433x |
        return(kappa);
32 |  |
      } else {
33 | 2x |
        return(list(kappa = kappa, above = kappa > kappaThreshold))
34 |  |
      }
35 |  |
    }

1 |  |
    ###
2 |  |
    # Get Recall
3 |  |
    #
4 |  |
    # This gets the recall from a pair
5 |  |
    # @param kappa This is the kappa of the pair
6 |  |
    # @param BR This is the baserate of the pair
7 |  |
    # @param P this is the precision of the pair
8 |  |
    # @keywords recall, get, kappa, baserate, precision
9 |  |
    # @return returns the recall given the parameters
10 |  |
    ###
11 |  |
    getR = function(kappa, BR, P){
12 |  |
      #gets recall that corresponds with the kappa, baserate, and precision
13 | 28355x |
      top = kappa*P;
14 | 28355x |
      R = top / (2*P-2*BR-kappa+2*BR*kappa);
15 | 28355x |
      return (R);
16 |  |
    }

1 |  |
    ###
2 |  |
    # Check Baserate, Precision, Kappa Combo
3 |  |
    #
4 |  |
    # This function checks to make sure the combination of BR, P, and K will create a valid set
5 |  |
    # @param BR This is the supplied baserate
6 |  |
    # @param P This is the supplied precision
7 |  |
    # @param K This is the supplied kappa
8 |  |
    # @keywords check, combo, baserate, precision, kappa
9 |  |
    # @return returns TRUE if the combination is valid, FALSE otherwise
10 |  |
    ###
11 |  |
    checkBRPKcombo = function(BR, P, K) {
12 |  |
      #if right is less than P, then it is a valid combination. If not, then the pair of kappa, precision, and baserate does not correspond to a valid set
13 | 37625x |
      right = (2*BR*K - 2*BR - K)/(K - 2);
14 | 37625x |
      if(P > right){
15 | 34955x |
        return(TRUE);
16 |  |
      }else{
17 | 2670x |
        return(FALSE);
18 |  |
      }
19 |  |
    }

1 |  |
    #include <RcppArmadilloExtensions/sample.h>
2 |  |
    // [[Rcpp::depends(RcppArmadillo)]]
3 |  |
4 |  |
    #include <RcppArmadillo.h>
5 |  |
    #include <Rcpp.h>
6 |  |
    using namespace Rcpp;
7 |  |
    using namespace arma;
8 |  |
9 |  |
    // NumericMatrix toNumericMatrix(arma::imat x) {
10 |  |
    //   int nRows=x.n_rows;
11 |  |
    //
12 |  |
    //   NumericMatrix y(nRows, x.n_cols);
13 |  |
    //   for (int j=0; j<x.n_rows; j++) {
14 |  |
    //     for (int i=0; i<x.n_cols; i++) {
15 |  |
    //       y(j,i) = x(j, i);
16 |  |
    //     }
17 |  |
    //   }
18 |  |
    //
19 |  |
    //   return y;
20 |  |
    // }
21 |  |
22 |  |
    // // @name sample_weighted
23 |  |
    // // @title sample_weighted
24 |  |
    // // @param x vector
25 |  |
    // // @param n integer
26 |  |
    // // @export
27 |  |
    // //' @export
28 |  |
    // // [[Rcpp::export]]
29 |  |
    // arma::vec sample_weighted(arma::vec xx, arma::vec probs, int n) {
30 |  |
    //   // NumericVector ret(n);
31 |  |
    //   // NumericVector xx(x);
32 |  |
    //   arma::vec ret(n);
33 |  |
    //
34 |  |
    //   // IntegerVector rng = seq(0, xx.size()-1);
35 |  |
    //   arma::ivec rng = regspace<ivec>(0, xx.size() - 1);
36 |  |
    //   for(int w = 0; w < n; w++) {
37 |  |
    //     // int tot = sum(xx);
38 |  |
    //     // NumericVector probs = xx / tot;
39 |  |
    //     // arma::vec probs2 = toArmaVec(probs);
40 |  |
    //   //   IntegerVector s = Rcpp::sugar::sample(rng, 1, false,  probs);
41 |  |
    //     arma::ivec s = Rcpp::RcppArmadillo::sample(rng, 1, false, probs);
42 |  |
    //     xx[s[0]] -= 1;
43 |  |
    //     ret[w] = s[0];
44 |  |
    //   }
45 |  |
    //   // return(xx);
46 |  |
    //   return(ret + 1);
47 |  |
    // }
48 |  |
49 |  |
    // // ' @export
50 |  |
    // // [[Rcpp::export]]
51 |  |
    // arma::ivec sample_weighted2(arma::imat xx, int n, bool forR = true) {
52 |  |
    //   arma::ivec ret = arma::ivec(n);
53 |  |
    //   NumericMatrix x = toNumericMatrix(xx);
54 |  |
    //   IntegerVector rng = seq(0, x.size() - 1);
55 |  |
    //
56 |  |
    //   for(int w = 0; w < n; w++) {
57 |  |
    //     double tot = accu(xx);
58 |  |
    //     NumericVector probs = x / tot;
59 |  |
    //     IntegerVector s = sample(rng, 1, false,  probs);
60 |  |
    //     x[s[0]] -= 1;
61 |  |
    //     ret[w] = s[0];
62 |  |
    //   }
63 |  |
    //
64 |  |
    //   if(forR == true) {
65 |  |
    //     return(ret + 1);
66 |  |
    //   } else {
67 |  |
    //     return(ret);
68 |  |
    //   }
69 |  |
    // }
70 |  |
71 |  |
    //' @title sample_contingency_table
72 |  |
    //' @name sample_contingency_table
73 |  |
    //'
74 |  |
    //' @param xx contingency table matrix
75 |  |
    //' @param n int size of the contingency table
76 |  |
    //' @param forR bool if true, add 1 to the results accounting for R indices starting at 1
77 |  |
    //'
78 |  |
    //' @export
79 |  |
    // [[Rcpp::export]]
80 |  |
    arma::ivec sample_contingency_table(arma::imat xx, int n, bool forR = true) {
81 | 31000x |
      arma::ivec ret = arma::ivec(n);
82 |  |
83 |  |
      int maxInt;
84 |  |
      int s;
85 | 3182000x |
      for(int w = 0; w < n; w++) {
86 | 3120000x |
        maxInt = accu(xx);
87 | 4680000x |
        s = randi<int>(distr_param(0, maxInt));
88 | 16678539x |
        ret[w] = (s <= xx(0) ? 0 : (s <= (xx(0)+xx(1))) ? 1 : (s <= (xx(0)+xx(1)+xx(2))) ? 2 : 3);
89 | 4680000x |
        xx(ret[w])--;
90 |  |
      }
91 |  |
92 | 31000x |
      if(forR == true) return(ret + 1);
93 | 62000x |
      else return(ret);
94 |  |
    }
95 |  |
96 |  |
    //' @name getBootPvalue_c
97 |  |
    //' @title getBootPvalue_c
98 |  |
    //'
99 |  |
    //' @param distribution vector of calculated kappas
100 |  |
    //' @param result double calculated kappa to compare against
101 |  |
    //'
102 |  |
    //' @description
103 |  |
    //' returns the percentage of the time that the distribution was greater or equal to the observed kappa
104 |  |
    //' if the result is less than the mean of the distribution, than the p value is 1
105 |  |
    //' else return the number of times that the distribution is greater than the result as a percentage of the total number
106 |  |
    //' of items in the distribution
107 |  |
    //'
108 |  |
    //' @return double calculated p-value
109 |  |
    //'
110 |  |
    //' @export
111 |  |
    // [[Rcpp::export]]
112 |  |
    double getBootPvalue_c(arma::vec distribution, double result) {
113 | 82x |
      if(result < mean(distribution)) {
114 | 1x |
        return(1.0);
115 |  |
      }
116 |  |
117 |  |
      else {
118 | 243x |
        arma::uvec matched = find(distribution >= result);
119 |  |
120 | 243x |
        double rho = (matched.size() * 1.0) / distribution.size();
121 | 81x |
        return( rho );
122 |  |
      }
123 |  |
    }
124 |  |
125 |  |
    // Validate a Baserate/Kappa combo against the supplied Precision.
126 |  |
    // Return true if Precision is greater than the calculated right value
127 |  |
    // [[Rcpp::export]]
128 |  |
    bool check_BRK_combo(double BR, double P, double K) {
129 | 37079x |
      double right = ((2 * BR * K) - (2 * BR) - K) / (K - 2);
130 |  |
131 | 37079x |
      return(P > right);
132 |  |
    }
133 |  |
134 |  |
    //' @name recall
135 |  |
    //' @title recall
136 |  |
    //'
137 |  |
    //' @param kappa double
138 |  |
    //' @param BR double
139 |  |
    //' @param P double
140 |  |
    //'
141 |  |
    //' @return Recall calculated from provided kappa, BR, and P
142 |  |
    //'
143 |  |
    //' @export
144 |  |
    // [[Rcpp::export]]
145 |  |
    double recall(double kappa, double BR, double P) {
146 | 19502x |
      double top = kappa * P;
147 | 19502x |
      double R = top / (2*P-2*BR-kappa+2*BR*kappa);
148 | 19502x |
      return(R);
149 |  |
    }
150 |  |
151 |  |
    // Find a valid precision and kappa combo from a distribution of Kappas and Precisions given a baserate
152 |  |
    //
153 |  |
    // [[Rcpp::export]]
154 |  |
    NumericVector find_valid_pk(
155 |  |
      arma::colvec kappaDistribution, arma::vec kappaProbability,
156 |  |
      arma::colvec precisionDistribution, arma::vec precisionProbability,
157 |  |
      double baserate
158 |  |
    ) {
159 | 37079x |
      arma::ivec kappaRng = regspace<arma::ivec>(0, kappaDistribution.size() - 1);
160 |  |
161 | 37079x |
      arma::ivec whKappa;
162 | 37079x |
      if (kappaProbability.n_elem > 0) {
163 | 32934x |
        whKappa = Rcpp::RcppArmadillo::sample(kappaRng, 1, false, kappaProbability);
164 |  |
      } else {
165 | 78303x |
        whKappa = Rcpp::RcppArmadillo::sample(kappaRng, 1, false);
166 |  |
      }
167 | 111237x |
      double currKappa = kappaDistribution[whKappa.at(0)];
168 |  |
169 | 74158x |
      arma::ivec precRng = regspace<arma::ivec>(0, precisionDistribution.size() - 1);
170 | 37079x |
      arma::ivec whPrec = Rcpp::RcppArmadillo::sample(precRng, 1, false);
171 | 111237x |
      double currPrec = precisionDistribution[whPrec.at(0)];
172 |  |
173 | 74158x |
      if (!check_BRK_combo(baserate, currPrec, currKappa)) {
174 | 978x |
        double precisionMin = (2 * baserate * currKappa - 2 * baserate - currKappa) / (currKappa - 2);
175 | 3912x |
        arma::uvec indicies = find(precisionDistribution > precisionMin);
176 |  |
177 | 1956x |
        if (indicies.size() == 0) {
178 | 4880x |
          return(find_valid_pk(kappaDistribution, kappaProbability, precisionDistribution, precisionProbability, baserate));
179 |  |
        }
180 |  |
181 | 6x |
        precisionDistribution = precisionDistribution.elem(indicies);
182 |  |
183 | 4x |
        if (precisionProbability.size() > 0) {
184 | 3x |
          precisionProbability = precisionProbability.elem(indicies);
185 |  |
        }
186 |  |
187 | 8x |
        precRng = regspace<ivec>(0, precisionDistribution.size() - 1);
188 | 6x |
        whPrec = Rcpp::RcppArmadillo::sample(precRng, 1, false);
189 | 6x |
        currPrec = precisionDistribution[whPrec.at(0)];
190 |  |
      }
191 |  |
192 | 72206x |
      return(Rcpp::NumericVector::create(currPrec, currKappa));
193 |  |
    }
194 |  |
195 |  |
196 |  |
    //' @name generateKPs_c
197 |  |
    //' @title generate_kp_list
198 |  |
    //'
199 |  |
    //' @param numNeeded int
200 |  |
    //' @param baserate double
201 |  |
    //' @param kappaMin double
202 |  |
    //' @param kappaMax double
203 |  |
    //' @param precisionMin double
204 |  |
    //' @param precisionMax double
205 |  |
    //' @param distributionType int 0 - normal (default), 1 - bell
206 |  |
    //' @param distributionLength long
207 |  |
    //'
208 |  |
    //' @return matrix of kappa and precision values (column 1 as precision)
209 |  |
    //'
210 |  |
    //' @export
211 |  |
    // [[Rcpp::export]]
212 |  |
    Rcpp::NumericMatrix generate_kp_list(
213 |  |
        int numNeeded, double baserate,
214 |  |
        double kappaMin, double kappaMax,
215 |  |
        double precisionMin, double precisionMax,
216 |  |
        int distributionType = 0, long distributionLength = 10000
217 |  |
    ) {
218 | 24x |
      double kappaStep = ((kappaMax - kappaMin) / (distributionLength - 1));
219 |  |
220 | 24x |
      colvec kappaDistribution = regspace(kappaMin, kappaStep, kappaMax);
221 |  |
222 | 24x |
      arma::vec kappaProbability;
223 | 24x |
      if(distributionType == 1) {
224 | 3x |
        kappaProbability = normpdf(kappaDistribution, 0.65, 0.1);
225 |  |
      }
226 |  |
227 | 24x |
      double precStep = (precisionMax - precisionMin) / (10000 - 1);
228 | 24x |
      colvec precisionDistribution = regspace(precisionMin, precStep, precisionMax);
229 |  |
230 | 24x |
      arma::vec precisionProbability;
231 | 24x |
      Rcpp::NumericMatrix KPs(numNeeded,2);
232 |  |
233 | 72250x |
      for(int i=0; i < numNeeded; i++) {
234 | 144404x |
        KPs(i, _) = find_valid_pk(
235 | 36101x |
          kappaDistribution,
236 | 36101x |
          kappaProbability,
237 | 36101x |
          precisionDistribution,
238 | 36101x |
          precisionProbability,
239 | 36101x |
          baserate
240 |  |
        );
241 |  |
      }
242 |  |
243 | 24x |
      return (KPs);
244 |  |
    }
245 |  |
246 |  |
    //' @name contingency_table
247 |  |
    //' @title contingency_table
248 |  |
    //' @description Create a contingency table using the provied precision, recall, baserate, and length.
249 |  |
    //'
250 |  |
    //' @param precision double
251 |  |
    //' @param rec double
252 |  |
    //' @param length int
253 |  |
    //' @param baserate double
254 |  |
    //'
255 |  |
    //' @export
256 |  |
    // [[Rcpp::export]]
257 |  |
    arma::imat contingency_table(double precision, double rec, int length, double baserate) {
258 | 58509x |
      int gold1s = max(NumericVector::create(round(baserate * length), 1));
259 | 19503x |
      int gold0s = length - gold1s;
260 | 58509x |
      int TP = max(NumericVector::create(round(gold1s * rec), 1));
261 | 97515x |
      int FP = min(NumericVector::create(gold0s, max(NumericVector::create(round(TP * (1 - precision) / precision),  1))));
262 |  |
263 | 19503x |
      arma::imat ct = arma::imat(2,2);
264 | 39006x |
      ct(0,0) = TP;
265 | 39006x |
      ct(0,1) = gold1s - TP;
266 | 39006x |
      ct(1,0) = FP;
267 | 39006x |
      ct(1,1) = gold0s - FP;
268 |  |
269 | 19503x |
      return(ct);
270 |  |
    }
271 |  |
272 |  |
    //' @name random_contingency_table
273 |  |
    //' @title random_contingency_table
274 |  |
    //' @param setLength [TBD]
275 |  |
    //' @param baserate [TBD]
276 |  |
    //' @param kappaMin [TBD]
277 |  |
    //' @param kappaMax [TBD]
278 |  |
    //' @param minPrecision [TBD]
279 |  |
    //' @param maxPrecision [TBD]
280 |  |
    //'
281 |  |
    //' @export
282 |  |
    // [[Rcpp::export]]
283 |  |
    arma::imat random_contingency_table(
284 |  |
        int setLength,double baserate,
285 |  |
        double kappaMin, double kappaMax,
286 |  |
        double minPrecision = 0, double maxPrecision = 1
287 |  |
    ) {
288 | 1x |
      NumericMatrix KP = generate_kp_list(1, baserate, kappaMin, kappaMax, minPrecision, maxPrecision);
289 | 2x |
      double kappa = KP(0,1);
290 | 2x |
      double precision = KP(0,0);
291 | 2x |
      double rec = recall(kappa, baserate, precision);
292 |  |
293 | 1x |
      arma::imat ct = contingency_table(precision, rec, setLength, baserate);
294 |  |
295 | 1x |
      return(ct);
296 |  |
    }
297 |  |
298 |  |
    //' @title kappa_ct
299 |  |
    //' @description Calculate kappa from a contingency table
300 |  |
    //' @param ct [TBD]
301 |  |
    //'
302 |  |
    //' @export
303 |  |
    // [[Rcpp::export]]
304 |  |
    double kappa_ct(arma::imat ct) {
305 | 39006x |
      double a = ct(0,0); //gold 1 & silver 1
306 | 19503x |
      double c = ct(0,1); //gold 1 & silver 0
307 | 19503x |
      double b = ct(1,0); //gold 0 & silver 1
308 | 19503x |
      double d = ct(1,1); //gold 0 & silver 0
309 | 19503x |
      double size = accu(ct);
310 |  |
311 | 19503x |
      double pZero = (a+d) / size;
312 | 19503x |
      double pYes = ((a + b) / size) * ((a + c) / size);
313 | 19503x |
      double pNo =  ((c + d) / size) * ((b + d) / size);
314 | 19503x |
      double pE = (pYes + pNo);
315 | 19503x |
      double k = (pZero - pE) / (1 - pE);
316 |  |
317 | 19503x |
      return(k);
318 |  |
    }
319 |  |
320 |  |
    //' @title Get Handset from Conntingency Table
321 |  |
    //' @description This function returns a Handset from a larger Contingency Table
322 |  |
    //' @name getHand_ct
323 |  |
    //'
324 |  |
    //' @param ct [TBD]
325 |  |
    //' @param handSetLength [TBD]
326 |  |
    //' @param handSetBaserate [TBD]
327 |  |
    //'
328 |  |
    //' @return Contingency Table matrix
329 |  |
    //'
330 |  |
    //' @export
331 |  |
    // [[Rcpp::export]]
332 |  |
    arma::imat getHand_ct(arma::imat ct, int handSetLength, double handSetBaserate) {
333 | 19500x |
      int positives = ceil(handSetLength * handSetBaserate);
334 | 19500x |
      arma::ivec positiveIndices;
335 | 19500x |
      arma::imat otherCT(ct);
336 |  |
337 | 19500x |
      if(positives > 0) {
338 | 34500x |
        arma::irowvec gold1s = ct.row(0);
339 | 46000x |
        positiveIndices = sample_contingency_table(gold1s, positives, false);
340 |  |
341 | 34500x |
        double sumOnes = sum(positiveIndices == 0);
342 | 34500x |
        double sumTwos = sum(positiveIndices == 1);
343 | 11500x |
        positiveIndices *= 2;
344 | 34500x |
        otherCT(0,0) = otherCT(0,0) - sumOnes;
345 | 34500x |
        otherCT(0,1) = otherCT(0,1) - sumTwos;
346 |  |
      }
347 |  |
348 | 58500x |
      arma::ivec otherAsVector = otherCT.as_col();
349 | 58500x |
      arma::ivec otherIndices = sample_contingency_table(otherCT, handSetLength - positives, false);
350 | 58500x |
      arma::ivec allIndices = join_cols<ivec>(positiveIndices, otherIndices);
351 | 19500x |
      arma::imat newCT = arma::imat(2,2);
352 |  |
353 | 78000x |
      newCT(0,0) = sum(allIndices == 0);
354 | 78000x |
      newCT(1,0) = sum(allIndices == 1);
355 | 78000x |
      newCT(0,1) = sum(allIndices == 2);
356 | 78000x |
      newCT(1,1) = sum(allIndices == 3);
357 |  |
358 | 19500x |
      return(newCT);
359 |  |
    }
360 |  |
361 |  |
    //' @name getHand_kappa
362 |  |
    //' @title getHand_kappa
363 |  |
    //'
364 |  |
    //' @description This function returns kappa calculated from a Handset taken from  a larger Contingency Table
365 |  |
    //'
366 |  |
    //' @param ct KPs matrix of kappa (column 1) and precision (column 2) values
367 |  |
    //' @param handSetLength The length of the \code{\link[=getTestSet]{testSet}} (ignored unless \emph{data} is an observed kappa value)
368 |  |
    //' @param handSetBaserate baserate to inflate the sampled contingency table to
369 |  |
    //'
370 |  |
    //' @return Kappa as double
371 |  |
    //'
372 |  |
    //' @export
373 |  |
    // [[Rcpp::export]]
374 |  |
    double getHand_kappa(arma::imat ct, int handSetLength, double handSetBaserate) {
375 | 39000x |
      arma::imat newCT = getHand_ct(ct, handSetLength, handSetBaserate);
376 |  |
377 | 58500x |
      double k = kappa_ct(newCT);
378 | 19500x |
      return(k);
379 |  |
    }
380 |  |
381 |  |
    // This function calculates rho given the parameters of the initial set and
382 |  |
    // the kappa of the observed handset and is called by rhoK()
383 |  |
    //
384 |  |
    // [[Rcpp::export]]
385 |  |
    double calcRho_c(
386 |  |
      double x,
387 |  |
      double OcSBaserate,
388 |  |
      int testSetLength,
389 |  |
      double testSetBaserateInflation = 0,
390 |  |
      int OcSLength = 10000,
391 |  |
      int replicates = 800,
392 |  |
      double ScSKappaThreshold = 0.65,
393 |  |
      double ScSKappaMin = 0.40,
394 |  |
      double ScSPrecisionMin = 0.6,
395 |  |
      double ScSPrecisionMax = 1.0,
396 |  |
      NumericMatrix KPs = NumericMatrix(0)
397 |  |
    ) {
398 | 100x |
      if (KPs.size() == 1 && KPs[0] == 0) {
399 | 40x |
        KPs = generate_kp_list(replicates, OcSBaserate, ScSKappaMin, ScSKappaThreshold, ScSPrecisionMin, ScSPrecisionMax);
400 |  |
      }
401 |  |
402 | 80x |
      if(KPs.nrow() < replicates) {
403 | 10x |
        replicates = KPs.nrow();
404 |  |
      }
405 |  |
406 | 80x |
      arma::vec savedKappas = arma::vec(replicates);
407 | 39160x |
      for (int KProw = 0; KProw < replicates; KProw++) {
408 | 58500x |
        NumericVector KP = KPs(KProw, _);
409 | 39000x |
        double precision = KP[0];
410 | 39000x |
        double kappa = KP[1];
411 | 39000x |
        double rec = recall(kappa, OcSBaserate, precision);
412 |  |
413 | 19500x |
        arma::imat fullCT = contingency_table(precision, rec, OcSLength, OcSBaserate);
414 | 58500x |
        double kk = getHand_kappa(fullCT, testSetLength, testSetBaserateInflation);
415 | 39000x |
        savedKappas[KProw] = kk;
416 |  |
      };
417 |  |
418 | 240x |
      return(getBootPvalue_c(savedKappas, x));
419 |  |
    }
420 |  |
421 |  |
    // // @name rhoCT_c
422 |  |
    // // @title rhoCT_ct
423 |  |
    // // @param contingencyTable [TBD]
424 |  |
    // // @param testSetBaserateInflation [TBD]
425 |  |
    // // @param replicates [TBD]
426 |  |
    // // @param OcSLength [TBD]
427 |  |
    // // @param OcSBaserate [TBD]
428 |  |
    // // @param ScSKappaThreshold [TBD]
429 |  |
    // // @param ScSKappaMin [TBD]
430 |  |
    // // @param ScSPrecisionMin [TBD]
431 |  |
    // // @param ScSPrecisionMax [TBD]
432 |  |
    // //' @export
433 |  |
    // // [[Rcpp::export]]
434 |  |
    // arma::vec rhoCT_c(
435 |  |
    //   arma::imat contingencyTable,
436 |  |
    //   double testSetBaserateInflation = 0.2,
437 |  |
    //   int replicates = 800,
438 |  |
    //   int OcSLength = 10000,
439 |  |
    //   double OcSBaserate = -1.0,
440 |  |
    //   double ScSKappaThreshold = 0.65,
441 |  |
    //   double ScSKappaMin = 0.40,
442 |  |
    //   double ScSPrecisionMin = 0.6,
443 |  |
    //   double ScSPrecisionMax = 1
444 |  |
    // ){
445 |  |
    //   // Observed kappa
446 |  |
    //   double kappa = kappa_ct(contingencyTable);
447 |  |
    //   double size = accu(contingencyTable);
448 |  |
    //   double gold1;
449 |  |
    //
450 |  |
    //   if(OcSBaserate == -1.0) {
451 |  |
    //     gold1 = accu(contingencyTable.row(0));
452 |  |
    //     OcSBaserate = gold1 / size;
453 |  |
    //   }
454 |  |
    //
455 |  |
    //   double rho = calcRho_c(
456 |  |
    //     kappa, OcSBaserate,
457 |  |
    //     size, testSetBaserateInflation,
458 |  |
    //     OcSLength, replicates,
459 |  |
    //     ScSKappaThreshold, ScSKappaMin, ScSPrecisionMin, ScSPrecisionMax
460 |  |
    //   );
461 |  |
    //
462 |  |
    //   return(NumericVector::create(kappa, rho));
463 |  |
    // }
464 |  |
465 |  |
466 |  |
    // //' @export
467 |  |
    // // [[Rcpp::export]]
468 |  |
    // bool run_rho_in_c(
469 |  |
    //     double threshold,  double infl,
470 |  |
    //     int handsetLength = 20,
471 |  |
    //     double handsetBaserate = 0.2,
472 |  |
    //     int replicates = 10000
473 |  |
    // ) {
474 |  |
    //   arma::imat ct = random_contingency_table(10000, infl, 0.4, threshold, 0.2, 0.8);
475 |  |
    //   // handset --> testset
476 |  |
    //   arma::imat handset = getHand_ct(ct, handsetLength, handsetBaserate);
477 |  |
    //   arma::vec res = rhoCT_c(handset, handsetBaserate, replicates, 10000, -1, 0.65, 0.4, 0.6, 1);
478 |  |
    //
479 |  |
    //   return(res[0]>threshold);
480 |  |
    // }
481 |  |