8.0 パッケージの準備

下記のパッケージを使用

library(tidyverse)
library(pwr)
library(epiR)

8.1

発展途上国では、正確なメカニズムは明らかにされていないが、ビタミンAの補給が感染症に対する耐性を付与することが知られている。このメカニズムを解明するために、健康な子どもたちにビタミンAのサプリメントを長期にわたって投与する。ベースライン時とサプリメント投与期間終了時に、すべての子供たちから尿を採取し、ネオプテリンの排泄量を測定する。現在のネオプテリンのレベルはクレアチニン1mmolあたり600mmolと推定され、サプリメントによる変化の標準偏差は（小さなパイロットサンプルから）クレアチニン200mmol/mmolと推定されている。

※ネオプテリン：サイトカインのひとつであるIFN-γ刺激に伴いマクロファージによって作られる物質。細胞性免疫の活性指標としてよく使われている。要は、ここでは感染症があると上がる。

8.1.1

仮にネオプテリンの10%の低下が意味のある差であると考え、両側1%検定を用いて95%の検出力で見つけたい。何人の子供を募集すべきか？

effect sizeを算出する関数を作成する。

es_two <- function(H0, H1, SD){
  samplesize <- abs(H1-H0)/SD
  return(samplesize)
}
es_dif <- function(dif, SD){
  samplesize <- abs(dif)/SD
  return(samplesize)
}


まずはpwrパッケージを使用する。
one sample, t-testのサンプルサイズ計算。

es_811 <- es_two(600, 600*(1-0.1), 200)
ans_811 <- pwr.t.test(d = es_811, 
                      sig.level = 0.01, 
                      power = 0.95, 
                      type = "one.sample", 
                      alternative = "two.sided")
ans_811
## 
##      One-sample t test power calculation 
## 
##               n = 201.2665
##               d = 0.3
##       sig.level = 0.01
##           power = 0.95
##     alternative = two.sided

答え＝202人

8.1.2

150人の子供にしかサプリメントを与えられないとする。1%の両面テストでネオプテリンの10%の低下を検出する確率はどれくらいか？

ans_812 <- pwr.t.test(n = 150, 
                      d = es_811, 
                      sig.level = 0.01, 
                      type = "one.sample", 
                      alternative = "two.sided")
ans_812
## 
##      One-sample t test power calculation 
## 
##               n = 150
##               d = 0.3
##       sig.level = 0.01
##           power = 0.8548325
##     alternative = two.sided
round(100*ans_812$power, 1)
## [1] 85.5

答え＝85.5%

8.1.3

150人の子供が募集された場合、1%の両側検定を用いて95%の検出力で検出できると期待できる、ネオプテリンの最小の割合の減少は？

ans_813 <- pwr.t.test(n = 150, 
                      sig.level = 0.01, 
                      power = 0.95,
                      type = "one.sample", 
                      alternative = "two.sided")
ans_813
## 
##      One-sample t test power calculation 
## 
##               n = 150
##               d = 0.3485311
##       sig.level = 0.01
##           power = 0.95
##     alternative = two.sided
round(100* (200*ans_813$d / 600) , 1)
## [1] 11.6

答え＝11.6%

8.2

喫煙者と非喫煙者のBMIを比較する研究を行う。喫煙者のBMIが低いという対立仮説に対して、差がないという帰無仮説の片側5%の有意性検定を計画する。過去の経験から、BMIの標準偏差は約4kg/m^2であることが知られている。

8.2.1

喫煙者と非喫煙者が同数募集されたと仮定して、喫煙者のBMIが1kg/m2低い場合、95％の検出力で検出できるようにするには、どのようなサンプルサイズを用いるか？

es_821 <- es_dif(1, 4)
ans_821 <- pwr.t.test(d = -1*es_821, 
                      sig.level = 0.05, 
                      power = 0.95,
                      type = "two.sample", 
                      alternative = "less")
