看護応用top

市原清志先生「バイオサイエンスの統計学」をRで実装したものです。
市原先生は直線回帰・相関の順に説明しておられます。両者の違いが重要だと認識されていたのだと思います。

回帰 相関
誤　差 y方向のみ x,y両方が確率変数
信頼域 yの信頼区間 x,yの分布の中心からの等確率距離
目　的 どのように直線に関係づけられるか どの程度直線的か
単　位 単位に依存 単位に依存しない

相関係数.

\({r_{xy} = \dfrac{\sum_{i = 1}^n (x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\sum_{i = 1}^n (x_i - \bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i - \bar{y})^2}} =\dfrac{1}{n}\sum_{i=1}^n \dfrac{x_i - \bar{x}}{s_x} \cdot \dfrac{y_i - \bar{y}}{s_y} }\)

1 直線回帰.

回帰係数.

\(\displaystyle b=\frac{\sum ^{n}_{i=1}\left( x_{i}-\bar{x}\right) \left( y_{i}-\bar{y}\right) }{ \sum ^{n}_{i=1}\left( x_{i}-\bar{x}\right) ^{2}} = \dfrac{\sqrt{\sum_{i = 1}^n (y_i - \bar{y})^2}}{{\sqrt{\sum_{i = 1}^n (x_i - \bar{x})^2}}}\cdot\displaystyle \frac{\sum ^{n}_{i=1}\left( x_{i}-\bar{x}\right) \left( y_{i}-\bar{y}\right) }{ \sqrt{\sum_{i = 1}^n (x_i - \bar{x})^2}\sqrt{\sum_{i = 1}^n (y_i - \bar{y})^2}}=\dfrac{s_y}{s_x}r_{xy}\)
重要なのは分子の共分散で、分母は単位を調整するためのものです。標準化して無名数にしておけば、相関係数と回帰係数は一致します。

相関係数と回帰係数の違い.

  • 変数が予測値にどの程度影響するかを表すのが回帰係数.
    （yの予測値なので楕円の青い縦の弦の中点の軌跡）
  • 回帰直線にどの程度集まっているかを表すのが相関係数.
    （楕円の扁平率、つぶれの程度）

1.1 回帰係数を帰納.

1.1.1 

x <- c(-2,2)
y <- c(-2,2)
plot(x,y,asp = 1,xlim=c(-3,3),ylim=c(-3,3))
abline(h=0,v=0,col="lightblue")

(x-mean(x))*(y-mean(y))#偏差積
## [1] 4 4
(x-mean(x))*(x-mean(x))#偏差平方
## [1] 4 4
sum((x-mean(x))*(y-mean(y)))/sum((x-mean(x))*(x-mean(x)))
## [1] 1

1.1.2 

x <- c(-2,2)
y <- c(2,-2)
plot(x,y,asp = 1,xlim=c(-3,3),ylim=c(-3,3))
abline(h=0,v=0,col="lightblue")

(x-mean(x))*(y-mean(y))#偏差積
## [1] -4 -4
(x-mean(x))*(x-mean(x))#偏差平方
## [1] 4 4
sum((x-mean(x))*(y-mean(y)))/sum((x-mean(x))*(x-mean(x)))
## [1] -1

1.2 乱数による単回帰モデル.

S=c()
DATA=data.frame(X=c(),Y=c())
X <- c(0,1,2,3,4)
for (i in 1:100){
  e <- rnorm(5,0,0.5)
  Y <- 0.5*X+1+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  f <- fit$ fitted.values#予測値
  a <- fit$ coefficients;a#回帰係数
  e <- fit$ residuals;e#残差
  se <- sum(e^2)#残差平方和
  s <- sqrt(se/3)#se/2#残差分散
  S <- c(S,s)
  DATA <- rbind(DATA, data)
}
mean(S)
## [1] 0.4520193

1.2.1 回帰直線.

model = lm( formula = Y ~ X, data = DATA )
plot(DATA,asp=1,col="grey")
abline(model)

1.3 例7-1（p208）.

