ポリコリック相関係数のThresholdsの意味

院生からのメールで、Rのpolycorパッケージのpolychor()やpolyserial()のthreshold=TRUEで表示されるThresholdsの意味について問い合わせがあったので調べてみた。以下返信メールの該当部分を貼り付ける。

======（ここから）======
Thresholdsは計算途中で使われる値で、  

https://www.researchgate.net/profile/Fritz-Drasgow/publication/313196484_Polychoric_and_polyserial_correlations/links/5968d1d2458515e9afa7a0b8/Polychoric-and-polyserial-correlations.pdf 

や  

https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/procstat/procstat_corr_details16.htm 

を見ると、その順序カテゴリまでに含まれるサンプルの割合の正規分布のパーセント点です。講義テキスト（第８章、pp.112）の例だと、 

> library(MASS)
> library(polychor)
> survey$Smoke <- factor(survey$Smoke, levels=levels(survey$Smoke)[c(2, 3, 4, 1)])
> survey$Exer <- factor(survey$Exer, levels=levels(survey$Exer)[c(2, 3, 1)])
> polychor(survey$Smoke, survey$Exer, std.err=TRUE)

Polychoric Correlation, 2-step est. = 0.1269 (0.1011)
Test of bivariate normality: Chisquare = 4.274, df = 5, p = 0.5107

  Row Thresholds
1 0.8447
2 1.1820
3 1.6790

  Column Thresholds
1 -1.29600
2  0.03187
となりますが、

> print(dat <- xtabs(~Smoke+Exer, data=survey))
       Exer
Smoke   None Some Freq
  Never   18   84   87
  Occas    3    4   12
  Regul    1    7    9
  Heavy    1    3    7  

> print(RR <- rowSums(dat))
Never Occas Regul Heavy
  189    19    17    11  

> print(CC <- colSums(dat))
None Some Freq
  23   98  115 

のように、行ごと、列ごとの合計を計算しておいて、

> qnorm(cumsum(RR)/sum(RR))
    Never     Occas     Regul     Heavy
0.8446521 1.1817923 1.6786497       Inf 

> qnorm(cumsum(CC)/sum(CC))
       None        Some        Freq
-1.29617484  0.03186931         Inf 

として出てくる値がThresholdsです。計算途中で使われる値で、分布の様子を示しますが、とくに意味を解釈する必要はありません。

コメント

コメントを投稿

このブログの人気の投稿

統計教育サンプルデータ作成プロジェクト

このところ保健学研究共通特講IV, VIIIにおける統計教育の改善方法についてつらつら考えていて（というか結果的には無駄になってしまった某業務のため、強制的に内容の再整理を考える必要があったのだが）、自由に使える適当なデータが欲しいと思うようになった。作図の説明のためには、量的なデータと質的なデータが混在している何百人規模のものが適していて、これまではMASSパッケージのsurveyという、アデレード大学の統計関係の学生のデータを使ってきた。自分のデータを公開するのは、匿名化したとしても倫理審査関係のクリアが難しいなあと思って、そこで思考が停まっていたが、ふと、架空の人物について想像したデータをバザール形式で作ってしまえば良いのではないかと思いついた。 質問紙はウェブのフォームで作り、完全匿名で（何の登録も必要ないゲストとして）ボランティアを募って（botに入力されないような仕組みは必要だと思うが、メール認証までしなくても、良くある CAPTCHA 、例えばPerlのcgiなら このモジュールを使う とか、あるいは Turstile で十分だろう）、自分が好きな小説や漫画の登場人物について作品中で明らかになっている情報や、性格的には読者が想像した情報で、登場人物の気になって埋めてもらうのだ。そうすれば、個人情報にも抵触しないし、倫理的な問題も生じないはずだ。データやプロジェクトのタイトルを 読者が勝手に想像した架空人物の社会心理的特徴 、とでもするか。 たぶん架空人物も成長するので、何歳時点のデータというのを想像した場合は何度でも入力できるようにしたい。例えば20代のシャーロック・ホームズと、ライヘンバッハの滝から生還後の40代のシャーロック・ホームズでは別データとなる。高村薫『李歐』の李歐や吉田一彰は、22歳の時と再会後では別データとなる。名前、年齢、身長、体重（を入れるためには、人間に限るとしなくてはデータにならないが、明らかに定量的なデータとしては身長と体重くらいしか想像で入力できるものが思いつかない）、性別（架空人物の場合は、新井素子『二分割幽霊奇譚』の主人公のようなインターセックスの登場人物は珍しくないので、それ以外と不明というカテゴリは必須だな）、同居世帯構成（ワンピースのサニー号乗組員のような場合もあるので、血縁関係なしが何人、というカテゴリが実在す...
続きを読む

