医学情報処理演習第８回「相関と回帰」
中澤

港 (nminato@med.gunma-u.ac.jp)

2004 年 11 月 29 日¶
前回の課題の回答例³

アヤメの種（Species）によって萼片の幅（Sepal.Width）の平均値に差があるかどうかを調べるためには，まず Bartlett の検定によって，「萼片の幅の母分散がどの種でも同じ」という帰無仮説を検定しなくてはならない。bartlett.test(Sepal.Width~Species) とすると，Bartlett’s K-squared = 2.0911, df = 2, p-value = 0.3515
より，帰無仮説は棄却できないので，母分散は等しいとみて，「どの種でも萼片の幅の母平均は同じ」という帰無仮説を，一元配置分散分析で検定する。summary(aov(Sepal.Width~Species)) と打って得られる分散分析表をみると，Pr(>F) の値は< 2.2e-16 となっていて，ほとんどゼロなので帰無仮説は棄却される。つまり，萼片の幅の母平均が異なる種が少なくとも１つあると考えられる。
次 に ど の 種 と ど の 種 が 異 な る か を み る た め に 多 重 比 較 を す る 。

TukeyHSD(aov(Sepal.Width~Species)) と 入 力 す る と ，versicolor と setosa の 差 の 95%信 頼 区
間 は [-0.82,-0.50]，virginica と setosa の 差 の 95%信 頼 区 間 は [-0.61,-0.29]，virginica とversicolor の差の 95%信頼区間は [0.04,0.36] と，いずれも 0 を含まないので，有意水準 5%でテュー
キーの多重比較を行った結果，「どの２種の比較でも萼片の幅の母平均には差がないとはいえない」ことがわかる。
なお，多重比較は，pairwise.t.test(Sepal.Width,Species) として，ホルムの方法で調整した有意確率を得ても，いずれの組み合わせでも 5%より小さいので，テューキーの HSD でやった場合と結論は同じである。

µ´

相関と回帰の違い
相関と回帰は，どちらも２つの変数の関係を扱うので混同されやすいが，思想はまったく違う。相関は，変数間の関連の強さを表すものである。回帰は，ある変数の値のばらつきが，どの程度他の変数の値のばらつきによって説明されるかを示すものである。回帰の際に，説明される変数を従属変数または目的変数，説明するための変数を独立変数または説明変数と呼ぶ。２つの変数間の関係を予測に使うためには，回帰を用いる。

