No.14115　因子得点　　【SK】　2011/01/06(Thu) 01:24

因子分析の勉強のため，http://aoki2.si.gunma-u.ac.jp/R/pfa.html　を読んでいます。
因子得点を算出する部分ですが，引用しますと，
--
w <- solve(r0)%*%evec # 因子得点係数を計算する
scores <- (scale(dat)*sqrt(nr/(nr-1)))%*%w # 因子得点を計算する
--
というコードが書かれています。
この部分の sqrt(nr/nr-1)の部分は必要なのでしょうか？Rのヘルプを見ると,scale(dat)のみで，datの標準化が可能のようですが，この理解は間違いでしょうか？なにとぞよろしくお願いいたします。
No.14116　Re: 因子得点　　【青木繁伸】　2011/01/06(Thu) 01:36

不偏分散を1にするか，分散を1にするかのこだわりで，どうでもよいと思う人にはどうでもよいことでしょう。
No.14117　Re: 因子得点　　【sk】　2011/01/06(Thu) 02:10

ありがとうございます。寝てからもう一度，式を確認してみます。失礼しました．
No.14202　Re: 因子得点　　【新参者】　2011/01/23(Sun) 21:53

統計もRも初心者です。
近赤外スペクトルデータを主成分分析した結果からマハラノビス平方距離を計算して判別分析をしてみようと考えていて，*(n-1)/n で迷っていたところ，
この質問に出会いましたので，返信の形で質問させていただきます。

*(n-1)/n を使うか使わないかは，サンプル分散をとるか不偏分散をとるかということで，元のデータを集団全体とみなすか，母集団からサンプリングされたものとみなす かによって変わるのではないかと考えているのですが，青木先生の関数では，サンプル分散を使用されていて，Rの関数では不偏分散が使用されているように思 います。そもそも，主成分分析やマハラノビス平方距離の計算に不偏分散を使うのは不適切で，厳密には，サンプル統計量を使用すべきなのでしょうか？
少し齧っただけ口幅ったいことを申しますが，Rの関数の中には，この点に無頓着なものもあるような気がしてきましたので，質問させていただきました。
無頓着と感じました理由は，
例えば，パッケージ{timeDate}にある関数 skewnessとkurtosis でmethod="moment"を指定したときの計算が，
skewness = sum((x - mean(x))^3/sqrt(as.numeric(var(x)))^3)/length(x)
kurtosis = sum((x - mean(x))^4/as.numeric(var(x))^2)/length(x)
となっていたからです。これは，サンプルskewnessなので，不偏分散 var(x)を使うのは不適切ではないかと考えています。実際，青木先生の関数や dagoTestが計算する skewnessとkurtosisとは異なる結果を返します。

話はそれるのですが，method ="fisher"を指定した場合に，population skewnessとpopulation kurtosisを計算するのだとすると，式
skewness = ((sqrt(n * (n - 1))/(n - 2)) * (sum(x^3)/n))/((sum(x^2)/n)^(3/2))
kurtosis = ((n + 1) * (n - 1) * ((sum(x^4)/n)/(sum(x^2)/n)^2 - (3 * (n - 1))/(n + 1)))/((n - 2) * (n - 3))
は，分子，分母とも，- mean(x)が抜けているのではないでしょうか？

この点についてもご教示いただけましたら幸いです。
No.14204　Re: 因子得点　　【青木繁伸】　2011/01/24(Mon) 00:22

変動をサンプルサイズで割ったものをサンプル分散と呼ぶのは，ほかのソフトなどが「標本分散」と呼んでいる「標 本」を「サンプル」としたのでしょうか。いずれも不適切だと思いますけど？そもそも，標本データで計算したものが標本分散ですから，「サンプルサイズ」で 割ろうが「サンプルサイズ−1」で割ろうが標本分散でしょう。その二つの標本分散のうち，母分散を推定するときには「サンプルサイズ−1」で割った不偏分 散が適切ですよということですよね。また，「サンプルサイズ」で割ると「母分散」を計算しているというのは，「得られたデータを母集団と見なすから母分 散」などというのは本末転倒な定義だと思いますけどねえ。有限母集団，しかも小規模な有限母集団でなければ，母集団全部のデータを得るなどということは普 通できないことで，「サンプルサイズ」で割ったら「母分散」というのもおかしな話。このあたりは，母分散，不偏分散，標本分散が混乱して用いられていると 思いますね。