ans_821
## 
##      Two-sample t test power calculation 
## 
##               n = 346.9881
##               d = -0.25
##       sig.level = 0.05
##           power = 0.95
##     alternative = less
## 
## NOTE: n is number in *each* group
2*ans_821$n
## [1] 693.9763

H0>H1の片側検定なので、dをマイナスに、alternativeを“less”に指定する。

答え＝694人

8.2.2

喫煙者と非喫煙者が同数募集されたと仮定して、600人の被験者しか募集できない場合、1kg/m2の差を検出できると予想される検出力はどれくらいか？

ans_822 <- pwr.t.test(n = 600/2, 
                      d = -1*es_821,
                      sig.level = 0.05, 
                      type = "two.sample", 
                      alternative = "less")
ans_822
## 
##      Two-sample t test power calculation 
## 
##               n = 300
##               d = -0.25
##       sig.level = 0.05
##           power = 0.9212518
##     alternative = less
## 
## NOTE: n is number in *each* group
round(100*ans_822$power, 1)
## [1] 92.1

答え＝92.1%

8.2.3

200人の喫煙者と400人の非喫煙者を募集した場合、1kg/m2の差を検出できると予想される検出力はどれくらいか？

2群の人数が異なる場合は、別の関数を使う。

ans_823 <- pwr.t2n.test(n1 = 400, n2 = 200, 
                        d = -1*es_821,
                        sig.level = 0.05, 
                        alternative = "less")
ans_823
## 
##      t test power calculation 
## 
##              n1 = 400
##              n2 = 200
##               d = -0.25
##       sig.level = 0.05
##           power = 0.8922587
##     alternative = less
round(100*ans_823$power, 1)
## [1] 89.2

答え＝89.2%

8.2.4

ここで、別のタイプの研究を想定する。禁煙するつもりの被験者を募集し、禁煙前と数年後のBMIを記録する。その後のデータに適用する検定条件が8.2.1のような場合（5%の片側検定、1kg/m2の差、95%の検出力）、このデザインで、8.2.1で用いたデザインよりも少ないBMIの記録数となるのはどのような条件か？　被験者が時間の経過とともに「失われる」ことはないと仮定する。

（解答例）
一般的に、前後比較など対応のある2群の検定の場合は、そうでない場合と比較してサンプルサイズを節約できる。

※8.4.4 Paired dataを参照のこと
Provided that the correlation between the ‘before’ and ‘after’ measurements is high and positive (as is likely in practice), pairing will give a considerable saving in sample numbers.

8.3

（本文中の）例8.13と例8.14では、より一般的な両側検定ではなく、片側検定が使用された。それぞれのケースで両側検定が使用された場合，サンプルサイズの要件はどのように変わるか？

（解答例）
片側検定の有意水準がαであったのなら、両側検定ではα/2を代わりに用いる。Zα < Zα/2なのは自明なので、必要となるサンプルサイズも大きくなる。

8.4

Basnayakeら(1983)は、スリランカの女性は低用量経口避妊薬(OC)に適している(“more suited”)という仮説を検証するための介入研究を提案した。女性は、ノリニル（標準用量OC）またはブレビコン（低用量OC）のいずれかを使用するように無作為に割り付けられた。目的は、試験開始から12ヵ月後に割り付けられたOCを継続して使用している女性の割合を明らかにすることである。

8.4.1

研究参加者の適切な包含基準と除外基準を提案せよ。

（解答例）
包含基準：女性、初潮後かつ閉経前
除外基準：子宮摘出術後、著しい月経不順、ホルモン動態に影響を与える疾患（乳腺疾患、卵巣悪性腫瘍など？）の既往

OCの投与目的が避妊なのか、月経随伴症状の緩和なのかによっても変わる気がする。

8.4.2

