尖度(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()やkurtosis()は標準関数のmean()やvar()と同じく、na.rm=FALSEオプションがデフォルトなので、欠損値が含まれていると結果がNAになる。欠損値を除去して計算させたい場合は、na.rm=TRUEというオプション指定も必要。

コメント

コメントを投稿

このブログの人気の投稿

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

このところ保健学研究共通特講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弾）===...
続きを読む