6.6アプリケーション
1)表6.1は、スクリーニングツールとして運動ストレステスト(EST)測定を使用する利点を決定することを目的としたコホート研究の結果を示している。 このような検査では、冠状動脈疾患（ＣＡＤ）の診断中に合否結果を得ることができる。
テーブル6.1。 Screening test for physical stress
| CAD | ||||
|---|---|---|---|---|
| Healthy | Ill | Total | ||
| EST | Negative | 327 | 208 | 535 |
| Positive | 115 | 815 | 930 | |
| Total | 442 | 1 023 | 1 465 | |
検証バイアスがないと仮定されます。
a)
この混同行列に基づいて、感度と特異度、正と負の予測値（95％の信頼区間）を示します。
b)
現在のデータの曲線の下の領域の値は何ですか？データテーブルの作成は、以下のように、コマンドmatrix()で作成されたテーブルから行うことができます。
>tab<-as。table(matrix(c(815,115,208,327),nrow=2,byrow=TRUE,
dimnames=list(EST=c),nrow=2,byrow=TRUE,
dimnames=list(EST=c),nrow=2,byrow=TRUE(“+”,”−”),
CAD=c(“+”,”-“))))
>タブ
テーブルの最初の行（通常は暴露）と最初の列（通常は病気）に”陽性”イベントがあり、通常の表記に対応するようにテーブルが再構Rを使用した非常に少ない算術演算を使用して、必要なすべての数量を計算することが可能になります。
しかし、パッケージepiRには、2×2表の疫学的質問に対応するために必要なすべてのツールが含まれています。 したがって、コマンドepi。tests()は、有病率、感度/特異度値、および正および負の予測値を提供します。
>ライブラリ（epiR）
>epi。テスト（タブ）
曲線の下の領域については、外部コマンドrocを使用することもできます。から。パッケージepicalc内のテーブル（）：
>ライブラリ（epicalc）
>roc。から。表(tab,graph=FALSE)
2)表6.2は、40歳から59歳の男性で観察され、血圧のレベルとコレステロールの速度が測定された心筋梗塞の割合をまとめたもので、順序付けられたクラスの形で考慮されています。
表6.2。 /Th>
gt;284
The data are available in the file hdis.それぞれ、血圧（8つのカテゴリ、評価1-8）、コレステロールの割合（7つのカテゴリ、評価1-7）、心筋梗塞の数と個人の総数を示す4つの列を持つテーブルの形 血圧と心筋梗塞に罹患する可能性との関連は、関心のある点である：
a）
血圧の各レベルの心筋梗塞の割合を計算し、それらを表およびグラフ形式で表;
c)
ロジスティック回帰モデルに基づいて、血圧の閾値α=0.05で有意な関連があるかどうかを判断し、クラスセンターと心臓発作の可能性を考慮して定量的変数として扱われる。
d)
ロジット単位で、各血圧レベルのモデルによって予測される梗塞の確率を表現する。
e)
同じグラフに経験的割合と血圧値（クラスセンター）に応じたロジスティック回帰曲線を表示する。
データhdisをインポートします。datは以下のように達成されます：
>bp<-読み取り。テーブル(“hdis.dat”,header=TRUE)
>str(bp)
検証される可能性があるため、関心のある変数のレベルに関連付けられたラベルはありません（血圧 ラベルを生成して関連付けるには、次のコマンドを挿入できます:
> blab < – c(“<117”,”117-126”,”127-136”,”137-146”,
“147-156”,”157-166”,”167-186”,”>186”)
> clab < – c(“<200”,”200-209”,”210-219”,”220-244”,
“245-259”,”260-284”,”>284”)
> bp$bpress < – factor(bp$bpress, labels = blab)
> bp$chol <-factor(bp$chol,labels=clab)
最後のコマンドは、元の変数を質的変数に変換し、blabとclabで定義されたラベルをモダリティに関連付けます。 データベースが目的の形式になったことを確認するには、次のコマンドが役立ちます。
>str(bp)
>summary(bp)
文で公開され:p>
>round(xtabs(hdis/total bpress+chol,data=bp),2)
血圧と心筋梗塞の関係のみに焦点を当てるので、コレステロール値に関するデータを集計する必要があります。 言い換えれば、コレステロールレベルに関係なく、各血圧レベルの値を蓄積する必要があります。 また、クラス間隔の中心を使用して変数bpressのレベルの名前を変更する機会を取る必要があります。
>blab2<-c(111.5,121.5,131。5,141.5,151.5,161.5,176.5,191.5)
>bp$bpress<-rep(blab2,each=7)
>dfrm<-aggregate(bp,
list(bpress=bp),sum)
これは最後のコマンドであるaggregate()であり、データを集計することができます。chol変数のすべての値（sum）を合計します。chol分析のために保持される変数のリ 同時に、結果はdfrmという名前の新しいデータベースに格納されます。 この場合、式表記ではなく、aggregate()の代替構文（応答変数の後に1つ以上の説明変数からなるリストが続く）を使用することに注意してください。 集計されたデータの概要は、入力としてhead()を与えることによって得ることができます。
>head(dfrm,5)
比率pが利用可能な場合、単位がlogitsであるスケール上の値は、関係log(p/(1−p))によって与えられます。 そこから、次のコマンドを導出して、logit単位のhdis/totalのように計算された比率を変換することができます:
><機能(x)のlog(x/(1-x))
>dfrm$プロップ<-dfrm$hdis/dfrm$計
>dfrm$ロジット<-ロジット(dfrm$hdis/dfrm$合計)
で観察することに小さな機能で定義された変換値x、こちらを仮定しているのでプロポーションでそれに相当するlog(x/(1−x). 同様に、一つは書くことができます:
>log((dfrm$hdis/dfrm$total)/(1-dfrm$hdis/dfrm$total))
計算された値は、propとlogitという名前でデータフレームに直接追加されました。
ロジスティック回帰モデルは次のように書かれています。
>glm(cbind(hdis,total-hdis)bpress,data=dfrm,family=binomial)
使用されている定式化cbind(hdis,total-hdis)bpressは、個々の応答ではなくグループ化されたデータにアクセスしているという事実を考慮に入れています。 オプションfamily=binomialを指定したコマンドglm()はロジスティック回帰に対応し、詳細にはあまり触れずに、デフォルトでlogit関数を正規リンクとして使用します。
回帰係数は、summary()を使用して表示することができます。
>summary(glm(cbind(hdis,total-hdis)bpress,data=dfrm,
family=binomial))
前の結果には、血圧と心臓発作の可能性との関連の統計的有意性の質問に答えるために不可欠な情報が含まれています: (log-oddsスケールの)回帰係数は0.024に等しく、通常のしきい値5%で有意です(列Pr(>|z|)を参照)。
検討中の血圧の異なるレベルに応じて心臓発作を起こす可能性は、次のように得られます。
>glm。res<-glm(cbind(hdis,total-hdis)bpress,data=dfrm,
ファミリー=二項式)
>predict(glm.rによって生成された中間結果はglmという名前で格納されていることに注意する必要があります。コマンドpredict()を使用する前にres。 Rによって生成された予測はロジット形式で表され、観測されたロジットと予測されたロジットを比較することができます。
>cbind(dfrm,logit.predit=predict(glm.res))
実質的にすべての要素は、血圧のレベルに応じてグラフで観察され予測される心筋梗塞の割合を表すために利用可能である。 Logoddsスケールではなく、比率の形で表現されたモデルの予測は欠落しています。 したがって、予測曲線、すなわち、血圧に応じた心臓発作を有する可能性を、後者を変数ｂｐｒｅｓｓについて観察される８つの値に限定することなく、描くこと 可能な解決策は次のとおりです。
>dfrm prop prop。predit<-predit(glm.res,type=”response”)
>f<-function(x)
1/(1+exp(-(coef(glm.res) + coef(glm.res)*x)))
> xyplot(hdis/total ˜ bpress, data = dfrm, aspect = 1.2, cex = .8,
xlab = “Blood pressure”,ylab = “Infarction likelihood”,
panel = function(x, y, …) {
panel.xyplot(x, y, col = “gray30”, pch = 19, …)
panel.curve(f, lty = 3, col = “gray70”)
panel.points(x, dfrm$prop.predit,col=”gray70″,…)
})
3）ヒトにおけるアルコールの消費とタバコおよび食道癌の消費との関係に焦点を当てた症例対照調査（”Ille and Villaine”研究）。 症例のグループは、食道癌に罹患している200人の患者で構成され、1972年から1974年の間に診断された。 総選挙では、775人の支配者が選挙人名簿から選出された。 表6.図3は、80gを超える消費が危険因子であると考えられることを考慮して、アルコールの毎日の消費量に応じた被験者の全体の分布を示しています。
表6.3。 アルコール消費量と食道癌
| アルコール消費量（g/日） | |||||||||||||||||||||||||||||||||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ≥80 |
|
| アルコール消費量（g/日） |
≥80 |
|
|
|
|
|
|
|
|
|
|
|
|
|
| < |
ケース |
96 |
104 |
200 |
コントロール |
109 |
コントロール |
109 |
コントロール |
109 |
コントロール |
109 |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
コントロール |
05 770 |
975 |
|
|||
オッズ比の値とその95％信頼区間（Woolf法）？ それは相対的なリスクの良い推定値ですか？b)
リスクのある消費者の割合は、ケース間およびコントロール間で同じですか（α=0.05を考慮してください）？
c）
アルコール消費と個人の状態との関連性をテストすることを可能にするロジスティック回帰モデルを構築する。
c）
アルコール消費と個人の状 回帰係数は重要ですか？
d)
回帰分析の結果に基づいて、(b)で計算された観測オッズ比とその信頼区間の値を求めます。
d)
(b)で計算された観測オッズ比の値を求めます。
d)
（個々の）生データは利用できないため、問題文で提供されている頻度テーブルを直接操作する必要があります:
> alcohol < – matrix(c(666,104,109,96), nr = 2,
dimnames = list(c(“Control”,”Case”),
c(“<> = 80”)))
> alcohol
Regarding the odds ratio, the command oddsratio() from the package vcd will be used:
>library(vcd)
>oddsratio(alcohol,log=FALSE)
オプションlog=FALSEは、返された結果がオッズ比のログではなくオッズ比に正しく対応するこ 漸近信頼区間を取得するには、以下を使用します。
>confint(oddsratio(alcohol,log=FALSE))
同様に、summary(oddsratio(alcohol))を使用して、対数オッズ比(H0:log(λ)=0)の仮説検定を実行できます。
コマンドprop。test()は、毎日の消費量を持つ人の割合が≥80gであり、ケースとコントロールで同じであるという仮説を検定するために使用することができ、問題文で与えられたクロス集計から観察された値を示します。
>prop。test(c(96,109),c(200,775),correct=FALSE)
この検定は、データから推定された2つの比率の差を検定するZ検定とまったく同じです(連続性補正が使用されていない場合)。
ロジスティック回帰モデルに関しては、分割表は、二つの質的変数（疾患と曝露）が適切に表現されたデータ表、すなわちデータフレームに変換する必要が これは、パッケージreshape2のコマンドmelt()を使用して行うことができ、表形式のデータを「長い」形式に変換することができます。 この段階では、変数diseaseのレベルは0（コントロール）と1（ケース）で再コード化する必要がありますが、これは実際には必要ではなく、レベル0は参照カテゴリ（結果の解釈を容易にする）として考慮する必要があります。
>library(reshape2)
>アルコール。df<-メルト（アルコール）
>名前（アルコール。df）<-c（”病気”、”露出”、”n”）
>レベル（アルコール。df$disease)<-c(1,0)
>アルコール。df disease disease<-relevel(アルコール。次に、ロジスティック回帰モデルを次の方法で記述します。
>glm。res<-glm（疾患暴露、データ=アルコール。p>
>summary(glm.res)
関心のある結果は、expositionに関連付けられた行です>=80これは、Rによって推定された露出変数に関連付けられた回帰係数の値とその標準誤差、および検定統計量の値に関する情報を提供するためです。 ここで、回帰係数はオッズ比の対数として解釈されます。 前の定式化である暴露疾患における変数の役割を交換することによって、まったく同じ結果が得られることに注意してください。
上記で計算されたオッズ比は、関心因子（露出）に関連付けられた回帰係数とその95％信頼区間から、次の方法で求めることができます。
>exp(coef(glm.res))
>exp(confint(glm.res))
第二の解決策は、ケースの数と個人の総数を考慮することで構成されています：
>alcohol2<-データ。フレーム（expos=c（”<>=80″）、ケース=c（104,96）、
合計=c(104 + 666, 96 + 109))
>summary(glm(cbind(case,total-case)expos,data=alcohol2,
family=binomial))
4)このデータは、心筋梗塞の発生の予防に関する体内のクレアチンキナーゼ濃度の予後的妥当性を確
データはファイルsckで使用できます。dat: 最初の列は可変クレアチンキナーゼ（ck）に対応し、2番目の列は可変疾患の存在（pres）に対応し、最後の列は可変疾患の不在（abs）に対応する。
a）
被験者の総数は何ですか？
b)
相対的な病気/健康な頻度を計算し、散布図(ポイント+ポイントを接続するセグメント)を使用してクレアチンキナーゼ値に従って進化を表します。
c)
病気になる可能性を予測することを目的としたロジスティック回帰モデルに基づいて、このモデルが0.5のしきい値を考慮して人々が病気に罹患していると予測するCKの値を計算する(P(sick)≤0.5の場合、sick=1)。
d)
このモデルによってill予測される確率と、値ckによる経験的比率をグラフィカルに表します。
d)
D)
D)
d)
d)
d)
d)
e）
対応するROC曲線を確立し、曲線の下の領域の値を報告します。変数の名前はデータファイルにないため、データがインポートされた直後に変数を割り当てる必要があります。
>sck<-read。テーブル(“sck.p>
>names(sck)<-c(“ck”,”pres”,”abs”)
>summary(sck)
サブジェクトの総数は、合計に対応します2つの変数presとabsの頻度のうち、それは次のとおりです。
>sum(sck)
または、同等に: sum(sck$pres)+sum(sck$abs)(ただし、sck$pres+sck$absではありません)。
これら二つの変数の相対頻度を計算するには、変数ごとに合計量を知る必要があります。
これらの変数の相対頻度を計算するには、変数によ これらは、コマンドapplyとoperating by columnsを使用して取得できます。
>ni<-apply(sck,2,sum)
これに基づいて、変数presとabsの各値を以前に計算された値で除算すれば十分です。 得られた値は、同じデータテーブル内の2つの新しい変数に格納されます。
>sck$pres。prop<-sck$pres/ni
>sck$abs.prop<-sck$abs/ni
計算が正しいことを確認することができます：各変数の値の合計は1に等しくなければなりません：
>apply(sck,2,sum)
得られた比率は、変数ckの値を次のように考慮して、単純に単一のチャートで表されますp>
>xyplot(pres.プロップ+abs.prop ck,data=sck,type=c(“b”,”g”),
auto.命令type=c(“b”,”g”)は、グリッド（”g”）と同様に、線（”b”=”o”+”l”）で接続されたポイントを表示することを意味します。グループ化されたデータの回帰モデルは、次のように記述されます。
>glm。res<-glm(cbind(pres,abs)ck,data=sck,family=binomial)
>summary(glm.res)
log-oddsスケールではなく確率の形式で表される予測は、オプションtype=”response”を次のように指定することによってpredictコマンドによって取得されます。
>glm。pred<-predict(glm.res,type=”response”)
>names(glm.pred)<-sck≤ck
確率≤0.5が病気の個人を指定することを考慮することにより、以下の分布が得られます。
>glm。pred
このモデルによれば、値ckが80以上の場合、人々は病気とみなされると結論づけられています。 これは、変数ckの値に基づいて予測される確率が表されるチャートで簡単に検証されます。
>sck$sick<-sck$pres/(sck$pres+sck$abs)
>xyplot(glm.pred~sck$ck,type=”l”,
ylab=”Probability”,xlab=”ck”,
panel=function(…){
panel.xyplot(…)
パネル。xyplot(sck$ck,sck$sick,pch=19,col=”grey”)
})
最初のステップでは、グループ化されたデータを”解凍”し、二つの列を持つテーブルを作成する必要があります。 私たちは、以前に計算され、変数niで利用可能なサブグループ別のカウントを使用します：
>sck。展開<-データ。フレーム(ck=c(rep(sck≤ck,sck≤pres),
rep(sck≤ck,sck≤abs)),
sick=c(rep(1,ni),rep(0,ni)))
>テーブル(sck.expand sickを展開します）
>with(sck.expand,tapply(sick,ck,sum))
最後の二つのコマンドは、同じカウントが検出され、ckの値に応じた病気の人の分布が正しいことを確認することを目的としています。 ロジスティック回帰に関しても同じ結果が得られ、同時に前のデータテーブルで予測された値を追加できることも確認できます（以前と同じ手順に従res2<-glm(sick ck,data=sck.p>
>sckを展開します。expand予測を展開<-ifelse(predict(glm.res2,
type=”response”)>=0.5,
1,0)
>with(sck.expand,table(sick,prediction))
最後のコマンドは、実際の診断と回帰モデルによって予測された診断が交差する混同行列を表示することを可能にします。 参照しきい値を変更すると、正しい分類率を比較できます。
>classif。タブ<-with(sck.展開、テーブル（病気、予測））
>sum(diag(classif.タブ))/sum(classif.P>
ROC曲線を表示するには、パッケージROCRを使用します。
>ライブラリ（ROCR）
>pred<-予測（predict（glm。res2,type=”response”),
sck.expand$sick)
> perf < – performance(pred, “tpr”,”fpr”)
> plot(perf, ylab = “Sensitivity”, xlab = “1-Specificity”)
> grid()
> abline(0, 1, lty = 2)
The value of the area under the curve is obtained as follows:
> performance(pred, “auc”)@”y.values”