Skip to content

Instantly share code, notes, and snippets.

@kohske
Created May 9, 2018 09:29
Show Gist options
  • Save kohske/2da668635171993e768d56318124ada7 to your computer and use it in GitHub Desktop.
Save kohske/2da668635171993e768d56318124ada7 to your computer and use it in GitHub Desktop.
score付きpromax2
# zmatにはfactanalに渡したデータをもう一度渡す(因子得点計算用)
promax2 <- function (x, power = 4, kaiser = TRUE, zmat = NULL)
{
if (!is.matrix(x) & !is.data.frame(x)) {
if (!is.null(x$loadings))
x <- as.matrix(x$loadings)
}
else {
x <- x
}
if (ncol(x) < 2)
return(x)
dn <- dimnames(x)
xx <- varimax(x, normalize = kaiser)
temp <- list()
for(i in 1:ncol(xx$loading)){
temp[[i]] <- apply(abs(xx$loading^2),1,sum)
}
temp <- t(matrix(unlist(temp),nrow=ncol(xx$loading),byrow=nrow(xx$loading)))
xxx <- (xx$loadings/temp^0.5)
Q <- xxx * abs(xxx)^(power - 1)
U <- lm.fit(x, Q)$coefficients
d <- try(diag(solve(t(U) %*% U)), silent = TRUE)
if (class(d) == "try-error") {
warning("Factors are exactly uncorrelated and the model produces a singular matrix. An approximation is used")
ev <- eigen(t(U) %*% U)
ev$values[ev$values < .Machine$double.eps] <- 100 * .Machine$double.eps
UU <- ev$vectors %*% diag(ev$values) %*% t(ev$vectors)
diag(UU) <- 1
d <- diag(solve(UU))
}
U <- U %*% diag(sqrt(d))
dimnames(U) <- NULL
z <- x %*% U
U <- xx$rotmat%*%U
ui <- solve(U)
Phi <- ui %*% t(ui)
dimnames(z) <- dn
class(z) <- "loadings"
result <- list(loadings = z, rotmat = U, Phi = Phi)
# score
# 例外処理してないので要注意・・・
if (!is.null(zmat)) {
zmat = as.matrix(zmat)
scores = "Bartlett"
Lambda <- result$loadings
zz <- scale(zmat, TRUE, TRUE)
switch(scores, regression = {
covmat <- cov.wt(zmat)
cv <- covmat$cov
sc <- zz %*% solve(cv, Lambda)
if (!is.null(Phi <- attr(Lambda, "covariance"))) sc <- sc %*%
Phi
}, Bartlett = {
# uniqunessを計算
uniqueness = fa2un(result)
Lambda <- result$loadings
d <- 1/(uniqueness[[2]])
tmp <- t(Lambda * d)
sc <- t(solve(tmp %*% Lambda, tmp %*% t(zz)))
})
rownames(sc) <- rownames(zmat)
colnames(sc) <- colnames(Lambda)
# 欠損値処理は省略
# if (!is.null(na.act))
# sc <- napredict(na.act, sc)
result$scores <- sc
}
class(result) <- c("psych", "fa")
return(result)
}
# print.psych.faよりuniquness計算
fa2un = function(x, digits = 2, all = FALSE, cut = NULL, sort = FALSE,
suppress.warnings = TRUE, ...)
{
#if (!is.matrix(x) && !is.null(x$fa) && is.list(x$fa))
# x <- x$fa
load <- x$loadings
cut <- 0
nitems <- dim(load)[1]
nfactors <- dim(load)[2]
if (sum(x$uniqueness) + sum(x$communality) > nitems) {
covar <- TRUE
}
else {
covar <- FALSE
}
loads <- data.frame(item = seq(1:nitems), cluster = rep(0,
nitems), unclass(load))
u2.order <- 1:nitems
if (sort) {
loads$cluster <- apply(abs(load), 1, which.max)
ord <- sort(loads$cluster, index.return = TRUE)
loads[1:nitems, ] <- loads[ord$ix, ]
rownames(loads)[1:nitems] <- rownames(loads)[ord$ix]
items <- table(loads$cluster)
first <- 1
item <- loads$item
for (i in 1:length(items)) {
if (items[i] > 0) {
last <- first + items[i] - 1
ord <- sort(abs(loads[first:last, i + 2]), decreasing = TRUE,
index.return = TRUE)
u2.order[first:last] <- item[ord$ix + first -
1]
loads[first:last, 3:(nfactors + 2)] <- load[item[ord$ix +
first - 1], ]
loads[first:last, 1] <- item[ord$ix + first -
1]
rownames(loads)[first:last] <- rownames(loads)[ord$ix +
first - 1]
first <- first + items[i]
}
}
}
if (max(abs(load) > 1) && !covar)
cat("\n Warning: A Heywood case was detected. \n")
ncol <- dim(loads)[2] - 2
rloads <- round(loads, digits)
fx <- format(rloads, digits = digits)
nc <- nchar(fx[1, 3], type = "c")
fx.1 <- fx[, 1, drop = FALSE]
fx.2 <- fx[, 3:(2 + ncol), drop = FALSE]
load.2 <- as.matrix(loads[, 3:(ncol + 2)])
fx.2[abs(load.2) < cut] <- paste(rep(" ", nc), collapse = "")
if (sort) {
fx <- data.frame(V = fx.1, fx.2)
if (dim(fx)[2] < 3)
colnames(fx) <- c("V", colnames(x$loadings))
}
else {
fx <- data.frame(fx.2)
colnames(fx) <- colnames(x$loadings)
}
if (nfactors > 1) {
if (is.null(x$Phi)) {
h2 <- rowSums(load.2^2)
}
else {
h2 <- diag(load.2 %*% x$Phi %*% t(load.2))
}
}
else {
h2 <- load.2^2
}
if (!is.null(x$uniquenesses)) {
u2 <- x$uniquenesses[u2.order]
}
else {
u2 <- (1 - h2)
}
list(h2, u2)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment