単回帰分析

データの準備と要約

僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。

例題
PT単位数がFIM運動項目(mFIM)に与える影響について検証します
対象:9名の脳卒中患者
目的変数:FIM差 = PT後mFIM - PT前mFIM
説明変数:理学療法実施単位数(PT単位数)

データセットlm01をdatに格納します(ファイルの読み込み

dat <- read.csv("lm01.csv", header=T, fileEncoding = "UTF-8")
#View(dat)

散布図

plot(FIM差~PT単位数, data=dat)

Rで単回帰分析

lmはlinear model(linear regression model)の略です

fit <- lm(FIM差~PT単位数, data=dat)
print(fit)

これで以下のような回帰式が求められます

$y=2.86+0.50x$

切片と傾きは、有意な数値なのか有意水準5%で検定してみましょう

summary(fit)

傾きのp値(0.0199)は0.05未満なので有意です

PT単位数が1単位増えるとFIM差は約0.5増加することが予測できました

F値(9.012)が出力されています

結果のみが必要な場合は、ここまで大丈夫です

ここからはもう少し詳しくなります

分散分析

分散分析を実行することで、単回帰分析の結果に記載されているF値(9.012)が算出されます

anova(fit)

分散分析表

詳細は分散分析のページをご参照ください

回帰分析表と分散分析表について、重回帰分析のページに記載してます

グラフ

散布図に回帰直線を加えてみます

グラフの描き方

plot(FIM差~PT単位数, data=dat)
x0 <- range(dat$PT単位数) #x軸の範囲
b0 <- fit$coef[1] #切片
b1 <- fit$coef[2] #傾き
lines(x0, b0 + b1*x0) #回帰直線

p値の意味

おそらくFIM差22が1名いるので、切片が十分に予測できないようですね

試しに6行目のFIM差22を15に置き換えた以下のデータでやってみましょう

注意)あくまでも統計学の勉強のためのデータ操作です

dat2にデータを格納します

dat2 <- read.csv("lm02.csv", header=T, fileEncoding = "UTF-8")
#View(dat2)
fit2 <- lm(FIM差~PT単位数, data=dat2)
summary(fit2)

予測通りにこのデータだと切片も有意になりました(p値=0.017)

#グラフの描き方
plot(FIM差~PT単位数, data=dat2)
x0 <- range(dat2$PT単位数) #x軸の範囲
b0 <- fit$coef[1] #切片
b1 <- fit$coef[2] #傾き
lines(x0, b0 + b1*x0) 

ほとんどのリハビリテーション研究では、複数の因子が結果に影響を与えるため、単回帰分析では結果を十分に考察することが困難です。特に、原因から結果を予測する場合、複数の交絡する要因が存在する可能性があります。これらの因子を統計学では交絡因子(confounder、confounding factor)と呼び、独立変数と目的変数の両方に影響を与えることで結果の解釈を歪める可能性があります。交絡因子の影響を調整するためには、単変量解析だけでなく、層別解析、傾向スコア、多変量解析などの手法を用いた十分な吟味が必要です。

繰り返しになりますが、ここでは交絡因子は無視して「PT単位数でFIM差を予測する」という、あり得ない仮定で進めていますのでご注意ください。

外挿に注意

ただし、今回取得したデータの範囲外(外挿:Extrapolation)の予測は、誤った予測値になる可能性があるので注意が必要です

統計モデル

線形単回帰モデル (linear simple regression model)

単回帰分析は、説明変数 $x_i$​ と誤差項 $e_i$ を用いて目的変数 $y_i$ を説明する分析手法です

$y_i = \beta_0 + \beta_1x_i + e_i \quad e_i \thicksim N(0,\sigma^2)$

このモデルでは、$\beta_0$ は切片、$\beta_1$​ は傾き($x_1$​ の効果)を表します。誤差項 $e_i$ は、通常、「平均 0 と一定の分散を持ち、独立な正規分布に従う $(iid)$」と仮定されます。これらの仮定のもとで、最小二乗法などの手法を用いて $\beta_0​$ と $\beta_1​$ をデータから推定します。この分析により、説明変数が目的変数にどのように影響を与えるかを定量的に評価できます。

$iid$

独立同分布 (independent and identically distributed)

最小二乗法

パラメタ推定には主に最小二乗法と最尤法があります.ここでは最小二乗法によるパラメタ推定を紹介します.統計モデルによって与えられた予測式を以下のように書きます.

$y=b_0+b_1x$

これは見慣れた式だと思います

例題(dat)では上述したように

$y=2.86+0.50x$

となりました

この式の2.86と0.50をどのようにして求めるかというと・・・

$\sum(y-(b_0+b_1x))^2$

を$b_0$、$b_1$で偏微分して=0となるような連立方程式を解きます

$b_1$で微分した場合

$2\sum(b_1x^2-x_iy_i+b_0x_i)=0$

$b_0$で微分した場合

$-2\sum(y-(b_0+b_1x))=0$

これは正規方程式と呼ばれています

\(
\left(
\begin{array}{ccc}

\sum_{i=1}^nx_i^{2*1} & \sum_{i=1}^nx_i^{2*1-1} \\

\sum_{i=1}^nx_i^{2*1-1} & \sum_{i=1}^nx_i^{2*1-2} \\


\end{array}\right)
\left(\begin{array}{cc}
b_1 \\
b_0 \\
\end{array}\right)
=
\left(\begin{array}{cc}
\sum_{i=1}^nx_i^1y_i \\
\sum_{i=1}^nx_i^{1-1}y_i \\
\end{array}\right)
\)

上記の式は以下のような2式になります

$b_1\sum{x_i}^2+b_0\sum{x_i}=\sum{x_iy_i}$

