僕が作った架空のデータで練習しましょう。結果についてはご意見があると思いますが、統計学の練習のためのサンプルですのでご了承ください。
例題）FIM運動項目が60±5点の患者を対象に、理学療法実施単位数（PT単位数）とFIM運動項目の差（FIM差；PT後FIM-PT前のFIM）の関連性について検証した。この研究の目的は、PT単位数（原因）でFIM差（結果）を予測することです。
データの準備と要約
データセットlm01をdatに格納します（ファイルの読み込み）
head(dat)
散布図
plot(FIM差~PT単位数, data=dat)
Rで単回帰分析
lmはlinear model（linear regression model）の略です
(fit <- lm(FIM差~PT単位数, data=dat))
これで以下のような回帰式が求められます
$y=2.86+0.50x$
切片と傾きは、有意な数値なのか有意水準5％で検定してみましょう
summary(fit)
傾きのｐ値（0.0199）は0.05未満なので有意です
PT単位数が1単位増えるとFIM差は約0.5増加することが予測できました
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が１名いるので、切片が十分に予測できないようですね
試しにFIM差22を15に置き換えた以下のデータでやってみましょう
注意）あくまでも統計学の勉強のためのデータ操作です
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)
$y_i=\beta_0+\beta_1x_i+\varepsilon_{i} \qquad \varepsilon_{i} \thicksim N(0,\sigma^2)$
このモデルを利用してパラメタ（$\beta_0, \beta_1$）を予測する分析を単回帰分析と呼んでいます
最小二乗法
パラメタ推定には主に最小二乗法と最尤法があります
ここでは最小二乗法によるパラメタ推定を紹介します
統計モデルによって与えられた予測式を以下のように書きます
$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}$
この２式から$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}}}$
となります・・・つまり共分散を標準化した値が相関係数です
ｘとｙが共有する情報の割合を示します
回帰係数
$b_1=\dfrac{\sum(x_i-\bar{x})(y_i-\bar{y})}{\sum(x_i-\bar{x})^2}$
相関係数と似てますね・・・
ただしこれはｘが原因でｙが結果を想定した回帰係数です
ｙが原因、ｘが結果を想定した回帰係数は
$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^,}$
相関係数はｘとｙのどちらもお互いに影響を与え合っている程度を表す指標です
詳しい数式はいいから、「r=0.4～0.7で強い相関」と覚えなさい・・・というようなこともよく耳にしますが、理解できる人は理解した方がよいと思います。なんでもかんでも”0.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）
*証明は後日アップする予定です
コメント