相関関係とか関連とかいっても，その中身は多様である。例えば，pV = nRT のような物理法則は，測定誤差を別にすれば 100%成り立つ関係である。身長と体重の間の関係はそうではないが，無関係ではないことは直感的にも理解できるし，「身長の高い人は体重も概して重い傾向がある」ことは間違いない。一般に，２個以上1の変量が「かなりの程度の規則正しさをもって，増減をともにする関係」のことを相関関係 (correlation) という。相関には正の相関 (positive correlation) と負の相関 (negative correlation) があり，一方が増えれば他方も増える場合を正の相関，一方が増えると他方は減る場合を負の相関と呼ぶ。例えば，身長と体重の関係は正の相関である。
相関関係があっても，それが見かけ上のものである（それらの変量がともに，別の変量と真の相関関係をもっている）場合がある。具体例としては，血圧と所得の間に正の相関があるという命題は，データをとってみれば，多くの場合に成り立つであろう。しかしこれは，おそらくどちらも年齢や摂取エネルギー量との間に真の相関関係が存在するのであって，それらの影響を制御したら（例えば同年齢で同じような食生活をしている人だけについて見る，という層別化をしたら），血圧と所得の間の正の相関は消えてしまうだろう。この場
合，見かけ上の相関があることは，たまたまそのデータで成り立っているだけであって，科学的仮説としての意味に乏しい。
時系列データや地域相関のデータでは，擬似相関 (spurious correalation) が見られる場合もある。例えば，日本の砂糖輸入量と溺死・溺水者数の年次別データをプロットしてみると，負の相関関係があるように見えるのだが，両者の間には真の関係はない。ある年に日本で植えた木の幹の太さと，同じ年にイギリスで生れた少年の身長を 15 年分，毎年１回測ったデータをプロットすると，おそらくは正の相関関係があるように見えるのだが，両者の間には関係がないのは明らかである（どちらも年次と真の相関があるとはいえるだろうが）。
複数の種類の異なるデータをまとめて見ることで見かけの相関が生じてしまう場合もあるので，注意が必要である。
上で定義したように，相関関係は増減をともにする関係であればいいので，その関係が線形（一次式で表される，散布図で直線として表される）であろうと非線形（二次式以上または階段関数などで表される）であろうと問題ない。しかし，一般には，線形の関係があるという限定的な意味で使われる場合が多い。なぜなら，相関を表すための代表的な指標である相関係数∗1 r が，線形の関係を示すための指標だからである。もっといえば，r が意味をもつためには，厳密に言えば，２つの変量が二次元正規分布に従っていなければならない。
非線形の相関関係を捉えるには，２つのアプローチがある。１つは線形になるように対数変換などの変換をほどこすことで，もう１つはノンパラメトリックな相関係数（分布の形によらない，例えば順位の情報だけを使った相関係数）を使うことである。ノンパラメトリックな相関係数にはスピアマン (Spearman) の順位相関係数 ρ や，ケンドール (Kendall) の順位相関係数 τ がある。
ピアソンの積率相関係数とは，X と Y の共分散を X の分散と Y の分散の積の平方根で割った値である。
式で書けば，相関係数の推定値 r は，X の平均を X̄ ，Y の平均を Ȳ と書けば，Pni=1 (Xi − X̄)(Yi − Ȳ )Pnn22i=1 (Xi − X̄)
i=1 (Yi − Ȳ )r = qPとなる。母相関係数がゼロかどうかという両側検定のためには，それがゼロであるという帰無仮説の下で，検定統計量√r n−2t0 = √1 − r2が，自由度 n − 2 の t 分布に従うことを利用して検定すればよい。
なお R では，

r <- cov(X,Y)/sqrt(var(X)*var(Y))∗1 普通，ただ相関係数といえば，ピアソンの積率相関係数 (Pearson’s Product Moment Correlation Coefficient) を指し，通常，r という記号で表す。

2

n <- length(X)t0 <- r*sqrt(n-2)/sqrt(1-r^2)
として，2*(1-pt(t0,n-2)) で有意確率が得られるが，cor.test() 関数（下記）を使う方が簡単である∗2 。

cor.test() 関数を使った場合は，信頼区間も計算される。なお，信頼区間は，サンプルサイズがある程度大きければ（通常は 20 以上），正規近似を使って計算できる。すなわち，a=1 1+r1ln−√Z(α/2),
2 1−rn−3b=1 1+r1ln+√Z(α/2)
2 1−rn−3と書くことにすると∗3 ，母相関係数の 100 × (1 − α)%信頼区間の下限は (exp(2a) − 1)/(exp(2a) + 1)，上限は (exp(2b) − 1)/(exp(2b) + 1) である∗4 。
順位相関係数は，非線形の相関関係を捉えたい場合以外にも，分布が歪んでいたり，外れ値がある場合に使うと有効である。スピアマンの順位相関係数 ρ は∗5 ，値を順位で置き換えた（同順位には平均順位を与えた）ピアソンの積率相関係数になる。Xi の順位を Ri ，Yi の順位を Qi とかけば，

ρ=1−nX6(Ri − Qi )2n(n2 − 1) i=1となる。スピアマンの順位相関係数がゼロかどうかという両側検定は，サンプル数が 10 以上ならばピアソン√ρ n−2の場合と同様に，T = √が自由度 n − 2 の t 分布に従うことを利用して行うことができる。
21−ρ

ケンドールの順位相関係数 τ は，

τ=(A − B)n(n − 1)/2によって得られる。ここで A は順位の大小関係が一致する組の数，B は不一致数である。