par(family = "HiraKakuProN-W3") #日本語フォントの指定
data = matrix(0,10,2) 
data[,1]<-c(1,3,4,5,6,7,9,10,11,12)
data[,2]<-c(4,4,5,6,5,7,6,8,9,7)
colnames(data)<- c("X","Y") 
plot(data)

dat<- data.frame(data) 
dat
X Y
1 4
3 4
4 5
5 6
6 5
7 7
9 6
10 8
11 9
12 7 
glimpse(dat)
## Rows: 10
## Columns: 2
## $ X <dbl> 1, 3, 4, 5, 6, 7, 9, 10, 11, 12
## $ Y <dbl> 4, 4, 5, 6, 5, 7, 6, 8, 9, 7
#データフレーム形式(行・列はラベルを必ず持ち， ラベルによる操作が可能である
dat.lm<- lm(Y~X,data=dat) 
dat.lm
## 
## Call:
## lm(formula = Y ~ X, data = dat)
## 
## Coefficients:
## (Intercept)            X  
##      3.4164       0.3946
par(family = "HiraKakuProN-W3") #日本語フォントの指定
par(mfrow=c(2,2))
plot(dat.lm)

summary(dat.lm) 
## 
## Call:
## lm(formula = Y ~ X, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.15217 -0.73829  0.09699  0.63043  1.24247 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.41639    0.61769   5.531 0.000553 ***
## X            0.39465    0.08097   4.874 0.001233 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8855 on 8 degrees of freedom
## Multiple R-squared:  0.7481, Adjusted R-squared:  0.7166 
## F-statistic: 23.76 on 1 and 8 DF,  p-value: 0.001233
residuals(dat.lm)#残差
##            1            2            3            4            5            6 
##  0.188963211 -0.600334448  0.005016722  0.610367893 -0.784280936  0.821070234 
##            7            8            9           10 
## -0.968227425  0.637123746  1.242474916 -1.152173913
sqrt(sum(residuals(dat.lm)^2)/8)#標準偏差
## [1] 0.8854783
#taikei.lm$residuals でもよい
coefficients(dat.lm)#回帰係数
## (Intercept)           X 
##   3.4163880   0.3946488
plot(residuals(dat.lm) )
abline(h=0,col="lightblue")

1.4 残差平方和.

Sx <-var(dat$X)*9;Sx
## [1] 119.6
Sy <-var(dat$Y)*9;Sy
## [1] 24.9
Sxy <-cov(dat$X,dat$Y)*9
sqrt((Sy-coefficients(dat.lm)[2]*Sxy)/8)
##         X 
## 0.8854783

テキストの標準偏差と少しずれているが、四捨五入による誤差と考えられる。

1.5 注意.

p211の図でわかるように、回帰直線は.

  • 縦の弦の中点の軌跡.
  • 縦線と楕円の接点を結ぶ直線.

x,y の尺度が同じとき、
yを回帰する直線が　y=ax+b　であれば、
xを回帰する直線も　x=ay+b　である（同じ傾き）。

  • \(x=\dfrac{1}{a}y\)　とはならないことに注意しておこう。

2 回帰の有意性の検定.

回帰分析から分散分析へ.

1.1.2 aov による回帰分析 aov は分散分析の関数です。R は因子型（factor）の説明変数に対し、特性関数（indicator function）を用いて回帰分析を行い、分散分析として出力しているそうです。このため、水準数が４ならば、使う特性関数は３つ、もう一つの水準は x=0 のときの切片の値として使うので、自由度は 4-1=3 、分散分析の場合の自由度が 水準数−１であるのは、こういう理由によります。 説明変数 x を数値型（numeric）として aov を適応すると、回帰分析が行われ、説明変数の自由度は１になります。 （factor(x)として因子型にすると、自由度は因子数−１になります）

大森宏「回帰分析」

3 Advanced.

以下の内容は、2変量正規分布を扱うので高校のレベルを超える。
ここでは項目だけを列挙し、詳しい説明は次のページで行う。

4 回帰係数の差の検定.

共分散分析（ANCOVA）の第１段階に相当する。

5 偏相関係数.