#-----------------------------------------------------------------------
#
# Statistical shape analysis routines
# written by Ian Dryden in R (see http://cran.r-project.org)
# (c) Ian Dryden
# University of Nottingham. version 1.2.6
# 2003-2021
#
# Includes contributions by many other authors, including
# Mohammad Faghihi, Kwang-Rae Kim, Alfred Kume,
# Gregorio Quintana-Orti, Amelia Simo.
#
###########################################################################
tangentcoords.partial.inv = function(v, p, R)
{
return(matrix(sqrt(1 - sum(v^2)) * c(p) + v, nrow = nrow(p)) %*% t(R))
}
preshape2shape = function(z)
{
k = nrow(z) + 1
H = defh(k - 1)
return(t(H) %*% z)
}
plot3darcs<-function(x,pcno=1,c=1,nn=100,boundary.data=TRUE,view.theta=0,view.phi=0,type="pnss"){
# points along principal arcs
pns.out <- x
k<-pns.out$GPAout$k
m<-pns.out$GPAout$m
n.pc<-dim(pns.out$resmat)[1]
npts = 100
arc1 = t(get.prinarc(resmat = pns.out$resmat, PNS = pns.out$PNS, arc = 1, n = npts, boundary.data = boundary.data))
arc2 = t(get.prinarc(resmat = pns.out$resmat, PNS = pns.out$PNS, arc = 2, n = npts, boundary.data = boundary.data))
arc3 = t(get.prinarc(resmat = pns.out$resmat, PNS = pns.out$PNS, arc = 3, n = npts, boundary.data = boundary.data))
# PNS mean
PNSmean = pns.out$PNS$mean
GPAout = pns.out$GPAout
# c * sd : lower and upper along principal arcs. lower -> 0 -> upper
# total 2 * nn + 1 configurations
# : nn configurations from negative to 0
# 1 configuration at 0
# nn configurations from 0 to positive
#default c = 1
{
cat("stdev of PNS1 score:", round(sd(pns.out$resmat[ 1, ]), 4), "\n")
cat("stdev of PNS2 score:", round(sd(pns.out$resmat[ 2, ]), 4), "\n")
cat("stdev of PNS3 score:", round(sd(pns.out$resmat[ 3, ]), 4), "\n")
}
rng = c * sd(pns.out$resmat[1,])
val = c(seq(-rng, 0, length = nn + 1)[-(nn + 1)],
0,
seq(0, rng, length = nn + 1)[-1])
lu.arc1 = t(get.prinarc.value(PNS = pns.out$PNS, arc = 1, res = val))
rng = c * sd(pns.out$resmat[2,])
val = c(seq(-rng, 0, length = nn + 1)[-(nn + 1)],
0,
seq(0, rng, length = nn + 1)[-1])
lu.arc2 = t(get.prinarc.value(PNS = pns.out$PNS, arc = 2, res = val))
rng = c * sd(pns.out$resmat[3,])
val = c(seq(-rng, 0, length = nn + 1)[-(nn + 1)],
0,
seq(0, rng, length = nn + 1)[-1])
lu.arc3 = t(get.prinarc.value(PNS = pns.out$PNS, arc = 3, res = val))
# spherical data to PC scores
scores.arc1 = sphere2pcscore(x = arc1)
scores.arc2 = sphere2pcscore(x = arc2)
scores.arc3 = sphere2pcscore(x = arc3)
scores.PNSmean = sphere2pcscore(x = t(PNSmean))
scores.lu.arc1 = sphere2pcscore(x = lu.arc1)
scores.lu.arc2 = sphere2pcscore(x = lu.arc2)
scores.lu.arc3 = sphere2pcscore(x = lu.arc3)
# PC scores to tangent coordinates
U1 = matrix(0, npts, ncol(GPAout$pcar))
U2 = matrix(0, npts, ncol(GPAout$pcar))
U3 = matrix(0, npts, ncol(GPAout$pcar))
for (i in 1:npts)
{
for (j in 1:n.pc)
{
U1[i, ] = U1[i, ] + scores.arc1[i, j] * GPAout$pcar[ , j]
U2[i, ] = U2[i, ] + scores.arc2[i, j] * GPAout$pcar[ , j]
U3[i, ] = U3[i, ] + scores.arc3[i, j] * GPAout$pcar[ , j]
}
}
U.mean = matrix(0, 1, ncol(GPAout$pcar))
for (j in 1:n.pc)
{
U.mean = U.mean + scores.PNSmean[j] * GPAout$pcar[ , j]
}
tan.lu.arc1 = matrix(0, nrow(lu.arc1), ncol(GPAout$pcar))
tan.lu.arc2 = matrix(0, nrow(lu.arc2), ncol(GPAout$pcar))
tan.lu.arc3 = matrix(0, nrow(lu.arc3), ncol(GPAout$pcar))
for (i in 1:nrow(lu.arc1))
{
for (j in 1:n.pc)
{
tan.lu.arc1[i, ] = tan.lu.arc1[i, ] + scores.lu.arc1[i, j] * GPAout$pcar[ , j]
tan.lu.arc2[i, ] = tan.lu.arc2[i, ] + scores.lu.arc2[i, j] * GPAout$pcar[ , j]
tan.lu.arc3[i, ] = tan.lu.arc3[i, ] + scores.lu.arc3[i, j] * GPAout$pcar[ , j]
}
}
# tangent coordinates to preshapes
# preshapes to shapes
shapes.arc1 = array(NA, c(k, m, npts))
shapes.arc2 = array(NA, c(k, m, npts))
shapes.arc3 = array(NA, c(k, m, npts))
H = defh(k - 1)
for (i in 1:npts)
{
shapes.arc1[ , , i] = preshape2shape(tangentcoords.partial.inv(v = U1[i, ], p = H %*% GPAout$mshape, R = diag(m)))
shapes.arc2[ , , i] = preshape2shape(tangentcoords.partial.inv(v = U2[i, ], p = H %*% GPAout$mshape, R = diag(m)))
shapes.arc3[ , , i] = preshape2shape(tangentcoords.partial.inv(v = U3[i, ], p = H %*% GPAout$mshape, R = diag(m)))
}
shapes.PNSmean = preshape2shape(tangentcoords.partial.inv(v = U.mean, p = H %*% GPAout$mshape, R = diag(m)))
shapes.lu.arc1 = array(NA, c(k, m, nrow(lu.arc1)))
shapes.lu.arc2 = array(NA, c(k, m, nrow(lu.arc2)))
shapes.lu.arc3 = array(NA, c(k, m, nrow(lu.arc3)))
for (i in 1:nrow(lu.arc1))
{
shapes.lu.arc1[ , , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc1[i, ], p = H %*% GPAout$mshape, R = diag(m)))
shapes.lu.arc2[ , , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc2[i, ], p = H %*% GPAout$mshape, R = diag(m)))
shapes.lu.arc3[ , , i] = preshape2shape(tangentcoords.partial.inv(v = tan.lu.arc3[i, ], p = H %*% GPAout$mshape, R = diag(m)))
}
#convert PCs to configuration space if not already in place
h <- defh(k - 1)
zero <- matrix(0, k - 1, k)
H <- cbind(h, zero, zero)
H1 <- cbind(zero, h, zero)
H2 <- cbind(zero, zero, h)
H <- rbind(H, H1, H2)
if (dim(GPAout$pcar)[1] == (3 * (k - 1))) {
pcarot <- (t(H) %*% GPAout$pcar)
GPAout$pcar <- pcarot
}
#------------------------#
# #
# PLOT TYPE I #
# : full plot #
# #
#------------------------#
if (pcno==1){
shapes.lu.arc<-shapes.lu.arc1
}
if (pcno==2){
shapes.lu.arc<-shapes.lu.arc2
}
if (pcno==3){
shapes.lu.arc<-shapes.lu.arc3
}
if (type=="pca"){
open3d()
par3d(windowRect = c(20, 30, 800, 800))
view3d(view.theta,view.phi)
# GPA mean + PC
plot3d(GPAout$mshape, type = "s", col = rainbow(k), size = 1, add = TRUE)
lines3d(GPAout$mshape, col = rainbow(k), lwd = 5)
pcu <- GPAout$mshape + c * GPAout$pcasd[pcno] *
cbind(GPAout$pcar[1:k, pcno], GPAout$pcar[(k + 1):(2 *
k), pcno], GPAout$pcar[(2 * k + 1):(3 * k), pcno])
pcl <- GPAout$mshape - c * GPAout$pcasd[pcno] *
cbind(GPAout$pcar[1:k, pcno], GPAout$pcar[(k + 1):(2 *
k), pcno], GPAout$pcar[(2 * k + 1):(3 * k), pcno])
spheres3d( pcu, radius = 0.004, color = "black")
spheres3d( pcl, radius = 0.004, color = "grey")
for (j in 1:k) {
lines3d(rbind(pcl[j, ], pcu[j, ]), col = rainbow(k)[j] )
if (j>1){
lines3d(rbind(pcu[j-1, ],pcu[j, ]) ,col="black" )
lines3d(rbind(pcl[j-1, ],pcl[j, ]) ,col="grey" )
}
}
}
if (type=="pnss"){
open3d()
par3d(windowRect = c(20, 30, 800, 800))
view3d(view.theta,view.phi)
# principal arcs in landmark space
# PNS mean in landmark space
plot3d(shapes.PNSmean, type = "s", col = rainbow(k), size = 1.3, add = TRUE)
lines3d(shapes.PNSmean,lwd=5,col=rainbow(k))
# lower and upper
for (i in 1:k)
{
lines3d(t(shapes.lu.arc[i, , ]), col = rainbow(k)[i], lwd = 1,lty=2)
# farthest negative from 0
spheres3d(head(t(shapes.lu.arc[i, , ]), 1), radius = 0.004, color = "black")
if (i>1){
lines3d( (shapes.lu.arc[(i-1):i,,1]) ,col="black" )
}
# # farthest positive from 0
spheres3d(tail(t(shapes.lu.arc[i, , ]), 1), radius = 0.004, color = "grey")
if (i>1){
lines3d( (shapes.lu.arc[(i-1):i,,201]) ,col="grey" )
}
}
}
out<-list(PNSmean=0,lu.arc=0)
out$PNSmean<-shapes.PNSmean
out$lu.arc<-shapes.lu.arc
out
}
########
pnss3d<-
function (x, sphere.type = "seq.test", alpha = 0.1, R = 100,
nlast.small.sphere = 0, n.pc = 3)
{
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]
if (n.pc =="Full" ) {
n.pc=m*k-m*(m-1)/2-m
}
if (m==2){
tem1 <- array( 0, c(k,3,n) )
tem1[,1:2,]<-x
x<-tem1
m<-3
}
if (n < ((k - 1) * m)) {
print("Note: n < (k - 1) * m.")
jj<- round( (k-1)*m/n + 0.5)
print("Adding extra copies of the data")
tem<- array(0,c(k,m,jj*n))
tem[,,1:n]<-x
for (i in 2:jj){
for (j in 1:n){
tem[,,(i-1)*n+ j ]<-x[,,j] + 0*matrix( rnorm(k*m), k,m)
}
}
x<-tem
}
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]
out = pc2sphere2(x = x, n.pc = n.pc)
spheredata = t(out$spheredata)
GPAout = out$GPAout
pns.out = pns(x = spheredata, sphere.type = sphere.type,
alpha = alpha, R = R, nlast.small.sphere = nlast.small.sphere)
print("Radii of spheres")
print(pns.out$PNS$radii)
pns.out$percent = pns.out$percent * sum(GPAout$percent[1:n.pc])/100
pns.out$GPAout = GPAout
pns.out$spheredata = spheredata
return(pns.out)
}
pc2sphere2<-function (x, n.pc)
{
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]
GPAout = procGPA(x = x, scale = TRUE, reflect = FALSE, tol1=1e-5,tangentcoords = "partial",
distances = FALSE)
cat("First ", n.pc, " principal components explain ", round(sum(GPAout$percent[1:n.pc])),
"% of total variance. \n", sep = "")
H = defh(k - 1)
X.hat = H %*% GPAout$mshape
S = array(NA, c(k - 1, m, n))
for (i in 1:n) {
S[, , i] = H %*% GPAout$rotated[, , i]
}
T.c = GPAout$tan - apply(GPAout$tan, 1, mean)
out = pcscore2sphere(n.pc = n.pc, X.hat = X.hat, S = S, Tan = T.c,
V = GPAout$pcar)
return(list(spheredata = out, GPAout = GPAout))
}
#==================================================================================
# PNS The Principal Nested Spheres code (PNS) for spheres and shapes has
# been written by Kwang-Rae Kim, and builds closely on the original matlab
# code for PNS by Sungkyu Jung
#==================================================================================
#==================================================================================
pns = function(x,
sphere.type = "seq.test",
alpha = 0.1,
R = 100,
nlast.small.sphere = 0)
{
n = ncol(x)
k = nrow(x)
if (abs(sum(apply(x ^ 2, 2, sum)) - n) > 1e-8)
{
stop("Error: Each column of x should be a unit vector, ||x[ , i]|| = 1.")
}
svd.x = svd(x, nu = nrow(x))
uu = svd.x$u
maxd = which(svd.x$d < 1e-15)[1]
if (is.na(maxd) | k > n)
{
maxd = min(k, n) + 1
}
nullspdim = k - maxd + 1
d = k - 1
cat("Message from pns() : dataset is on ", d, "-sphere. \n", sep = "")
if (nullspdim > 0)
{
cat(" .. found null space of dimension ",
nullspdim,
", to be trivially reduced. \n",
sep = "")
}
resmat = matrix(NA, d, n)
orthaxis = list()
orthaxis[[d - 1]] = NA
dist = rep(NA, d - 1)
pvalues = matrix(NA, d - 1, 2)
ratio = rep(NA, d - 1)
currentSphere = x
if (nullspdim > 0)
{
for (i in 1:nullspdim)
{
oaxis = uu[, ncol(uu) - i + 1]
r = pi / 2
pvalues[i,] = c(NaN, NaN)
res = acos(t(oaxis) %*% currentSphere) - r
orthaxis[[i]] = oaxis
dist[i] = r
resmat[i,] = res
NestedSphere = rotMat(oaxis) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)
uu = rotMat(oaxis) %*% uu
uu = uu[1:(k - i),] / repmat(matrix(sqrt(1 - uu[nrow(uu),] ^ 2), nrow = 1), k - i, 1)
cat(d - i + 1,
"-sphere to ",
d - i,
"-sphere, by ",
"NULL space \n",
sep = "")
}
}
if (sphere.type == "seq.test")
{
cat(" .. sequential tests with significance level ",
alpha,
"\n",
sep = "")
isIsotropic = FALSE
for (i in (nullspdim + 1):(d - 1))
{
if (!isIsotropic)
{
sp = getSubSphere(x = currentSphere, geodesic = "small")
center.s = sp$center
r.s = sp$r
resSMALL = acos(t(center.s) %*% currentSphere) - r.s
sp = getSubSphere(x = currentSphere, geodesic = "great")
center.g = sp$center
r.g = sp$r
resGREAT = acos(t(center.g) %*% currentSphere) - r.g
pval1 = LRTpval(resGREAT, resSMALL, n)
pvalues[i, 1] = pval1
if (pval1 > alpha)
{
center = center.g
r = r.g
pvalues[i, 2] = NA
cat(
d - i + 1,
"-sphere to ",
d - i,
"-sphere, by GREAT sphere, p(LRT) = ",
pval1,
"\n",
sep = ""
)
} else {
pval2 = vMFtest(currentSphere, R)
pvalues[i, 2] = pval2
if (pval2 > alpha)
{
center = center.g
r = r.g
cat(
d - i + 1,
"-sphere to ",
d - i,
"-sphere, by GREAT sphere, p(LRT) = ",
pval1,
", p(vMF) = ",
pval2,
"\n",
sep = ""
)
isIsotropic = TRUE
} else {
center = center.s
r = r.s
cat(
d - i + 1,
"-sphere to ",
d - i,
"-sphere, by SMALL sphere, p(LRT) = ",
pval1,
", p(vMF) = ",
pval2,
"\n",
sep = ""
)
}
}
} else if (isIsotropic) {
sp = getSubSphere(x = currentSphere, geodesic = "great")
center = sp$center
r = sp$r
cat(
d - i + 1,
"-sphere to ",
d - i,
"-sphere, by GREAT sphere, restricted by testing vMF distn",
"\n",
sep = ""
)
pvalues[i, 1] = NA
pvalues[i, 2] = NA
}
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res
cur.proj = project.subsphere(x = currentSphere,
center = center,
r = r)
NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)
if (nrow(currentSphere) == 3)
{
plot3d(
x = t(currentSphere),
xlab = "",
ylab = "",
zlab = "",
xlim = c(-1, 1),
ylim = c(-1, 1),
zlim = c(-1, 1),
box = FALSE,
axes = FALSE,
aspect = "iso"
)
axis3d(
edge = 'x',
labels = TRUE,
tick = TRUE,
at = c(-1, 1),
pos = c(NA, 0, 0)
)
axis3d(
edge = 'y',
labels = TRUE,
tick = TRUE,
at = c(-1, 1),
pos = c(0, NA, 0)
)
axis3d(
edge = 'z',
labels = TRUE,
tick = TRUE,
at = c(-1, 1),
pos = c(0, 0, NA)
)
}
if (nrow(currentSphere) == 2)
{
points3d(t(cur.proj), col = "red", size = 2)
}
}
} else if (sphere.type == "BIC") {
cat(" .. with BIC \n")
for (i in (nullspdim + 1):(d - 1))
{
sp = getSubSphere(x = currentSphere, geodesic = "small")
center.s = sp$center
r.s = sp$r
resSMALL = acos(t(center.s) %*% currentSphere) - r.s
sp = getSubSphere(x = currentSphere, geodesic = "great")
center.g = sp$center
r.g = sp$r
resGREAT = acos(t(center.g) %*% currentSphere) - r.g
BICsmall = n * log(mean(resSMALL ^ 2)) + (d - i + 1 + 1) * log(n)
BICgreat = n * log(mean(resGREAT ^ 2)) + (d - i + 1) * log(n)
cat("BICsm: ", BICsmall, ", BICgr: ", BICgreat, "\n", sep = "")
if (BICsmall > BICgreat)
{
center = center.g
r = r.g
cat(d - i + 1,
"-sphere to ",
d - i,
"-sphere, by ",
"GREAT sphere, BIC \n",
sep = "")
} else {
center = center.s
r = r.s
cat(d - i + 1,
"-sphere to ",
d - i,
"-sphere, by ",
"SMALL sphere, BIC \n",
sep = "")
}
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res
NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)
}
} else if (sphere.type == "small" | sphere.type == "great") {
pvalues = NaN
for (i in (nullspdim + 1):(d - 1))
{
sp = getSubSphere(x = currentSphere, geodesic = sphere.type)
center = sp$center
r = sp$r
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res
NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)
}
} else if (sphere.type == "bi.sphere") {
if (nlast.small.sphere < 0)
{
cat("!!! Error from pns(): \n")
cat("!!! nlast.small.sphere should be >= 0. \n")
return(NULL)
}
mx = (d - 1) - nullspdim
if (nlast.small.sphere > mx)
{
cat("!!! Error from pns(): \n")
cat("!!! nlast.small.sphere should be <= ",
mx,
" for this data. \n",
sep = "")
return(NULL)
}
pvalues = NaN
if (nlast.small.sphere != mx)
{
for (i in (nullspdim + 1):(d - 1 - nlast.small.sphere))
{
sp = getSubSphere(x = currentSphere, geodesic = "great")
center = sp$center
r = sp$r
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res
NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)
}
}
if (nlast.small.sphere != 0)
{
for (i in (d - nlast.small.sphere):(d - 1))
{
sp = getSubSphere(x = currentSphere, geodesic = "small")
center = sp$center
r = sp$r
res = acos(t(center) %*% currentSphere) - r
orthaxis[[i]] = center
dist[i] = r
resmat[i,] = res
NestedSphere = rotMat(center) %*% currentSphere
currentSphere = NestedSphere[1:(k - i),] /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), k - i, 1)
}
}
} else {
print("!!! Error from pns():")
print("!!! sphere.type must be 'seq.test', 'small', 'great', 'BIC', or 'bi.sphere'")
print("!!! Terminating execution ")
return(NULL)
}
S1toRadian = atan2(currentSphere[2,], currentSphere[1,])
meantheta = geodmeanS1(S1toRadian)$geodmean
orthaxis[[d]] = meantheta
resmat[d,] = mod(S1toRadian - meantheta + pi, 2 * pi) - pi
par(
mfrow = c(1, 1),
mar = c(4, 4, 1, 1),
mgp = c(2.5, 1, 0),
cex = 0.8
)
plot(
currentSphere[1,],
currentSphere[2,],
xlab = "",
ylab = "",
xlim = c(-1, 1),
ylim = c(-1, 1),
asp = 1
)
abline(h = 0, v = 0)
points(
cos(meantheta),
sin(meantheta),
pch = 1,
cex = 3,
col = "black",
lwd = 5
)
abline(
a = 0,
b = sin(meantheta) / cos(meantheta),
lty = 3
)
l = mod(S1toRadian - meantheta + pi, 2 * pi) - pi
points(
cos(S1toRadian[which.max(l)]),
sin(S1toRadian[which.max(l)]),
pch = 4,
cex = 3,
col = "blue"
)
points(
cos(S1toRadian[which.min(l)]),
sin(S1toRadian[which.min(l)]),
pch = 4,
cex = 3,
col = "red"
)
legend(
"topright",
legend = c("Geodesic mean", "Max (+)ve from mean", "Min (-)ve from mean"),
col = c("black", "blue", "red"),
pch = c(1, 4, 4)
)
{
cat("\n")
cat(
"length of BLUE from geodesic mean : ",
max(l),
" (",
round(max(l) * 180 / pi),
" degree)",
"\n",
sep = ""
)
cat(
"length of RED from geodesic mean : ",
min(l),
" (",
round(min(l) * 180 / pi),
" degree)",
"\n",
sep = ""
)
cat("\n")
}
radii = 1
for (i in 1:(d - 1))
{
radii = c(radii, prod(sin(dist[1:i])))
}
resmat = flipud0(repmat(matrix(radii, ncol = 1), 1, n) * resmat)
PNS = list()
PNS$radii = radii
PNS$orthaxis = orthaxis
PNS$dist = dist
PNS$pvalues = pvalues
PNS$ratio = ratio
PNS$basisu = NULL
PNS$mean = c(PNSe2s(matrix(0, d, 1), PNS))
if (sphere.type == "seq.test")
{
PNS$sphere.type = "seq.test"
} else if (sphere.type == "small") {
PNS$sphere.type = "small"
} else if (sphere.type == "great") {
PNS$sphere.type = "great"
} else if (sphere.type == "BIC") {
PNS$sphere.type = "BIC"
}
varPNS = apply(abs(resmat) ^ 2, 1, sum) / n
total = sum(varPNS)
propPNS = varPNS / total * 100
return(list(
resmat = resmat,
PNS = PNS,
percent = propPNS
))
}
#==================================================================================
pns4pc = function(x,
sphere.type = "seq.test",
alpha = 0.1,
R = 100,
nlast.small.sphere = 0,
n.pc = 2)
{
if (n.pc < 2)
{
stop("Error: n.pc should be >= 2.")
}
out = pc2sphere(x = x, n.pc = n.pc)
spheredata = t(out$spheredata)
GPAout = out$GPAout
pns.out = pns(
x = spheredata,
sphere.type = sphere.type,
alpha = alpha,
R = R,
nlast.small.sphere = nlast.small.sphere
)
pns.out$percent = pns.out$percent * sum(GPAout$percent[1:n.pc]) / 100
pns.out$GPAout = GPAout
pns.out$spheredata = spheredata
return(pns.out)
}
pns.pc = function(x,
sphere.type = "seq.test",
alpha = 0.1,
R = 100,
nlast.small.sphere = 0,
n.pc = 0)
{
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]
if (n.pc == 0)
{
GPAout = procGPA(
x = x,
scale = TRUE,
reflect = FALSE,
tangentcoords = "partial",
distances = FALSE
)
spheredata = matrix(NA, k * m, n)
for (i in 1:n)
{
spheredata[, i] = c(GPAout$rotated[, , i])
}
pns.out = pns(
x = spheredata,
sphere.type = sphere.type,
alpha = alpha,
R = R,
nlast.small.sphere = nlast.small.sphere
)
resmat = pns.out$resmat
PNS = pns.out$PNS
npts = 200
prinarc1 = get.prinarc(
resmat,
PNS,
arc = 1,
n = npts,
boundary.data = FALSE
)
prinarc2 = get.prinarc(
resmat,
PNS,
arc = 2,
n = npts,
boundary.data = FALSE
)
prinarc1.ar = array(NA, c(k, m, npts))
prinarc2.ar = array(NA, c(k, m, npts))
for (i in 1:npts)
{
prinarc1.ar[, , i] = matrix(prinarc1[, i], nrow = k)
prinarc2.ar[, , i] = matrix(prinarc2[, i], nrow = k)
}
scores.prinarc1 = shape.pcscores.partial(PCAout = GPAout, x = prinarc1.ar)
scores.prinarc2 = shape.pcscores.partial(PCAout = GPAout, x = prinarc2.ar)
out = pns.out
out$GPAout = GPAout
out$scores.prinarc1 = scores.prinarc1
out$scores.prinarc2 = scores.prinarc2
} else {
pns.out = pns4pc(
x = x,
sphere.type = sphere.type,
alpha = alpha,
R = R,
nlast.small.sphere = nlast.small.sphere,
n.pc = n.pc
)
GPAout = pns.out$GPAout
resmat = pns.out$resmat
PNS = pns.out$PNS
npts = 200
prinarc1 = get.prinarc(
resmat,
PNS,
arc = 1,
n = npts,
boundary.data = FALSE
)
prinarc2 = get.prinarc(
resmat,
PNS,
arc = 2,
n = npts,
boundary.data = FALSE
)
scores.prinarc1 = matrix(NA, npts, n.pc)
scores.prinarc2 = matrix(NA, npts, n.pc)
for (g in 1:npts)
{
size1 = acos(prinarc1[1, g])
size2 = acos(prinarc2[1, g])
scores.prinarc1[g,] = prinarc1[2:(n.pc + 1), g] / (sin(size1) / size1)
scores.prinarc2[g,] = prinarc2[2:(n.pc + 1), g] / (sin(size2) / size2)
}
out = pns.out
out$scores.prinarc1 = scores.prinarc1
out$scores.prinarc2 = scores.prinarc2
}
return(out)
}
#==================================================================================
rotMat = function(b, a = NULL, alpha = NULL)
{
if (is.matrix(b))
{
if (min(dim(b)) == 1)
{
b = c(b)
} else {
stop("Error: b should be a unit vector.")
}
}
d = length(b)
b = b / norm(b, type = "2")
if (is.null(a) & is.null(alpha))
{
a = c(rep(0, d - 1), 1)
alpha = acos(sum(a * b))
} else if (!is.null(a) & is.null(alpha)) {
alpha = acos(sum(a * b))
} else if (is.null(a) & !is.null(alpha)) {
a = c(rep(0, d - 1), 1)
}
if (abs(sum(a * b) - 1) < 1e-15)
{
rot = diag(d)
return(rot)
}
if (abs(sum(a * b) + 1) < 1e-15)
{
rot = -diag(d)
return(rot)
}
c = b - a * sum(a * b)
c = c / norm(c, type = "2")
A = a %*% t(c) - c %*% t(a)
rot = diag(d) + sin(alpha) * A + (cos(alpha) - 1) * (a %*% t(a) + c %*% t(c))
return(rot)
}
#==================================================================================
ExpNPd = function(x)
{
if (is.vector(x))
{
x = as.matrix(x)
}
d = nrow(x)
nv = sqrt(apply(x ^ 2, 2, sum))
Exppx = rbind(matrix(rep(sin(nv) / nv, d), nrow = d, byrow = T) * x, cos(nv))
Exppx[, nv < 1e-16] = repmat(matrix(c(rep(0, d), 1)), 1, sum(nv < 1e-16))
return(Exppx)
}
#==================================================================================
LogNPd = function(x)
{
n = ncol(x)
d = nrow(x)
scale = acos(x[d,]) / sqrt(1 - x[d,] ^ 2)
scale[is.nan(scale)] = 1
Logpx = repmat(t(scale), d - 1, 1) * x[-d,]
return(Logpx)
}
#==================================================================================
objfn = function(center, r, x)
{
return(mean((acos(t(
center
) %*% x) - r) ^ 2))
}
#==================================================================================
getSubSphere = function(x, geodesic = "small")
{
svd.x = svd(x)
initialCenter = svd.x$u[, ncol(svd.x$u)]
c0 = initialCenter
TOL = 1e-10
cnt = 0
err = 1
n = ncol(x)
d = nrow(x)
Gnow = 1e+10
while (err > TOL)
{
c0 = c0 / norm(c0, type = "2")
rot = rotMat(c0)
TpX = LogNPd(rot %*% x)
fit = sphereFit(
x = TpX,
initialCenter = rep(0, d - 1),
geodesic = geodesic
)
newCenterTp = fit$center
r = fit$r
if (r > pi)
{
r = pi / 2
svd.TpX = svd(TpX)
newCenterTp = svd.TpX$u[, ncol(svd.TpX$u)] * pi / 2
}
newCenter = ExpNPd(newCenterTp)
center = solve(rot, newCenter)
Gnext = objfn(center, r, x)
err = abs(Gnow - Gnext)
Gnow = Gnext
c0 = center
cnt = cnt + 1
if (cnt > 30)
{
break
}
}
i1save = list()
i1save$Gnow = Gnow
i1save$center = center
i1save$r = r
U = princomp(t(x))$loadings[,]
initialCenter = U[, ncol(U)]
c0 = initialCenter
TOL = 1e-10
cnt = 0
err = 1
n = ncol(x)
d = nrow(x)
Gnow = 1e+10
while (err > TOL)
{
c0 = c0 / norm(c0, type = "2")
rot = rotMat(c0)
TpX = LogNPd(rot %*% x)
fit = sphereFit(
x = TpX,
initialCenter = rep(0, d - 1),
geodesic = geodesic
)
newCenterTp = fit$center
r = fit$r
if (r > pi)
{
r = pi / 2
svd.TpX = svd(TpX)
newCenterTp = svd.TpX$u[, ncol(svd.TpX$u)] * pi / 2
}
newCenter = ExpNPd(newCenterTp)
center = solve(rot, newCenter)
Gnext = objfn(center, r, x)
err = abs(Gnow - Gnext)
Gnow = Gnext
c0 = center
cnt = cnt + 1
if (cnt > 30)
{
break
}
}
if (i1save$Gnow == min(Gnow, i1save$Gnow))
{
center = i1save$center
r = i1save$r
}
if (r > pi / 2)
{
center = -center
r = pi - r
}
return(list(center = c(center), r = r))
}
#==================================================================================
LRTpval = function(resGREAT, resSMALL, n)
{
chi2 = max(n * log(sum(resGREAT ^ 2) / sum(resSMALL ^ 2)), 0)
return(pchisq(
q = chi2,
df = 1,
lower.tail = FALSE
))
}
#==================================================================================
vMFtest = function(x, R = 100)
{
d = nrow(x)
n = ncol(x)
sumx = apply(x, 1, sum)
rbar = norm(sumx, "2") / n
muMLE = sumx / norm(sumx, "2")
kappaMLE = (rbar * d - rbar ^ 3) / (1 - rbar ^ 2)
sp = getSubSphere(x = x, geodesic = "small")
center.s = sp$center
r.s = sp$r
radialdistances = acos(t(center.s) %*% x)
xi_sample = mean(radialdistances) / sd(radialdistances)
xi_vec = rep(0, R)
for (r in 1:R)
{
rdata = randvonMisesFisherm(d, n, kappaMLE)
sp = getSubSphere(x = rdata, geodesic = "small")
center.s = sp$center
r.s = sp$r
radialdistances = acos(t(center.s) %*% rdata)
xi_vec[r] = mean(radialdistances) / sd(radialdistances)
}
pvalue = mean(xi_vec > xi_sample)
return(pvalue)
}
#==================================================================================
geodmeanS1 = function(theta)
{
n = length(theta)
meancandi = mod(mean(theta) + 2 * pi * (0:(n - 1)) / n, 2 * pi)
theta = mod(theta, 2 * pi)
geodvar = rep(0, n)
for (i in 1:n)
{
v = meancandi[i]
dist2 = apply(cbind((theta - v) ^ 2, (theta - v + 2 * pi) ^ 2, (v - theta + 2 * pi) ^
2), 1, min)
geodvar[i] = sum(dist2)
}
m = min(geodvar)
ind = which.min(geodvar)
geodmean = mod(meancandi[ind], 2 * pi)
geodvar = geodvar[ind] / n
return(list(geodmean = geodmean, geodvar = geodvar))
}
#==================================================================================
PNSe2s = function(resmat, PNS)
{
dm = nrow(resmat)
n = ncol(resmat)
NSOrthaxis = rev(PNS$orthaxis[1:(dm - 1)])
NSradius = flipud0(matrix(PNS$dist, ncol = 1))
geodmean = PNS$orthaxis[[dm]]
res = resmat / repmat(flipud0(matrix(PNS$radii, ncol = 1)), 1, n)
T = t(rotMat(NSOrthaxis[[1]])) %*%
rbind(repmat(sin(NSradius[1] + matrix(res[2,], nrow = 1)), 2, 1) *
rbind(cos(geodmean + res[1,]), sin(geodmean + res[1,])),
cos(NSradius[1] + res[2,]))
if (dm > 2)
{
for (i in 1:(dm - 2))
{
T = t(rotMat(NSOrthaxis[[i + 1]])) %*%
rbind(repmat(sin(NSradius[i + 1] + matrix(
res[i + 2,], nrow = 1
)), 2 + i, 1) * T,
cos(NSradius[i + 1] + res[i + 2,]))
}
}
if (!is.null(PNS$basisu))
{
T = PNS$basisu %*% T
}
return(T)
}
#==================================================================================
PNSs2e = function(spheredata, PNS)
{
if (nrow(spheredata) != length(PNS$mean))
{
cat(" Error from PNSs2e() \n")
cat(" Dimensions of the sphere and PNS decomposition do not match")
return(NULL)
}
if (!is.null(PNS$basisu))
{
spheredata = t(PNS$basisu) %*% spheredata
}
kk = nrow(spheredata)
n = ncol(spheredata)
Res = matrix(0, kk - 1, n)
currentSphere = spheredata
for (i in 1:(kk - 2))
{
v = PNS$orthaxis[[i]]
r = PNS$dist[i]
res = acos(t(v) %*% currentSphere) - r
Res[i,] = res
NestedSphere = rotMat(v) %*% currentSphere
currentSphere = as.matrix(NestedSphere[1:(kk - i),]) /
repmat(matrix(sqrt(1 - NestedSphere[nrow(NestedSphere),] ^
2), nrow = 1), kk - i, 1)
}
S1toRadian = atan2(currentSphere[2,], currentSphere[1,])
devS1 = mod(S1toRadian - rev(PNS$orthaxis)[[1]] + pi, 2 * pi) - pi
Res[kk - 1,] = devS1
EuclidData = flipud0(repmat(PNS$radii, 1, n) * Res)
return(EuclidData)
}
#==================================================================================
randvonMisesFisherm = function(m, n, kappa, mu = NULL)
{
if (is.null(mu))
{
muflag = FALSE
} else {
muflag = TRUE
}
if (m < 2)
{
print("Message from randvonMisesFisherm(): dimension m must be > 2")
print("Message from randvonMisesFisherm(): Set m to be 2")
m = 2
}
if (kappa < 0)
{
print("Message from randvonMisesFisherm(): kappa must be >= 0")
print("Message from randvonMisesFisherm(): Set kappa to be 0")
kappa = 0
}
b = (-2 * kappa + sqrt(4 * kappa ^ 2 + (m - 1) ^ 2)) / (m - 1)
x0 = (1 - b) / (1 + b)
c = kappa * x0 + (m - 1) * log(1 - x0 ^ 2)
nnow = n
w = c()
while (TRUE)
{
ntrial = max(round(nnow * 1.2), nnow + 10)
Z = rbeta(n = ntrial,
shape1 = (m - 1) / 2,
shape2 = (m - 1) / 2)
U = runif(ntrial)
W = (1 - (1 + b) * Z) / (1 - (1 - b) * Z)
indicator = kappa * W + (m - 1) * log(1 - x0 * W) - c >= log(U)
if (sum(indicator) >= nnow)
{
w1 = W[indicator]
w = c(w, w1[1:nnow])
break
} else {
w = c(w, W[indicator])
nnow = nnow - sum(indicator)
}
}
V = UNIFORMdirections(m - 1, n)
X = rbind(repmat(sqrt(1 - matrix(w, nrow = 1) ^ 2), m - 1, 1) * V, matrix(w, nrow = 1))
if (muflag)
{
mu = mu / norm(mu, "2")
X = t(rotMat(mu)) %*% X
}
return(X)
}
#==================================================================================
UNIFORMdirections = function(m, n)
{
V = matrix(0, m, n)
nr = matrix(rnorm(m * n), nrow = m)
for (i in 1:n)
{
while (TRUE)
{
ni = sum(nr[, i] ^ 2)
if (ni < 1e-10)
{
nr[, i] = rnorm(m)
} else {
V[, i] = nr[, i] / sqrt(ni)
break
}
}
}
return(V)
}
#==================================================================================
trans.subsphere = function(x, center)
{
return(repmat(1 / sqrt(1 - (t(
center
) %*% x) ^ 2), length(center) - 1, 1) *
(rotMat(center)[-length(center),] %*% x))
}
#==================================================================================
get.prinarc.value = function(PNS, arc, res)
{
d = length(PNS$orthaxis)
n = length(res)
prinarc = matrix(NA, d + 1, n)
for (g in 1:n)
{
newres = matrix(0, d, 1)
newres[arc] = res[g]
T = PNSe2s(newres, PNS)
prinarc[, g] = T
}
return(prinarc)
}
#==================================================================================
get.prinarc = function(resmat, PNS, arc, n, boundary.data = FALSE)
{
d = nrow(resmat)
if (boundary.data)
{
mn = min(resmat[arc,])
mx = max(resmat[arc,])
} else {
mn = -pi * tail(PNS$radii, arc)[1]
mx = pi * tail(PNS$radii, arc)[1]
}
prinarcgrid = seq(mn, mx, length = n)
prinarc = matrix(NA, d + 1, n)
for (g in 1:n)
{
newres = matrix(0, d, 1)
newres[arc] = prinarcgrid[g]
T = PNSe2s(newres, PNS)
prinarc[, g] = T
}
return(prinarc)
}
#==================================================================================
get.prinarc.subsphere = function(resmat,
PNS,
arc,
n,
subsphere = arc,
boundary.data = FALSE)
{
if (subsphere < arc)
{
stop("Error: subsphere >= arc.")
}
if (subsphere < 1)
{
stop("Error: subsphere >= 1.")
}
prinarc = get.prinarc(
resmat = resmat,
PNS = PNS,
arc = arc,
n = n,
boundary.data = boundary.data
)
d = nrow(resmat)
prinarc.sub = prinarc
if (subsphere < d)
{
for (i in 1:(d - subsphere))
{
prinarc.sub = trans.subsphere(x = prinarc.sub, center = PNS$orthaxis[[i]])
}
}
return(prinarc.sub)
}
#==================================================================================
get.data.subsphere = function(resmat, PNS, x, subsphere)
{
if (subsphere < 1)
{
stop("Error: subsphere >= 1.")
}
d = nrow(resmat)
x.sub = x
if (subsphere < d)
{
for (i in 1:(d - subsphere))
{
x.sub = trans.subsphere(x = x.sub, center = PNS$orthaxis[[i]])
}
}
return(x.sub)
}
#==================================================================================
mod = function(x, y)
{
return(x %% y)
}
#==================================================================================
repmat = function(x, m, n)
{
return(kronecker(matrix(1, m, n), x))
}
#==================================================================================
flipud0 = function(x)
{
return(apply(x, 2, rev))
}
#==================================================================================
sphere.obj = function(center, x, is.greatCircle)
{
di = sqrt(apply((x - repmat(
matrix(center, ncol = 1), 1, ncol(x)
)) ^ 2, 2, sum))
if (is.greatCircle)
{
r = pi / 2
} else {
r = mean(di)
}
sum((di - r) ^ 2)
}
#==================================================================================
sphere.res = function(center, x, is.greatCircle)
{
center = c(center)
xmc = x - center
di = sqrt(apply(xmc ^ 2, 2, sum))
if (is.greatCircle)
{
r = pi / 2
} else {
r = mean(di)
}
(di - r)
}
#==================================================================================
sphere.jac = function(center, x, is.greatCircle)
{
center = c(center)
xmc = x - center
di = sqrt(apply(xmc ^ 2, 2, sum))
di.vj = -xmc / repmat(matrix(di, nrow = 1), length(center), 1)
if (is.greatCircle)
{
c(t(di.vj))
} else {
r.vj = apply(di.vj, 1, mean)
c(t(di.vj - repmat(matrix(r.vj, ncol = 1), 1, ncol(x))))
}
}
#==================================================================================
sphereFit = function(x,
initialCenter = NULL,
geodesic = "small")
{
if (is.null(initialCenter))
{
initialCenter = apply(x, 1, mean)
}
op = nls.lm(
par = initialCenter,
fn = sphere.res,
jac = sphere.jac,
x = x,
is.greatCircle = ifelse(geodesic == "great", TRUE, FALSE),
control = nls.lm.control(maxiter = 300)
)
center = coef(op)
di = sqrt(apply((x - repmat(
matrix(center, ncol = 1), 1, ncol(x)
)) ^ 2, 2, sum))
if (geodesic == "great")
{
r = pi / 2
} else {
r = mean(di)
}
list(center = center, r = r)
}
#==================================================================================
tr = function(x)
{
return(sum(diag(x)))
}
#==================================================================================
Enormalize = function(x)
{
return(x / Enorm(x))
}
#==================================================================================
sphere2pcscore = function(x)
{
n = nrow(x)
p = ncol(x)
scores = matrix(NA, n, p - 1)
for (i in 1:n)
{
size = acos(x[i, 1])
scores[i,] = (size / sin(size)) * x[i, 2:p]
}
return(scores)
}
#==================================================================================
pcscore2sphere = function(n.pc, X.hat, S, Tan, V)
{
d = nrow(Tan)
n = ncol(Tan)
W = matrix(NA, d, n)
for (i in 1:n)
{
W[, i] = acos(tr(S[, , i] %*% t(X.hat))) * Tan[, i] / sqrt(sum(Tan[, i] ^
2))
}
lambda = matrix(NA, n, d)
for (i in 1:n)
{
for (j in 1:d)
{
lambda[i, j] = sum(W[, i] * V[, j])
}
}
U = matrix(0, n, d)
for (i in 1:n)
{
for (j in 1:n.pc)
{
U[i,] = U[i,] + lambda[i, j] * V[, j]
}
}
S.star = matrix(NA, n, n.pc + 1)
for (i in 1:n)
{
U.norm = sqrt(sum(U[i,] ^ 2))
S.star[i,] = c(cos(U.norm),
sin(U.norm) / U.norm * lambda[i, 1:n.pc])
}
return(S.star)
}
#==================================================================================
pc2sphere = function(x, n.pc)
{
k = dim(x)[1]
m = dim(x)[2]
n = dim(x)[3]
if (n < ((k - 1) * m))
{
stop("Error: n must be >= (k - 1) * m.")
}
GPAout = procGPA(
x = x,
scale = TRUE,
reflect = FALSE,
tangentcoords = "partial",
distances = FALSE
)
cat(
"First ",
n.pc,
" principal components explain ",
round(sum(GPAout$percent[1:n.pc])),
"% of total variance. \n",
sep = ""
)
H = defh(k - 1)
X.hat = H %*% GPAout$mshape
S = array(NA, c(k - 1, m, n))
for (i in 1:n)
{
S[, , i] = H %*% GPAout$rotated[, , i]
}
T.c = GPAout$tan - apply(GPAout$tan, 1, mean)
out = pcscore2sphere(
n.pc = n.pc,
X.hat = X.hat,
S = S,
Tan = T.c,
V = GPAout$pcar
)
return(list(spheredata = out, GPAout = GPAout))
}
#==================================================================================
rot.mat = function(Y,
X,
reflect = FALSE,
center = TRUE)
{
svd.out = svd(t(X) %*% Y)
R = svd.out$u %*% t(svd.out$v)
if (!reflect)
{
if (det(R) < 0)
{
u = svd.out$u
v = svd.out$v
if (det(u) < 0)
{
u[, dim(u)[2]] = -u[, dim(u)[2]]
} else if (det(v) < 0) {
v[, dim(v)[2]] = -v[, dim(v)[2]]
}
R = u %*% t(v)
}
}
return(R)
}
#==================================================================================
Procrustes.dist.full = function(x1, x2)
{
m = ncol(x1)
z1 = preshape(x1)
z2 = preshape(x2)
Q = t(z1) %*% z2 %*% t(z2) %*% z1
ev = eigen(Q)$values
sign = ifelse(det(t(z1) %*% z2) >= 0, 1,-1)
dF = sqrt(abs(1 - sum(sqrt(abs(
ev[1:(m - 1)]
)), sign * sqrt(abs(
ev[m]
))) ^ 2))
R = rot.mat(
Y = z2,
X = z1,
reflect = FALSE,
center = FALSE
)
scale = sum(svd(t(z1) %*% z2)$d)
return(list(dF = dF, R = R, scale = scale))
}
#==================================================================================
tangent.coords.partial = function(x, p)
{
k = nrow(x)
m = ncol(x)
if (abs(norm(p, "F") - 1) > 1e-15)
{
print("||p|| is not 1. Normalised one is used.")
p = Enormalize(p)
}
tmp = Procrustes.dist.full(x, p)
R = tmp$R
scale = tmp$scale
pre.p = preshape(p)
pre.x = preshape(x)
ident = diag(k * m - m)
tan = (ident - matrix(pre.p) %*% t(c(pre.p))) %*% c(pre.x %*% R)
tan.scale = (ident - matrix(pre.p) %*% t(c(pre.p))) %*% c(pre.x %*% R * scale)
return(list(
tan = c(tan),
tan.scale = c(tan.scale),
R = R,
scale = scale
))
}
#==================================================================================
shape.pcscores = function(PCAout, x, tangentcoords = "partial")
{
if (tangentcoords == "partial")
{
if (abs(norm(PCAout$mshape, "F") - 1) > 1e-15)
{
print("||PCAout$mshape|| is not 1. Normalised one is used.")
mshape = Enormalize(PCAout$mshape)
} else {
mshape = PCAout$mshape
}
if (abs(norm(x, "F") - 1) > 1e-15)
{
print("||x|| is not 1. Normalised one is used.")
x = Enormalize(x)
}
opa.out = procOPA(mshape, x, scale = FALSE)
matched = opa.out$Bhat
tan.out = tangent.coords.partial(matched, mshape)
mean.tan = apply(PCAout$tan, 1, mean)
scores = t(tan.out$tan - mean.tan) %*% PCAout$pcar
scores.scale = t(tan.out$tan.scale - mean.tan) %*% PCAout$pcar
return(
list(
rotated = matched,
tan = tan.out$tan,
tan.scale = tan.out$tan.scale,
scores = c(scores),
scores.scale = c(scores.scale)
)
)
}
}
#==================================================================================
shape.pcscores.partial = function(PCAout, x)
{
n = dim(x)[3]
scores = c()
for (i in 1:n)
{
s = shape.pcscores(PCAout, x[, , i], tangentcoords = "partial")
scores = rbind(scores, s$scores)
}
return(scores)
}
#==================================================================================
plotshapes3d.pns = function(x,
type = "p",
col = "black",
size = 5,
aspect = "iso",
joinline = TRUE,
col.joinline = "#d4d2d2",
lwd.joinline = 0.5,
tick = FALSE,
labels.tick = FALSE,
xlab = "",
ylab = "",
zlab = "")
{
k = dim(x)[1]
n = dim(x)[3]
aa = c()
bb = c()
cc = c()
for (i in 1:n)
{
aa = c(aa, x[, 1, i])
bb = c(bb, x[, 2, i])
cc = c(cc, x[, 3, i])
}
xlim = range(aa)
ylim = range(bb)
zlim = range(cc)
plot3d(
x[, , 1],
type = "n",
xlab = "",
ylab = "",
zlab = "",
box = FALSE,
axes = FALSE,
aspect = aspect,
xlim = xlim,
ylim = ylim,
zlim = zlim
)
for (i in 1:n)
{
plot3d(
x[, , i],
type = type,
col = col,
size = size,
add = TRUE
)
}
if (tick)
{
axis3d(
edge = 'x',
labels = labels.tick,
tick = TRUE,
pos = c(NA, 0, 0),
cex = 0.6,
lwd = 0.5
)
axis3d(
edge = 'y',
labels = labels.tick,
tick = TRUE,
pos = c(0, NA, 0),
cex = 0.6,
lwd = 0.5
)
axis3d(
edge = 'z',
labels = labels.tick,
tick = TRUE,
pos = c(0, 0, NA),
cex = 0.6,
lwd = 0.5
)
} else {
}
r = cbind(xlim, ylim, zlim)
pos = r[2,] + apply(r, 2, diff) / 20
text3d(pos[1], 0, 0, texts = xlab, cex = 0.8)
text3d(0, pos[2], 0, texts = ylab, cex = 0.8)
text3d(0, 0, pos[3], texts = zlab, cex = 0.8)
if (joinline)
{
for (i in 1:n)
{
lines3d(x[, , i], col = col.joinline, lwd = lwd.joinline)
}
}
}
#==================================================================================
Plot3D = function(x,
type = "s",
col = "black",
size = 1.2,
aspect = "iso",
joinline = FALSE,
col.joinline = "#d4d2d2",
lwd.joinline = 0.5,
tick = TRUE,
tick.boundary = FALSE,
labels.tick = TRUE,
xlab = "",
ylab = "",
zlab = "")
{
n = nrow(x)
plot3d(
x,
type = "n",
xlab = "",
ylab = "",
zlab = "",
box = FALSE,
axes = FALSE,
aspect = aspect
)
plot3d(
x,
type = type,
col = col,
size = size,
add = TRUE
)
if (tick)
{
axis3d(
edge = 'x',
labels = labels.tick,
tick = TRUE,
pos = c(NA, 0, 0),
cex = 0.6,
lwd = 0.5
)
axis3d(
edge = 'y',
labels = labels.tick,
tick = TRUE,
pos = c(0, NA, 0),
cex = 0.6,
lwd = 0.5
)
axis3d(
edge = 'z',
labels = labels.tick,
tick = TRUE,
pos = c(0, 0, NA),
cex = 0.6,
lwd = 0.5
)
}
if (tick.boundary)
{
tks = pretty(x[, 1], n = 10)
axis3d(
edge = 'x',
labels = labels.tick,
tick = TRUE,
at = c(tks[1], tks[length(tks)]),
pos = c(NA, 0, 0),
cex = 0.6,
lwd = 0.5
)
tks = pretty(x[, 2], n = 10)
axis3d(
edge = 'y',
labels = labels.tick,
tick = TRUE,
at = c(tks[1], tks[length(tks)]),
pos = c(0, NA, 0),
cex = 0.6,
lwd = 0.5
)
tks = pretty(x[, 3], n = 10)
axis3d(
edge = 'z',
labels = labels.tick,
tick = TRUE,
at = c(tks[1], tks[length(tks)]),
pos = c(0, 0, NA),
cex = 0.6,
lwd = 0.5
)
}
r = apply(x, 2, range)
pos = r[2,] + apply(r, 2, diff) / 20
text3d(pos[1], 0, 0, texts = xlab, cex = 0.8)
text3d(0, pos[2], 0, texts = ylab, cex = 0.8)
text3d(0, 0, pos[3], texts = zlab, cex = 0.8)
if (joinline)
{
lines3d(x, col = col.joinline, lwd = lwd.joinline)
}
}
#==================================================================================
col2RGB = function(col, alpha = 255)
{
n = length(col)
out = c()
for (i in 1:n)
{
out[i] = rgb(
red = col2rgb(col[i])[1],
green = col2rgb(col[i])[2],
blue = col2rgb(col[i])[3],
alpha = alpha,
maxColorValue = 255
)
}
return(out)
}
#==================================================================================
project.subsphere = function(x, center, r)
{
n = ncol(x)
d = nrow(x)
x.proj = matrix(NA, d, n)
for (i in 1:n)
{
rho = acos(sum(x[, i] * center))
x.proj[, i] = (sin(r) * x[, i] + sin(rho - r) * center) / sin(rho)
}
return(x.proj)
}
##############################################################end of PNS###########
##### Penalised Euclidean Distance Regression
#==================================================================================
ped <- function(X, Y, method = c("AIC")) {
if (method == "AIC") {
aicmin <- 999999999
for (lam in c(0.2, 0.5, 1.0)) {
for (cofp in c(0.75, 1, 1.35, 1.5)) {
out <-
pedreg(
X,
Y,
nlambda = 1,
constc0 = 1.1,
constc1 = cofp,
lambdainit = lam
)
if (out$aic < aicmin) {
minout <- out
mincofp <- cofp
aicmin <- out$aic
}
}
}
out <- minout
}
if (method == "BIC") {
bicmin <- 999999999
for (lam in c(0.2, 0.5, 1.0)) {
for (cofp in c(0.75, 1, 1.35, 1.5)) {
out <-
pedreg(
X,
Y,
nlambda = 1,
constc0 = 1.1,
constc1 = cofp,
lambdainit = lam
)
if (out$bic < bicmin) {
minout <- out
mincofp <- cofp
bicmin <- out$bic
}
}
}
out <- minout
}
if (method == "khat") {
aicmin <- 999999999
for (lam in c(0.2, 0.5, 1.0)) {
for (cofp in c(0.75, 1, 1.25, 1.5)) {
out <-
pedreg(
X,
Y,
nlambda = 1,
constc0 = 1.1,
constc1 = cofp,
lambdainit = lam
)
if (-out$khat < aicmin) {
minout <- out
mincofp <- cofp
aicmin <- -out$khat
}
}
}
out <- minout
}
if (method == "CV") {
n <- length(Y)
cvmin <- 999999999
for (lam in c(0.2, 0.5, 1.0)) {
for (cofp in c(0.75, 1, 1.25, 1.5)) {
cverr <- 0
for (jj in 1:10) {
subsample <- ((jj - 1) * 10 + 1):((jj - 1) * 10 + 10)
out <-
pedreg(
X[-subsample, ],
Y[-subsample],
nlambda = 1,
constc0 = 1.1,
constc1 = cofp,
lambdainit = lam
)
cverr <-
cverr + Enorm(Y[subsample] - out$intercept - X[subsample, ] %*% out$betahat) **
2
}
if (cverr < cvmin) {
minout <- out
minlam <- lam
mincofp <- cofp
cvmin <- cverr
}
}
}
out <-
pedreg(
X,
Y,
nlambda = 1,
constc0 = 1.1,
constc1 = mincofp,
lambdainit = minlam
)
}
out1 <- list(
betahat = 0,
yhat = 0,
lambda = 0,
coef = 0,
resid = 0
)
out1$intercept <- out$intercept
out1$coef <- c(out$intercept, out$betahat)
out1$betahat <- out$betahat
out1$lambda <- out$lambda
out1$delta <- mincofp
out1$yhat <- out$yhat
out1$resid <- Y - out$yhat
out1
}
###########################function for PED#####################
#==================================================================================
pedreg <-
function(X,
Y,
constc0 = 1.1,
constc1 = 1.35,
alpha = 0.05,
LMM = 50,
MIT = 10000,
NUM_METHOD = 1,
nlambda = 1,
lambdamax = 1,
PLOT = TRUE,
BIC = FALSE,
lambdainit = 1) {
# NUM_METHOD = 1 = L-BFGS-B
# LMM = Parameter M in L-BFGS method 1
# MIT = Max iterations for optimization
p <- dim(X)[2]
n <- dim(X)[1]
constc <- constc0
Ymean <- mean(Y)
Ysd <- sd(Y)
Yinit <- Y
pinit <- p
ans0 <- rep(0, times = pinit)
Xorig <- X
Yorig <- Y
vm <- rep(0, times = ncol(X))
vsd <- rep(0, times = ncol(X))
for (i in 1:ncol(X)) {
vm[i] <- mean(X[, i])
}
for (i in 1:ncol(X)) {
vsd[i] <- sd(X[, i])
}
#standardize to sphere
X <- scale(X) / sqrt(n - 1)
Y <- scale(Y) / sqrt(n - 1)
X0 <- X
Y0 <- Y
lambdainit1 <- constc / sqrt(n - 1) * sqrt(sqrt(p)) / n * qnorm(1 - alpha /
(2 * p))
if (nlambda == 1) {
lambdainit1 <- lambdainit
}
METHOD1 <- "L-BFGS-B"
xi <- -9999999
c1 <- 1
nlam <- nlambda
betamat <- matrix(0, p, nlam)
betamat.sparse <- betamat
lambdamat <- rep(0, times = nlam)
ximat <- rep(0, times = nlam)
aic <- rep(0, times = nlam)
bic <- rep(0, times = nlam)
npar <- rep(0, times = nlam)
selectmat <- betamat
#cat(c("Lambda iteration (out of ",nlam,"):"))
for (ilam in (nlam:1)) {
#cat(c(ilam," "))
if (nlam == 1) {
lambda <- lambdainit1
}
if (nlam > 1) {
c1 <- sqrt(n) + (ilam - 1) / (nlam - 1) * 1 / lambdainit1
c1 <-
sqrt(n) + ((ilam - 1) / (nlam - 1)) * 1 / lambdainit1 * (lambdamax - lambdainit1 *
sqrt(n))
lambda <- lambdainit1 * c1
}
if (ilam == nlam) {
x0 <- rep(1 / sqrt(p), times = p)
}
if (ilam != nlam) {
x0 <- betahat + rnorm(p) / sqrt(p)
}
#x0<-rnorm(p)/sqrt(p)
pedfun <- function(pars, Y = 0, X = 0) {
p <- length(pars)
pars <- matrix(pars, p, 1)
ped <- Enorm(Y - X %*% pars) + lambda * sqrt(Enorm(pars) * sum(abs(pars)))
ped
}
pedgrad <- function(pars, X = 0, Y = 0) {
GM <- sqrt(Enorm(pars) * sum(abs(pars)))
gradL <- rep(0, times = p)
gradL <- -t(X) %*% (Y - X %*% pars) / Enorm(Y - X %*% pars)
gradL <- gradL + matrix(
lambda / 2 * pars / Enorm(pars) * sum(abs(pars)) / GM +
lambda / 2 * sign(pars) * Enorm(pars) / GM ,
p,
1
)
c(gradL)
}
if (NUM_METHOD == 1) {
repeat {
#,ndeps=1e-3,factr=1e-5,pgtol=1e-5
res2 <-
optim(
par = x0,
fn = pedfun,
gr = pedgrad,
method = METHOD1,
control = list(lmm = LMM, maxit = MIT),
X = X,
Y = Y
)
betahat <- res2$par
if (res2$convergence == 0) {
break
}
x0 <- rnorm(p) / sqrt(p)
}
}
oldxi <- xi
xi <-
sqrt(Enorm(betahat) / sum(abs(betahat))) - sqrt(n) / (constc * c1 * p ^
(1 / 4))
dif <- (xi - oldxi)
REGC <- 0.0001
betamat[, ilam] <- betahat / (Enorm(betahat) + REGC)
lambdamat[ilam] <- lambda
ximat[ilam] <- xi
ximat[ilam] <- sqrt(Enorm(betahat) / sum(abs(betahat)))
MM <- constc1 / (sqrt(n))
select <- (abs(betahat) / (Enorm(betahat) + REGC) > MM)
selectmat[, ilam] <- select
betamat.sparse[, ilam] <- betamat[, ilam]
betamat.sparse[select == FALSE, ilam] <- 0 * betamat[select == FALSE, ilam]
pp <- sum(select)
npar[ilam] <- pp
bic[ilam] <- log(Enorm(Y) ** 2 / n) * (n) + log(n) * (1)
aic[ilam] <- log(Enorm(Y) ** 2 / n) * (n) + 2 * (1)
if (sum(select) > 0) {
aa <- lm(Y ~ X[, c(1:p)[select]] - 1)
pred <- predict(aa)
# Use AIC with finite sample correction
aic[ilam] <-
log(Enorm(Y - pred) ** 2 / n) * (n) + 2 * (pp + 1) + 2 * (pp + 1) * (pp +
2) / (n - pp - 2)
# Use AIC/BIC
bic[ilam] <- log(Enorm(Y - pred) ** 2 / n) * (n) + log(n) * (pp + 1)
aic[ilam] <- log(Enorm(Y - pred) ** 2 / n) * (n) + 2 * (pp + 1)
}
}
best <- 1
if (nlam > 1) {
###########################################choose best via AIC
best <- c(1:nlam)[aic == min(aic)][1]
select <- as.logical(selectmat[, best])
lambdaaic <- lambdamat[best]
###########################################choose best via Corollary 1 with khat
best2 <- nlam
if ((sum(ximat > 0.25)) > 0) {
best2 <- c(1:nlam)[(ximat) > 0.25][1]
}
#################### biggest xi from Corollary 1
xism <- (ximat)
if (sum(diff(xism) < 0.01) > 0) {
best2 <- c(2:nlam)[diff(xism) < 0.01][1] - 1
}
############## biggest sqrt( Enorm(beta) / norm(beta)_1 )
best2 <- c(1:nlam)[ximat == max(ximat)]
selectcor <- as.logical(selectmat[, best2])
lambdacor1 <- lambdamat[best2]
if (BIC == FALSE) {
best <- best2
select <- selectcor
}
}
################last part######estimate with reduced p###########
if (sum(select) > 0) {
X <- as.matrix(X[, select])
p <- sum(select)
p14 <- sqrt(sqrt(p))
final <- c(1:pinit)[select]
}
#######
lambda <- constc / sqrt(sqrt(p)) / sqrt(n) * qnorm(1 - alpha / (2 * p))
if (sum(select) > 0) {
x0 <- rep(1 / p, times = p)
if (NUM_METHOD == 1) {
repeat {
res2 <-
optim(
par = x0,
fn = pedfun,
gr = pedgrad,
method = METHOD1,
control = list(lmm = LMM, maxit = MIT),
X = X,
Y = Y
)
betahat <- res2$par
if (res2$convergence == 0) {
break
}
x0 <- rnorm(p) / p
}
}
ans0[final] <- betahat
}
#ind<-which(abs(ans0/Enorm(ans0))<10^(-5))
#ans0[ind]<-0
out <- list(betahat = 0,
yhat = 0,
lambda = 0)
out$betahatscale <- ans0
out$yhatscale <- c(X0 %*% ans0)
if (nlam > 1) {
out$lambdacor1 <- lambdacor1
out$lambdaaic <- lambdaaic
out$betamat.sparse <- betamat.sparse
out$betamat.rescale <- betamat
out$betamat <- betamat
for (i in 1:nlam) {
out$betamat.rescale[, i] <- c(out$betamat[, i] / vsd) * sd(Yorig)
}
out$lambdamat <- lambdamat
out$ximat <- ximat
out$MM <- MM
out$fmax <- res2$value
out$npar <- npar
out$selectmat <- selectmat
}
#use AIC
out$aic <- aic
#use BIC
out$bic <- bic
#use khat
out$khat <- ximat
out$lambdath3 <- lambdainit1 * sqrt(n)
out$lambda <- out$lambdacor1
if (BIC == TRUE) {
out$lambda <- out$lambdaaic
}
if (nlam == 1) {
out$lambda <- lambdainit1
out$constc1 <- constc1
}
sol <- sd(Yorig) * c(ans0 / vsd)
inter <- drop(mean(Yorig) - sd(Yorig) * (vm / vsd) %*% ans0)
out$intercept <- drop(mean(Yorig) - sd(Yorig) * (vm / vsd) %*% ans0)
out$betahat <- sd(Yorig) * c(ans0 / vsd)
out$best <- best
out$Yinit <- Yinit
out$yhat <- Xorig %*% sol + inter
out
}
############################################################
#
# FUNCTIONS FOR CALCULATING NON-EUCLIDEAN MEANS AND DISTANCES
# OF COVARIANCE MATRICES
#
############################################################
# Log Euclidean mean: Sigma_L
#==================================================================================
estLogEuclid <- function(S, weights = 1) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
sum <- S[, , 1] * 0
for (j in 1:M) {
eS <- eigen(S[, , j], symmetric = TRUE)
sum <- sum +
weights[j] * eS$vectors %*% diag(log(eS$values)) %*% t(eS$vectors) /
sum(weights)
}
ans <- sum
eL <- eigen(ans, symmetric = TRUE)
eL$vectors %*% diag(exp(eL$values)) %*% t(eL$vectors)
}
#==================================================================================
estPowerEuclid <- function(S, weights = 1, alpha = 0.5) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
sum <- S[, , 1] * 0
for (j in 1:M) {
eS <- eigen(S[, , j], symmetric = TRUE)
sum <- sum +
weights[j] * eS$vectors %*% diag(abs(eS$values) ** alpha) %*% t(eS$vectors) /
sum(weights)
}
ans <- sum
eL <- eigen(ans, symmetric = TRUE)
eL$vectors %*% diag(abs(eL$values) ** (1 / alpha)) %*% t(eL$vectors)
}
# Riemannian (weighted mean) : Sigma_R
#==================================================================================
estLogRiem2 <- function(S, weights = 1) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
check <- 9
tau <- 1
Hold <- 99999
mu <- estLogEuclid(S, weights)
while (check > 0.0000000001) {
ev <- eigen(mu, symmetric = TRUE)
logmu <- ev$vectors %*% diag(log(ev$values)) %*% t(ev$vectors)
Hnew <- Re(Hessian2(S, mu, weights))
logmunew <- logmu + tau * Hnew
ev <- eigen(logmunew, symmetric = TRUE)
mu <- ev$vectors %*% diag(exp(ev$values)) %*% t(ev$vectors)
check <- Re(Enorm(Hold) - Enorm(Hnew))
if (check < 0) {
tau <- tau / 2
check <- 999999
}
Hold <- Hnew
}
mu
}
# Hessian used in calculating Sigma_R
#==================================================================================
Hessian2 <- function(S, Sigma, weights = 1) {
M <- dim(S)[3]
k <- dim(S)[1]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
ev0 <- eigen(Sigma, symmetric = TRUE)
shalf <-
ev0$vectors %*% (diag(sqrt((ev0$values)))) %*% t(ev0$vectors)
sumit <- matrix(0, k, k)
for (i in 1:(M)) {
ev2 <- eigen(shalf %*% solve(S[, , i]) %*% shalf, symmetric = TRUE)
sumit <-
sumit + weights[i] * ev2$vectors %*% diag (log((ev2$values))) %*% t(ev2$vectors) /
sum(weights)
}
- sumit
}
# Euclidean : Sigma_E
#==================================================================================
estEuclid <- function(S, weights = 1) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
sum <- S[, , 1] * 0
for (j in 1:M) {
sum <- sum + S[, , j] * weights[j] / sum(weights)
}
sum
}
# Cholesky mean : Sigma_C
#==================================================================================
estCholesky <- function(S, weights = 1) {
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
sum <- S[, , 1] * 0
for (j in 1:M) {
sum <- sum + t(chol(S[, , j])) * weights[j] / sum(weights)
}
cc <- sum
cc %*% t(cc)
}
#==================================================================================
ild_estSS <- function(S, weights = 1) {
M <- dim(S)[3]
k <- dim(S)[1]
H <- defh(k)
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
Q <- array(0, c(k + 1, k, M))
for (j in 1:M) {
Q[, , j] <- t(H) %*% (rootmat(S[, , j]))
}
ans <- procWGPA(
Q,
fixcovmatrix = diag(k + 1),
scale = FALSE,
reflect = TRUE,
sampleweights = weights
)
H %*% ans$mshape %*% t(H %*% ans$mshape)
}
#==================================================================================
ild_estShape <- function(S, weights = 1) {
M <- dim(S)[3]
k <- dim(S)[1]
H <- defh(k)
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
Q <- array(0, c(k + 1, k, M))
for (j in 1:M) {
Q[, , j] <- t(H) %*% (rootmat(S[, , j]))
}
ans <- procWGPA(
Q,
fixcovmatrix = diag(k + 1),
scale = TRUE,
reflect = TRUE,
sampleweights = weights
)
H %*% ans$mshape %*% t(H %*% ans$mshape)
}
#==================================================================================
estRiemLe <- function(S, weights) {
M <- dim(S)[3]
k <- dim(S)[1]
if (M != 2)
print("Sorry - Calculation not implemented for M>2 yet")
if (M == 2) {
P1 <- S[, , 1]
P2 <- S[, , 2]
detP1 <- prod(eigen(P1)$values)
detP2 <- prod(eigen(P2)$values)
P1 <- P1 / (detP1) ^ (1 / k)
P2 <- P2 / (detP2) ^ (1 / k)
P1inv <- solve(P1)
P12sq <- P1inv %*% P2 %*% P2 %*% P1inv
tem <- eigen(P12sq, symmetric = TRUE)
A2 <- tem$vectors %*% diag(log(tem$values)) %*% t(tem$vectors)
logPs2 <- weights[2] * A2
tem2 <- eigen(logPs2, symmetric = TRUE)
Ps2 <- tem2$vectors %*% diag(exp(tem2$values)) %*% t(tem2$vectors)
P12s <- P1 %*% Ps2 %*% P1
tem3 <- eigen(P12s, symmetric = TRUE)
P12sA <-
tem3$vectors %*% diag(sqrt(tem3$values)) %*% t(tem3$vectors)
Ptildes <- (detP1 * (detP2 / detP1) ^ weights[2]) ^ (1 / k) * P12sA
Ptildes
}
}
##########distances#################################
#==================================================================================
distRiemPennec <- function(P1, P2) {
eig <- eigen(P1, symmetric = TRUE)
P1half <- eig$vectors %*% diag(sqrt(eig$values)) %*% t(eig$vectors)
P1halfinv <- solve(P1half)
AA <- P1halfinv %*% P2 %*% P1halfinv
tem <- eigen(AA, symmetric = TRUE)
A2 <- tem$vectors %*% diag(log(tem$values)) %*% t(tem$vectors)
dd <- Enorm(A2)
dd
}
#==================================================================================
distLogEuclidean <- function(P1, P2) {
eig <- eigen(P1, symmetric = TRUE)
logP1 <- eig$vectors %*% diag(log(eig$values)) %*% t(eig$vectors)
tem <- eigen(P2, symmetric = TRUE)
logP2 <- tem$vectors %*% diag(log(tem$values)) %*% t(tem$vectors)
dd <- Enorm(logP1 - logP2)
dd
}
#==================================================================================
distRiemannianLe <-
function(P1, P2) {
dd <- distRiemPennec(P1 %*% t(P1), P2 %*% t(P2)) / 2
dd
}
#==================================================================================
ild_distProcrustesSizeShape <-
function (P1, P2)
{
H <- defh(dim(P1)[1])
Q1 <- t(H) %*% rootmat(P1)
Q2 <- t(H) %*% rootmat(P2)
ans <- sqrt(
centroid.size(Q1) ^ 2 + centroid.size(Q2) ^ 2 - 2 *
centroid.size(Q1) * centroid.size(Q2) * cos(riemdist(Q1,
Q2, reflect = TRUE))
)
ans
}
#==================================================================================
ild_distProcrustesFull <- function(P1, P2) {
H <- defh(dim(P1)[1])
Q1 <- t(H) %*% rootmat(P1)
Q2 <- t(H) %*% rootmat(P2)
ans <- riemdist(Q1, Q2, reflect = TRUE)
ans
}
#==================================================================================
distPowerEuclidean <- function(P1, P2, alpha = 1 / 2) {
if (alpha != 0) {
eS <- eigen(P1, symmetric = TRUE)
Q1 <- eS$vectors %*% diag(abs(eS$values) ^ alpha) %*% t(eS$vectors)
eS <- eigen(P2, symmetric = TRUE)
Q2 <- eS$vectors %*% diag(abs(eS$values) ^ alpha) %*% t(eS$vectors)
dd <- Enorm(Q1 - Q2) / abs(alpha)
}
if (alpha == 0) {
dd <- distLogEuclidean(P1, P2)
}
dd
}
#==================================================================================
ild_distCholesky <- function(P1, P2) {
H <- defh(dim(P1)[1])
Q1 <- t(H) %*% t(chol(P1))
Q2 <- t(H) %*% t(chol(P2))
ans <- Enorm(Q1 - Q2)
ans
}
#==================================================================================
distEuclidean <- function(P1, P2) {
ans <- Enorm(P1 - P2)
ans
}
##################
#==================================================================================
distcov <- function(S1,
S2 ,
method = "Riemannian",
alpha = 1 / 2) {
if (method == "Procrustes") {
dd <- distProcrustesSizeShape(S1, S2)
}
if (method == "ProcrustesShape") {
dd <- distProcrustesFull(S1, S2)
}
if (method == "Riemannian") {
dd <- distRiemPennec(S1, S2)
}
if (method == "Cholesky") {
dd <- distCholesky(S1, S2)
}
if (method == "Power") {
dd <- distPowerEuclidean(S1, S2, alpha)
}
if (method == "Euclidean") {
dd <- distEuclidean(S1, S2)
}
if (method == "LogEuclidean") {
dd <- distLogEuclidean(S1, S2)
}
if (method == "RiemannianLe") {
dd <- distRiemannianLe(S1, S2)
}
dd
}
#==================================================================================
estcov <-
function (S,
method = "Riemannian",
weights = 1,
alpha = 1 / 2,
MDSk = 2)
{
out <- list(
mean = 0,
sd = 0,
pco = 0,
eig = 0,
dist = 0
)
M <- dim(S)[3]
if (length(weights) == 1) {
weights <- rep(1, times = M)
}
if (method == "Procrustes") {
dd <- estSS(S, weights)
}
if (method == "ProcrustesShape") {
dd <- estShape(S, weights)
}
if (method == "Riemannian") {
dd <- estLogRiem2(S, weights)
}
if (method == "Cholesky") {
dd <- estCholesky(S, weights)
}
if (method == "Power") {
dd <- estPowerEuclid(S, weights, alpha)
}
if (method == "Euclidean") {
dd <- estEuclid(S, weights)
}
if (method == "LogEuclidean") {
dd <- estLogEuclid(S, weights)
}
if (method == "RiemannianLe") {
dd <- estRiemLe(S, weights)
}
out$mean <- dd
sum <- 0
for (i in 1:M) {
sum <-
sum + weights[i] * distcov(S[, , i], dd, method = method) ^ 2 / sum(weights)
}
out$sd <- sqrt(sum)
dist <- matrix(0, M, M)
for (i in 2:M) {
for (j in 1:(i - 1)) {
dist[i, j] <- distcov(S[, , i], S[, , j], method = method)
dist[j, i] <- dist[i, j]
}
}
out$dist <- dist
if (M > MDSk) {
ans <-
cmdscale(
dist,
k = MDSk,
eig = TRUE,
add = TRUE,
x.ret = TRUE
)
out$pco <- ans$points
out$eig <- ans$eig
if (MDSk > 2) {
shapes3d(out$pco[, 1:min(MDSk, 3)], axes3 = TRUE)
}
if (MDSk == 2) {
plot(out$pco,
type = "n",
xlab = "MDS1",
ylab = "MDS2")
text(out$pco[, 1], out$pco[, 2], 1:length(out$pco[, 1]))
}
}
out
}
rootmat <- function(P1) {
eS <- eigen(P1, symmetric = TRUE)
if (min(eS$values) < -0.001) {
print("Not positive-semi definite")
}
else{
Q1 <- eS$vectors %*% diag(sqrt(abs(eS$values))) %*% t(eS$vectors)
Q1
}
}
##########################
#==================================================================================
shapes.cva <- function(X ,
groups ,
scale = TRUE,
ncv = 2) {
g <- dim(table (groups))
ans <- procGPA(X , scale = scale)
if (scale == TRUE)
pp <- (ans$k - 1) * ans$m - (ans$m * (ans$m - 1) / 2) - 1
if (scale == FALSE)
pp <- (ans$k - 1) * ans$m - (ans$m * (ans$m - 1) / 2)
pracdim <- min(pp, ans$n - g)
out <- lda(ans$scores[, 1:pracdim] , groups)
print((out))
cv <- predict(out, dimen = 3)$x
if (dim(cv)[2] == 1) {
cv <- cbind(cv, rnorm(dim(cv)[1]) / 1000)
}
if (ncv == 2) {
eqscplot(cv,
type = "n",
xlab = "CV1",
ylab = "CV2")
text(cv, labels = groups)
}
if (ncv == 3) {
shapes3d(cv, color = groups, axes3 = TRUE)
}
cv
}
#==================================================================================
groupstack <- function(A1,
A2,
A3 = 0,
A4 = 0,
A5 = 0,
A6 = 0,
A7 = 0,
A8 = 0) {
out <- list(x = 0, groups = "")
dat <- abind(A1, A2)
group <- c(rep(1, times = dim(A1)[3]), rep(2, times = dim(A2)[3]))
if (is.array(A3)) {
dat <- abind(dat, A3)
group <- c(group, rep(3, times = dim(A3)[3]))
if (is.array(A4)) {
dat <- abind(dat, A4)
group <- c(group, rep(4, times = dim(A4)[3]))
if (is.array(A5)) {
dat <- abind(dat, A5)
group <- c(group, rep(5, times = dim(A5)[3]))
if (is.array(A6)) {
dat <- abind(dat, A6)
group <- c(group, rep(6, times = dim(A6)[3]))
if (is.array(A7)) {
dat <- abind(dat, A7)
group <- c(group, rep(7, times = dim(A7)[3]))
if (is.array(A8)) {
dat <- abind(dat, A8)
group <- c(group, rep(8, times = dim(A8)[3]))
}
}
}
}
}
}
out$x <- dat
out$groups <- group
out
}
###########################
#==================================================================================
procdist <- function(x, y, type = "full", reflect = FALSE) {
if (type == "full") {
out <- sin(riemdist(x, y, reflect = reflect))
}
if (type == "partial") {
out <- sqrt(2) * sqrt(abs(1 - cos(riemdist(x, y, reflect = reflect))))
}
if (type == "Riemannian") {
out <- riemdist(x, y, reflect = reflect)
}
if (type == "sizeandshape") {
out <- ssriemdist(x, y, reflect = reflect)
}
out
}
#==================================================================================
transformations <- function(Xrotated, Xoriginal) {
# outputs the translations, rotations and
# scalings for ordinary Procrustes rotation
# of each individual in Xoriginal to the
# Procrustes rotated individuals in Xrotated
X1 <- Xrotated
X2 <- Xoriginal
n <- dim(X1)[3]
m <- dim(X1)[2]
translation <- matrix(0, m, n)
scale <- rep(0, times = n)
rotation <- array(0, c(m, m, n))
for (i in 1:n) {
translation[, i] <- -apply(X2[, , i] - X1[, , i], 2, mean)
ans <- procOPA(X1[, , i], X2[, , i])
scale[i] <- ans$s
rotation[, , i] <- ans$R
}
out <- list(translation = 0,
scale = 0,
rotation = 0)
out$translation <- translation
out$scale <- scale
out$rotation <- rotation
out
}
#==================================================================================
iglogl <- function(x , lam, nlam) {
gamma <- abs(x[1])
alpha <- gamma / mean(1 / lam[1:nlam])
ll <-
-(gamma + 1) * sum(log(lam[1:nlam])) - alpha * sum (1 / lam[1:nlam]) +
nlam * gamma * log(alpha) - nlam * lgamma (gamma)
- ll
}
#==================================================================================
procWGPA <-
function(x,
fixcovmatrix = FALSE,
initial = "Identity",
maxiterations = 10,
scale = TRUE,
reflect = FALSE,
prior = "Exponential",
diagonal = TRUE,
sampleweights = "Equal") {
X <- x
priorargument <- prior
alpha <- "not estimated"
gamma <- "not estimated"
k <- dim(X)[1]
n <- dim(X)[3]
m <- dim(X)[2]
if (initial[1] == "Identity") {
Sigmak <- diag(k)
}
else{
if (initial[1] == "Rawdata") {
tol <- 0.0000000001
if (m == 2) {
Sigmak <- diag(diag(var(t(X[, 1, ]))) + diag(var(t(X[, 2, ])))) / 2 + tol
}
if (m == 3) {
Sigmak <-
diag(diag(var(t(X[, 1, ]))) + diag(var(t(X[, 2, ]))) + diag(var(t(X[, 3, ])))) /
3 + tol
}
}
else
{
Sigmak <- initial
}
}
mu <- procGPA(X, scale = scale)$mshape
#cat("Iteration 1 \n")
if (fixcovmatrix[1] != FALSE) {
Sigmak <- fixcovmatrix
}
ans <-
procWGPA1(
X,
mu,
metric = Sigmak,
scale = scale,
reflect = reflect,
sampleweights = sampleweights
)
if ((maxiterations > 1) && (fixcovmatrix[1] == FALSE)) {
ans0 <- ans
dif <- 999999
it <- 1
while ((dif > 0.00001) && (it < maxiterations)) {
it <- it + 1
if (it == 2) {
cat("Differences in norm of Sigma estimates... \n ")
}
if (prior[1] == "Identity") {
prior <- diag(k)
}
if (prior[1] == "Inversegamma") {
lam <- eigen(ans$Sigmak)$values
nlam <- min(c(n * m - m - 3, k - 3))
mu <- mean(1 / lam[1:(nlam)])
alpha <- 1 / mu
out <- nlm(iglogl,
p = c(1) ,
lam = lam,
nlam = nlam)
#print(out)
gamma <- abs(out$estimate[1])
alpha <- gamma / mean(1 / lam[1:nlam])
newmetric <-
n * m / (n * m + 2 * (1 + gamma)) * (ans$Sigmak + (2 * alpha / (n * m)) *
diag(k))
#dif2<-999999
#while (dif2> 0.000001){
#old<-alpha
#lam <- eigen(newmetric)$values
#out <- nlm( iglogl, p=c(1) ,lam=lam, nlam=nlam)
#gamma <- abs(out$estimate[1])
#alpha<- gamma/ mean(1/lam[1:nlam])
#newmetric <- n*m/(n*m+2*(1+gamma))*( ans$Sigmak + (2*alpha/(n*m))*diag(k) )
#dif2<- abs(alpha- old)
#print(dif2)
#}
}
if (prior[1] == "Exponential") {
lam <- eigen(ans$Sigmak)$values
nlam <- min(c(n * m - m - 2, k - 2))
mu <- mean(1 / lam[1:(nlam)])
alpha <- 1 / mu
gamma <- 1
newmetric <-
n * m / (n * m + 2 * (2)) * (ans$Sigmak + (2 * alpha / (n * m)) * diag(k))
#dif2<-999999
#while (dif2> 0.000001){
#old<-alpha
#newmetric <- n*m/(n*m+2*(2))*( ans$Sigmak + (2*alpha/(n*m))*diag(k) )
#lam <- eigen(newmetric)$values
#mu <- mean(1/lam[1:( nlam)])
#alpha <- 1/mu
#newmetric <- n*m/(n*m+2*(2))*( ans$Sigmak + (2*alpha/(n*m))*diag(k) )
#dif2<- abs(alpha- old)
#}
}
if (is.double(prior[1])) {
newmetric <- (ans$Sigmak + prior)
}
if (diagonal == TRUE) {
newmetric <- diag(diag(newmetric))
}
if (fixcovmatrix[1] != FALSE) {
newmetric <- fixcovmatrix
}
ans2 <-
procWGPA1(
X,
ans$mshape,
metric = newmetric ,
scale = scale,
sampleweights = sampleweights
)
plotshapes(ans2$rotated)
dif <- Enorm((ans$Sigmak - ans2$Sigmak))
ans <- ans2
cat(c(it, " ", dif, " \n"))
}
}
if ((priorargument[1] == "Exponential") ||
(priorargument[1] == "Inversegamma")) {
ans$alpha <- alpha
ans$gamma <- gamma
}
cat(" \n")
ans
}
#==================================================================================
procWGPA1 <- function(X,
mu,
metric = "Identity",
scale = TRUE,
reflect = FALSE,
sampleweights = "Equal") {
k <- dim(X)[1]
n <- dim(X)[3]
m <- dim(X)[2]
sum <- 0
for (i in 1:n) {
sum <- sum + centroid.size(X[, , i]) ** 2
}
size1 <- sqrt(sum)
if (sampleweights[1] == "Equal") {
sampleweights <- rep(1 / n, times = n)
}
if (length(sampleweights) != n) {
cat("Sample weight vector not of correct length \n")
}
if (metric[1] == "Identity") {
Sigmak <- diag(k)
}
else{
Sigmak <- metric
}
eig <- eigen(Sigmak, symmetric = TRUE)
Sighalf <-
eig$vectors %*% diag (sqrt(abs(eig$values))) %*% t(eig$vectors)
Siginvhalf <-
eig$vectors %*% diag(1 / sqrt(abs(eig$values))) %*% t(eig$vectors)
Siginv <- eig$vectors %*% diag (1 / (eig$values)) %*% t(eig$vectors)
one <- matrix(rep(1, times = k), k, 1)
Xstar <- X
for (i in 1:n) {
Xstar[, , i] <- Xstar[, , i] - one %*% t(one) %*% Siginv %*% Xstar[, , i] /
c(t(one) %*% Siginv %*% one)
Xstar[, , i] <- Siginvhalf %*% Xstar[, , i]
}
mu <- mu - one %*% t(one) %*% Siginv %*% mu / c(t(one) %*% Siginv %*% one)
ans <- procGPA(Xstar, eigen2d = FALSE)
ans2 <- ans
dif3 <- 99999999
while (dif3 > 0.00001) {
for (i in 1:n) {
old <- mu
tem <-
procOPA(Siginvhalf %*% mu ,
Xstar[, , i],
scale = scale,
reflect = reflect)
Gammai <- tem$R
betai <- tem$s
#ci <- t(one)%*% Siginvhalf %*% X[,,i] %*% Gammai*betai/k
#Yi <- Sighalf%*% ans$rotated[,,i] + Sighalf%*%one%*% ci
#Zi <- Yi - one %*% t(one)%*% Siginv %*% Yi / c( t(one)%*%Siginv%*%one )
Zi <- Sighalf %*% Xstar[, , i] %*% Gammai * betai
ans2$rotated[, , i] <- Zi
}
sum2 <- 0
for (i in 1:n) {
sum2 <- sum2 + centroid.size(ans2$rotated[, , i]) ** 2
}
size2 <- sqrt(sum2)
tem <- ans2$mshape * 0
for (i in 1:n) {
ans2$rotated[, , i] <- ans2$rotated[, , i] * size1 / size2
tem <- tem + ans2$rotated[, , i] * sampleweights[i] / sum(sampleweights)
}
mu <- tem
dif3 <- riemdist(old, mu)
}
z <- ans2
z$mshape <- tem
tan <- z$rotated[, 1,] - z$mshape[, 1]
for (i in 2:m) {
tan <- rbind(tan, z$rotated[, i,] - z$mshape[, i])
}
pca <- prcomp1(t(tan))
z$tan <- tan
npc <- 0
for (i in 1:length(pca$sdev)) {
if (pca$sdev[i] > 1e-07) {
npc <- npc + 1
}
}
z$scores <- pca$x
z$rawscores <- pca$x
z$stdscores <- pca$x
for (i in 1:npc) {
z$stdscores[, i] <- pca$x[, i] / pca$sdev[i]
}
z$pcar <- pca$rotation
z$pcasd <- pca$sdev
z$percent <- z$pcasd ^ 2 / sum(z$pcasd ^ 2) * 100
size <- rep(0, times = n)
rho <- rep(0, times = n)
x <- X
size <- apply(x, 3, centroid.size)
rho <- apply(x, 3, y <- function(x) {
riemdist(x, z$mshape)
})
z$rho <- rho
z$size <- size
z$rmsrho <- sqrt(mean(rho ^ 2))
z$rmsd1 <- sqrt(mean(sin(rho) ^ 2))
z$k <- k
z$m <- m
z$n <- n
tem <- matrix(0, k, k)
for (i in 1:n) {
tem <-
tem + (z$rotated[, , i] - z$mshape) %*% t((z$rotated[, , i] - z$mshape))
}
tem <- tem / (n * m)
z$Sigmak <- tem
return(z)
}
#==================================================================================
testshapes <-
function(A,
B,
resamples = 1000,
replace = TRUE,
scale = TRUE) {
if (replace == TRUE) {
out <- bootstraptest(A, B, resamples = resamples, scale = scale)
}
if (replace == FALSE) {
out <- permutationtest(A, B, nperms = resamples, scale = scale)
}
out
}
#==================================================================================
testmeanshapes <-
function(A,
B,
resamples = 1000,
replace = FALSE,
scale = TRUE) {
if (replace == TRUE) {
out <- bootstraptest(A, B, resamples = resamples, scale = scale)
}
if (replace == FALSE) {
out <- permutationtest(A, B, nperms = resamples, scale = scale)
}
if (resamples > 0) {
aa <- list(
H = 0,
H.pvalue = 0,
H.table.pvalue = 0,
G = 0,
G.pvalue = 0,
G.table.pvalue = 0,
J = 0,
J.pvalue = 0,
J.table.pvalue = 0
)
aa$H <- out$H
aa$H.pvalue <- out$H.pvalue
aa$H.table.pvalue <- out$H.table.pvalue
aa$G <- out$G
aa$G.pvalue <- out$G.pvalue
aa$G.table.pvalue <- out$G.table.pvalue
aa$J <- out$J
aa$J.pvalue <- out$J.pvalue
aa$J.table.pvalue <- out$J.table.pvalue
}
if (resamples == 0) {
aa <- list(
H = 0,
H.table.pvalue = 0,
G = 0,
G.table.pvalue = 0,
J = 0,
J.table.pvalue = 0
)
aa$H <- out$H
aa$H.table.pvalue <- out$H.table.pvalue
aa$G <- out$G
aa$G.table.pvalue <- out$G.table.pvalue
aa$J <- out$J
aa$J.table.pvalue <- out$J.table.pvalue
}
aa
}
#==================================================================================
permutationtest2 <- function (A, B, nperms = 1000, scale = scale)
{
A1 <- A
A2 <- B
mdim <- dim(A1)[2]
B <- nperms
nsam1 <- dim(A1)[3]
nsam2 <- dim(A2)[3]
pool <-
procGPA(
abind (A1, A2) ,
scale = scale,
tangentcoords = "partial",
pcaoutput = FALSE
)
tempool <- pool
for (i in 1:(nsam1 + nsam2)) {
tempool$tan[, i] <- pool$tan[, i] / Enorm(pool$tan[, i]) * pool$rho[i]
}
pool <- tempool
permpool <- pool
Gtem <- Goodall(pool, nsam1, nsam2)
Htem <- Hotelling(pool, nsam1, nsam2)
Jtem <- James(pool, nsam1, nsam2, table = TRUE)
Ltem <- Lambdamin(pool, nsam1, nsam2)
Gumc <- Gtem$F
Humc <- Htem$F
Jumc <- Jtem$Tsq
Lumc <- Ltem$lambdamin
Gtabpval <- Gtem$pval
Htabpval <- Htem$pval
Jtabpval <- Jtem$pval
Ltabpval <- Ltem$pval
if (B > 0) {
Apool <- array(0, c(dim(A1)[1], dim(A1)[2], dim(A1)[3] +
dim(A2)[3]))
Apool[, , 1:nsam1] <- A1
Apool[, , (nsam1 + 1):(nsam1 + nsam2)] <- A2
out <-
list(
H = 0,
H.pvalue = 0,
H.table.pvalue = 0,
J = 0,
J.pvalue = 0,
J.table.pvalue = 0,
G = 0,
G.pvalue = 0,
G.table.pvalue = 0
)
Gu <- rep(0, times = B)
Hu <- rep(0, times = B)
Ju <- rep(0, times = B)
Lu <- rep(0, times = B)
cat("Permutations - sampling without replacement: ")
cat(c("No of permutations = ", B, "\n"))
for (i in 1:B) {
if (i / 100 == trunc(i / 100)) {
cat(c(i, " "))
}
select <- sample(1:(nsam1 + nsam2))
permpool$tan <- pool$tan[, select]
Gu[i] <- Goodall(permpool, nsam1, nsam2)$F
Hu[i] <- Hotelling(permpool, nsam1, nsam2)$F
Ju[i] <- James(permpool, nsam1, nsam2)$Tsq
Lu[i] <-
Lambdamin(permpool, nsam1, nsam2)$lambdamin
}
Gu <- sort(Gu)
numbig <- length(Gu[Gumc < Gu])
pvalG <- (1 + numbig) / (B + 1)
Lu <- sort(Lu)
numbig <- length(Lu[Lumc < Lu])
pvalL <- (1 + numbig) / (B + 1)
Ju <- sort(Ju)
numbig <- length(Ju[Jumc < Ju])
pvalJ <- (1 + numbig) / (B + 1)
Hu <- sort(Hu)
numbig <- length(Hu[Humc < Hu])
pvalH <- (1 + numbig) / (B + 1)
cat(" \n")
out$Hu <- Hu
out$Ju <- Ju
out$Gu <- Gu
out$Lu <- Lu
out$H <- Humc
out$H.pvalue <- pvalH
out$H.table.pvalue <- Htabpval
out$J <- Jumc
out$J.pvalue <- pvalJ
out$J.table.pvalue <- Jtabpval
out$G <- Gumc
out$G.pvalue <- pvalG
out$G.table.pvalue <- Gtabpval
out$L <- Lumc
out$L.pvalue <- pvalL
out$L.table.pvalue <- Ltabpval
}
if (B == 0) {
out <- list(
H = 0,
H.table.pvalue = 0,
G = 0,
G.table.pvalue = 0
)
out$H <- Humc
out$H.table.pvalue <- Htabpval
out$J <- Jumc
out$J.table.pvalue <- Jtabpval
out$G <- Gumc
out$G.table.pvalue <- Gtabpval
out$L <- Lumc
out$L.table.pvalue <- Ltabpval
}
out
}
#==================================================================================
bootstraptest <- function (A,
B,
resamples = 200,
scale = TRUE)
{
A1 <- A
A2 <- B
mdim <- dim(A1)[2]
B <- resamples
nsam1 <- dim(A1)[3]
nsam2 <- dim(A2)[3]
pool <-
procGPA(
abind (A1, A2) ,
scale = scale ,
tangentcoords = "partial",
pcaoutput = FALSE
)
tempool <- pool
for (i in 1:(nsam1 + nsam2)) {
tempool$tan[, i] <- pool$tan[, i] / Enorm(pool$tan[, i]) * pool$rho[i]
}
pool <- tempool
bootpool <- pool
Gtem <- Goodall(pool, nsam1, nsam2)
Htem <- Hotelling(pool, nsam1, nsam2)
Jtem <- James(pool, nsam1, nsam2, table = TRUE)
Ltem <- Lambdamin(pool, nsam1, nsam2)
Gumc <- Gtem$F
Humc <- Htem$F
Jumc <- Jtem$Tsq
Lumc <- Ltem$lambdamin
Gtabpval <- Gtem$pval
Htabpval <- Htem$pval
Jtabpval <- Jtem$pval
Ltabpval <- Ltem$pval
if (B > 0) {
Apool <- array(0, c(dim(A1)[1], dim(A1)[2], dim(A1)[3] +
dim(A2)[3]))
Apool[, , 1:nsam1] <- A1
Apool[, , (nsam1 + 1):(nsam1 + nsam2)] <- A2
out <-
list(
H = 0,
H.pvalue = 0,
H.table.pvalue = 0,
J = 0,
J.pvalue = 0,
J.table.pvalue = 0,
G = 0,
G.pvalue = 0,
G.table.pvalue = 0
)
Gu <- rep(0, times = B)
Hu <- rep(0, times = B)
Ju <- rep(0, times = B)
Lu <- rep(0, times = B)
pool2 <- pool
pool2$tan[, 1:nsam1] <-
pool$tan[, 1:nsam1] - apply(pool$tan[, 1:nsam1], 1, mean)
pool2$tan[, (nsam1 + 1):(nsam1 + nsam2)] <-
pool$tan[, (nsam1 + 1):(nsam1 + nsam2)] -
apply(pool$tan[, (nsam1 + 1):(nsam1 + nsam2)], 1, mean)
cat("Bootstrap - sampling with replacement within each group under H0: ")
cat(c("No of resamples = ", B, "\n"))
for (i in 1:B) {
if (i / 100 == trunc(i / 100)) {
cat(c(i, " "))
}
select1 <- sample(1:nsam1, replace = TRUE)
select2 <- sample((nsam1 + 1):(nsam1 + nsam2), replace = TRUE)
bootpool$tan <- pool2$tan[, c(select1, select2)]
Gu[i] <- Goodall(bootpool, nsam1, nsam2)$F
Hu[i] <- Hotelling(bootpool, nsam1, nsam2)$F
Ju[i] <- James(bootpool, nsam1, nsam2)$Tsq
Lu[i] <- Lambdamin(bootpool, nsam1, nsam2)$lambdamin
}
Gu <- sort(Gu)
numbig <- length(Gu[Gumc < Gu])
pvalG <- (1 + numbig) / (B + 1)
Ju <- sort(Ju)
numbig <- length(Ju[Jumc < Ju])
pvalJ <- (1 + numbig) / (B + 1)
Hu <- sort(Hu)
numbig <- length(Hu[Humc < Hu])
pvalH <- (1 + numbig) / (B + 1)
numbig <- length(Lu[Lumc < Lu])
pvalL <- (1 + numbig) / (B + 1)
cat(" \n")
out$Hu <- Hu
out$Ju <- Ju
out$Gu <- Gu
out$Lu <- Lu
out$H <- Humc
out$H.pvalue <- pvalH
out$H.table.pvalue <- Htabpval
out$J <- Jumc
out$J.pvalue <- pvalJ
out$J.table.pvalue <- Jtabpval
out$G <- Gumc
out$G.pvalue <- pvalG
out$G.table.pvalue <- Gtabpval
out$L <- Lumc
out$L.pvalue <- pvalL
out$L.table.pvalue <- Ltabpval
}
if (B == 0) {
out <-
list(
H = 0,
H.table.pvalue = 0,
G = 0,
G.table.pvalue = 0,
J = 0,
J.table.pvalue = 0
)
out$H <- Humc
out$H.table.pvalue <- Htabpval
out$J <- Jumc
out$J.table.pvalue <- Jtabpval
out$G <- Gumc
out$G.table.pvalue <- Gtabpval
out$L <- Lumc
out$L.table.pvalue <- Ltabpval
}
out
}
#==================================================================================
Lambdamin <-
function (pool, n1, n2, p = 0)
{
censiz <- centroid.size(pool$mshape)
tan1 <- pool$tan[, 1:n1]
tan2 <- pool$tan[, (n1 + 1):(n1 + n2)]
kt <- dim(tan1)[1]
n <- n1 + n2
k <- pool$k
m <- pool$m
if (p == 0) {
p <- min(k * m - (m * (m - 1)) / 2 - 1 - m, n1 + n2 - 2)
}
HH <- diag(k)
mu1 <- pool$mshape
if (dim(tan1)[1] == k * m - m) {
HH <- defh(k - 1)
mu1 <- preshape(pool$mshape)
}
if (m == 2) {
mu <- c(mu1[, 1], mu1[, 2])
}
if (m == 3) {
mu <- c(mu1[, 1], mu1[, 2], mu1[,
3])
}
dd <- kt
X1 <- tan1 * 0
X2 <- tan2 * 0
S1 <- matrix(0, dd, dd)
S2 <- matrix(0, dd, dd)
for (i in 1:n1) {
X1[, i] <- (mu + tan1[, i]) / Enorm(mu + tan1[, i])
S1 <- S1 + X1[, i] %*% t(X1[, i])
}
for (i in 1:n2) {
X2[, i] <- (mu + tan2[, i]) / Enorm(mu + tan2[, i])
S2 <- S2 + X2[, i] %*% t(X2[, i])
}
sumx1 <- 0
sumx2 <- 0
for (i in 1:n1) {
sumx1 <- sumx1 + X1[, i]
}
for (i in 1:n2) {
sumx2 <- sumx2 + X2[, i]
}
sum1 <- apply(X1, 1, sum)
sum2 <- apply(X2, 1, sum)
mean1 <- sum1 / Enorm(sum1)
mean2 <- sum2 / Enorm(sum2)
bb1 <- mean1[1:(dd - 1)]
cc1 <- mean1[dd]
bb2 <- mean2[1:(dd - 1)]
cc2 <- mean2[dd]
A1 <- cc1 / abs(cc1) * diag(dd - 1) - cc1 / (abs(cc1) + cc1 ^ 2) *
bb1 %*% t(bb1)
M1 <- cbind(A1,-bb1)
A1 <- cc2 / abs(cc2) * diag(dd - 1) - cc2 / (abs(cc2) + cc2 ^ 2) *
bb2 %*% t(bb2)
M2 <- cbind(A1,-bb2)
G1 <- matrix(0, dd - 1, dd - 1)
G2 <- matrix(0, dd - 1, dd - 1)
for (iu in 1:(dd - 1)) {
for (iv in iu:(dd - 1)) {
G1[iu, iv] <- G1[iu, iv] + t((t(M1))[, iu]) %*% S1 %*%
(t(M1))[, iv]
G1[iv, iu] <- G1[iu, iv]
G2[iu, iv] <- G2[iu, iv] + t((t(M2))[, iu]) %*% S2 %*%
(t(M2))[, iv]
G2[iv, iu] <- G2[iu, iv]
}
}
G1 <- G1 / n1 / Enorm(sumx1 / n1) ^ 2
G2 <- G2 / n2 / Enorm(sumx2 / n2) ^ 2
# eva1 <- eigen(G1, symmetric = TRUE, EISPACK = TRUE)
eva1 <- eigen(G1, symmetric = TRUE)
pcar1 <- eva1$vectors[, 1:p]
pcasd1 <- sqrt(abs(eva1$values[1:p]))
# eva2 <- eigen(G2, symmetric = TRUE, EISPACK = TRUE)
eva2 <- eigen(G2, symmetric = TRUE)
pcar2 <- eva2$vectors[, 1:p]
pcasd2 <- sqrt(abs(eva2$values[1:p]))
if ((pcasd1[p] < 1e-06) || (pcasd2[p] < 1e-06)) {
offset <- 1e-06
cat("*")
pcasd1 <- sqrt(pcasd1 ^ 2 + offset)
pcasd2 <- sqrt(pcasd2 ^ 2 + offset)
}
Ahat1 <-
n1 * t(M1) %*% (pcar1 %*% diag(1 / pcasd1 ^ 2) %*% t(pcar1)) %*%
M1
Ahat2 <-
n2 * t(M2) %*% (pcar2 %*% diag(1 / pcasd2 ^ 2) %*% t(pcar2)) %*%
M2
Ahat <- (Ahat1 + Ahat2)
# eva <- eigen(Ahat, symmetric = TRUE, EISPACK = TRUE)
eva <- eigen(Ahat, symmetric = TRUE)
lambdamin <- eva$values[p + 1]
pval <- 1 - pchisq(lambdamin, p)
#print(lambdamin)
#print(pval)
z <- list()
z$pval <- pval
z$df <- p
z$lambdamin <- lambdamin
return(z)
}
#==================================================================================
Goodall <- function(pool , n1, n2, p = 0) {
tan1 <- pool$tan[, 1:n1]
tan2 <- pool$tan[, (n1 + 1):(n1 + n2)]
kt <- dim(tan1)[1]
n <- n1 + n2
k <- pool$k
m <- pool$m
if (p == 0) {
p <- min(k * m - (m * (m - 1)) / 2 - 1 - m, n1 + n2 - 2)
}
top <- Enorm(apply(tan1, 1, mean) - apply(tan2, 1, mean)) ** 2
bot <-
sum(diag(var(t(tan1)))) * (n1 - 1) + sum(diag(var(t(tan2)))) * (n2 - 1)
Fstat <- ((n1 + n2 - 2) / (1 / n1 + 1 / n2) * top) / bot
pval <- 1 - pf(Fstat, p, (n1 + n2 - 2) * p)
z <- list()
z$F <- Fstat
z$pval <- pval
z$df1 <- p
z$df2 <- (n1 + n2 - 2) * p