R では cor.test(X, Y, method="pearson") とすれば（あるいは method=オプションをつけなければ）ピアソンの相関係数が，cor.test(X, Y, method="spearman") でスピアマンの順位相関係数が，cor.test(X, Y, method="kendall") でケンドールの順位相関係数が得られる。同時に，alternative を指定しないときは，「相関係数がゼロである」を帰無仮説として両側検定した有意確率と 95%信頼区間が表示される。なお，例えば cor.test(X, Y, alternative="g") とすれば，ピアソンの相関係数が計算され，対立仮説を「正の相関がある」とした片側検定の結果が得られる。なお，ケンドールに関しては並べ換えによる正確な確率も求めることができ，その場合は exact=T というオプションを指定する。

¶例題³ToothGrowth は各群 10 匹ずつのモルモットに 3 段階の用量のビタミン C をアスコルビン酸としてある
いはオレンジジュースとして投与したときの象牙芽細胞（歯）の長さを比較するデータである。変数 lenが長さ，supp が投与方法，dose が用量を示す。投与方法の違いを無視して用量と長さの相関関係を調べよ。

µ´

まず，以下のようにして ToothGrowth に含まれている変数が使えるようにする。

∗2 また，cov(X,Y)/sqrt(var(X)*var(Y)) と同値な関数として cor(X,Y) がある。
∗3 Z(α/2)は ，標 準 正 規 分 布 の 100 × (1 − α/2) パ ー セ ン ト 点 ，つ ま り R で は（α を alpha と 書 く こ と に し て ），qnorm(1-alpha/2,0,1) である。例えば有意水準 5%，すなわち α = 0.05 なら，qnorm(0.975,0,1) である。
∗4 なお，ln は自然対数，exp は指数関数を表す。この式から明らかなように，母相関係数の信頼区間はどんなに広がっても下限は−1 以下にはならず，上限は 1 以上にならない。
∗5 ピアソンの相関係数の母相関係数を ρ と書き，スピアマンの順位相関係数を r と書く流儀もある。
s3¶³

data(ToothGrowth)
attach(ToothGrowth)µ´

作図によって用量と長さの関係を概観するためには，散布図を描けばよいので，横軸を dose，縦軸を lenとしたプロットをするために，plot(dose,len) とする。なんとなく dose が増すにつれて len が長くなっていくような，正の相関関係があるように見えるだろう。そこで，この相関係数を算出し，「相関係数がゼロと差がない」という帰無仮説を検定してみるために，cor.test(dose,len) と打てば，以下の出力が得られる。

¶³

Pearson’s product-moment correlationdata:

dose and len

t = 10.2501, df = 58, p-value = 1.243e-14
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.6892521 0.8777169
sample estimates:
cor0.8026913µ´

この結果をみると，ピアソンの相関係数の推定値は 0.80，95%信頼区間が [0.69, 0.88]（ここでは四捨五入で示しているが，真の区間を含む最小の幅に丸める方がよいという意見もあり，それなら [0.68, 0.88] と記載する）となる。95%信頼区間がゼロを含んでいないので，帰無仮説が有意水準 5%で棄却されるのは明らかだが，有意確率をみても，p-value=1.243e-14 とほとんどゼロであることが確認できる。method="spearman"とか method="kendall"でも試してみよう。

回帰
実験によって，あるサンプルの濃度を求めるやり方の１つに，検量線の利用がある。検量線とは，予め濃度がわかっている標準物質を測ったときの吸光度のばらつきが，その濃度によってほぼ完全に（通常 98%以上）説明されるときに（そういう場合は，散布図を描くと，点々がだいたい直線上に乗るように見える），その関係を利用して，サンプルを測ったときの吸光度からサンプルの濃度を逆算するための回帰直線である（曲線の場合もあるが，通常は何らかの変換をほどこし，線形回帰にして利用する）。検量線の計算には，(A) 試薬ブ
ランクでゼロ点調整をした場合の原点を通る回帰直線を用いる場合と，(B) 純水でゼロ点調整をした場合の切片のある回帰直線を用いる場合がある。
いずれも，量がわかっているもの（この場合は濃度）を x，誤差を含んでいる可能性がある測定値（この場合は吸光度）を y として y = bx + a という形の回帰式の係数 a と b を最小二乗法で推定し，サンプルを測定した値 y から x = (y − a)/b によってサンプルの濃度 x を求める。測定値から濃度を推定するときには，回帰式をそのまま使うのではなく，逆算する形になるので注意が必要である。
回帰直線の適合度の目安としては，学生実習でも相関係数の２乗が 0.98 以上あることが望ましい。また，データ点の最小，最大より外で直線関係が成立する保証はない。従って，サンプル測定値が標準物質の測定値の最小より低いか，最大より高いときは，限界を超えていることになってしまう∗6 。
∗6 余談だが，このような場合はサンプルを希釈するか濃縮して測定するのが普通である。

