R/trimmedoid.R

In Anthropometry: Statistical Methods for Anthropometric Data

trimmedoid <- function(D, numClust, alpha, niter, algSteps = 7, verbose){
  n <- dim(D)[1]
  
  no.trim <- floor(n * (1 - alpha))
  
  ll <- 1:numClust
  
  ind <- 1:n
  
  dist <- ind
  
  # Initialize the objective function by a large enough value:
  vopt <- 1e+08
  
  # Random restarts:
  for (iter in 1:niter) {
    if (verbose) { 
      cat("new iteration")
      print(iter)
    }
    
    # Randomly choose the numClust initial centers:
    cini <- sample(1:n, size = numClust, replace = FALSE)
    
    # C-steps.
    for (t in 1:algSteps) {
      disti = c()
      ind = c()
      
      # Distances of each data point to its closest medoid:
      for (h in 1:n) {
        for (k in 1:numClust) {
          ll[k] <- D[h, cini[k]]
        }
        
        disti[h] <- min(ll)
        ind[h] <- which.min(ll)
      }
      
      # Modified data (Dmod) with the non-trimmed points and a vector indqq equal to the clusters allocations:
      qq <- (1:n)[disti <= sort(disti)[no.trim]]
      
      Dmod <- D[qq, qq]
      
      indqq = ind[qq]
      
      if (length(unique(indqq)) < numClust) {
        break
      }else{
        # Calculus of the new k centers:
        for (k in 1:numClust) {
          ni <- sum(indqq == k)
          
          if (ni > 1) {
            # Convert Dmod in dist:
            Dmodkv = as.dist(Dmod[indqq == k, indqq == k])
            
            rpam = pam(Dmodkv, k = 1) # here initial medoid, obtained by build.
            
            medk = rpam$id.med
          }
          
          if (ni == 1) {
            medk = 1
          }
          
          if (ni == 0) {
            medk = 1
          }
          
          aux = qq[indqq == k]
          
          cini[k] = aux[medk]
        }
        
        obj <- 0
    
        for (l in qq) {
          for (k in 1:numClust) {
            ll[k] <- D[l, cini[k]]
          }
          obj <- obj + min(ll)
        }
        
        if (iter < 10) {
          if (verbose) { 
            print(obj / no.trim)
          }  
        }
        
        rm(obj)
      }		
    }
    
    # Calculus of the trimmed k-variance:
    obj <- 0
    
    for (l in qq) { 
      for (k in 1:numClust) {
        ll[k] <- D[l, cini[k]]
      }
      obj <- obj + min(ll)
    }
    obj <- obj / no.trim
    
    # Change the optimal value and the optimal centers (copt) if a reduction in the objective function happens:
    if (obj < vopt) {
      vopt <- obj
      # Improvements in the objective functions are printed:
      if (verbose) {
        cat("optimal")
        print(vopt)
      }
      
      copt <- cini
    } 
  }
  
  # Obtain the final cluster allocations: 
  # This is necessary, because a final cluster assignment movement is possible.
  asig <- rep(0, n)
  
  # Distances of each data point to its closest medoid:
  for (h in 1:n) {
    for (k in 1:numClust) {
      ll[k] <- D[h, copt[k]] 
    }
    
    disti[h] <- min(ll)
    ind[h] <- which.min(ll)
  }
  
  # a vector indqq equal to the clusters allocations:
  qq <- (1:n)[disti <= sort(disti)[no.trim]]
  
  Dmod <- D[qq, qq]
  
  indqq = ind[qq]
  
  #Assign every observation to each cluster and 0 for the trimmed observations:
  asig[qq] = indqq
  
  # Between clusters sum of distances:
  b = 0
  for (k in 1:(numClust - 1)) {
    for (j in (k + 1):numClust) {
      b = b + D[copt[k], copt[j]] 
    }
  }
  
  # ch goodness index:
  ch = b * (no.trim - numClust) / (vopt * no.trim * (numClust - 1))
  
  rt = list(vopt = vopt, copt = copt, asig = asig, ch = ch, Dmod = Dmod, qq = qq) 
  
  return(rt)
}