ある患者群の主観的幸福感（QOL；VASで取得）、栄養状態（MNA）、日中の身体活動（Mets）、痛み（VASで取得）を想定したサンプルデータです。ここからQOLが、栄養状態、身体活動、痛みで説明できるかを検証します。
|
QOL |
mna |
mets |
pain |
|
83 |
24 |
2.4 |
36 |
|
84 |
26 |
2.2 |
40 |
|
中略 |
中略 |
中略 |
中略 |
|
79 |
24 |
2.2 |
40 |
|
93 |
26 |
2.5 |
55 |
コードと結果は以下の通りです。
> model1<-lm(qol ~ mna + mets + pain, data = data1) > model1 |> summary() Call: lm(formula = qol ~ mna + mets + pain, data = data1) Residuals: Min 1Q Median 3Q Max -6.2957 -1.5591 0.0094 1.5762 5.7335 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 30.95741 4.67349 6.624 1.44e-08 *** mna 0.96400 0.12507 7.708 2.34e-10 *** mets 6.53627 1.83639 3.559 0.000766 *** pain 0.34312 0.03929 8.734 4.87e-12 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.757 on 56 degrees of freedom Multiple R-squared: 0.7574, Adjusted R-squared: 0.7444 F-statistic: 58.28 on 3 and 56 DF, p-value: < 2.2e-16
重回帰のモデルにmodel1という名前を付けてコードを記載しました。lm()関数を使うことで、単回帰だけでなく重回帰のモデルを作成することができます。今回はQOLが目的変数で、MNA、Mets、painが説明変数です。目的変数と説明変数は ~ で結び、説明変数は+で加えていき、モデル式をつくります。今回はqol ~ mna + mets + painとなります。その後、 , で区切ってデータの名前を指定します。そして、model1をパイプ演算子でsummary()関数とつなぎます。
結果です。Coefficients:を確認します。Estimateは回帰係数、Std. Errorは標準誤差、t valueがt値で、Prがp値です。ただし、InterceptのEstimateは回帰係数ではなく、切片の値です。ここから以下のようなモデル式が導かれます。
qol = 0.97329×mna + 16.90631×mets + 0.37702×pain + 5.35034
また、5%を有意水準と定めると、すべての回帰係数は有意です。ここから、回帰係数は０ではないと結論づけられ、すべての説明変数は、QOLに正の影響を及ぼしていると言えそうです（もし係数が負の値なら負の影響を及ぼすことになります）。注意が必要ですが、回帰係数だけみると約16.9のmetsの影響が抜き出て大きいようにも見えます。Metsの値は、0.1刻みで変化する値で、１変化することは考えにくいと言えます。一方で、VASで測定している痛みの測定値は、最大で100となるため1の変化はあまり大きいと言えません。このように、説明変数の値の意味が異なるため、単純な比較はできないことに注意しましょう。どの回帰係数がどれだけの影響があるのか、を考えるときには、その説明変数の値が持っている意味をドメイン知識に照らし、判断していくのが大切です。
また、モデル式のデータに対する当てはまりをチェックするための決定係数は、Multiple R-squaredで、Adjusted R-squaredは調整済みの決定係数になります。調整済み決定係数は0.7444と非常に高い値となっています。
多重共線性を確認します。各説明変数の相関をチェックすることも行われますが、ここではVIFを確認します。carパッケージのvif()関数が必要になります。carパッケージは、rstatixパッケージをインストールした際に同時にインストールされているパッケージです。ですので改めてインストールする必要はありません。コードと結果は以下のようになります。
> model1 |> vif()
mna mets pain
1.094704 1.097706 1.003791 library()関数でcarパッケージを読み込みます。それからmodel1とvif()関数つなぐ、もしくはvif()関数の()内にモデル名を指定します。結果は10よりも小さくすべての説明変数を投入しても問題なさそうです。
回帰モデルのモデル選択について解説します。複数のモデルから予測性能の高いモデルを選択するために赤池情報量基準を用います。model1以外のモデルを作成します。mna、mets、painを一つずつ減らした説明変数が2つだけのモデルを作りました。コードは以下の通りで、全部で4つになります。
> model2<-lm(qol ~ mna + mets, data = data1) > model3<-lm(qol ~ mna + pain, data = data1) > model4<-lm(qol ~ mets + pain, data = data1) > > model1 |> AIC() [1] 297.8261 > model2 |> AIC() [1] 347.4006 > model3 |> AIC() [1] 308.0626 > model4 |> AIC() [1] 339.2155
コメント