Juliaで統計解析プログラミング入門：１０万個の子宮

Julia で統計解析をやるためのプログラミング入門用の例を示す。解析対象は村中璃子著「10万個の子宮」(2018)から、26の身体症状についてワクチン接種の有無と症状の有無を調べたデータで、各症状ごとにカイ二乗検定（およびフィッシャーの正確確率検定）をやってｐ値を計算する。

まずはデータファイル。CSVファイル。

```Vaccine_Symptom.csv
#n=30279,非接種無症状,非接種有症状,接種無症状,接種有症状,不詳
月経不順,6812,2330,15354,5515,268
月経量の異常,8569,565,19205,1638,302
関節や体が痛む,8412,729,19324,1522,292
ひどく頭が痛い,8232,928,18714,2168,237
身体がだるい,8116,1047,18587,2291,238
すぐ疲れる,8163,996,18578,2297,245
集中できない,8433,728,19407,1448,263
視野の異常,8986,171,20470,388,264
光を異常にまぶしく感じる,8802,359,19964,915,239
視力が急に低下した,8358,799,19466,1400,256
めまいがする,8060,1095,18564,2299,261
足が冷たい,8004,1155,18317,2536,267
なかなか眠れない,8454,698,19379,1492,256
異常に長く寝てしまう,8080,1073,18357,2488,281
皮膚が荒れてきた,8076,1076,18789,2081,257
過呼吸,8834,333,20187,704,225
物覚えが悪くなった,8944,220,20257,632,228
簡単な計算ができなくなった,9082,81,20697,189,230
簡単な漢字が思い出せなくなった,8986,185,20471,417,220
身体が自分の意思に反して動く,9107,58,20689,200,225
普通に歩けなくなった,9135,22,20811,73,238
杖や車椅子が必要になった,9139,16,20853,33,238
突然力が抜ける,9054,100,20586,284,255
手や足に力が入らない,9007,124,20461,357,330
その他１,2641,118,5539,528,21453
その他２,2467,28,5201,89,22494
```
まずは必要なパッケージを読み込み、データファイルを読み込んでデータフレームオブジェクトに入れ、データの概要を見てみる。以下、`Julia>` は Julia の実行環境である REPL が自動的に出力するプロンプトである。

```Julia
Julia> using CSV, DataFrames, HypothesisTests, Printf

Julia> dat = CSV.read("Vaccine_Symptom.csv")  # データファイル -> データフレーム
Julia> describe(dat)                          # データフレームの内容の概要の表示
6×8 DataFrame. Omitted printing of 2 columns
│ Row │ variable     │ mean    │ min        │ median  │ max          │ nunique │
│     │ Symbol       │ Union…  │ Any        │ Union…  │ Any          │ Union…  │
├─────┼──────────────┼─────────┼────────────┼─────────┼──────────────┼─────────┤
│ 1   │ #n=30279     │         │ すぐ疲れる │         │ 集中できない │ 26      │
│ 2   │ 非接種無症状 │ 8075.12 │ 2467       │ 8443.5  │ 9139         │         │
│ 3   │ 非接種有症状 │ 578.231 │ 16         │ 462.0   │ 2330         │         │
│ 4   │ 接種無症状   │ 18393.3 │ 5201       │ 19393.0 │ 20853        │         │
│ 5   │ 接種有症状   │ 1307.08 │ 33         │ 1157.5  │ 5515         │         │
│ 6   │ 不詳         │ 1925.46 │ 220        │ 256.0   │ 22494        │         │

```

まずは基本的な処理の例。カイ二乗検定を行う ChisqTest() 関数は、データの型が整数のベクトルか行列でないと受け取ってくれない（データフレームそのままだと拒否されてしまう）ので、整数ベクトルを作って、それを２ｘ２の行列の形にしてから ChisqTest() に渡す。以下、REPL のプロンプト `Julia>` は省略する。

```Julia
N = length(dat[:, 1]) # データ行数 -> N
V = zeros(Int64, 4)   # ChisqTest() にはデータフレームではなく整数ベクトル/行列が必要
for i in 1:N
    for j in 1:4
        V[j] = dat[i, j+1]             # データフレームのi行j+1列 -> V[j]
    end
    M = reshape(V, 2, 2)               # ベクトルを行列に整形
    res = ChisqTest(M)                 # 行列に対してカイ二乗検定
    print(pvalue(res), dat[i,1], "\n") # ｐ値と症状の表示
end
```

これで、以下のような出力が出てくる。

```
0.08806870795706069月経不順
3.227130107483255e-7月経量の異常
0.04147764676986717関節や体が痛む
0.5098118395881822ひどく頭が痛い
0.24991670604733196身体がだるい
0.741755519047255すぐ疲れる
...
```

もうちょっと見やすく出力させたいし、カイ二乗検定以外の方法との比較もできた方がよい｡また各症状のｐ値が有意水準と比較してどうかも表示させたい。ので、そうさせてみる。カイ二乗とフィッシャーのそれぞれで、ｐ値が 0.05 以下だったらその症状のところに `!` を表示させてみる｡

```Julia
@printf "\np.chi p.fis C F Symptom\n----------------------------------------\n"
for i in 1:N
    for j in 1:4 V[j] = dat[i, j+1] end          # for 文は一行でもかける
    R1 = ChisqTest(reshape(V, 2, 2))             # カイ二乗検定
    R2 = FisherExactTest(V[1], V[3], V[2], V[4]) # フィッシャーの正確確率検定
    S1 = pvalue(R1) < 0.05 ? "!" : " "           # カイ二乗検定のｐ値を有意水準と比較
    S2 = pvalue(R2) < 0.05 ? "!" : " "           # フィッシャーのｐ値を有意水準と比較
    @printf "%.3f %.3f %s %s %s\n" pvalue(R1) pvalue(R2) S1 S2 dat[i,1]
end
```
これによって、以下のような出力が得られる。だいぶ見やすい。有意水準による判断はどの症状でも、カイ二乗検定とフィッシャーの正確確率検定で同じになった。

```
p.chi p.fis C F Symptom
----------------------------------------
0.088 0.090     月経不順
0.000 0.000 ! ! 月経量の異常
0.041 0.045 ! ! 関節や体が痛む
0.510 0.524     ひどく頭が痛い
0.250 0.258     身体がだるい
0.742 0.758     すぐ疲れる
0.002 0.002 ! ! 集中できない
0.966 0.998     視野の異常
...
```