|
1 |
| -#' H+ Estimate |
| 1 | +#' @title H+ discordance metric |
2 | 2 | #'
|
3 |
| -#' Estimates H+ for two sets (vectors) A and B, or |
4 |
| -#' a dissimilarity matrix D and a label vector L. |
| 3 | +#' @description Estimates the H+ discordance metric |
| 4 | +#' for either (1) two sets (vectors) A and B, or |
| 5 | +#' (2) a dissimilarity matrix D and a label vector L. |
5 | 6 | #' Approximation is calculated using p+1 percentiles,
|
6 | 7 | #' with an accuracy bound of 1/p.
|
7 | 8 | #'
|
8 |
| -#' @param A (numeric vector) |
9 |
| -#' @param B (numeric vector) |
10 |
| -#' @param D (dist, nxn) |
11 |
| -#' @param L (numeric or character vector, nx1) |
12 |
| -#' @param p (integer) |
13 |
| -#' @param alg (integer), Algorithm 1 or 2 |
14 |
| -#' @return Estimated H+ (numeric) |
| 9 | +#' @param A numeric vector containing a set of length n |
| 10 | +#' @param B numeric vector containing a set of length n |
| 11 | +#' @param D distance matrix of dimension nxn |
| 12 | +#' @param L numeric or character vector of length n |
| 13 | +#' @param p integer representing the number of percentiles |
| 14 | +#' @param alg integer (1 or 2) representing the choice of algorithm used to estimate H+ (Algorithm 1 or 2) |
| 15 | +#' |
| 16 | +#' @return A numeric returning the estimated value of H+. |
15 | 17 | #' @export
|
| 18 | +#' |
16 | 19 | #' @examples
|
17 | 20 | #' a <- rnorm(n=500,mean=0)
|
18 |
| -#' b <- rnorm(n=500,mean=1 |
19 |
| -#' h <- hpe(a=a,b=b,p=101,alg=1) |
| 21 | +#' b <- rnorm(n=500,mean=1) |
| 22 | +#' h <- hpe(A=a,B=b,p=101,alg=1) |
20 | 23 | #'
|
21 | 24 | #' a <- sapply(1:500, function(i) rnorm(n=50,mean=0))
|
22 | 25 | #' b <- sapply(1:500, function(i) rnorm(n=50,mean=0))
|
23 | 26 | #' x <- cbind(a,b)
|
24 | 27 | #' d <- dist(t(x))
|
25 | 28 | #' l <- c(rep(0,500),rep(1,500))
|
26 | 29 | #' h <- hpe(D=d,L=l,p=101,alg=1)
|
27 |
| -hpe <- function(A,B,D,L,p=101,alg=1){ |
28 |
| - abflg<- missing(A) & missing(B) |
| 30 | +#' |
| 31 | +hpe <- function(A, B, D, L, p = 101, alg = 1) { |
| 32 | + abflg <- missing(A) & missing(B) |
29 | 33 | dlflg <- missing(D) & missing(L)
|
30 |
| - if( abflg & dlflg){ |
| 34 | + if (abflg & dlflg) { |
31 | 35 | stop("please provide either A B or D L")
|
32 |
| - } else if (!abflg){ |
| 36 | + } else if (!abflg) { |
33 | 37 | nmflg <- (!is.numeric(A) & !is.numeric(A))
|
34 |
| - if(nmflg){ |
| 38 | + if (nmflg) { |
35 | 39 | stop("please ensure A B are numeric")
|
36 | 40 | } else {
|
37 | 41 | print("Estimating H+ using A B formulation...")
|
38 | 42 | }
|
39 |
| - } else if(!dlflg){ |
40 |
| - # tyflg <- !(class(D)=='dist' | (is.matrix(D) & isSymmetric(D) & is.numeric(D)) ) |
41 |
| - tyflg <- !( class(D)=='dist' ) |
42 |
| - if(tyflg){ |
| 43 | + } else if (!dlflg) { |
| 44 | + # tyflg <- !(class(D)=='dist' | (is.matrix(D) & isSymmetric(D) & is.numeric(D)) ) |
| 45 | + tyflg <- !(class(D) == 'dist') |
| 46 | + if (tyflg) { |
43 | 47 | stop("please ensure D is a dist objcect")
|
44 | 48 | }
|
45 | 49 | tyflg <- !(is.numeric(L) | is.character(L))
|
46 |
| - if(tyflg){ |
| 50 | + if (tyflg) { |
47 | 51 | stop("please ensure L is numeric or character vector")
|
48 | 52 | }
|
49 | 53 | D <- as.matrix(D)
|
50 |
| - dmflg <- !( ncol(D) == nrow(D) & ncol(D) == length(L)) |
51 |
| - if(dmflg){ |
| 54 | + dmflg <- !(ncol(D) == nrow(D) & ncol(D) == length(L)) |
| 55 | + if (dmflg) { |
52 | 56 | stop("Dimension mismatch for D and L")
|
53 | 57 | } else {
|
54 | 58 | print("Estimating H+ using D L formulation...")
|
55 | 59 | }
|
56 |
| - ind <- sapply(L, function(x) sapply(L, function(y) x==y)) |
| 60 | + ind <- sapply(L, function(x) |
| 61 | + sapply(L, function(y) |
| 62 | + x == y)) |
57 | 63 | ind <- ind[upper.tri(ind)]
|
58 | 64 | iw <- which(ind)
|
59 | 65 | ib <- which(!ind)
|
60 | 66 | D <- D[upper.tri(D)]
|
61 | 67 | A <- D[iw]
|
62 | 68 | B <- D[ib]
|
63 | 69 | }
|
64 |
| - |
65 |
| - |
66 |
| - ps <- seq(0,1,length.out=p) |
| 70 | + |
| 71 | + |
| 72 | + ps <- seq(0, 1, length.out = p) |
67 | 73 | qA <- quantile(A, probs = ps)
|
68 | 74 | qB <- quantile(B, probs = ps)
|
69 |
| - |
| 75 | + |
70 | 76 | #call hp alg1 or 2
|
71 |
| - if(alg==1){ |
72 |
| - he <- alg1(qA,qB,p) |
73 |
| - }else if(alg==2){ |
74 |
| - he <- alg2(qA,qB,p) |
75 |
| - }else{ |
| 77 | + if (alg == 1) { |
| 78 | + he <- alg1(qA, qB, p) |
| 79 | + } else if (alg == 2) { |
| 80 | + he <- alg2(qA, qB, p) |
| 81 | + } else{ |
76 | 82 | stop("please specify a vaid algorithm (alg=1 or 2)")
|
77 | 83 | }
|
78 |
| - |
| 84 | + |
79 | 85 | #return estimate
|
80 | 86 | return(he)
|
81 | 87 | }
|
0 commit comments