score付きpromax2

promax2.r

# 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)
}