整理すると，母分散を推定する，あるいは，推定された母分散を用いてさらなる計算を行うという目的があるならば不偏分散を使 うべきでしょう（実際，検定では不偏分散が使われる）。また，得られたデータについてのみ何かをいう場合（例えば，データを正規化する）などの場合は分散 （不偏でない方）を使う。因子得点などを正規化する場合も同様。

> 分子，分母とも，- mean(x)が抜けているのではないでしょうか？

そ の定義式を確認したわけではないですが，その計算式では x は平均値が引かれた平均偏差データにされているのでは？正規化されていれば mean(x) = 0 は不要ですよね。それとも，fisher を指定すると誤った答えが出るのですか？もし，誤った答えが出るのなら，パッケージのメインテナーに連絡してあげた方がよいですよ。
No.14209　Re: 因子得点　　【青木繁伸】　2011/01/24(Mon) 11:59

http://core.ecu.edu/psyc/wuenschk/docs30/Skew-Kurt.doc
R のマニュアルは，注意を喚起していないなあ。間違えて使っている人が多いのかも。もっとも，method="fisher" を選択するような人は分かっているか。

No.14211　Re: 因子得点　　【青木繁伸】　2011/01/24(Mon) 15:46

平行移動したデータの skewness, kurtosis は同じであるはず。以下のようにやると，method="fisher" を指定して元データをそのまま与えると，異なった答えになってしまうことがわかる。
結論:　method="fisher" を指定する場合は，元データをそのまま与えるのではなく，平均偏差データか正規化（標準化）データを与えなければならない。
> library(timeDate)
> set.seed(123456)
> x1 <- rnorm(20)
> x2 <- x1+10

# method="moment" の場合
> skewness(x1)
[1] 0.1576622
attr(,"method")
[1] "moment"
> skewness(x2)
[1] 0.1576622
attr(,"method")
[1] "moment"

method="fisher" の場合

元データ　　===>  いずれも間違い
> skewness(x1, method="fisher")
[1] 1.420123
attr(,"method")
[1] "fisher"
> skewness(x2, method="fisher")
[1] 1.097535
attr(,"method")
[1] "fisher"

平均偏差データ
> skewness(x1-mean(x1), method="fisher")
[1] 0.1844003
attr(,"method")
[1] "fisher"
> skewness(x2-mean(x2), method="fisher")
[1] 0.1844003
attr(,"method")
[1] "fisher"

標準化データ（正規化データ）
> skewness((x1-mean(x1))/sd(x1), method="fisher")
[1] 0.1844003
attr(,"method")
[1] "fisher"
> skewness((x2-mean(x2))/sd(x2), method="fisher")
[1] 0.1844003
attr(,"method")
[1] "fisher"
No.14213　Re: 因子得点　　【新参者】　2011/01/24(Mon) 21:36

青木繁伸先生

日曜日の夜遅い時刻の投稿であったにもかかわりませず，早々にご回答ありがとうございました。
言葉の使い方が不適切で失礼いたしました。
Wikipedia の skewness の説明　http://en.wikipedia.org/wiki/Skewness　にありました sample skewness, sample variance という言葉から，
ちゃんと確かめもせず捏造してしまいました。

不偏分散と分散（不偏でない方）の使い分けについて，先生の解説を読ませていただき，少し理解が進んだように思います。
あるデータセットを学習用データとして主成分分析した結果を用いて，新しく測定したデータの判別分析に使用すること考えておりましたので，
この場合は不偏分散を使うべしと理解いたしました。未だ，勘違いしておりますようでしたら，ご指摘をお願いします。

また，skewness と kurtosisについても詳細に検討くださいましたこと，厚くお礼申し上げます。
Rの help もご確認下さったかと存じますが，method の説明は， These are either "moment" or "fisher"
The "moment" method is based on the definitions of skewness for distributions; these forms should be used when resampling (bootstrap or jackknife).
The "fisher" method correspond to the usual "unbiased" definition of sample variance, although in the case of skewness exact unbiasedness is not possible.

