私はもともとこれをクロス検証に投稿しましたが、純粋にソフトウェア構文に関するものであるため、SO の方が適切かもしれないと思います。
これは、この投稿のフォローアップの質問 です。多項ロジスティック回帰を実行して、回答者がさまざまな病状 ( pain、sleep、メンタルヘルス/物質使用 ( mhsu) および他のすべての状態 ( allOther)) を医療用大麻licitまたはillicit医療用大麻で治療したことを示す対数オッズの差を調べました。
グッズデータはこちら
df <- tibble(mcType = factor(rep(c("licit", "illicit"),
times = c(534,1207))),
cond = factor(c(rep(c("pain","mhsu","allOther","sleep"),
times = c(280,141,82,31)),
rep(c("pain","mhsu","allOther","sleep"),
times = c(491,360,208,148))),
levels = c("pain","sleep","mhsu","allOther")))
そして、大麻の種類ごとに報告された各種類の状態の割合
mcType cond n tot perc
<fct> <fct> <int> <int> <dbl>
1 illicit pain 491 1207 40.7
2 illicit sleep 148 1207 12.3
3 illicit mhsu 360 1207 29.8
4 illicit allOther 208 1207 17.2
5 licit pain 280 534 52.4
6 licit sleep 31 534 5.81
7 licit mhsu 141 534 26.4
8 licit allOther 82 534 15.4
使用した大麻の種類に基づいて、各タイプの状態を示す回答者の相対的な割合に違いがあるかどうかを確認するためmultinom()に、nnetパッケージで使用して多項ロジスティック回帰を実行しました。以下の出力、
library(nnet)
summary(mm <- multinom(cond ~ mcType,
data = df))
# output
Coefficients:
(Intercept) mcTypelicit
sleep -1.1992431 -1.0014884
mhsu -0.3103369 -0.3756443
allOther -0.8589398 -0.3691759
Std. Errors:
(Intercept) mcTypelicit
sleep 0.09377333 0.2112368
mhsu 0.06938587 0.1244098
allOther 0.08273132 0.1503720
Residual Deviance: 4327.814
AIC: 4339.814
emmeansパッケージを使用して、簡単な効果のテストを実行しました。このブログ投稿で著者は、emmeans パッケージがデフォルトでエラー修正を適用することを提案していますが、これはadjust = 引数を介して制御できます。
# testing effect of mc type at each level of condition. first create emmeans object
library(emmeans)
(em_mcTypeByCond <- emmeans(object = mm,
specs = ~mcType|cond,
adjust = "bonferroni"))
# output
cond = pain:
mcType prob SE df lower.CL upper.CL
illicit 0.4068 0.01414 6 0.3648 0.4488
licit 0.5243 0.02161 6 0.4602 0.5885
cond = sleep:
mcType prob SE df lower.CL upper.CL
illicit 0.1226 0.00944 6 0.0946 0.1506
licit 0.0581 0.01012 6 0.0280 0.0881
cond = mhsu:
mcType prob SE df lower.CL upper.CL
illicit 0.2983 0.01317 6 0.2592 0.3374
licit 0.2641 0.01908 6 0.2074 0.3207
cond = allOther:
mcType prob SE df lower.CL upper.CL
illicit 0.1723 0.01087 6 0.1401 0.2046
licit 0.1535 0.01560 6 0.1072 0.1999
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 2 estimates
問題は、エラー訂正の他の方法 (たとえば、"BH"、"fdr"、"westfall"、"holm") を選択できないように見えることです。テストを適用する前など、間違ったステップで修正を適用しているためかどうかはわかりません。
そこで、pairs()関数内で調整引数を適用してみました（2種類の大麻間の各条件の確率の違いをテストします）
(mcTypeByCond_test <- pairs(em_mcTypeByCond,
adjust = "bonferroni"))
cond = pain:
contrast estimate SE df t.ratio p.value
illicit - licit -0.1175 0.0258 6 -4.551 0.0039
cond = sleep:
contrast estimate SE df t.ratio p.value
illicit - licit 0.0646 0.0138 6 4.665 0.0034
cond = mhsu:
contrast estimate SE df t.ratio p.value
illicit - licit 0.0342 0.0232 6 1.476 0.1905
cond = allOther:
contrast estimate SE df t.ratio p.value
illicit - licit 0.0188 0.0190 6 0.987 0.3616
しかし、ご覧のとおり、どのタイプのエラー修正が適用されたかを示すメッセージは表示されません (何もないと仮定して、いくつかの異なる方法を試しました)。また、4 つのペアごとの比較すべてでエラーを制御したいと考えています。
したがって、p 値の調整を指定する引数をどのように、どの段階で作成する必要があるかを知る必要があります。
どんな助けでも大歓迎です