投稿

尖度(kurtosis)と歪度(skewness)についての、R/EZR/jamoviの実装

尖度(kurtosis)と歪度(skewness)について質問があり、R/EZR/jamoviでの実装を調べてみた。 R本体には関数として実装されていないが、e1071パッケージに入っていて、helpを見たら、 Joanes and Gill (1988) に基づいて、それぞれ3通りの計算方法が実装されているとわかった。 古典的定義（歪度を3次のモーメントを2次のモーメントの平方根の3乗で割った値、尖度を4次のモーメントを2次のモーメントの二乗で割って３を引いた値とする定義。ちなみにr次のモーメントとは平均からの偏差のr乗和をサンプルサイズで割った値である）に基づくのが type=1 で、SASやSPSSと同じ実装が type=2 （e1071パッケージのヘルプでは表示されないが、Joanes and Gill, 1988によると、EXCELもこの実装を採用しているとのこと）、MINITABやBMDPと同じ実装が type=3 で、どちらも デフォルトは type=3 とのこと。 正規分布の下で不偏推定量であることがわかっているのは、歪度の方は３つともだが、尖度は type=2 だけとのことなので（なぜ type=2 をデフォルトにしなかったのか謎だが）、R本体で library(e1071) から skewness() や kurtosis() を使う場合は、特別な理由がない限り、 type=2 を指定すべきと思う。 EZRでは統計解析＞連続変数の解析＞正規性の検定でアウトプットウィンドウの中に表示されるが、 type=2 の結果と同じようだ。 jamoviは探索(Exploration)の記述統計(Descriptives)で変数を選んで計算させた後で統計量(Statistics)の分布(Distribution)の下にある歪度(skewness)と尖度(kurtosis)にチェックを入れると type=2 と同じ結果が表示される。 医学統計ではSASやSPSSと同じ結果になるのが無難な選択であろうから、EZRの方針は合理的だと思う。元々心理統計のために開発されたらしいjamoviも、最近は医学論文でも使われているから、妥当な選択であろう。 なお、EZRやjamoviでは欠損値があっても自動的に除去して計算してくれるが、e1071パッケージの skewness(...
コメントを投稿
続きを読む

対応のある多群間のノンパラメトリックな比較のための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弾）===...
1 件のコメント
続きを読む

差がゼロという帰無仮説以外の仮説の下での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, ...
1 件のコメント
続きを読む

baseグラフィクスの方が楽

イメージ
gridグラフィクスが流行している理由が、やっぱり良くわからない。 例えば、 bob3bob3さんの『データ可視化学入門』をPythonからRに翻訳した話 で、日本の人口変化をChatGPT4ベースのCopilotを使ってPythonからRに翻訳させるとggplot2を使ったコードが出てくるようだが、この程度ならデータ定義部分だけ元コードからコピペして X <- data.frame(year = 1990:2022, pop = c(1.23611, 1.24101, 1.24567, 1.24938, 1.25265, 1.2557, 1.25859, 1.26157, 1.26472, 1.26667, 1.26926, 1.27316, 1.27486, 1.27694, 1.27787, 1.27768, 1.27901, 1.28033, 1.28084, 1.28032, 1.28057, 1.27834, 1.27593, 1.27414, 1.27237, 1.27095, 1.27042, 1.26919, 1.26749, 1.26555, 1.26146, 1.25502, 1.24947))   とデータフレーム X を定義しさえすれば、コードとしては、 plot(pop ~ year, data=X, xlab="西暦", ylab="日本の総人口[億人]", type="b", col="blue", pch=16, las=1); grid() で済む。パイプを使うよりも、データはいったんオブジェクトとして定義（付値）しておく方がRの思想的にもわかりやすいし、すべて標準のbaseグラフィクスで書けば、間違いなくどのR環境でも動くわけだから、Copilotすらggplot2を使うコードを出してくるのが謎でならない。念の為、bob3bob3さんのコード移植にはAI活用が便利という発表主旨には賛成で、スライドも素晴らしいと思うが、世の中で広まっているやり方に引きずられるのは欠点だと思う。  
コメントを投稿
続きを読む

分位数によるカテゴリ化の是非

そもそも連続量が正規分布していない場合に分位数で区切ってカテゴリ化することの是非はどうなのかも書いておく。 2012年に出た Bennette C, Vickers A: "Against quantiles: categorization of continuous variables in epidemiologic research, and its discontents", BMC Med Res Methodol , 12: 21 はデメリットを強く訴えていて、著者の一人は2020年にも、本人のアカウントかどうか知らないがGreenland大先生への＠tweetとして 分位数を巡るカテゴリ化は昔から論争の的であった と書いていて、それについているコメントで知った このスライド など見ると、怖くて連続量のカテゴリ化などできなくなりそうではある。もっと新しいものをみても、 Busch EL (2021) Cut points and contexts. Cancer, 127: 4348-4355. というコメンタリーは、連続量をカテゴリ化する方法は一長一短で、どれか１つのカテゴリ化が正しいというようなものは存在しないと論じている。 一方、 このSASのブログ記事 は、（解決策の一案として分割するカテゴリ数を減らすことを提案しているだけで、根本的な問題解決にはなっていないが）まさに丸めによって整数値になっていたり、元々整数値だったものを分位数によってカテゴリ化する際に同順位が多くなることで生じる問題を指摘している。結局、同順位のデータがちょうど分割したい分位数になっている場合に、その数値を上の区間に含めても下の区間に含めても歪んでしまうので、truequantile()のような方法で実数の分位数を見つけて、ちょうど分位数になっている同順位の人数を上と下に按分して、どの人を上の区間に割り当て、どの人を下の区間に割り当てるかを乱数で決めるという動作を1000回くらい繰り返して仮想データセットを作り、必要な統計解析をした後でそれを統合するという方法を試した先人は見つからなかった。欠損値への多重代入法で似たようなことは許されているので、これも理屈の上ではできると思うが、仮想データで試してみたいところではある。そのためにも、小説や漫画の登場人...
コメントを投稿
続きを読む

ポリコリック相関係数の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   Colu...
コメントを投稿
続きを読む

truequantile( )というRの関数を書いてみた

院生から相談があって、通常連続量とみなして分析する整数のスケールを分位数で区切ってカテゴリ化する場合に（そういう手口の論文は多数ある）、サンプルサイズがある程度大きいと、ちょうど分位数で同順位になる人が複数いて、その分位数を下と上のどちらの範囲カテゴリに入れるかで結果が大きく変わってしまうので、どちらに入れたら良いのか？　という問題を考えてみた。このカテゴリ化の意味は、順位によってほぼ同数のサンプルサイズのサブグループに分けるということだから、分位数ちょうどのサンプルをどちらの範囲カテゴリに含めてもサイズのアンバランスは不可避であって、ベストな解はないと思う（そもそもカテゴリ化によって情報量は落ちるが、分布によっては悪くない）。先行研究があるなら、それに従うべきであろう。中央値より上か下かで２つのサブグループにするだけなら、 fmsbパッケージのtruemedian()関数 （リンク先pdfのpp.39-41）を使えば、たぶん同順位があっても、同順位の測定値が測定単位より高い精度で測ったら等間隔に分布すると仮定した場合の真の中央値が出るので、下の区間にも上の区間にもその値を含まないことにすれば、一応問題は解決する。例えば、8という値をもつ人が多数いて、それが中央値だった場合、真の中央値が7.5以上8未満であれば8は上の区間に含まれることになり、8以上8.5未満であれば下の区間に含まれることになる（真の中央値もちょうど8になってしまった場合にどうするかという問題は残るが、タイが多いのにそうなってしまう可能性は無視できるくらい小さいと思うので、とりあえず気にしないことにする）。もちろん、等間隔に分布すると仮定すれば、同順位となる分位数ちょうどのデータの一部は下の区間、残りは上の区間に振り分けられるべきだが、実際に同順位となった値のどの個体が下の区間、どの個体が上の区間に振り分けられるべきかを判定する方法が原理的に存在しないので、それは諦めるしかない。 しかしそれさえ我慢すれば、かつて 奥村さんが議論されていた ように、中央値だけではなく、任意の分位数について、同じ手口が使えるはずである。fmsbパッケージのtruemedian()関数は、等間隔の分布を仮定した真の値を生成してからmedian()をとっているので、returnされるmedian(YY)をquantile...
1 件のコメント
続きを読む