研究者は、12ヶ月後にOCを継続する女性の割合がノリニル群よりもブレビコン群の方が10％高いことを90％確実に検出したいと考えている。
以前の研究では、ノリニル使用者の12ヶ月継続率は55%であった。本研究では、各治療群で同数の女性を対象とした両側5％有意差検定を使用する。
追跡調査中の脱落がないと仮定して、何人の女性がこの研究に登録されるべきか？　5%が12ヵ月以内に追跡調査で脱落すると予想される場合、何人の女性を研究に組入れるべきか？

比の検定のサンプルサイズ計算には、pwr.2p.testを使用する。
effect size計算には、ES.hが使用できる。
“10% higher”が55%+10%=65%なのか、55%×1.1=60.5%なのか微妙に迷ったが、前者と解釈すると8.5で算出するpowerが100%になるので違うと思い、後者の解釈とした。

ans_842 <- pwr.2p.test(h = ES.h(p1 = 0.55, 
                                p2 = 0.55*1.1),
                       sig.level = 0.05,
                       power = 0.90,
                       alternative = "two.sided")
2*ans_842$n
## [1] 3386.311

答え＝3388人

5%が脱落するとすると、

2*ans_842$n/(1-0.05)
## [1] 3564.538

答え＝3566人

8.4.3

2つのOCの特性を比較するために12ヶ月継続率のみを使用することに、批判的な点はあるか？　他に検討すべき方策を述べよ。

（解答例）
月経随伴症状の緩和を目的としたOCの使用を検討したいのであれば、月経随伴症状をPROとして評価してもよいかもしれない。
避妊が適切に行われているかの評価であれば、現実的には薬剤継続率で評価することになるだろう（望まない妊娠の発生率を評価することは現実的でない）。
primary endpointではないかもしれないが、有害事象についても客観的所見、PROいずれも評価すべきと思われる。

8.5

健康増進ユニットの研究者たちは、利尿薬の服用と高齢者の転倒との関係についての研究を計画している。利尿薬は血圧を下げるため、転倒の危険因子になる可能性があると考えられている。
研究では、薬局の薬剤師に、高齢者に処方された薬を処方する際に、過去1年間に転倒したことがあるかどうかを尋ねてもらうことを計画している。利尿剤をもらったことがあるかどうか、過去1年間に重大な転倒を経験したことがあるかどうか、という2つの変数を高齢者の顧客ごとに記録する。
合計2000人の被験者を募集する。

8.5.1

利尿薬を服用していない高齢者の3分の1は、過去1年以内に重度の転倒を経験したことがあると研究者は推定している。
利尿薬を服用している人の中で重度の転倒の可能性が20％高い場合を検出するために、各群で同数であると仮定して、両側5％の有意性検定を用いて、検定の検出力を推定せよ。

ans_851 <- pwr.2p.test(n = 2000/2, 
                       h = ES.h(p1 = 0.33, 
                                p2 = 0.33*1.2),
                       sig.level = 0.05,
                       alternative = "two.sided")
round(100*ans_851$power, 1)
## [1] 86.7

答え＝86.7%

8.5.2

提案された研究で予測される問題点について述べよ。

（解答例）
cross-sectional studyなので、両者の因果関係について明確な結論は得られない。

8.6

コーヒーを飲むことによる健康への影響を調べるために、研究チームは疫学調査を計画している。冠動脈性心疾患(CHD)は調査対象となる人口の中で死亡原因の第一位であるため、サンプルサイズの計算の基礎としてこれを選択した。
コーヒーを飲まない人のCHD発生率は100人あたり2.3人で、人口の40％がコーヒーを飲んでいることが知られている（ここでは、1日に4杯以上のコーヒーを定期的に飲むことを意味する）。
コーヒーを飲むとCHDのリスクが（少なくとも約）2倍になることを、95％の検出力で検出したいと考えている。コーヒーを飲むこととCHDの間に関連性がないという帰無仮説の、両側5％有意性検定を実施するために必要なサンプルサイズ（1:1のグループを仮定して）を、以下の研究デザインの下で算出せよ。

