https://github.com/cran/robCompositions
Tip revision: b0cf7e2c31b281933351acf2fd497c40e4e77996 authored by Matthias Templ on 08 February 2016, 15:39:39 UTC
version 2.0.0
version 2.0.0
Tip revision: b0cf7e2
factanal.fit.principal1.R
# PF, 2008-09-18
# is used in pfa1.R which
# computes principal factor analysis for compositional data
# Uniquenesses are nor longer of diagonal form
factanal.fit.principal1 <-
function (cmat, factors, p = ncol(cmat), start = NULL, iter.max = 10,
unique.tol = 1e-04)
{
dof <- 0.5 * ((p - factors)^2 - p - factors)
if (dof < 0)
warning("negative degrees of freedom")
if (any(abs(diag(cmat) - 1) > .Machine$single.eps))
stop("must have correlation matrix")
if (length(start)) {
if (length(start) != p)
stop("start is the wrong length")
if (any(start < 0 | start >= 1))
stop("all values in start must be between 0 and 1")
oldcomm <- 1 - start
}
else {
diag(cmat) <- NA
oldcomm <- apply(abs(cmat), 1, max, na.rm = TRUE)
}
# PF 10.9.2008
H <- diag(p)-matrix(1,p,p)/p
psi <- 1-oldcomm
psistar <- H%*%diag(psi)%*%H
cmatstar <- cmat-psistar
if (iter.max < 0)
stop("bad value for iter.max")
ones <- rep(1, factors)
if (iter.max == 0) {
z <- eigen(cmatstar, symmetric = TRUE)
kvals <- z$values[1:factors]
if (any(kvals <= 0))
stop("impermissible estimate reached")
Lambda <- z$vectors[, 1:factors, drop = FALSE] * rep(kvals^0.5,
rep.int(p, factors))
psinew <- diag(cmat) - Lambda^2 %*% ones
}
if (iter.max > 0) {
for (i in 1:iter.max) {
z <- eigen(cmatstar, symmetric = TRUE)
kvals <- z$values[1:factors]
if (any(kvals <= 0))
stop("impermissible estimate reached")
Lambda <- z$vectors[, 1:factors, drop = FALSE] *
rep(kvals^0.5, rep.int(p, factors))
psinew <- drop(diag(cmat) - Lambda^2 %*% ones)
psinewstar <- H%*%diag(psinew)%*%H
if (all(abs(psinew - psi) < unique.tol)) {
iter.max <- i
break
}
psistar <- psinewstar
cmatstar <- cmat-psistar
}
}
dn <- dimnames(cmat)[[1]]
dimnames(Lambda) <- list(dn, paste("Factor", 1:factors, sep = ""))
diag(cmat) <- 1
uniq <- diag(psistar)
names(uniq) <- dn
ans <- list(loadings = Lambda, uniquenesses = uniq, correlation = cmat,
criteria = c(iterations = iter.max), factors = factors,
dof = dof, method = "principal",psi = psistar)
class(ans) <- "factanal"
ans
}