4

測定点 (x1 , y1 ), (x2 , y2 ), ..., (xn , yn ) が得られたときに，検量線 y = bx + a を推定するには，図に示した線分の二乗和が最小になるように a と b を設定すればよい，というのが最小二乗法の考え方である。つまり，

f (a, b) =nX{yi − (bxi + a)}2i=1= b2nXx2i − 2bi=1nXxi yi + 2abi=1nXxi − 2ai=1nX

yi + na2 +i=1nXyi2i=1が最小になるような a と b を推定すればよい。通常，a と b で偏微分した値がそれぞれ 0 となることを利用して計算すると簡単である。つまり，
nnXX∂f (a, b)= 2na + 2(bxi −yi ) = 0∂ai=1i=1i.e. na =n
Xyi − bi=1i.e.

nXxii=1a = (y の平均) − (x の平均) ∗ bnnn
XXX
∂f (a, b)= 2bx2i + 2(axi −xi yi ) = 0∂bi=1i=1i=1i.e. bnX
i=1x2i =n
Xxi yi − ai=1n
Xxii=1を連立方程式として a と b について解けばよい。これを解くと，

b=n

Pn
Pn
Pn
xi yi − i=1 xi i=1 yii=1Pn
Pn
n i=1 x2i − ( i=1 xi )2が得られる∗7 。b の値を上の式に代入すれば a も得られる。検量線に限らず，一般の回帰直線でも，計算方法は原則として同じである。名称の説明をしておくと，一般に，y = bx + a という回帰直線について，b を回帰係数 (regression coefficient)，a を切片 (intercept) と呼ぶ。

R では，線形回帰を行うための関数は lm() である。例えば，lm(Y~X) のように用いれば，回帰直線の推定値が得られる。散布図の上に回帰直線を重ね描きするには，plot(Y~X) としてから，abline(lm(Y~X)) とすればよい∗8 。決定係数や回帰係数と切片の検定結果は，summary(lm(Y~X)) とすれば出力される。

∗7 分母分子を n2 で割れば，b は x y の平均から x の平均と y の平均の積を引いて，x の二乗の平均から x の平均の二乗を引いi iiiiiた値で割った形になる。
∗8 なお，回帰直線の 95%予測区間と 95%信頼区間を点線と破線で重ね描きするには，XX <- data.frame(X=seq(range(X)[1],range(X)[2],length=20))
plim <- predict(res, XX, interval="prediction")clim <- predict(res, XX, interval="confidence")lines(XX$X, plim[,2], lty=2)
lines(XX$X, plim[,3], lty=2)
lines(XX$X, clim[,2], lty=3)
lines(XX$X, clim[,3], lty=3)とすればよい。

5¶
例題

³

血清鉄濃度を Fe-Test Wako というキットで測定するため，鉄の標準希釈系列を 0, 0.5, 1, 2 (mg/L) として作成し，それをこのキットで処理して発色させた溶液の波長 562 nm の吸光度を測った結果が，0.012,0.058, 0.104, 0.193 として得られた。これから検量線を求めて，測定に使えるかどうか評価せよ。次に，6 人の血清サンプルを同じ方法で処理して発色させた溶液の吸光度が 0.107, 0.075, 0.077, 0.099, 0.096,0.108 だったときに，この 6 人の血清鉄濃度を求めよ。
µ´