8.6.1

コホート研究

コホート研究やケースコントロール研究のサンプルサイズ計算には、epiRパッケージの関数を使用する。最新バージョン(2.0.17)でないとエラーが出るので注意する。

細かい引数の条件はmanualかvignetteを参照のこと（章末にリンクあり）。ここでは不要な引数は省略している。

ans_861 <- epi.sscohortc(irexp1 = 2*2.3/100, 
                         irexp0 = 2.3/100, 
                         pexp = 0.4, 
                         n = NA, 
                         power = 0.95, 
                         r = 1, 
                         sided.test = 2, 
                         conf.level = 0.95)
ans_861$n.total
## [1] 3268

答え＝3268人（1634人ずつ）

8.6.2

非マッチケースコントロール研究

ここではORはRRに近似できると考え、OR=2で求める。

ans_862 <- epi.sscc(OR = 2, 
                    p0 = 0.4, 
                    n = NA, 
                    power = 0.95, 
                    r = 1, 
                    sided.test = 2, 
                    conf.level = 0.95, 
                    method = "unmatched")
ans_862$n.total
## [1] 436

答え＝436人（218人ずつ）

8.6.3

コーヒーを飲まない人の特定の種類の癌のリスクは、1000人あたり8人であることも知られている。
CHDではなく癌をサンプルサイズの計算の基礎として選んだ場合（他はすべて変更なし）、8.6.1と8.6.2で必要なサンプルサイズはどのように変わるか？

①コホート研究

ans_863 <- epi.sscohortc(irexp1 = 2*8/1000, 
                         irexp0 = 8/1000, 
                         pexp = 0.4, 
                         n = NA, 
                         power = 0.95, 
                         r = 1, 
                         sided.test = 2, 
                         conf.level = 0.95)
ans_863$n.total
## [1] 9624

答え＝9624人（4812人ずつ）

②非マッチケースコントロール研究

本文やepi.ssccの引数からもわかるように、ケースコントロール研究のサンプルサイズは、アウトカムのリスク（頻度）に影響を受けない。そのため、8.6.2と同じ。

8.7

職業が精巣がんの発生率に及ぼす影響を確認するために、非マッチの症例対照研究を行う。症例は、特定の地域の病院に新たに紹介された患者である。対照者は、同じ病院に入院している患者の中から、事前に指定されたいくつかの許容可能な診断のうちのどれか1つを持っている患者から募集する。各被験者は、現在または最後の雇用を明らかにするために、面接を受けるか、または記録を参照する。記録されたさまざまな雇用形態は、その後、相互に排他的ないくつかのグループに結合される。
オッズ比が計算され、各職業グループを残りのすべての組み合わせと比較する。オッズ比（おおよその相対リスク）が2であることは非常に重要であると考えられているので、両側1%の有意性検定を使用する際には見逃すべきではない。
研究者らは、利用可能な時間内に400例の症例を募集することを予想している。必要であれば、1600人の対照者を募集できると推定している。

全国の男性成人人口の、
(8.7.1) 5％
(8.7.2) 10％
(8.7.3) 20％
が属する職業群について、症例と対照群を同数(各400人)用いた場合に、オッズ比2を検出する検出力を求めよ。

対照群1600人（症例400人と合わせて）をすべて用いた場合の、
(8.7.4) 5％
(8.7.5) 10％
(8.7.6) 20％
の職業群で、同様に求めよ。

8.7.1 ~ 8.7.3

ans_871 <- epi.sscc(OR = 2, 
                    p0 = 0.05, 
                    n = 400+400, 
                    power = NA, 
                    r = 1, 
                    sided.test = 2, 
                    conf.level = 0.95, 
                    method = "unmatched")
ans_872 <- epi.sscc(OR = 2, 
                    p0 = 0.1, 
                    n = 400+400, 
                    power = NA, 
                    r = 1, 
                    sided.test = 2, 
                    conf.level = 0.95, 
                    method = "unmatched")