差がゼロという帰無仮説以外の仮説の下でのp値の計算

イメージ
『宇宙怪人しまりす　統計よりも重要なことを学ぶ』朝倉書店（以下しまりす本３）のpp.12「回復割合の差は20%」のp値を計算しようとしたが合わない。何かぼくのミスだと思うが原因がわからない。 x <- c(58, 62); n <- c(80, 100); X <- cbind(x, n-x) colnames(X) <- c("回復", "未回復") rownames(X) <- c("ヨクナール", "経過観察") print(X) s <- c(colSums(X)[1]/sum(X), colSums(X)[2]/sum(X)) # 帰無仮説：回復割合が処置と独立の場合の期待度数E print(E <- rbind(ヨクナール=rowSums(X)[1]*s, 経過観察=rowSums(X)[2]*s)) E/rowSums(E) # 確認：回復割合にヨクナールと経過観察で差はない 1-pchisq(sum((X-E)^2/E), 1) # 帰無仮説の下でのp値 Z <- c(0.1, -0.1) # ヨクナール群の回復割合が経過観察群より20%良い場合の期待度数EE print(EE <- rbind(ヨクナール=rowSums(X)[1]*(s+Z), 経過観察=rowSums(X)[2]*(s-Z))) EE/rowSums(EE) # 確認：回復割合がヨクナールの方が20％良い 1-pchisq(sum((X-EE)^2/EE), 1) # ヨクナールの回復割合が20%高い仮説の下でのp値 というRコードを実行すると、 差がないという帰無仮説（しまりす本３では「ゼロ仮説」）の下でのp値は本と同じく13.8%になるのだけれども、最後に表示される「回復割合の差が20%」という仮説の下でのp値は16.4%になってしまって、しまりす本３に書かれている17.2%とは合わない。しまりす本３ではpp.13-14に17.2%の算出方法が書かれていて、上のコードで使ったカイ二乗適合度検定とは違う方法（下記）なのだけれども、どちらでも一致するはずなんだがなあ。自分の頭の回転が鈍くて嫌になる。 rdp <- function(a, ...
続きを読む

対応のある多群間のノンパラメトリックな比較のためのFriedmanの検定後の対比較について

Friedmanの検定後のpost hocな多重比較の結果がEZRとSPSSとjamoviで食い違うのは何故かという問い合わせメールをもらったので調べてみた。以下メールの返事を貼り付けて若干改変。 （返信第1弾）============================= jamoviは、 https://bookdown.org/sbtseiji/jamovi_complete_guide/sec-ANOVA-friedman.html に書かれている通り、Friedmanの後の対比較ではbonferroni補正なしのDurbin-Conover法を使っています。  このDurbin-Conover法というのは検定の多重性の調整をしていません。内部的にはPMCMRplusというパッケージに含まれている、frdAllPairsConoverTest()という関数を使っているようで、 この関数自体にはオプションでp.adj="bon"を指定できる のですが、jamoviではそこをいじれない実装になっています。 もっとも、Bonferroniの補正という手法は、対比較で出てきたp値に比較回数を掛けるだけなので、5群間の2群ずつの総当たり比較なら10回の比較がされているので、得られたp値を10倍すれば良いわけです。 （データを使った検証部分は省略） EZRが何故合わないかですが、EZRは対比較にWilcoxonの 符号順位検定を使っています（とメッセージにでてきます）。EZRでは正規近似でp値を出しているので、同順位の値が含まれていると、 wilcox.test.default(pairwise.response1, pairwise.response2, alternative = "two.sided", で警告がありました: タイがあるため、正確な p 値を計算することができません という警告が出ます。 ここで表示されているp値は、wilcox.test(..., paired=TRUE)で得られるp値の10倍になっています。同順位があっても正確なp値が欲しい場合は、exactRankTestsパッケージのwilcox.exact()関数を使えばOKです（もちろん個々のp値は10倍する必要があります）。 （返信第2弾）===...
続きを読む