医学基礎技術演習・実験基本技術（医学統計学）テキスト (2)中澤

港（公衆衛生学准教授）nminato@med.gunma-u.ac.jp

2009 年 6 月 3 日

参考文献の続き• http://socserv.mcmaster.ca/jfox/Misc/Rcmdr/Getting-Started-with-the-Rcmdr.pdf（作成者である John Fox 自身による R Commander の入門テキスト）• http://www.ec.kansai-u.ac.jp/user/arakit/documents/Getting-Started-with-the-Rcmdr-ja.pdf（日本語版メニュー作成者である荒木孝治さんによる邦訳）

1 相関と回帰
相関も回帰も２つの量的な変数間の関係を調べる点は共通である。そのため，まず散布図を描く。

¶
例題 1

³

http://phi.med.gunma-u.ac.jp/msb/data/p01.txt は，男女合わせて 100 人の集団の身長 (HT) と体重 (WT) のデータ
（欠損値を含む）である。身長を横軸，体重を縦軸とした散布図を描け。記号は性別 (SEX) ごとに変えること。

µ´

まず，先週やったように R のアイコンをダブルクリックして R を起動後，プロンプトに library(Rcmdr) として

Rcmdr を起動する。
次 い で デ ー タ を 読 み 込 む 。Rcmdr の メ ニ ュ ー の「 デ ー タ 」か ら「 デ ー タ の イ ン ポ ー ト 」の「 テ キ ス ト フ ァ
イ ル ま た は ク リ ッ プ ボ ー ド か ら 」を 選 ん で ，表 示 さ れ る ダ イ ア ロ グ で「 区 切 り 」を「 タ ブ 」に 変 え ，OK ボ
タ ン を ク リ ッ ク し て ，フ ァ イ ル を 選 択 す る ウ ィ ン ド ウ が 出 て き た ら ，フ ァ イ ル 名 を 入 力 す る 枠 に デ ー タ のURL(http://phi.med.gunma-u.ac.jp/msb/data/p01.txt) を打って OK するとネットワーク経由でデータファ
イルをデータフレームとして読み込むことができる。とくに変えていなければ，データフレーム名は Dataset となっているはずである（ここまでは前回の復習）。
次いで，散布図を描く。Rcmdr のメニューの「グラフ」から「散布図」を選ぶ。表示されるウィンドウの中で，「x変数」として HT を選び，「y 変数」として WT を選ぶ。下の方は，周辺箱ヒゲ図と最小２乗直線と平滑線の右側のボックスにチェックが入っているが，相関をみる場合は最小２乗直線のチェックは外す（平滑線もない方がいい）。周辺箱ヒゲ図は横軸の変数，縦軸の変数別々に箱ヒゲ図を描いてくれるので，チェックが入ったままでよい。
この例題では性別にプロット記号を変えることとなっているので，下の方の「層別のプロット」というボタンをクリックして，出てくるウィンドウの中で層別変数として SEX を選ぶ。その下の層別して線を描くというボックスにチェックが入っているが，最小２乗直線のチェックを外してあれば，このボックスの指定は無効である。後は OK ボタ
ンをクリックしていけば，次の散布図ができる。

1

1.1 相関と回帰の違い
大雑把に言えば，相関が変数間の関連の強さを表すのに対して，回帰はある変数の値のばらつきがどの程度他の変数の値のばらつきによって説明されるかを示す。回帰の際に，説明される変数を（従属変数または）目的変数，説明するための変数を（独立変数または）説明変数と呼ぶ。２つの変数間の関係を予測に使うためには，回帰を用いる。

1.2 相関関係とは一般に，2 個以上の変量が「かなりの程度の規則正しさをもって，増減をともにする関係」のことを相関関係(correlation) という。相関には正の相関 (positive correlation) と負の相関 (negative correlation) があり，一方が増えれば他方も増える場合を正の相関，一方が増えると他方は減る場合を負の相関と呼ぶ。例えば，身長と体重の関係は正の相関である。

1.3 見かけの相関，擬似相関
相関関係があっても，それが見かけ上の関係に過ぎない場合がある。具体例としては，血圧と所得の間に正の相関があるという命題は，データをとってみれば，多くの場合に成り立つであろう。しかし，おそらくどちらも年齢や摂取エ
ネルギー量との間に真の相関関係があって，それらの影響を制御したら（例えば同年齢で同じような食生活をしている
人だけについて見る，という限定をしたら），相関関係は消えてしまうだろう。この場合，見かけ上の相関があることは科学的仮説としての意味に乏しい。
時系列データや地域相関のデータでは，擬似相関 (spurious correalation) が見られる場合もある。例えば，ある年に日本で植えた木の幹の太さと，同じ年にイギリスで生れた少年の身長を 15 年分，毎年１回測ったデータをプロットすると，おそらくは正の相関関係があるように見えるのだが，両者の間に直接関係がないのは明らかである（どちらも時間が経つにつれて大きくなっているだけである）。

1.4 直線的な相関，直線に載らない相関
相関関係は増減をともにすればいいので，直線的な関係である必要はなく，二次式でも指数関数でもシグモイドでもよいが，通常，直線的な関係をいうことが多い（指標はピアソンの積率相関係数）。曲線的な関係の場合，直線的になるように変換したり，順位の情報だけを使った相関の指標（順位相関係数）を計算する。
普通，ただ相関係数といえば，ピアソンの積率相関係数 (Pearson’s Product Moment Correlation Coefficient) を指し，r という記号で表すが，この r は直線的な関係の強さの指標である。X と Y の共分散を X の分散と Y の分散の2

積の平方根で割った値であり，範囲は [−1, 1] である。最も強い負の相関があるとき r = −1，最も強い正の相関があるとき r = 1，まったく相関がないとき（2 つの変数が独立なとき），r = 0 となることが期待される。X の平均を X̄ ，Yの平均を Ȳ と書けば，

Pn

i=1 (Xi − X̄)(Yi − Ȳ )Pnn22
i=1 (Xi − X̄)
i=1 (Yi − Ȳ )r = qPである。

相関係数の有意性の検定においては，母相関係数がゼロ（＝相関が無い）という帰無仮説の下で，実際に得られている相関係数よりも絶対値が大きな相関係数が偶然得られる確率（これを「有意確率」という。通常，記号 p で表すので，「p 値」とも呼ばれる）の値を調べる。偶然ではありえないほど珍しいことが起こったと考えて，帰無仮説が間違っていたと判断するのは有意確率がいくつ以下のときか，という水準を有意水準といい，検定の際には予め有意水準を（例えば 5% と）決めておく必要がある。例えば p = 0.034 であれば，有意水準 5% で有意な相関があるという意思決定を行なうことができる。p 値は，検定統計量√
r n−2
t0 = √
1 − r2が自由度 n − 2 の t 分布に従うことを利用して求めるられる。

¶
例題 2

³

例題 1 のデータで身長と体重のピアソンの積率相関係数を計算し，有意性を検定せよ。

µ¤¡
Rcmdr では，「統計量」の「要約」の「相関の検定」を選び，変数として WT と HT を選ぶ（£Ctrl ¢キーを押しなが´

ら変数名をクリックすれば複数選べる）。相関のタイプとして「ピアソンの積率相関」と「スピアマンの順位」と「ケンドールのタウ」が選べるようになっている。この例題ではピアソンの積率相関係数を求めるので，初期設定のまま「ピアソンの積率相関」にしておけばよい。検定についても「対立仮説」の下に「両側」「相関＜０」「相関＞０」の３つから選べるようになっているが，通常は「両側」でよい。OK をクリックすると，Rcmdr の出力ウィンドウに次の内容が表示される。

¶³

Pearson’s product-moment correlation
data: Dataset$HT and Dataset$WT
t = 16.4519, df = 95, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.7977988 0.9045751
sample estimates:
cor
0.860348µ´

これより，身長と体重の関係について求めたピアソンの積率相関係数は，r = 0.86（95% 信頼区間が [0.798, 0.905]）であり，p-value < 2.2e-16（有意確率が 2.2 × 10−16 より小さいという意味）より，「相関が無い」可能性はほとん
どゼロなので，有意な相関があるといえる。なお，相関の強さは相関係数の絶対値の大きさによって判定し，伝統的に0.7 より大きければ「強い相関」，0.4〜0.7 で「中程度の相関」，0.2〜0.4 で「弱い相関」とみなすのが目安である。
せっかく男女別にプロットしたので，相関係数の検定も男女別に実行したいところだが，残念ながらボタン１つというわけにはいかない。男女別に相関係数の検定を実行するには，「データ」の「アクティブデータセット」の「アクティブデータセットの部分集合を抽出」を使って男女別のデータフレームを作成しなくてはならない。表示されるウィンドウで，「すべての変数を含む」はチェックが入ったままでよく，「部分集合の表現」のボックスに SEX=="M"と入力
し，「新しいデータセットの名前」に Males（既にある名前と重複しなければ何でもよい）と入力して OK ボタンをクリックすると男性だけのデータフレーム Males ができてアクティブになる。ここで先ほどと同じ「統計量」「要約」「相関の検定」をすれば男性の身長と体重についてピアソンの積率相関係数を求めて有意性の検定をすることができる。

3

¶

³

Pearson’s product-moment correlation
data: Males$HT and Males$WT
t = 8.636, df = 46, p-value = 3.476e-11
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6468709 0.8750522
sample estimates:cor
0.7864562µ´

女性について同じことをするには，まず「データ」の「アクティブデータセット」の「アクティブデータセットの選択」で Dataset を選び直し，「アクティブデータセットの部分集合を抽出」の「部分集合の表現」で SEX=="F"，「新しいデータセットの名前」で Females として OK ボタンをクリックしてから「統計量」→「要約」→「相関の検定」を実行すればよい。

¶³

Pearson’s product-moment correlation
data: Females$HT and Females$WT
t = 7.1667, df = 47, p-value = 4.569e-09
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5539837 0.8342857
sample estimates:cor
0.7226128µ´

以上より，身長と体重の相関係数は，男性で 0.786（95% 信頼区間が [0.647,0.875]），女性で 0.723（95% 信頼区間が [0.554,0.834]）とわかった。男女とも統計学的に有意な強い正の相関があったといえる。
このデータでは必要なかったが，相関関係が直線的でなかったり，外れ値があったりする場合は，順位相関係数を使うのが適切な場合もある。

¶
例題 3

³

組み込みデータ airquality は，1973 年 5 月 1 日から 9 月 30 日まで 154 日間のニューヨーク市の大気環境データである。
含まれている変数は，Ozone（ppb 単位でのオゾン濃度），Solar.R（セントラルパークでの 8:00 から 12:00 までの 4000 から 7700 オングストロームの周波数帯の太陽放射の強さを Langley 単位で表した値），Wind（LaGuardia 空港での 7:00 から10:00 までの平均風速，マイル／時），Temp（華氏での日最高気温），Month（月），Day（日）である。太陽放射の強さとオゾン濃度の相関関係を検討せよ。

µ´

まずデータをアクティブにして散布図を描くのはいつも同じである。「データ」の「パッケージ内のデータ」の「ア
タッチされたパッケージからデータセットを読み込む」を選び，「データセット名を入力」のボックスに airqualityと打って OK ボタンをクリックしてから，「グラフ」の「散布図」で「x 変数」として「Solar.R」を選び，「y 変数」として「Ozone」を選び，「最小２乗直線」と「平滑線」のチェックを外してから OK ボタンをクリックする。

4

どう見ても直線的な関係とは言いがたいので，スピアマンの順位相関係数を計算して，その有意性の検定をしてみ
る。
「統計量」
「要約」
「相関の検定」で変数として Solar.R と Ozone を選び，相関のタイプを「スピアマンの順位」にして OK ボタンをクリックすると，次の結果が得られる。弱いけれども有意な正の相関があるといえる。

¶³

Spearman’s rank correlation rho
data: airquality$Ozone and airquality$Solar.R
S = 148561.3, p-value = 0.0001806
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.3481865µ´

もっとも，このデータの場合はピアソンの積率相関係数でも似たような結果が得られ（各自確かめよ），直線的な相関でないことの影響はあまりクリアでない。

¶
順位相関係数の定義³

なお，スピアマンの順位相関係数 ρ はa ，値を順位で置き換えた（同順位には平均順位を与えた）ピアソンの積率相関係数と同じである。Xi の順位を Ri ，Yi の順位を Qi とかけば，

ρ=1−

n
X

6
n(n2 − 1)(Ri − Qi )2i=1となる。スピアマンの順位相関係数がゼロと差がないことを帰無仮説とする両側検定は，サンプル数が 10 以上ならばピアソンの場合と同様に，

√
ρ n−2
T = p
1 − ρ2が自由度 n − 2 の t 分布に従うことを利用して行うことができる。ケンドールの順位相関係数 τ は，

τ =(A − B)n(n − 1)/2

によって得られる。ここで A は順位の大小関係が一致する組の数，B は不一致数である。
a ピアソンの相関係数の母相関係数を ρ と書き，スピアマンの順位相関係数を r と書く流儀もある。
s

µ

´

5

1.5 回帰モデルの数理
既に述べたとおり，回帰は，従属変数のばらつきを独立変数のばらつきで説明するというモデルの当てはめである。
十分な説明ができるモデルであれば，そのモデルに独立変数の値を代入することによって，対応する従属変数の値が予測あるいは推定できるし，従属変数の値を代入すると，対応する独立変数の値が逆算できる。こうした回帰モデルの実用例の最たるものが検量線である。検量線とは，実験において予め濃度がわかっている標準物質を測ったときの吸光度のばらつきが，その濃度によってほぼ完全に（通常 98% 以上）説明されるときに（そういう場合は，散布図を描くと，点々がだいたい直線上に乗るように見える），その関係を利用して，サンプルを測ったときの吸光度からサンプルの濃度を逆算するための回帰直線である（曲線の場合もあるが，通常は何らかの変換をほどこし，線形回帰にして利用する）。
検量線の計算には，(A) 試薬ブランクでゼロ点調整をした場合の原点を通る回帰直線を用いる場合と，(B) 純水でゼロ点調整をした場合の切片のある回帰直線を用いる場合がある。いずれも，量がわかっているもの（この場合は濃度）を x，誤差を含んでいる可能性がある測定値（この場合は吸光度）を y として y = bx + a という形の回帰式の係数 aと b（a を切片，b を回帰係数と呼ぶ）を最小二乗法で推定し，サンプルを測定した値 y から x = (y − a)/b によってサンプルの濃度 x を求める。ただし，データ点の最小，最大より外で直線関係が成立する保証はない。従って，サンプル
測定値が標準物質の測定値の最小より低いか，最大より高いときは，限界を超えていることになってしまう*1 。

1.6 回帰モデルの当てはまりデータから得た回帰直線は，完璧にデータに乗ることはない。そこで，回帰直線の当てはまりのよさを評価する。aと b が決まったとして，zi = a + bxi とおいたとき，ei = yi − zi を残差 (residual) と呼ぶ。残差は，yi のばらつきのうち，回帰直線では説明できなかった残りに該当する。つまり，残差が大きいほど，回帰直線の当てはまりは悪いと考えられる。残差にはプラスもマイナスもあるので二乗和をとり，

Q=

n
X

e2i =

n
X

i=1

=

n
X
i=1

n
X

yi2 − (yi )2 /n −(yi − zi )2

i=1

(n

Pn
Pn
x
yi )2
xy −
i=1
Pni i 2 i=1
Pn i i=1/n2n i=1 xi − ( i=1 xi )

Pn

i=1

として得られる Q は，回帰直線の当てはまりの悪さを示す尺度となる。
この Q を「残差平方和」と呼び，それを n で割った Q/n を残差分散という。残差分散（var(e) と書くことにする）と Y の分散 var(Y ) とピアソンの相関係数 r の間には，var(e) = var(Y )(1 − r2 ) という関係が常に成り立つので，

r2 = 1 − var(e)/var(Y ) となる。このことから r2 が 1 に近いほど回帰直線の当てはまりがよいことになる。その意味
で，r2 を「決定係数」と呼ぶ。また，決定係数は，Y のばらつきがどの程度 X のばらつきによって説明されるかを意味するので，X の「寄与率」と呼ぶこともある。
データによっては，何通りもの回帰直線の残差平方和が大差ないという状況がありうる。例えば，目的変数と説明変数が実はまったく無関係であった場合は，データの重心を通るどのような傾きの線を引いても残差平方和はほとんど同じになってしまう。言い換えれば，傾きや切片の推定値が不安定になる。
回帰直線のパラメータ（回帰係数 b と切片 a）の推定値の安定性を評価するためには，t 値が使われる。いま，Y とX の関係が Y = a0 + b0 X + e というモデルで表されるとして，誤差項 e が平均 0，分散 σ 2 の正規分布に従うものとすれば，回帰係数の推定値 a も，平均 a0 ，分散 (σ 2 /n)(1 + M 2 /V )（ただし M と V は x の平均と分散）の正規分布に従い，残差平方和 Q を誤差分散 σ 2 で割った Q/σ 2 が自由度 (n − 2) のカイ二乗分布に従うことから，

pn(n − 2)(a − a0 )t0 (a0 ) = p(1 + M 2 /V )Q
が自由度 (n − 2) の t 分布に従うことになる。
*1 このような場合はサンプルを希釈するか濃縮して測定するのが普通である。

6

しかしこの値は a0 がわからないと計算できない。a0 が 0 に近ければこの式で a0 = 0 と置いた値（つまり t0 (0)。これを切片に関する t 値と呼ぶ）を観測データから計算した値が t0 (a0 ) とほぼ一致し，自由度 (n − 2) の t 分布に従うはずなので，その絶対値は 95% の確率で t 分布の 97.5% 点（サンプルサイズが大きければ約 2 である）よりも小さくなる。つまり，データから計算された t 値がそれより大きければ，切片は 0 でない可能性が高いことになるし，t 分布の分布関数を使えば，「切片が 0 である」という帰無仮説に対する有意確率が計算できる。
回帰係数についても同様に，

pt0 (b) =n(n − 2)V b√
Qが自由度 (n − 2) の t 分布に従うことを利用して，「回帰係数が 0」であるという帰無仮説に対する有意確率が計算できる。有意確率が充分小さければ，切片や回帰係数がゼロでない何かの値をとるといえるので，これらの推定値は安定していることになる。

1.7 回帰分析の向きと予測の際の留意点
身長と体重のように，どちらも誤差を含んでいる可能性がある測定値である場合には，一方を説明変数，他方を目的変数とすることは妥当でないかもしれない（一般には，身長によって体重が決まるなど方向性が仮定できれば，身長を説明変数にしてもよいことになっている）。また，最小二乗推定の説明から自明なように，回帰式の両辺を入れ替えた
回帰直線は一致しない。従って，どちらを目的変数とみなし，どちらを説明変数とみなすか，因果関係の方向性に基づいて（先行研究や臨床的知見を参照し）きちんと決めるべきである。
回帰を使って予測をするとき，外挿には注意が必要である。とくに検量線は外挿してはいけない。実際に測った濃度より濃かったり薄かったりするサンプルに対して，同じ関係が成り立つという保証はどこにもないからである（吸光度を y とする場合は，濃度が高くなると分子の重なりが増えるので飽和 (saturate) してしまい，吸光度の相対的な上がり方が小さくなっていき，直線から外れていく）。サンプルを希釈したり濃縮したりして，検量線の範囲内で定量しなくてはいけない。

¶
例題 4

³

例題 3 のニューヨーク大気環境データについて，日照の強さを説明変数，オゾン濃度を目的変数として回帰分析せよ。

µ´

既に散布図は描いたが，回帰分析の場合は最小 2 乗直線も描くのが普通なのでやり直す。次に，「統計量」の「モデルへの適合」の「線形回帰」を選ぶ。目的変数として Ozone を，説明変数として Solar.R を選んで OK ボタンをクリッ
クすると，「出力ウィンドウ」に次の結果が得られる。

¶

³

Call:
lm(formula = Ozone ~ Solar.R, data = airquality)
Residuals:
Min
1Q-48.292 -21.361

Median-8.864

3Q
Max
16.373 119.136

Coefficients:
Estimate Std. Error t value Pr(>|t|)(Intercept) 18.59873
6.747902.756 0.006856 **Solar.R
0.127170.032783.880 0.000179 ***--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Residual standard error: 31.33 on 109 degrees of freedom(42 observations deleted due to missingness)
Multiple R-Squared: 0.1213,Adjusted R-squared: 0.1133
F-statistic: 15.05 on 1 and 109 DF, p-value: 0.0001793µ´

得られた回帰式は Ozone = 18.599 + 0.127 · Solar.R であり，最下行をみると F 検定の結果の p 値が 0.0001793 と7きわめて小さいので，モデルの当てはまりは有意である。しかし，その上の行の Adjusted R-squared の値が 0.11 ということは，このモデルではオゾン濃度のばらつきの 10% 余りしか説明されないことになり，あまりいい回帰モデルではない。

1.8 共分散分析
複数のグループがあって，どのグループに属するサンプルについても，同じ説明変数と目的変数が調べられているとき，それらの関係がグループによって異なるかどうか調べたい場合がある。共分散分析は，このような場合に用いられる。
典型的には，Y = β0 + β1 C + β2 X + β12 CX + ε というモデルになる。２値変数 C によって示される２群間で，量的変数 Y の平均値に差があるかどうかを比べるのだが，Y が量的変数 X と相関がある場合に（このとき X を共変量と呼ぶ），X と Y の回帰直線の傾き (slope) が C の示す２群間で差がないときに，X による影響を調整した Y の修正平均（adjusted mean; 調整平均ともいう）を，C の２群間で比べる。修正平均は C の各変数についての係数（２群の場合，基準にする変数の係数はゼロ）に，共変量の平均に共変量の係数を掛けたものを加え，さらに切片を加えることによって計算できる。
ただし，２本の回帰直線がともに充分な説明力をもっていて，かつ２本の回帰直線の間で傾きに差がない場合でないと，修正平均の比較には意味がない。そもそも回帰直線の説明力が低ければその変数は共変量として考慮する必要がな
いし，傾きが違っていれば群分け変数と独立変数の交互作用が従属変数に関して有意に影響しているということなので，２群を層別して別々に解釈する方が良い。
いま，C で群分けされる２つの母集団における，(X, Y) の間の母回帰直線を，y = α1 + β1 x，y = α2 + β2 x とすれ
ば，共分散分析は次の手順で進める。

H0 : β1 = β2 ，H1 : β1 6= β2 を検定する。各群について，X と Y の平均と変傾きに差がないという帰無仮説の検定

動と共変動を出しておけば*2 ，仮説 H1 のもとでの残差平方和

d1 = SSY 1 − (SSXY 1 )2 /SSX1 + SSY 2 − (SSXY 2 )2 /SSX2と仮説 H0 のもとでの残差平方和

d2 = SSY 1 + SSY 2 − (SSXY 1 + SSXY 2 )2 /(SSX1 + SSX2 )を計算して F = (d2 − d1 )/(d1 /(N − 4)) が H0 のもとで第 1 自由度 1，第 2 自由度 N − 4 の F 分布に従うことを使って傾きが等しいかどうかの検定ができる。
傾きに差がないとき，y 切片に差がない帰無仮説の検定 β1 = β2 のもとで（即ち，共通の傾き β を，β = (SSXY 1 +SSXY 2 )/(SSX1 + SSX2 ) として推定し），H00 : α1 = α2 ，H10 : α1 6= α2 を検定する。帰無仮説 H00 のもとで全部のデータを使った残差平方和 d3 = SSY − (SSXY )2 /SSX を計算して，F = (d3 − d2 )/(d2 /(N − 3)) が第 1
自由度 1，第 2 自由度 N − 3 の F 分布に従うことを使って検定できる。H00 が棄却された場合は各群の平均を共通の傾きに代入すれば各群の切片が求められるし，棄却されない場合は２群を一緒にして普通の単回帰分析をすることになる。
傾きに有意差があるとき，層別解析

β1 = SSXY 1 /SSX1 ，β2 = SSXY 2 /SSX2 として別々に傾きを推定し，y 切片 α

もそれぞれの式に各群の平均値を入れて計算する。

¶
例題 5

³

組み込みデータ swiss（1888 年頃のスイスのフランス語を話す 47 州についての，標準化された出生力水準 Fertility，農業
就業割合 Agriculture，陸軍の試験で最高ランクを記録した人の割合 Examination，初等教育を超える教育を受けた人の割合 Education，カソリック信者割合 Catholic，乳児死亡割合 Infant.Mortality からなるデータ）を使って，教育水準が高いほど出生力が低いけれども，それがカソリック信者割合に影響を受ける（カソリック信者の方がプロテスタント信者よりも一般に出生力が高い）という仮説を検討してみよう。

µ´*2 サンプルサイズ N 1 の第１群に属する xi , yi について，EX1 =P(yi − EY 1 )2 ，EXY 1 =P

xi yi /N 1，SSXY 1 =

P

P

P

xi /N 1，SSX1 =(xi − EX1 )2 ，EY 1 =(xi yi − EXY 1 )2 。第２群も同様。

8P

yi /N 1，SSY 1 =「データ」→「パッケージ内のデータ」→「アタッチされたパッケージからデータセットを読み込む」として，パッケージとして datasets，データとして swiss を選択する。次に，「データ」→「アクティブデータセット内の変数の管理」→「新しい変数を計算」として，「新しい変数名」を MoreCatholic，計算式を Catholic>=50 とすれば，カソリック信者割合が 50% 以上の州で TRUE，そうでない州で FALSE となる新しい論理型の変数 MoreCatholic ができる。しかし，Rcmdr では論理型の変数がカテゴリ変数とみなされないので，さらに「データ」→「アクティブデータ
セット内の変数の管理」→「変数の再コード化」で「再コード化の変数」として MoreCatholic を選び，「New variablename or prefix for multiple recodes」に CatholicShare，計算式を FALSE="Less";TRUE="More"として OK すると，やっと層別に使える因子型の変数 CatholicShare が得られる。
¶³（注）このプロセスは，RConsole では 1 行でできる。次のどちらかで良い（上の場合は水準名が”Less”と”More”ではな
く，”[0,50)”と”[50,101)”になるが）。

swiss$CatholicShare <- cut(swiss$Catholic,c(0,50,101),right=F)
swiss$CatholicShare <- as.factor(ifelse(swiss$Catholic>=50,"More","Less"))µ´

次に CatholicShare で層別して散布図を描かせ，最小 2 乗直線も層別に描かせる。つまり，「グラフ」→「散布図」
から，「x 変数」として Education，「y 変数」として Fertility を選び，平滑線のチェックを外し，「層別のプロット」をクリックして CatholicShare を選んでから元のウィンドウで OK すれば目的のグラフが得られる。グラフを見ると最小 2 乗直線が交差しているので交互作用がありそうにも見える。
次に層別に回帰分析をして線型回帰モデルが有意に当てはまるか調べる。「統計量」→「モデルへの適合」→「線形回帰」で「目的変数」として Fertility，「説明変数」として Education を選び，「部分集合の表現」としてCatholicShare=="Less"として OK したときと CatholicShare=="More"として OK したときの出力を両方見ると，どちらでも Education の係数は 5% 水準で有意であることがわかる。
そこで，「統計量」→「モデルへの適合」→「線型モデル」で，モデルとして（R を起動後初めてモデルへの適合を呼び出す場合，モデル名 LinearModel.1 となっている）左辺に Fertility を，右辺に Education*CatholicShare を指定すれば交互作用項（Education:CatholicShare）により傾きの差を検討することができる。この場合，交互作用項
が 5% 水準で有意なので，カソリック信者の多い郡と少ない郡の間では，教育水準と出生力の関係が異なっているといえ，層別回帰の結果を採用することになる。
ちなみに，もし傾きの差が有意でなければ，もう一度「線型モデル」を呼び出して，モデル名 LinearModel.2 とし，右辺を Education+CatholicShare として，CatholicShare の係数がゼロと有意差があるかどうかをみれば，カソリックが多い郡と少ない郡の間で教育水準の影響を調整した出生力の修正平均に差があるかどうかがわかる。修正平均
そのものは Rcmdr では得られないので，R Console で以下を打つ。

¶³

cfs <- dummy.coef(LinearModel.2)
cfs[[1]] + cfs$Education * mean(swiss$Education) + cfs$CatholicShare

µ
¶
例題 6

´
³

http://phi.med.gunma-u.ac.jp/grad/sample3.dat は，都道府県別のタブ区切りテキストデータファイルである。変数としては，都道府県名 (PREF)，日本の東西（REGION），1990 年の 100 世帯当たり乗用車台数（CAR1990)，1989 年の人
口 10 万人当たり交通事故死者数 (TA1989)，1985 年の国勢調査による人口集中地区居住割合 (DIDP1985) が含まれている（REGION の 1 は東日本，2 は西日本を意味する）。
このデータについて，東日本と西日本で，100 世帯当たり乗用車台数で調整した人口 10 万人当たり交通事故死者数に差がある
か，共分散分析によって検討せよa 。
a （注）実は乗用車台数の影響を調整しなければ人口当たり交通事故死者数は東西で有意な差はない。

µ´

データセット名 sample3 として web からデータを読み込む。次に REGION で層別した散布図を描く。「x 変数」を

CAR1990，「y 変数」を TA1989 とすると，2 本の回帰直線はほぼ平行に見える。次に東西日本別々に層別して，CAR1990
によって TA1989 が説明されるかをみるため，単回帰分析を行う（やり方は例題 5 と同じなので省略する）。CAR1990の係数は東西どちらでも有意にゼロと異なる。したがって，その影響を調整することに意味はあると思われる。

9

そこで，次に，傾きに差があるかを解析する。「統計量」→「モデルへの適合」→「線型モデル」でモデル名

LinearModel.3 として左辺に TA1989 を，右辺に CAR1990*REGION を指定する。結果をみると，交互作用効果は有意でないので，2 本の回帰直線の傾きに有意差はないことがわかる。
そこで今度は，乗用車所有台数で調整した交通事故死者数の修正平均に差があるかどうかをみるため，交互作用項を除いて回帰を行う。モデル名を LinearModel.4 とし，右辺を CAR1990+REGION として線型モデルの当てはめを実行する。この結果，REGION の効果は有意なので，乗用車保有台数で調整すると西日本の方が東日本よりも人口当たり交通事故死者数が多いことがわかる。修正平均は以下の枠内を R Console に打てば得られる。単純な平均値は東日本が10.5，西日本が 9.87 であるが，乗用車保有台数の影響を調整した修正平均は，東日本が 9.44，西日本が 11.0 と逆転し，かつ有意差があることがわかる。

¶³

cfs <- dummy.coef(LinearModel.4)
cfs[[1]] + cfs$CAR1990 * mean(sample3$CAR1990) + cfs$REGIONµ´

1.9 ロジスティック回帰分析
ロジスティック回帰分析は，従属変数（目的変数。ロジスティック回帰分析では反応変数，あるいは応答変数と呼ぶこともある）が２値変数であり，二項分布に従うので lm() ではなく，glm() を使う。
ロジスティック回帰分析の思想としては，例えば疾病の有無を，複数のカテゴリ変数によって表される要因の有無と年齢のような交絡因子によって説明するモデルをデータに当てはめようとする。量的な変数によって表される交絡を調整しながらオッズ比を計算できるのが利点であり，医学統計ではもっともよく使われる手法の一つである。
疾病の有無は 0/1 で表され，データとしては有病割合（総数のうち疾病有りの人数の割合）となるので，そのままではモデルの左辺は 0 から 1 の範囲しかとらないが，右辺は複数のカテゴリ変数と量的変数（多くは交絡因子）からなるので実数のすべての範囲をとる。そのため，左辺をロジット変換（自身を 1 から引いた値で割って自然対数をとる）
する。
つまり，疾病の有病割合を P とすると，ln(P/(1 − P )) = b0 + b1 X1 + ...bk Xk と定式化できる。X1 が要因の有無を示す２値変数で，X2 , ...Xk が交絡であるとき，X1 = 0 の場合を X1 = 1 の場合から引けば，

b1 = ln(P1 /(1 − P1 )) − ln(P0 /(1 − P0 )) = ln(P1 ∗ (1 − P0 )/(P0 ∗ (1 − P1 )))となるので，b1 が他の変数の影響を調整したオッズ比の対数になる。対数オッズ比が正規分布するとすれば，オッズ比の 95% 信頼区間がexp(b1 ± 1.96 × SE(b1 ))として得られる。

¶
例題 7

³

library(MASS) の data(birthwt) は，Springfield の Baystate 医療センターの 189 の出生について，低体重出生とそのリス
ク因子の関連を調べたデータで，次の変数を含んでいる。低体重出生の有無を反応変数としたロジスティック回帰分析をせよ。

low
低体重出生の有無を示す２値変数（児の出生時体重 2.5 kg 未満が 1）
age
年齢
lwt
最終月経時体重（ポンド単位。略号 lb. で，1 lb. は 0.454 kg に当たる）
race 人種（１＝白人，２＝黒人，３＝その他）
smoke 喫煙の有無（１＝あり）
ptl
非熟練労働経験数
ht
高血圧の既往（１＝あり）
ui
子宮神経過敏の有無（１＝あり）
ftv
妊娠の最初の３ヶ月の受診回数
bwt
児の出生時体重 (g)µ´

データには多くの変数が含まれているが，本来，ロジスティック回帰分析では，反応変数に対する効果を見たい変数と交絡因子となっている変数はすべて説明変数としてモデルに投入するべきである（説明変数と反応変数の両方と10

有意な相関があれば交絡因子となっている可能性がある）。ここでは，丁寧な考察を経て，独立変数が人種，喫煙の有無，高血圧既往の有無，子宮神経過敏の有無，最終月経時体重，非熟練労働経験数となったとしよう。ロジスティック
回帰分析の前に，数値型で入っているカテゴリ変数を要因型に変換しておく。「データ」「アクティブデータセット内の変数の管理」「数値変数を因子に変換」を選び，まず変数として low を選び，新しい変数名を lowc として OK ボタ
ンをクリックする。数値 0 が水準 1 となり（NBW と名付ける），数値 1 が水準 2 となる（LBW と名付ける）。次に

race を選び，新変数名を racec として OK ボタンをクリックし，出てくるウィンドウで第 1 水準に”white”，第 2 水準に”black”，第 3 水準に”others”とカテゴリ名を指定し，OK ボタンをクリックする。smoke，ht，ui についても同様にカテゴリ変数 smokec，htc，uic に変換する。
ロジスティック回帰分析は，「統計量」「モデルへの適合」「一般化線型モデル」で，式の左辺に lowc (因子) をクリックして代入し（たんに lowc と入る），右辺に racec+smokec+htc+uic+lwt+ptl と打つ（またはクリックして選ぶ）。リンク関数族を binomial にして，リンク関数を logit にして OK すると，出力ウィンドウに次の枠内が表示さ
れる。

¶

³

Call:
glm(formula = lowc ~ racec + smokec + htc + uic + lwt + ptl, family = binomial(logit),
data = birthwt)
Deviance Residuals:
Min
1Q
Median-1.9049 -0.8124 -0.5241

3Q
0.9483Max
2.1812

Coefficients:
Estimate Std. Error z value Pr(>|z|)(Intercept)-0.086550
0.951760 -0.091 0.92754
racec[T.black]1.325719
0.522243
2.539 0.01113 *racec[T.other]
0.897078
0.433881
2.068 0.03868 *smokec[T.smoke]
0.938727
0.398717
2.354 0.01855 *htc[T.hypertension]
1.855042
0.695118
2.669 0.00762 **uic[T.uterine-hyper] 0.785698
0.456441
1.721 0.08519 .
lwt-0.015905
0.006855 -2.320 0.02033 *ptl
0.503215
0.341231
1.475 0.14029--Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1(Dispersion parameter for binomial family taken to be 1)
Null deviance: 234.67
Residual deviance: 201.99
AIC: 217.99

on 188on 181degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 4µ´

この出力からロジスティック回帰分析の結果は，下表のようにまとめられる。このように，量的な変数は表の下に共変量として調整したと書くのが普通である。また，係数は対数オッズ比でなく指数をとってオッズ比に直し，95% 信
頼区間も表示する。この操作は残念ながら Rcmdr ではまだできないので，分析結果が付値されている変数（回帰分析のモデルを指定するウィンドウの中で「モデル名」として指定したもの）が GLM.1 だったとすると，R Console で，exp(coef(GLM.1)) とすればオッズ比の点推定量が得られるし，exp(confint(GLM.1)) とすれば 95% 信頼区間が得
られる。

11

表. Baystate 医療センターにおける低体重出生リスクのロジスティック回帰分析結果
95% 信頼区間
独立変数 ∗

オッズ比

下限

p値

上限

人種（白人）

3.765 1.355 10.68 0.011
2.452 1.062 5.878 0.039
喫煙あり（なし）
2.557 1.185 5.710 0.019
高血圧既往あり（なし）6.392 1.693
27.3 0.008
子宮神経過敏あり（なし）2.194 0.888 5.388 0.085
AIC: 217.99，Dnull : 234.67（自由度 188），D: 201.99（自由度 181）
黒人

他の有色人種∗

カッコ内はリファレンスカテゴリ。これらの変数の他，最終月経時体重と非熟練労働経験数を共変量としてロジス

ティック回帰モデルに含んでいる。

¶
例題 8

³

前回 2 群の平均値の差の検定で用いたデータセット infert を用いて，不妊かどうか (case) を従属変数，年齢 (age)，既往出
生児数 (parity)，教育歴 (education)，人工妊娠中絶経験数 (induced)，自然流産経験数 (spontaneous) を独立変数とするロ
ジスティック回帰分析を実行せよ。

µ´
Rcmdr で「データ」→「パッケージ内のデータ」→「アタッチされたパッケージからデータセットを読み込む」として，パッケージとして datasets，データとして infert を選択するところまでは前回と同じである。
次 い で ，「 統 計 量 」→「 モ デ ル へ の 適 合 」→「 一 般 化 線 型 モ デ ル 」と し て ，出 て く る ウ ィ ン ド ウ 上 で ，「 モデ ル 名 を 入 力 」は GLM.2 の ま ま（ 別 の 名 前 を つ け て も 構 わ な い ），モ デ ル 式 と し て 左 辺 に case，右 辺 を

age + parity + education + induced + spontaneous とし，リンク関数族（ダブルクリックで選択）のところを binomial，リンク関数を logit として OK ボタンをクリックすれば，下記の結果が得られる。

¶

Deviance Residuals:
Min
1Q
Median-1.7603 -0.8162 -0.4956

3Q
0.8349³Max
2.6536

Coefficients:
Estimate Std. Error z value Pr(>|z|)(Intercept)-1.14924
1.41220 -0.814
0.4158age
0.03958
0.031201.269
0.2046
parity-0.82828
0.19649 -4.215 2.49e-05 ***education[T.6-11yrs] -1.04424
0.79255 -1.318
0.1876
education[T.12+ yrs] -1.40321
0.83416 -1.682
0.0925 .
induced1.288760.301464.275 1.91e-05 ***spontaneous
2.04591
0.31016
6.596 4.21e-11 ***--Signif. codes: 0 *** 0.001 ** 0.01 * 0.05 . 0.1

1(Dispersion parameter for binomial family taken to be 1)
Null deviance: 316.17
Residual deviance: 257.80
AIC: 271.8

on 247
on 241degrees of freedom
degrees of freedom

Number of Fisher Scoring iterations: 4µ´

この結果から，既往出生児数が多いほど不妊になりにくく，人工妊娠中絶や自然流産が多いほど不妊になりやすいことがわかる。R Console で exp(coef(GLM.2)) とすると，既往出生児数が 1 増えると不妊リスクは 0.44 倍になり，人12

工妊娠中絶経験数が 1 増えると不妊リスクが 3.6 倍，自然流産が 1 増えると不妊リスクが 7.7 倍になることがわかる。
なお，ペアごとに不妊リスクが異なることを仮定した条件付ロジスティック回帰分析をするには survival ライブラリの clogit() 関数を用いる必要があるが，Rcmdr ではできないので説明は省略する。使い方は example(infert) を実行すればわかる。

2 2 つのカテゴリ変数間の関係
カテゴリ変数間の関係を調べるにはどうしたらよいだろうか？ もちろん，カテゴリ変数についても関連の強さをみる指標はあって，ファイ係数（記号は ρ を用いるのが普通）と呼ばれる指標は，要因の有無，発症の有無を 1,0 で表した場合のピアソンの積率相関係数と同じ計算式で得られる。θ1 , θ2 を発症者中の要因あり割合，非発症者中の要因あり
割合として，ρ =p(π1 − π2 )(θ1 − θ2 ) である。また，疫学研究では，人数あるいは人年の比を取ることによって，要因があった群が，要因がなかった群に比べて，どれくらい発症しやすいかを調べることが多い（オッズ比やリスク比やハザード比を求め，その 95% 信頼区間が 1 を含まないかどうかで，要因の有無が発症の有無に有意に影響しているかどうかを判定することが慣例的に行われる）。
しかし，2 つのカテゴリ変数の関係を考えるとき，一般に，もっともよく行われるのは，それらが独立であるという帰無仮説を立てて検定することである。
カテゴリ変数のもつ統計的な情報は，カテゴリごとの度数だけである。そこで，２つのカテゴリ変数の間に関係について検討したいときには，まずそれらの組み合わせの度数を調べた表を作成する（R では table() という関数が使える）。これをクロス集計表と呼ぶ。とくに，２つのカテゴリ変数が，ともに２値変数のとき，そのクロス集計は２×２
クロス集計表 [2 by 2 cross tabulation]（２×２分割表 [2 by 2 contingency table]）と呼ばれ，その統計的性質が良く調べられている。

2.1 独立性のカイ二乗検定
独立性の検定としては，２つのカテゴリ変数の間に関連がないと仮定した場合に推定される期待度数を求めて，それに観測度数が適合するかを検定するカイ二乗検定が最も有名である*3 。

B
B̄

A
a人
c人

Ā
b人
d人

２つのカテゴリ変数 A と B が，それぞれ「あり」「なし」の２つのカテゴリ値しかとらないとき，これら２つのカテ
ゴリ変数の組み合わせは「A も B もあり（A ∩ B ）」「A なし B あり（Ā ∩ B ）」「A あり B なし（A ∩ B̄ ）」「A も B もなし（Ā ∩ B̄ ）」の４通りしかない。それぞれの度数を数えあげた結果が，上記の表として得られたときに，母集団の確
率構造が，

B
B̄

A
π11
π21

Ā
π12
π22

であるとわかっていれば，N = a + b + c + d として，期待される度数は，BB̄であるから，

χ2 =

A
N π11
N π21

Ā
N π12
N π22

(a − N π11 )2(b − N π12 )2(c − N π21 )2(d − N π22 )2
+
+
+
N π11
N π12
N π21
N π22

*3 もちろん，ある種の関連が仮定できれば，その仮定の元に推定される期待度数と観測度数との適合を調べてもいいが，一般に，２つのカテゴ

リ変数の間にどれくらいの関連がありそうかという仮定はできないことが多い。そこで，関連がない場合の期待度数を推定し，それが観測値に適合しなければ関連がないとはいえない，と推論するのである。

13

として，自由度３のカイ二乗検定をすればよいことになる。しかし，普通，π は未知である。そこで，P r(Ā) = 1−P r(A)かつ，この２つのカテゴリ変数が独立ならば P r(A ∩ B) = P r(A)P r(B) と考えれば良い*4 ことを使って，P r(A) とP r(B) を母数として推定する*5 。P r(A) の点推定量は，B を無視して A の割合と考えれば (a + c)/N であることは自明である。同様に，P r(B) の点推定量は，(a + b)/N となる。したがって，π11 = P r(A ∩ B) = P r(A)P r(B) =(a + c)(a + b)/(N 2 ) となる。
同様に考えれば，母集団の各セルの確率は下式で得られる。

π12 = (b + d)(a + b)/(N 2 )π21 = (a + c)(c + d)/(N 2 )π22 = (b + d)(c + d)/(N 2 )これらの値を使えば，
2

χ2 =

2

2

{b − (b + d)(a + b)/N }{c − (a + c)(c + d)/N }{d − (b + d)(c + d)/N }{a − (a + c)(a + b)/N }+++{(a + c)(a + b)/N }{(b + d)(a + b)/N }{(a + c)(c + d)/N }{(b + d)(c + d)/N }=2(ad − bc)2 {(b + d)(c + d) + (a + c)(c + d) + (b + d)(a + b) + (a + c)(a + b)}(a + c)(b + d)(a + b)(c + d)N分子の中括弧の中は N 2 なので，結局，

χ2 =N (ad − bc)2
(a + c)(b + d)(a + b)(c + d)

ただし通常は，イェーツの連続性の補正を行う。カイ二乗分布は連続分布なので，各度数に 0.5 を足したり引いたりしてやると，より近似が良くなるという発想である。この場合，

χc 2 =

N (|ad − bc| − N/2)2(a + c)(b + d)(a + b)(c + d)

が自由度１のカイ二乗分布に従うと考えて検定する。なお，|ad − bc| が N/2 より小さいときは補正の意味がないので，

χ2 = 0 とするのが普通である。|ad − bc| < N/2 のとき，R の chisq.test() では Yates の元論文の主旨に従うということで補正されてしまうけれども，prop.test() では補正されない。
実際の検定はクロス集計表が既に得られているとき，例えば a=12, b=8, c=9, d=10 などとわかっていれば，「統計量」「分割表」「２元表の入力と分析」で表示される表の各セルに直接数字を入力し，必要な統計量のチェックボックスにチェックを入れて OK ボタンをクリックするだけで非常に簡単である。

¶
例題 9

³

肺ガンの患者 100 人に対して，１人ずつ性・年齢が同じ健康な人を対照として 100 人選び（この操作をペアマッ
チサンプリングという），それぞれについて過去の喫煙の有無を尋ねた結果，患者群では過去に喫煙を経験した人が 80 人，対照群では過去に喫煙を経験した人が 55 人だった。肺ガンと喫煙は無関係といえるか？

独立性のカ

イ二乗検定をせよ。

µ´

帰無仮説は，肺ガンと喫煙が無関係（独立）ということである。クロス集計表を作ってみると，肺ガン患者群

健康な対照群

合計

80
20
100

55
45
100

135
65
200

過去の喫煙経験あり
過去の喫煙経験なし
合計*4 この帰無仮説は，合計に比例する割合で人数配分が行われていることに相当するので，B あり群と B なし群のそれぞれについて，A ありの割合に差がないという，比率の差の検定の帰無仮説と数学的に等価である。
*5 P r(X) はカテゴリ X の出現確率を示す記号である。また，２つの母数をデータから推定するので，得られるカイ二乗統計量が従う分布の自由度は３より２少なくなり，自由度１のカイ二乗分布となる。

14

となる。肺ガンと喫煙が無関係だという帰無仮説の下で期待される各カテゴリの人数は，

喫煙あり
喫煙なし

肺ガンあり

肺ガンなし

135 × 100/200 = 67.5
65 × 100/200 = 32.5

135 × 100/200 = 67.5
65 × 100/200 = 32.5となる。従って，連続性の補正を行なったカイ二乗統計量は，

χ2c = (80 − 68)2 /67.5 + (55 − 67)2 /67.5 + (20 − 32)2 /32.5 + (45 − 33)2 /32.5 = 13.128...
となり，自由度１のカイ二乗分布で検定すると 1-pchisq(13.128,1) より有意確率は 0.00029... となり，有意水準

5% で帰無仮説は棄却される。つまり，肺ガンの有無と過去の喫煙の有無には有意な関連があるといえる。生データがカテゴリ変数としてある場合は，「統計量」「分割表」「２元表」で２つのカテゴリ変数を選べばよい。
例題 7 のデータで low と smoke の独立性の検定をして試してみよ。

2.2 フィッシャーの直接確率（正確な確率）
期待度数が低い組み合わせがあるときには，カテゴリを併合して変数を作り直す方法もあるが，もっといい方法が考案されている。
ここで調べたいのは組み合わせの数なので，周辺度数を固定して（各々の変数については母比率が決まっていると仮定して）すべての組み合わせを考え，それらが起こる確率（超幾何分布に従う）を１つずつ計算し，得られている集計表が得られる確率よりも低い確率になるような表が得られる確率をすべて足し合わせてしまえば，２つのカテゴリ変数の間に関連がないという帰無仮説の下でそういう表が偶然得られる確率がどれほど低いのかを，直接計算することができる。こうして計算される確率を，フィッシャーの直接確率，あるいは，フィッシャーの正確な確率（検定）という。
これなら，近似ではないので，期待度数が低い組み合わせがあっても問題ない。
もう少し丁寧に言うと，サイズ N の有限母集団があって，そのうち変数 A の値が 1 である個体数が m1 ，1 でない個
体数が m2 あるときに，変数 B の値が 1 である個体数が n1 個（1 でない個体数が n2 = N − n1 個）あるという状況を考え，この n1 個のうち変数 A の値が 1 である個体数がちょうど a である確率を求めることになる。これは，m1 個か
ら a 個を取り出す組み合わせの数と m2 個から n1 − a 個を取り出す組み合わせの数を掛けて，N 個から n1 個を取り出す組み合わせの数で割った値になる。これと同じ周辺度数をもつ２×２分割表のうち，確率がこれと同じかこれよりも小さい表の確率をすべて足し合わせたものが，「変数 A と変数 B が独立」という帰無仮説が成り立つ確率になる*6 。
クロス集計表を使って２つのカテゴリ変数間の独立性の検定をするときは，コンピュータが使えるならば，サン
プルサイズがよほど大きくない限り，常にカイ二乗検定ではなく，フィッシャーの正確な確率を求めるべきである。

Rcmdr では，２元表の分析で「フィッシャーの正確検定」にチェックを入れておけば計算してくれる。

3 生存時間解析
生存時間解析は，まだ Rcmdr ではできないが，重要なのでごく簡単に説明を書いておく。

3.1 生存時間解析とは実験においては，化学物質などへの１回の曝露の影響を時間を追ってみていくことが良く行われる。時間ごとに何らかの量の変化を追うほかに，エンドポイントを死亡とした場合，死ぬまでの時間を分析することで毒性の強さを評価す*6 有限母集団からの非復元抽出になるので，平均 E(a) と分散 V (a) は，E(a) = n1 m1 /N
V (a) = {(N − n1 )/(N − 1)}n1 (m1 /N )(m2 /N ) = (m1 m2 n1 n2 )/{N 2 (N − 1)}となる。実際には組み合わせ計算が多いので，手計算で実行することはまずありえず，統計ソフトにやらせることになる。また，個々の２×２分割表の確率は離散値をとるので，同じ確率の表がありうる場合に，それを足し算に含めるのかどうかは難しい点である。これを乱数によって決める「ランダム検定」という手法もあるが，あまり一般的ではない。

15

ることができる。このような期間データを扱う方法としては，一般に生存時間解析 (Survival Analysis または Event

History Analysis) と呼ばれるものがある。なかでもよく知られているものが Kaplan-Meier の積・極限推定量である（現在では一般に，カプラン・マイヤ推定量と呼ばれている）。カプラン・マイヤ推定量は，イベントが起こった各時点
での，イベントが起こる可能性がある人口（リスク集合）あたりのイベント発生数を１から引いたものを掛け合わせて得られる，ノンパラメトリックな最尤推定量である。カプラン・マイヤ推定量は非常によく使われるので，具体例で説明しておく（後述）。複数の期間データ列があったときに，それらの差を検定したい場合は，ログランク検定や一般化
ウィルコクソン検定が使われる。
それらのノンパラメトリックな方法とは別に，イベントが起こるまでの時間が何らかのパラメトリックな分布に当て
はまるかどうかを調べる方法もある。当てはめる分布としては指数分布やワイブル分布がある。イベントが起こるまでの期間に何らかの別の要因が与える効果を調べたいときはコックス回帰（それらが基準となる個体のハザードに対してP
exp( βi zi ) という比例定数の形で掛かるとする比例ハザード性を仮定する方法）と，パラメトリックなモデルに対数線形モデルの独立変数項として入れてしまう加速モデルがある。R では生存時間解析をするための関数は survival
パッケージで提供されており，library(survival) とすれば使えるようになる。カプラン・マイヤ法は survfit()関数，コックス回帰は coxph() 関数，加速モデルは survreg() 関数で実行できる。なお，生存時間解析について，より詳しく知りたい方は，大橋，浜田 (1995) などを参照されたい。

3.2 カプラン・マイヤ法
まず，カプラン・マイヤ推定量についての一般論を示す。イベントが起こる可能性がある状態になってから，イベン
トが起こった時点を t1 , t2 , ... とし，t1 時点でのイベント発生数を d1 ，t2 時点でのイベント発生数を d2 ，以下同様であるとする。また，時点 t1 , t2 , ... の直前でのリスク集合の大きさを n1 , n2 , ... で示す。リスク集合の大きさとは，その直前でまだイベントが起きていない個体数である。観察途中で死亡や転居などによって打ち切りが生じるために，リスク集合の大きさはイベント発生によってだけではなく，打ち切りによっても減少する。従って ni は，時点 ti より前にイベント発生または打ち切りを起こした個体数を n1 から除いた残りの数となる。なお，イベント発生と打ち切りが同時点で起きている場合は，打ち切りをイベント発生直後に起きたと見なして処理するのが慣例である。このとき，カプラ
ン・マイヤ推定量 Ŝ(t) は，Ŝ(t) = (1 − d1 /n1 )(1 − d2 /n2 )... =Y(1 − di /ni )i<tとして得られる。その標準誤差はグリーンウッドの公式により，var(Ŝ) = Ŝ 2 ×

X

di
n
(n
− di )iii<tで得られる。なお，カプラン・マイヤ推定量を計算するときは，階段状のプロットを同時に行うのが普通である。

R では，library(survival) としてパッケージを呼び出し，dat <- Surv(生存時間，打ち切りフラグ) 関数で生存時間データを作り（打ち切りフラグは 1 でイベント発生，0 が打ち切り。ただし区間打ち切りの場合は 2 とか 3 も使う），res <- survfit(dat~1) でカプラン・マイヤ法によるメディアン生存時間が得られ（群分け変数 C によって群ごとにカプラン・マイヤ推定をしたい場合は，res <- survfit(dat~C) とすればよい），plot(res) とすれば階段状の生存曲線が描かれる。イベント発生時点ごとの値を見るには，summary(res) とすればよい*7 。

*7 参考までに書いておくと，生データがイベント発生の日付を示している場合，間隔を計算するには difftime() 関数や ISOdate() 関数を使う。例えば 1964 年 8 月 21 日生まれの人の今日の年齢は integer(difftime(ISOdate(2007,6,13),ISOdate(1964,8,21)))/365.24とすれば得られるし，さらに 12 を掛ければ月単位になる。

16

¶
例題 10

³

survival ライブラリに含まれているデータ aml は，急性骨髄性白血病 (acute myelogenous leukemia) 患者が化学療法によって寛解した後，ランダムに 2 群に分けられ，1 群は維持化学療法を受け（維持群），もう 1 群は維持化学療法を受けずに（非維持群），経過観察を続けて，維持化学療法が再発までの時間を延ばすかどうかを調べたデータであるa 。以下の３つの変数が含まれている。

time 生存時間あるいは観察打ち切りまでの時間（週）status 打ち切り情報（0 が観察打ち切り，1 がイベント発生）
x
維持化学療法が行われたかどうか（Maintained が維持群，Nonmaintained が非維持群）
薬物維持化学療法の維持群と非維持群で別々に，再発までの時間の中央値をカプラン＝マイヤ推定し，生存曲線をプロット
せよ。
a 出典:

Miller RG: Survival Analysis. John Wiley and Sons, 1981. 元々は，Embury SH, Elias L, Heller PH, Hood CE,
Greenberg PL, Schrier SL: Remission maintenance therapy in acute myelogenous leukaemia. Western Journal of Medicine,
126, 267-272, 1977. のデータ。研究デザインは Gehan と似ている。

µ
¶

´
³

library(survival)
print(res <- survfit(Surv(time,status)~x, data=aml))plot(res,xlab="(Weeks)",lty=c(1,2),main="Periods until remission of acute myelogenous leukaemia")legend("right",lty=c(1,2),legend=c("Maintained","Nonmaintained"))µ´

上枠内のように入力すると，2 行目でカプラン＝マイヤ法での計算がなされ，下枠内が表示される。

¶³

records n.max n.start events median 0.95LCL 0.95UCL
x=Maintained
11
11
11
7
31
18
NA
x=Nonmaintained
12
12
12
11
23
8
NA

µ

´

この表は，維持群が 11 人，非維持群が 12 人，そのうち再発が観察された人がそれぞれ 7 人と 11 人いて，維持群の再発までの時間の中央値が 31 週，非維持群の再発までの時間の中央値が 23 週で，95% 信頼区間の下限はそれぞれ 18
週と 8 週，上限はどちらも無限大であると読む。3 行目で 2 群別々の再発していない人の割合の変化が生存曲線として描かれ，4 行目で凡例が右端に描かれる。

3.3 ログランク検定
次に，ログランク検定を簡単な例で説明する。

8 匹のラットを 4 匹ずつ 2 群に分け，第 1 群には毒物 A を投与し，第 2 群には毒物 B を投与して，生存期間を追跡したときに，第 1 群のラットが 4,6,8,9 日目に死亡し，第 2 群のラットが 5,7,12,14 日目に死亡したとする。この場合，観察期間内にすべてのラットが死亡し，正確な生存時間がわかっているので，観察打ち切りがないデータとなっていて計算しやすい。
ログランク検定の思想は，大雑把に言えば，死亡イベントが起こったすべての時点で，群と生存／死亡個体数の 2 ×2 クロス集計表を作り，それをコクラン＝マンテル＝ヘンツェル流のやり方で併合するということである。上記の例では，死亡イベントが起こった時点 1〜8 において各群の期待死亡数を計算し，各群の実際の死亡数との差をとって，それに時点の重みを掛けたものを，各時点における各群のスコアとして，群ごとのスコアの合計を求める。２群しかない
ので，各時点において群１と群２のスコアの絶対値は同じで符号が反対になる。２群の生存時間に差がないという帰無仮説を検定するためには，群１のスコアの２乗を分散で割った値をカイ二乗統計量とし，帰無仮説の下でこれが自由度
１のカイ二乗分布に従うことを使って検定する。なお，重みについては，ログランク検定ではすべて 1 である。一般化
ウィルコクソン検定では，重みを，2 群を合わせたリスク集合の大きさとする（そうした場合，もし打ち切りがなければ，検定結果は，ウィルコクソンの順位和検定の結果と一致する）。つまり，ログランク検定でも一般化ウィルコクソ
ン検定でも，実は期間の情報はまったく使われず，死亡順位の情報だけが使われている。
記号で書けば次の通りである。第 i 時点の第 j 群の期待死亡数 eij は，時点 i における死亡数の合計を di ，時点 i に17

おける j 群のリスク集合の大きさを nij ，時点 i における全体のリスク集合の大きさを ni とすると，

eij = di · nij /niと表される。上の例では，e11 = 1 · n11 /n1 = 4/8 = 0.5 となる。時点 i における第 j 群の死亡数を dij ，時点の重みをwi と表せば，時点 i における群 j のスコア uij は，
uij = wi (dij − eij )となり，ログランク検定の場合（以下，重みは省略してログランク検定の場合のみ示す）の群 1 の合計スコアは

u1 =

X

di1 − ei1iとなる。上の例では，

u1 = (1 − 4/8) + (0 − 3/7) + (1 − 3/6) + (0 − 2/5) + (1 − 2/4) + (1 − 1/3) + (0 − 0/2) + (0 − 0/1)である。これを計算すると約 1.338 となる。分散は，分散共分散行列の対角成分を考えればいいので，

V = Vjj =X (ni − nij )nij di (ni − di )ni 2 (ni − 1)iとなる。この例の数値を当てはめると，

V =(8 − 4) × 4 (7 − 3) × 3 (6 − 3) × 3 (5 − 2) × 2 (4 − 2) × 2 (3 − 1) × 1
+
+
+
+
+
82
72
62
52
42
32

となり，計算すると，約 1.568 となる。したがって，χ2 = 1.3382 /1.568 = 1.14 となり，この値は自由度１のカイ二乗
分布の 95% 点である 3.84 よりずっと小さいので，有意水準 5% で帰無仮説は棄却されない。つまりこれだけのデータ
では，差があるとはいえないことになる（もちろん，サンプルサイズを大きくすれば違う結果になる可能性もある）。

R では，Surv(time,event) と group（注：ここで time は生存時間，event は 1 がイベント観察，0 が観察打ち切りを示すフラグ，group がグループを示す）を，survdiff() 関数に与えることによってログランク検定が実行できる。
打ち切りレコードがない場合は，event は省略できる。なお，生存時間解析の関数は survival パッケージに入っているので，まず library(survival) とすることは必須である。
上で説明した例では，¶³

library(survival)time <- c(4,6,8,9,5,7,12,14)
event <- c(1,1,1,1,1,1,1,1)
group <- c(1,1,1,1,2,2,2,2)
dat <- Surv(time,event)survfit(dat~group)
survdiff(dat~group)µ´2とすると解析結果が得られる。一番下を見ると，χ = 1.2，自由度 1，p = 0.268 となっているので，有意水準 5% で，2 群には差がないことがわかる。

18

4 レポート課題
http://phi.med.gunma-u.ac.jp/grad/sample4.dat は，web サイトで zip 形式に圧縮して公開されている「Theworld factbook 2007」(CIA)*8 から作ったタブ区切りテキスト形式のデータである。
このデータに含まれている変数は，次の通りである。
COUNTRY 国名GINI 国ごとの家族の所得の不平等度を示す Gini の集中係数（完全な平等のとき 0，完全な不平等で 100）
YEAR.x Gini の集中係数の報告年
LIFEEXP ゼロ歳平均余命
YEAR.y ゼロ歳平均余命の報告年
GDPPCUSD 米ドル換算購買力平価ベース 1 人当たり GDP
YEAR GDP 報告年
RP
GDP が 10000 以上だと RICH，10000 未満だと POOR となるカテゴリ変数
EU
GINI の集中係数が 35 以上だと Uneven（不平等），35 未満だと Egaritarian（平等主義者）となるカテゴリ変数
近年，社会疫学という研究分野において，健康が社会のありようによって影響を受けることが指摘されており，ゼロ
歳平均余命が所得の不平等度や国内総生産から受ける影響を調べることも行われているが，このデータをその視点で統計解析せよ。図示や記述統計をしてから，豊かな国と貧しい国の間で Gini の集中係数とゼロ歳平均余命の関係の違いを調べよ。違いがない場合には Gini の集中係数を共変量としてその影響を調整したゼロ歳平均余命の修正平均を，豊かな国と貧しい国で比較せよ。
提出期限は 6 月末日までとする。Microsoft Word，OpenOffice.org Writer，または pdf 形式でＡ４サイズ 3 枚以内にまとめて，中澤准教授 (nminato@med.gunma-u.ac.jp) までメールの添付ファイルとして提出すること。

¶
Report theme for foreign students³

The data http://phi.med.gunma-u.ac.jp/grad/sample4.dat is a tab-separated text file, which is originally published as “The world factbook 2007” (CIA) at the following url.
https://www.cia.gov/library/publications/the-world-factbook/index.html
The variables included there are:
COUNTRY The name of the country
GINI Gini index measuring the degree of inequality in the distribution of family income in a country (perfect equality
showing 0, perfect inequality showing 100)
YEAR.x Reported year of the above Gini index
LIFEEXP Life expectancy (years) at birth
YEAR.y Reported year of the above life expectancy
GDPPCUSD Gross Domestic Products per capita based on purchasing power parity in US dollar
YEAR Reported year of above GDP per capita
RP
Categorical variable as the RICH for “GDPPCUSD equal to or more than 10000” and the POOR for the rest
EU
Categorical variable as the Uneven for “GINI equal to or more than 35” and Egaritarian for the rest
Recently, many researchers of social epidemiology indicated that health levels of the country are affected by its social
settings; for instance, life expectancy at birth is affected by income inequality and by gross domestic products per
capita.
Let’s analyze this data from a such viewpoint of social epidemiology. At first, you must make graphs for each variable
and calculate descripitve statistics. After that, please test the difference of the relationships between Gini coefficientsand life expectancy at birth, between RICH and POOR countries.
If the difference was not significant, please compare adjusted means of life expectancy at birth with controlling Gini
index, between RICH and POOR countries, using the analysis of covariance.
Please submit the report to Minato Nakazawa via e-mail (nminato@med.gunma-u.ac.jp) until 30 June. The report
should be attached as a document file within 3 pages in A4, made by Microsoft Word or OpenOffice.org Writer, or
as a PDF file.

µ´*8 https://www.cia.gov/library/publications/the-world-factbook/index.html

19