この説明からは，"moment" が先生の関数で method="ordinary"を，"fisher"が先生の関数でmethod="Excel"を指定した場合に相当するものと理解されます。
http://aoki2.si.gunma-u.ac.jp/R/univariate.html　

また，引数として渡すデータ x の説明には，a numeric vector or object. としか書かれていません。

コードは以下のようになっており，
function (x, na.rm = FALSE, method = c("moment", "fisher"), ...) 
{
    method = match.arg(method)
    ＜中略＞
    n = length(x)
    if (is.integer(x)) 
        x = as.numeric(x)
    if (method == "moment") {
        skewness = sum((x - mean(x))^3/sqrt(as.numeric(var(x)))^3)/length(x)
    }
    if (method == "fisher") {
        if (n < 3) 
            skewness = NA
        else skewness = ((sqrt(n * (n - 1))/(n - 2)) * (sum(x^3)/n))/((sum(x^2)/n)^(3/2))
    }
    attr(skewness, "method") <- method
    skewness
}
method == "moment"のときは，x - mean(x)のようにちゃんと平均が引かれているのに，method == "fisher"のときは，- mean(x)　がありません。やはり間違って使ってしまうように思います。

method == "moment"　につきましても，分母　sqrt(as.numeric(var(x)))^3)　は， sqrt(as.numeric(var(x)*(n-1)/n))^3)　とするか，(sum((x-mean(x))^2)/n)^(3/2))のよ うにすべきではないかと思います。

実際，以下の計算をしますと，
> x<-rnorm(20,10,1)
> skewness(x)
[1] 0.4649626
attr(,"method")
[1] "moment"
> 
> skew(x,method="Excel")　　 #青木先生の関数
[1] 0.5438159
> skew(x,method="ordinary")　#青木先生の関数
[1] 0.502149
のようになり，青木先生の関数のいずれのmethodの結果とも一致しませんでした。

一方，Rのskewnessの式を下記のように書き換えますと，
method == "moment"のとき
        skewness = sum((x - mean(x))^3/sqrt(as.numeric(var(x)*(n-1)/n))^3)/length(x)
method == "fisher"のとき
　　　　skewness = ((sqrt(n * (n - 1))/(n - 2)) * (sum((x-mean(x))^3)/n))/((sum((x-mean(x))^2)/n)^(3/2))

> mySkew(x,method="fisher")
[1] 0.5438159
attr(,"method")
[1] "fisher"
> mySkew(x,method="moment")
[1] 0.502149
attr(,"method")
[1] "moment"
となり，先生の関数の結果と一致しました。
やはり，Rの関数は変だと思います。私にはかなり敷居が高いのですが，パッケージのメインテナーに連絡してみようと考えています。

初歩的な質問に，たいへん親切にご対応くださいましてありがとうございました。
また，返信が遅くなってしまいましたこと，お詫びいたします。
No.14214　Re: 因子得点　　【青木繁伸】　2011/01/24(Mon) 21:57

method="fisher" のときは，skewness, kurtosis を計算する式の前に，x <- scale(x) を挿入すればよいですね。
skewness

        else {
            x <- scale(x)
            skewness = ((sqrt(n * (n - 1))/(n - 2)) * (sum(x^3)/n))/((sum(x^2)/n)^(3/2))
        }

kurtosis

    if (method == "fisher") {
        x <- scale(x)
        kurtosis = ((n + 1) * (n - 1) * ((sum(x^4)/n)/(sum(x^2)/n)^2 - 
            (3 * (n - 1))/(n + 1)))/((n - 2) * (n - 3))
    }
No.14215　Re: 因子得点　　【新参者】　2011/01/24(Mon) 23:55

青木繁伸先生

再度，深夜にご返信とスマートな解決策，ありがとうございました。
つい今しがた，パッケージ timeDate のメインテナーに連絡のメイルを送ってみました。通じるか，心もとないかぎりですが…
重ねて，厚くお礼申し上げます。
No.14222　Re: 因子得点　　【新参者】　2011/01/25(Tue) 21:13

昨晩，メインテナーにメールを送ったところ，本日，早速返信が届いていました。
thanks I will care about this
という簡単なものでしたが，レスポンスの速さに感激しています。
ご報告まで。

● 「統計学関連なんでもあり」の過去ログ--- 044 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る