# R/HartiganShapes.R In Anthropometry: Statistical Methods for Anthropometric Data

#### Documented in HartiganShapes

```HartiganShapes <- function(array3D,numClust,algSteps=10,niter=10,stopCr=0.0001,simul,initLl,initials,
verbose){
#,computCost

time_iter <- list()
comp_time <- c()
list_ic1 <- list()
list_ic1_step <- list()
vect_all_rate <- c()

ll <- 1 : numClust
dist <- matrix(0, dim(array3D), numClust)

if(verbose){
print(Sys.time())
}
time_ini <- Sys.time()

#Initialize the objective function by a large enough value:
vopt <- 1e+08

#Ramdon restarts:
for(iter in 1 : niter){

wss_step <- list()

if(verbose){
cat("New iteration")
print(iter)

cat("Optimal value with which this iteration starts:")
print(vopt)
}
#STEP 1: For each point I, find its two closest centers, IC1(I) and IC2(I). Assign the point to IC1(I):
meanshapes <- 0 ; mean_sh <- list()
ic1 <- c() ; ic2 <- c() ; dt <- c()  ; nc <- c() #number of points in each cluster.
an1 <- c()   ; an2 <- c() ; itran <- c() ; ncp <- c()
indx <- c() ; d <- c() ; live <- c() ; wss <- c()

n <- dim(array3D)

initials_hart <- list()
if(initLl){
initials_hart[[iter]] <- initials[[iter]]
}else{
initials_hart[[iter]] <- sample(1:n, numClust, replace = FALSE)
}
if(verbose){
cat("Initial values of this iteration:")
print(initials_hart[[iter]])
}
meanshapes <- array3D[,,initials_hart[[iter]]]
#meanshapes_aux <- array3D[, , initials[[iter]]]

#if(computCost){
#time_ini_dist <- Sys.time()
#dist_aux = riemdist(array3D[,,1], y = meanshapes[,,1])
#time_end_dist <- Sys.time()
#cat("Computational cost of the Procrustes distance:")
#print(time_end_dist - time_ini_dist)
#}

for(i in 1 : n){
ic1[i] = 1
ic2[i] = 2

for(il in 1 : 2){
dt[il] = (riemdist(array3D[,,i], meanshapes[,,il]))^2
}

if(dt < dt){
ic1[i] = 2
ic2[i] = 1
temp = dt
dt = dt
dt = temp
}

if(simul == FALSE){
for(l in 3 : numClust){
db = (riemdist(array3D[,,i], meanshapes[,,l]))^2
if(db < dt){
if(dt <= db){
dt = db
ic2[i] = l
}else{
dt = dt
ic2[i] = ic1[i]
dt = db
ic1[i] = l
}
}
}
}
}

#if(computCost){
#time_ini_mean <- Sys.time()
#meanshapes_aux[,,1] = procGPA(array3D[, , ic1 == 1], distances = TRUE, pcaoutput = TRUE)\$mshape
#time_end_mean <- Sys.time()
#cat("Computational cost of the Procrustes mean:")
#print(time_end_mean - time_ini_mean)
#}

#STEP 2: Update the cluster centres to be the averages of points contained within them.
#Check to see if there is any empty cluster at this stage:
for(l in 1 : numClust){
nc[l] <- table(ic1)[l]
#print("Loop numClust")
#print(nc[l])
if(nc[l] <= 3){ # nc[l] == 0
#stop("At least one cluster is empty or has a very few elements after the initial assignment.
#     A better set of initial cluster centers is needed.")
#break
return(cat("At least one cluster is empty or has a very few elements after the initial assignment.
A better set of initial cluster centers is needed. No solution provided."))
}
}
#print("Loop checked successfully")

for(l in 1 : numClust){
aa = nc[l]
#print("This is array3D")
#print(array3D)
#print("This is ic1")
#print(ic1)
#print("This is l")
#print(l)
x <- array3D[, , ic1 == l]
#print("This is x")
#print(dim(x))
if (length(dim(x)) != 3) {
#stop("Please ensure that array3D has 3 dimensions.")
#break
return(cat("Please ensure that array3D has 3 dimensions.")) # This is not actually needed
# anymore because the previous return already stops the execution since there are some
# very small clusters.
}else{
meanshapes[,,l] = shapes::procGPA(x, distances = TRUE, pcaoutput = TRUE)\$mshape
}

#Initialize AN1, AN2, ITRAN and NCP.
#AN1(L) = NC(L) / (NC(L) - 1)
#AN2(L) = NC(L) / (NC(L) + 1)
#ITRAN(L) = 1 if cluster L is updated in the quick-transfer stage,
#         = 0 otherwise

#In the optimal-transfer stage, NCP(L) stores the step at which cluster L is last updated.
#In the quick-transfer stage, NCP(L) stores the step at which cluster L is last updated plus M:
an2[l] = aa / (aa + 1)
if(1 < aa){
an1[l] = aa / (aa - 1)
}else{
an1[l] = Inf
}
itran[l] = 1
ncp[l] = -1
}

indx <- 0
d[1:n] = 0
live[1:numClust] = 0

for(step in 1 : algSteps){
#In this stage, there is only one pass through the data. Each point is re-allocated, if necessary, to the
#cluster that will induce the maximum reduction in within-cluster sum of squares:
lis <- optraShapes(array3D,n,meanshapes,numClust,ic1,ic2,nc,an1,an2,ncp,d,itran,live,indx)

meanshapes <- lis[] ; ic1 <- lis[] ; ic2 <- lis[] ; nc <- lis[] ; an1 <- lis[] ; an2 <- lis[] ; ncp <- lis[]
d <- lis[] ; itran <- lis[] ; live <- lis[] ; indx <- lis[]

#Each point is tested in turn to see if it should be re-allocated to the cluster to which it is most likely
#to be transferred, IC2(I), from its present cluster, IC1(I). Loop through the data until no further change
#is to take place:
lis1 <- qtranShapes(array3D,n,meanshapes,ic1,ic2,nc,an1,an2,ncp,d,itran,indx)

meanshapes <- lis1[] ; ic1 <- lis1[] ; ic2 <- lis1[] ; nc <- lis1[] ; an1 <- lis1[] ; an2 <- lis1[] ; ncp <- lis1[]
d <- lis1[] ; itran <- lis1[]  ;  indx <- lis1[] ; icoun <- lis1[]

mean_sh[[step]] <- meanshapes

#NCP has to be set to 0 before entering OPTRA:
for( l in 1 : numClust ){
ncp[l] = 0
}

#Compute the within-cluster sum of squares for each cluster:
wss <- vector("list", numClust)
for(num_cl in 1 : numClust){
wss[[num_cl]] <- 0
array3D_cl <- array(0, dim = c(n, 3, table(ic1)[num_cl])) #table(ic1)[num_cl] is the number of observations that
#belong to each cluster.
array3D_cl <- array3D[,,ic1 == num_cl]

distances <- c()
for(num_mujs_cl in 1:table(ic1)[num_cl]){
distances[num_mujs_cl] <- riemdist(array3D_cl[,,num_mujs_cl], meanshapes[,,num_cl])^2
}

wss[[num_cl]] <-  sum(distances) / n
}

#Total within-cluster sum of squares:
wss_step[[step]] <- sum(unlist(wss))
list_ic1_step[[step]] <- ic1

if(verbose){
paste(cat("Clustering of the Nstep", step, ":\n"))
print(table(list_ic1_step[[step]]))
}
if(verbose){
if(iter <= 10){
paste(cat("Objective function of the Nstep", step))
print(wss_step[[step]])
}
}

if(step > 1){
aux <- wss_step[[step]]
aux1 <- wss_step[[step-1]]
if( ((aux1 - aux) / aux1) < stopCr ){
break
}
}
}#The algSteps loop ends here.

#Calculus of the objective function (the total within-cluster sum of squares):
wss1 <- vector("list", numClust)
for(num_cl in 1 : numClust){
wss1[[num_cl]] <- 0
array3D_cl1 <- array(0, dim = c(n, 3, table(ic1)[num_cl]))
array3D_cl1 <- array3D[,,ic1 == num_cl]

distances1 <- c()

for(num_mujs_cl in 1:table(ic1)[num_cl]){
distances1[num_mujs_cl] <- riemdist(array3D_cl1[,,num_mujs_cl], meanshapes[,,num_cl])^2
}

wss1[[num_cl]] <-  sum(distances1) / n
}
#Total within-cluster sum of squares:
wss_step1 <- 0
wss_step1 <- sum(unlist(wss1))

#Change the optimal value and the optimal centers (copt) if a reduction in the objective function happens:
if(wss_step1 > min(unlist(wss_step))){
if(min(unlist(wss_step)) < vopt){
vopt <- min(unlist(wss_step))
if(verbose){
#Improvements in the objective functions are printed:
cat("optimal")
print(vopt)
}
optim_wss <- which.min(unlist(wss_step))
copt <- mean_sh[[optim_wss]] #optimal centers.
ic1_opt <- list_ic1_step[[optim_wss]]
}
}else if(wss_step1 < vopt){
vopt <- wss_step1
if(verbose){
#Improvements in the objective functions are printed:
cat("optimal")
print(vopt)
}
optim_wss <- which.min(unlist(wss_step))
copt <- mean_sh[[optim_wss]] #optimal centers.
ic1_opt <- list_ic1_step[[optim_wss]]
}

time_iter[[iter]] <- Sys.time()

if(iter == 1){
comp_time <- difftime(time_iter[[iter]], time_ini, units = "mins")
if(verbose){
cat("Computational time of this iteration: \n")
print(time_iter[[iter]] - time_ini)
}
}else{
comp_time[iter] <- difftime(time_iter[[iter]], time_iter[[iter-1]], units = "mins")
if(verbose){
cat("Computational time of this iteration: \n")
print(time_iter[[iter]] - time_iter[[iter - 1]])
}
}
if(verbose){
cat("Optimal clustering of this iteration: \n")
}
optim_wss <- which.min(unlist(wss_step))
list_ic1[[iter]] <- list_ic1_step[[optim_wss]]
if(verbose){
print(table(list_ic1[[iter]]))
}

if(simul){
#Allocation rate:
as1 <- table(list_ic1[[iter]][1:(n/2)])
as2 <- table(list_ic1[[iter]][seq(n/2 + 1,n)])
if( max(as1) != n/2 & max(as2) != n/2 ){
suma <- min(as1) + min(as2)
all_rate <- 1 - suma / n
}else if( (max(as1) == n/2 & max(as2) != n/2) || (max(as1) != n/2 & max(as2) == n/2) ){
minim <- min(min(as1),min(as2))
all_rate <- 1 - minim / n
}else if( max(as1) == n/2 & max(as2) == n/2 ){
all_rate <- 1
}
vect_all_rate[iter] <- all_rate
if(verbose){
cat("Optimal allocation rate in this iteration:")
print(all_rate)
}
}
}#The niter loop ends here.

if(simul){
dimnames(copt) <- NULL
return(list(ic1=ic1_opt,cases=copt,vopt=vopt,compTime=comp_time,
AllRate=vect_all_rate))
}else{
return(list(ic1=ic1_opt,cases=copt,vopt=vopt))
}
}
```

## Try the Anthropometry package in your browser

Any scripts or data that you put into this service are public.

Anthropometry documentation built on Jan. 14, 2023, 1:16 a.m.