Santiago Gallón Gómez

Home » Code

Code

Template curve estimation based on manifold embedding (Dimeglio, C., S. Gallón, J-M. Loubes, and E. Maza (2014). A Robust Algorithm for Template Curve Estimation Based on Manifold Embedding. Computational Statistics & Data Analysis 70, 373–386.)

This R code requires the igraph R-package.

</pre>
# Function to calculate the Minimum Spanning Tree of a data set X
mst <- function(X) {
   D <- as.matrix(dist(X))
   G <- graph.adjacency(D, "undirected", weighted = TRUE)
   T <- minimum.spanning.tree(G)
   A <- get.adjacency(T)
   return(A)
}

# Function to obtain an extended graph
manif <- function(X) {
   D <- as.matrix(dist(X))
   A <- mst(X)
   DA <- D*A
   R <- apply(DA, 1, max)
   n <- dim(X)[1]
   d <- dim(X)[2]
   A <- matrix(0, nrow=n, ncol=n)
   L <- 100
   f <- function(l,k) {
      (1-l)*X[i,k]+l*X[j,k]
   }
   for (i in 1:(n-1)) {
      for (j in (i+1):n) {
         P <- outer(1:L/(L+1), 1:d, f)
         E <- as.matrix(dist(rbind(P, X)))
         F <- E[(L+1):(L+n), 1:L]
         if (all(apply(F<R, 2, any))) A[i,j]=1
      }
   }
   A <- A + t(A)
   return(A)
}

# Function to calculate the estimated intrinsic mean
imean <- function(X, moment=1) {
   n <- dim(X)[1]
   D <- as.matrix(dist(X))
   A <- manif(X)
   DA <- D*A
   G <- graph.adjacency(DA, "undirected", weighted = TRUE)
   GD <- shortest.paths(G)
   GD <- GD^moment
   d <- apply(GD, 2, sum)
   im <- which.min(d)
   return(im)
}
<span>

Manifold normalization (Gallón, S., J-M. Loubes, and E. Maza (2013). Statistical Properties of the Quantile Normalization Method for Density Curve Alignment. Mathematical Biosciences 242(2), 129–142.)

This R code is based on a slight modication of the normalizeQuantiles function in the package limma intended to normalize single channel microarray intensities between arrays. The code also depends on the function imean in the code for template curve estimation based on manifold embedding.

</pre>
# Manifold normalization
normanif <- function(X, ties = TRUE) {
   n <- dim(X)
   if (is.null(n))
      return(X)
   if (n[2] == 1)
      return(X)
   O <- S <- array(, n)
   nobs <- rep(n[1], n[2])
   i <- (0:(n[1] - 1))/(n[1] - 1)
   for (j in 1:n[2]) {
      Si <- sort(X[, j], method = "quick", index.return = TRUE)
      nobsj <- length(Si$x)
      if (nobsj < n[1]) {
         nobs[j] <- nobsj
         isna <- is.na(X[, j])
         S[, j] <- approx((0:(nobsj - 1))/(nobsj - 1), Si$x, i, ties = "ordered")$y
         O[!isna, j] <- ((1:n[1])[!isna])[Si$ix]
      }
      else {
         S[, j] <- Si$x
         O[, j] <- Si$ix
      }
   }
   m <- S[, imean(t(S), 1)]
      for (j in 1:n[2]) {
         if (ties)
            r <- rank(X[, j])
         if (nobs[j] < n[1]) {
            isna <- is.na(X[, j])
            if (ties)
               X[!isna, j] <- approx(i, m, (r[!isna] - 1)/(nobs[j] - 1), ties="ordered")$y
            else
               X[O[!isna,j], j] <- approx(i, m, (0:(nobs[j]-1))/(nobs[j]-1), ties="ordered")$y
         }
         else {
            if (ties)
               X[, j] <- approx(i, m, (r - 1)/(n[1] - 1), ties = "ordered")$y
            else
               X[O[, j], j] <- m
      }
   }
   X
}
Advertisements