まず濃度を conc，吸光度を abs としてデータを入力し，conc を横軸，abs を縦軸にして散布図を描く（描画命令は，ここでは plot(abs~conc) としたが，plot(conc,abs) でも同じことである）。

¶³

conc <- c(0, 0.5, 1, 2)
abs <- c(0.012, 0.058, 0.104, 0.193)
plot(abs~conc)µ´

だいたい直線に乗っているように見えるので，回帰分析をしてみる。以下のように，回帰分析の結果をいったん res に保存しておくと，描画や表示やその後の計算に便利である。

¶³

res <- lm(abs~conc)
abline(res)summary(res)µ´

出力される結果から，Adjusted R-squared: 0.9999 なので検量線として使ってよいと判断できる。回帰係数は 0.0904571，切片は 0.0126000 として得られているが，これらの値はそれぞれ res$coef[2]，res$coef[1] として参照できるので（注：res$coef の部分は coef(res) でも同じ意味），サンプルの吸光度から濃度を計算するときには，以下のように変数名のまま参照した方が間違えない。

¶³

dat <- c(0.107, 0.075, 0.077, 0.099, 0.096, 0.108)(dat-res$coef[1])/res$coef[2]µ´

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

Q=nX

e2i =nX(yi − zi )2i=1=i=1PnPnPnnX(n i=1 xi yi − i=1 xi i=1 yi )2PnPyi2 − (yi )2 /n −/nnn i=1 x2i − ( i=1 xi )2i=1i=1nXとすれば，この Q が回帰直線の当てはまりの悪さを示す尺度となる。この Q を「残差平方和」と呼び，それを n で割った Q/n を残差分散という。この残差分散 var(e) と Y の分散 var(Y ) とピアソンの相関係数 r の6

間には，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 分布に従うことになる。しかしこの値は 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」であるという帰無仮説に対する有意確率が計算できる。
以上の説明からすると，身長と体重のように，どちらも誤差を含んでいる可能性がある測定値である場合には，どちらかを独立変数，どちらかを従属変数，という風にみなしてよいのかということが問題になってくる。一般には，身長によって体重が決まってくるというように方向性が仮定できれば，身長を独立変数と見なしてもよいことになっているが，回帰分析をしてしまうと，独立変数に測定誤差がある可能性が排除されてしまうことには注意しておくべきである。つまり，測定誤差が大きい可能性がある変数を独立変数とした回帰分析は，できれば避けたほうが良い（が，そうもいかないのが実情である）。また，最小二乗推定の説明から自明なように，独立変数と従属変数を入れ替えた回帰直線は一致しないので，どちらを従属変数とみなし，どちらを独立変数とみなすか，ということは，因果関係の方向性に基づいてきちんと決めるべきである。
回帰を使って予測をするとき，外挿には注意が必要である。前述の通り，検量線は，原則として外挿してはいけない。実際に測った濃度より濃かったり薄かったりするサンプルに対して，同じ関係が成り立つという保証はどこにもないからである（吸光度を y とする場合は，濃度が高くなると分子の重なりが増えるので飽和(saturate) してしまい，吸光度の相対的な上がり方が小さくなっていき，直線から外れていく）。しかし，外挿による予測は，実際にはかなり行われている。例えば世界人口の将来予測とか，河川工学における基本高水計算式とか，感染症の発症数の将来予測は，回帰の外挿による場合が多い（この場合は逆算ではなく，実際にデータを得られていない横軸の値を代入したときに縦軸の値がいくつになるかを予測することになる）。このやり方が妥当性をもつためには，その回帰関係が（１）かなり説明力が大きく，（２）因果関係がある程度認め

7

られ，（３）それぞれの変数の分布が端の切れた分布でない（truncated distribution でない）という条件を満たす必要がある。そうでない場合は，その予測結果が正しい保証はどこにもない。

課題
data(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（日）である。
このデータから，オゾン濃度と風速の関係について検討せよ。また，もし風速からオゾン濃度を予想できるとしたら，風速 25 マイル／時の日のオゾン濃度はどのように予測されるか計算せよ。
結果は配布する紙に学籍番号，氏名と共に自筆して提出すること。結果の提出をもって出席確認とする。

8