統計相談に乗ったときのメモ

『宇宙怪人しまりす　統計よりも重要なことを学ぶ』朝倉書店（以下しまりす本３）の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,

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回くらい繰り返して仮想データセットを作り、必要な統計解析をした後でそれを統合するという方法を試した先人は見つからなかった。欠損値への多重代入法で似たようなことは許されているので、これも理屈の上ではできると思うが、仮想データで試してみたいところではある。そのためにも、小説や漫画の登場人

院生からのメールで、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.0318

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

小説や漫画の登場人物のスペックを読者が勝手に想像して入力することで、統計教育に使えるオープンデータを作ろうという思いつきについて、AIでできるかも？　と思いついたことから、試してみたが、失敗したという話である。詳細は以下。 Microsoft Edgeの右上のアイコンからCopilotを呼び出して、 小説や漫画の登場人物300人分の社会心理的特徴のデータをcsv形式で作ってください。 変数としては、名前、年齢、身長、体重、性別（それ以外と不明というカテゴリは必須）、同居世帯構成、親しい友人の数、自己効力感、Big5性格特性、PHQ-9、PHQ-15、GAD-7、ソーシャルサポート、ソーシャル・キャピタルの尺度を含めてください。 と質問を打ってみた。すると以下の回答が表示された。 もちろんです。小説や漫画の登場人物に関する社会心理的特徴のデータをCSV形式で作成します。以下にリストされた変数を含めて、300人分のデータを生成します。 **名前**: 登場人物の名前 **年齢**: 登場人物の年齢 **身長**: 登場人物の身長（単位はメートル） **体重**: 登場人物の体重（単位はキログラム） **性別**: 登場人物の性別（男性、女性、それ以外、不明のいずれか） **同居世帯構成**: 登場人物がどのような世帯で暮らしているか（単身、家族、ルームメイトなど） **親しい友人の数**: 登場人物が親しい友人とどれだけ交流しているか **自己効力感**: 登場人物の自己効力感の尺度 **Big5性格特性**: 登場人物の性格特性をBig Fiveモデルで評価 **PHQ-9**: 登場人物のうつ症状の程度を評価するための尺度 **PHQ-15**: 登場人物の身体的な不調の程度を評価するための尺度 **GAD-7**: 登場人物の一般的な不安症状の程度を評価するための尺度 **ソーシャルサポート**: 登場人物がどれだけソーシャルサポートを受けているか **ソーシャル・キャピタル**: 登場人物のソーシャル・キャピタルの尺度 これらの変数を含むCSVファイルを作成します。お待ちください。 ソース: Bing との会話 2024/2/3 (1) 物語はいかにして心を動かすのか： 物語説得研究の現状と態度 .... https://www.jstage.js