$b_1\sum{x_i}+nb_0=\sum{y_i}$

この2式から$b_0$と$b_1$を求めます

$b_1=\dfrac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}$

$b_0=\bar{y}-\bar{x}b_1$

回帰係数と相関係数の関係

相関係数

$r=\dfrac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum(x_i-\bar{x})^2\sum(y_i-\bar{y})^2}}$

もっと正確に書くと

$r=\dfrac{\dfrac{\sum(x_i-\bar{x})(y_i-\bar{y})}{n-1}}{\sqrt{\dfrac{\sum(x_i-\bar{x})^2}{n-1}\dfrac{\sum(y_i-\bar{y})^2}{n-1}}}$

となります・・・つまり共分散を標準化した値が相関係数です

xとyが共有する情報の割合を示します

回帰係数

$b_1=\dfrac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}$

相関係数と似てますね・・・

以下のように書くこともできます

$b_1=\dfrac{\dfrac{\sum(x_i-\bar{x})(y_i-\bar{y})}{n-1}}{\dfrac{\sum(x_i-\bar{x})^2}{n-1}}=\dfrac{cov(x, y)}{V[x]}$

yが原因、xが結果を想定した回帰係数は以下のようになります

$b_1^,=\dfrac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(y_i-\bar{y})^2}$

決定係数 (相関係数$^2$)

ここで$b_1$と$b_1^,~,$を掛けた値が決定係数になることが分かります

$b_1b_1^,=\dfrac{(\sum(x_i-\bar{x})(y_i-\bar{y}))^2}{\sum(x_i-\bar{x})^2\sum(y_i-\bar{y})^2}$

回帰係数と相関係数の関係

$b_1=\displaystyle r\frac{\sqrt{\sum(y_i-\bar{y})^2}}{\sqrt{\sum(x_i-\bar{x})^2}}$

つまり相関係数と回帰係数の関係は

$r=\sqrt{b_1b_1^,}$

相関係数はxとyのどちらもお互いに影響を与え合っている程度を表す指標です

詳しい数式はいいから、「r=0.4~0.7で強い相関」と覚えなさい・・・というようなこともよく耳にしますが、理解できる人は理解した方がよいと思います。なんでもかんでも”0.5=強”というわけにはいかないときもありますので・・・。統計ソフトがはじき出す結果をちゃんと理解できるようになりましょう!

統計学入門−第5章

おまけ1:線形代数を利用した計算方法

重回帰のときにとても便利です

表記の定義:$x_{ij}$の場合、$i$はデータセットの行、$j$はデータセットの列を示す

$y_{i}=b_0+b_1x_{ij}+\dotsc+b_jx_{ij}+\varepsilon_{i} \quad i={1,2,\dotsc,n} \quad j={1,2,\dotsc,k}$

\(
\mathbf{X}=
\left(\begin{array}{ccc}
1 & x_{11} & \dotsc & x_{1k} \\
1 & x_{21} & \dotsc & x_{2k} \\
\vdots & \vdots & \ddots & \vdots \\
1 & x_{n1} & \dotsc & x_{nk}
\end{array}\right)
\)

\(
\mathbf{y}=
\left(\begin{array}{cc}
y_1 \\
y_2 \\
\vdots \\
y_n
\end{array}\right)
\quad\)\(
\mathbf{b}=
\left(\begin{array}{cc}
b_0 \\
b_1 \\
\vdots \\
b_m
\end{array}\right)
\quad\)\(
\mathbf{0}=
\left(\begin{array}{cc}
0 \\
0 \\
\vdots \\
0
\end{array}\right)
\)

とした場合、パラメタ推定の公式は以下のようになります

$\mathbf{b}=(\mathbf{X^T}\mathbf{X})^{-1}\mathbf{X^T}\mathbf{y}$

もちろん単回帰分析でも使用可能です

今回のデータ(lm01)を当てはめて解いてみましょう

dat$b <- c(rep(1, dim(dat)[1]))
df <- dat[ ,c("b",  "PT単位数")]
mat <- as.matrix(df)

solve(
    t(mat) %*% mat
    ) %*% 
    t(mat) %*% 
    dat$FIM差

Rの関数で算出した係数と同じ値になりましたね!

おまけ2:係数検定のための標準誤差

$Y = \beta_0 + \beta_1x_i + e_i$

$\beta_0$ の推定値を $b_0$ 、$\beta_1$ の推定値を $b_1$ とすると

$E[b_0]=E[\bar{y}-\bar{x}b_1]=\beta_0$

$E[b_1]=E[\dfrac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}]=\beta_1$

$V[b_0]=\sigma^2(\dfrac{1}{n}+\dfrac{\bar{x}^2}{\sum(x_i-\bar{x})^2})$

$V[b_1]=\dfrac{\sigma^2}{\sum(x_i-\bar{x})^2}$

$\sigma^2$ の推定 $s^2$

$s^2=\dfrac{\sum{(y_i-\hat{y_i})^2}}{n-2} \qquad$ 自由度=$n-$パラメタ数(単回帰の場合は2)

$\beta_0$ の検定

Std.Error = $\sqrt{V[b_0]}=\sqrt{s^2(\dfrac{1}{n}+\dfrac{\bar{x}^2}{\sum(x_i-\bar{x})^2})}$

$\beta_1$ の検定

Std.Error = $\sqrt{V[b_1]}=\sqrt{\dfrac{s^2}{\sum(x_i-\bar{x})^2}}$

t統計量

$\dfrac{b}{Std.Error} \qquad$ 自由度=$n-$パラメタ数(単回帰の場合は2)

コメント欄 『間違い』や『分かりにくい部分』などのご意見もお寄せください

タイトルとURLをコピーしました