ans_873 <- epi.sscc(OR = 2, 
                    p0 = 0.2, 
                    n = 400+400, 
                    power = NA, 
                    r = 1, 
                    sided.test = 2, 
                    conf.level = 0.95, 
                    method = "unmatched")

round(100*ans_871$power, 1)
## [1] 69.4
round(100*ans_872$power, 1)
## [1] 91.5
round(100*ans_873$power, 1)
## [1] 99

答え＝順に69.4%, 91.5%, 99.0%

8.7.4 ~ 8.7.6

ans_874 <- epi.sscc(OR = 2, 
                    p0 = 0.05, 
                    n = 400+1600, 
                    power = NA, 
                    r = 4, 
                    sided.test = 2, 
                    conf.level = 0.95, 
                    method = "unmatched")
ans_875 <- epi.sscc(OR = 2, 
                    p0 = 0.1, 
                    n = 400+1600, 
                    power = NA, 
                    r = 4, 
                    sided.test = 2, 
                    conf.level = 0.95, 
                    method = "unmatched")
ans_876 <- epi.sscc(OR = 2, 
                    p0 = 0.2, 
                    n = 400+1600, 
                    power = NA, 
                    r = 4, 
                    sided.test = 2, 
                    conf.level = 0.95, 
                    method = "unmatched")

round(100*ans_874$power, 1)
## [1] 86.5
round(100*ans_875$power, 1)
## [1] 98.9
round(100*ans_876$power, 1)
## [1] 100

答え＝順に86.5%, 98.9%, 100%

8.7.7

結果から、研究のデザインについてどのような推奨事項があるか？

（解答例）
・ソース集団の曝露頻度が高い方が検出力は上がる。
・コントロールが多い方が検出力は上がるが、r=4~5以上はあまり変化ない。本文Figure 6.2を参照のこと。

8.8

Schlesselman (1974)は、新生児の円錐動脈幹奇形と母親の経口避妊薬の使用に関する、非マッチの症例対照研究を行った。
全女性の30%が経口避妊薬を使用していると仮定し、2つのグループで同数の5%の両側検定を行うと仮定する。パラメータの値は本文の例8.17と同じなので、90%の検出力でオッズ比(近似相対リスク)2を検出するには、376人の被験者(症例188人、対照188人)が必要であることを示している。

8.8.1

4倍の被験者（つまり1504人）を用いて、90%の検出力で検出できると期待できる、1以上の最小のオッズ比を求めよ。

ans_881 <- epi.sscc(OR = NA, 
                    p0 = 0.3, 
                    n = 1504, 
                    power = 0.90, 
                    r = 1, 
                    sided.test = 2, 
                    conf.level = 0.95, 
                    method = "unmatched")
ans_881$OR
## [1] 0.682810 1.425049

答え＝1.43

8.8.2

被験者の総数を1504人とし、検出可能な最小のオッズ比を1.30から1.50の間で変化させて、それに対応する検出力をプロットせよ。

OR <- seq(from = 1.3, to = 1.5, length = 100)
dat <- expand.grid(n.total = 1504, p0 = 0.3, OR = OR)

for(i in 1:nrow(dat)){
  power_tmp <- epi.sscc(OR = dat$OR[i], 
                        p0 = dat$p0[i], 
                        n = dat$n.total[i],
                        power = NA, 
                        r = 1, 
                        sided.test = 2,
                        conf.level = 0.95, 
                        method = "unmatched")
  dat$power[i] <- power_tmp$power
}

fig <- ggplot(data = dat, aes(x = OR, y = power)) +
  geom_point() +
  scale_x_continuous(limits = c(1.30, 1.50)) +
  scale_y_continuous(limits = c(0.6, 1.0)) + 
  labs(x = "Odd Ratio", 
       y = "Power",
       title = "Relationship between Odd Ratio and Power") +
  theme_bw()

fig


おわり。