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> 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> は省略する。
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 以下だったらその症状のところに ! を表示させてみる｡
@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 視野の異常
...