myscratchpad

psychパッケージのfa関数が不思議な挙動を見せるようなので
自分で調べてRjpwikiにレスしておいた。
せっかく調べたのでもう少し詳しく書いてみる

 library(psych)
 data(bfi)
 x <- factor.pa(bfi, nfactors=3, rotate="promax")

 # psych:::print.psych.faの中身を途中まで実行
 # パラメータを設定しておく
 digits = 2; all = FALSE; cut = NULL; sort = TRUE

#> psych:::print.psych.fa
#function (x, digits = 2, all = FALSE, cut = NULL, sort = FALSE,
#    ...) 
#{ ここまでは定義なので実行しないこと

    if (!is.null(x$fn)) {
        if (x$fn == "principal") {
            cat("Principal Components Analysis")
        }
        else {
            cat("Factor Analysis using method = ", x$fm)
        }
    }
    cat("\nCall: ")
    print(x$Call)
    load <- x$loadings
    if (is.null(cut))
        cut <- 0
    nitems <- dim(load)[1]
    nfactors <- dim(load)[2]
    loads <- data.frame(item = seq(1:nitems), cluster = rep(0,
        nitems), unclass(load))
    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)
                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 (cut > 0) {
        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)]
        load.2 <- as.matrix(loads[, 3:(ncol + 2)])
        fx.2[abs(load.2) < cut] <- paste(rep(" ", nc), collapse = "")
        fx <- data.frame(V = fx.1, fx.2)
        if (dim(fx)[2] < 3)
            colnames(fx) <- c("V", 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
        }
        else {
            u2 <- (1 - h2)
        }
        vtotal <- sum(h2 + u2)
        if (isTRUE(all.equal(vtotal, nitems))) {
            cat("Standardized loadings based upon correlation matrix\n")
            print(cbind(fx, h2, u2), quote = "FALSE", digits = digits)
        }
        else {
            cat("Unstandardized loadings based upon covariance matrix\n")
            print(cbind(fx, h2, u2, H2 = h2/(h2 + u2), U2 = u2/(h2 +
                u2)), quote = "FALSE", digits = digits)
        }
    } else {
        ncol <- dim(loads)[2] - 2
        load.2 <- as.matrix(loads[, 3:(ncol + 2)])
        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
        }
        else {
            u2 <- (1 - h2)
        }
        vtotal <- sum(h2 + u2)
        if (isTRUE(all.equal(vtotal, nitems))) {
            cat("Standardized loadings based upon correlation matrix\n")
            print(cbind(round(loads[-(1:2)], digits), h2, u2),
                quote = "FALSE", digits = digits)
        }
        else {
            cat("Unstandardized loadings based upon covariance matrix\n")
            print(cbind(round(loads[-c(1:2)], digits), h2, u2,
                H2 = h2/(h2 + u2), U2 = u2/(h2 + u2)), quote = "FALSE",
                digits = digits)
        }
    }
    if (is.null(x$Phi)) {
        if (nfactors > 1) {
            vx <- colSums(load.2^2)
        }
        else {
            vx <- sum(load.2^2)
        }
    } else {
        vx <- diag(x$Phi %*% t(load) %*% load)
    }
    names(vx) <- colnames(x$loadings)
    varex <- rbind(`SS loadings` = vx)
    varex <- rbind(varex, `Proportion Var` = vx/vtotal)
    if (nfactors > 1)
        varex <- rbind(varex, `Cumulative Var` = cumsum(vx/vtotal))
    cat("\n")
    print(round(varex, digits))

## この時点でh2はソートされているが、u2はソートされていない
h2
u2

# この後のコードでこの両者を使って標準化している (fx.2 ... sqrt(h2 + u2)) ため、誤りが生じます
    #if (!isTRUE(all.equal(vtotal, nitems))) {
        cat("\n Standardized loadings\n")
        fx <- format(loads, digits = digits)
        nc <- nchar(fx[1, 3], type = "c")
        fx.1 <- fx[, 1, drop = FALSE]
        fx.2 <- round(loads[, 3:(2 + ncol)]/sqrt(h2 + u2), digits)
         ########
         # 誤った負荷量
         fx.2
         # 修正コード
         u2 <- 1 - h2
         fx.2 <- round(loads[, 3:(2 + ncol)]/sqrt(h2 + u2), digits)
         fx.2
         ########
        load.2 <- loads[, 3:(ncol + 2)]/sqrt(h2 + u2)
        fx.2[abs(load.2) < cut] <- paste(rep(" ", nc), collapse = "")
        fx <- data.frame(V = fx.1, fx.2)

## 以下本題
一エンドユーザーとしてみると、RがSPSSやSASより優れている点はフリーであること以外にない。Rを利用しているえらい人は (自分の本を売るためか) Rの使用をすすめる人が多いが、個人的には止めたほうがいいと思う。
それは使用に慣れるまでが大変で、容易に間違った結果を出力するからである。巷にあふれるRの入門書も誤った使用法を紹介しているものが山ほどある (厳密にいえば誤っているわけではなく、教科書の例題のようなごく限られた場合にのみ適用可能なコードを紹介している) …ような気がする。

特に、「世界中の人が使ってチェックするので，計算の信頼性もかなり高い。」といった宣伝文句はどうかと思う。baseに含まれている関数には当てはまるが、あまたあるパッケージの関数には当てはまらない。というか、下手に使うと間違った結果しか出力されず、信頼性はかなり低い。なぜなら、上記のようにフリーであるので誰にでも使用できるが、そのため誤ったコード例が量産されている (本サイトもかなり誤りがあるんです) 。エンドユーザーはコピペのように手軽に実行でき、しかも正しい結果の得られるものを求める。Rで自らのデータに合わせて適切にオプションを選択し、分析を実行していくのは容易ではない。

コピペで間違ったら自己責任でしょ、フリーなんだし、自分の研究なんだから…というのであれば、やはりSPSSやSASの方が優れているとしか思えない。
「Rで分析しました」という記述があったら、その結果を疑わしく見てしまうのは自分だけかな。「下手に使って、どっか間違ってんじゃないの？」と思ってしまう。