Using a user-specified distance function

Note that we could also write our own distance function, and then use it with the discriminability driver. For a function to be a valid distance function for discriminability, it must accept an [n, d] array as the leading argument. Since our distance function returns the distance matrix as the only argument and requires no further parameters, we specify dist.return=NULL and dist.params=NULL. Note that the default for discriminability is the two-norm (euclidian distance) so we should expect the same result as previously:

# two norm between pairs of points
dist.fxn <- function(X) {
  n <- nrow(X)
  D <- array(0, dim=c(n, n))
  for (i in 1:(n - 1)) {
    for (j in i:n) {
      D[i,j] <- sum(abs(X[i,] - X[j,])^2)
    }
  }
  D <- D + t(D)
  return(D)
}

discr.stat(X, Y, dist.xfm=dist.fxn, dist.params=NULL, dist.return=NULL)$discr

## [1] 0.9005

Additionally, we can specify leading arguments as required to the distance function, and grab the appropriate return argument from the distance function that corresponds to the distance matrix. We can do this by specifying dist.params=list(method="2") and dist.return="Distance" in the below example:

# two norm between pairs of points
dist.fxn <- function(X, method="2") {
  if (method == "2") {
    n <- nrow(X)
    D <- array(0, dim=c(n, n))
    for (i in 1:(n - 1)) {
      for (j in i:n) {
        D[i,j] <- sum(abs(X[i,] - X[j,])^2)
      }
    }
    D <- D + t(D)
  } else {
    stop("Mistakes were made.")
  }
  return(list(Distance=D))
}

discr.stat(X, Y, dist.xfm=dist.fxn, dist.params=list(method="2"), dist.return="Distance")$discr

## [1] 0.9005

To investigate whether this discriminability exceeds that observed by random chance, we can use the one sample test of goodness-of-fit:

discr.test.one_sample(X, Y)$p.value

## [1] 0.001996008

The low p-value for a single test indicates that the discriminability exceeds that which would be observed via random chance.

Had we a second set of data, we could consider whether the discriminability of one dataset exceeds that of another. Consider a second dataset:

X2 <- t(sapply(Y, function(lab) rnorm(2, mean=lab, sd=2)))
X.dat.both <- rbind(X.dat, data.frame(x1=X2[,1], x2=X2[,2], Individual=factor(Y), Dataset="Second Dataset"))
ggplot(X.dat.both, aes(x=x1, y=x2, color=Individual)) +
  geom_point() +
  xlab("First Dimension") +
  ylab("Second Dimension") +
  ggtitle("Plot of Simulated Data") +
  theme_bw() +
  facet_grid(. ~ Dataset)

Note that the second dataset has a variance of 4 instead of 1, and therefore, appears less homogeneous than the first dataset. The “proximal balls” of measurements from each individual are now less pronounced. Therefore, we would expect the discriminability to be lower, as there is higher variance between measurements for a single subject:

discr.stat(X2, Y)$discr

## [1] 0.6408889

To formalize this notion, we can consider a test of whether the first dataset has a higher discriminability than the second dataset. This test can be performed using the two sample test of equality:

discr.test.two_sample(X, X2, Y, alt="greater")$p.value

## [1] 0.001996008

The low p-value indicates that we can conclude that the first dataset is more discriminable than the second dataset.