1 単回帰分析

回帰分析を、散布図の直線への当てはめだと考えることは初歩の理解としては有効ですが、統計的理解に進むには乱数によるシミュレーションなどが必要になります。
　　　「誤差と残差の違いは？」
この問いに、どう応えればよいでしょうか。
誤差は確率変数（確率モデル）であり、残差はその実現値ですが、このことを理解するにはモデルを用いたシミュレーションが有効です。

1.1 単回帰モデル

回帰分析では、説明変数は定数として固定して考えます。 真の回帰式として、
　　　\(y=1+0.5 x\)　を考え、誤差変動として \(\varepsilon_i \sim N(0,~0.5^2) ~:~ i.i.d.\) を加えた
　　　\(y_i=1+0.5 x_i+\varepsilon_i\)　　を単回帰モデルとします。

\(X=(x_i)=(0,1,2,3,4)\) としてYを出力させ、Rの線形回帰モデル lm（linear　model）によって 回帰係数を求めてみます。

X <- c(0,1,2,3,4)
e <- rnorm(5,0,0.5)
Y <- 0.5*X+1+e
data <- data.frame(X=X,Y=Y);data
X Y
0 1.373818
1 1.823634
2 2.153791
3 2.677286
4 3.408355 
plot(data,ylim=c(0,4))

誤差項があるため plot は直線から少しずれています。
この生成されたデータに当てはまる回帰直線を関数 lm によって 求めます。

fit = lm( formula = Y ~ X, data = data )
a <- fit$ coefficients#回帰係数
e <- fit$ residuals#残差
se <- sum(e^2)#残差平方和
#se/2#残差分散
sqrt(se/3)#残差の標準誤差
## [1] 0.132653
a[2]#回帰係数の推定値
##         X 
## 0.4922726
summary(fit)
## 
## Call:
## lm(formula = Y ~ X, data = data)
## 
## Residuals:
##        1        2        3        4        5 
##  0.07099  0.02853 -0.13359 -0.10236  0.13643 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  1.30283    0.10275   12.68  0.00106 **
## X            0.49227    0.04195   11.73  0.00133 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1327 on 3 degrees of freedom
## Multiple R-squared:  0.9787, Adjusted R-squared:  0.9716 
## F-statistic: 137.7 on 1 and 3 DF,  p-value: 0.00133

残差から求めた標準誤差と、summaryの”Residual standard error”は一致しています。

1.2 回帰の有意性

回帰が有意であるとは、どういうことでしょうか。
Rで回帰分析を行うと t値とF値が出力されます。
回帰の有意性を検定するための統計量で、次の帰無仮説・対立仮説のもとで算出されています。

  • H0：説明変数Xは被説明変数Yに影響を与えていない（回帰係数は０）
  • H1：説明変数Xは被説明変数Yに影響を与えている

帰無仮説のもとで生成されたデータでも、回帰係数の推定値は０になりません。
誤差変動の範囲を超えるほどの回帰係数かどうかが問題です。
つまり、今手にしたデータ \(Y=(y_1,y_2,\cdots,y_n)\) は、帰無仮説のもとでの無限回の実験で得られるデータの１つであり、統計量t、Fが有意差の範囲まで行っているかどうか判断するのが回帰の有意性の検定なのです。

「回帰係数は０」という帰無仮説のもとで１回シミュレーションします。

X <- c(0,1,2,3,4)
e <- rnorm(5,0,0.5)
Y <- 0*X+1+e
data <- data.frame(X=X,Y=Y);data
X Y
0 0.7249239
1 1.1889097
2 0.6470360
3 1.3867583
4 0.7227802 
plot(data,ylim=c(0,4))

fit = lm( formula = Y ~ X, data = data )
a <- fit$ coefficients;a#回帰係数
## (Intercept)           X 
##  0.89536936  0.01935613
e <- fit$ residuals#残差
f <- fit$ fitted.values#予測値
sr <- var(f)*4#回帰変動
se <- sum(e^2)#残差平方和
sqrt(se/3)#残差の標準誤差
## [1] 0.3816182
a[2]#回帰係数の推定値
##          X 
## 0.01935613
F <- var(f)*4/(se/3);F#F値
## [1] 0.0257264
t2 <- a[2]/sqrt(se/3)*sqrt(10);t2#t値
##         X 
## 0.1603945
t1 <- (a[1])/(sqrt(se/3)*sqrt(0.6));t1#t値
## (Intercept) 
##    3.028988
summary(fit)
## 
## Call:
## lm(formula = Y ~ X, data = data)
## 
## Residuals:
##       1       2       3       4       5 
## -0.1704  0.2742 -0.2870  0.4333 -0.2500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.89537    0.29560   3.029   0.0564 .
## X            0.01936    0.12068   0.160   0.8828  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3816 on 3 degrees of freedom
## Multiple R-squared:  0.008503,   Adjusted R-squared:  -0.322 
## F-statistic: 0.02573 on 1 and 3 DF,  p-value: 0.8828

F値とt値が、summaryの値と一致していることがわかります。
t1,t2は回帰係数が０を帰無仮説としたときのt値です。

1000回シミュレーションを繰り返してみます。

X <- c(0,1,2,3,4)
Se = c();Sr = c()
Tf = c();Tt = c()
for (i in 1:1000) {
  e <- rnorm(5,0,0.5)
  Y <- 0*X+1+e
  data <- data.frame(X=X,Y=Y)
  fit = lm( formula = Y ~ X, data = data )
  a <- fit$ coefficients#回帰係数
  e <- fit$ residuals#残差
  f <- fit$ fitted.values#予測値
  sr <- var(f)*4#回帰変動
  se <- sum(e^2)#残差平方和
  sqrt(se/3)#残差の標準誤差
  #a[2]#回帰係数の推定値
  f <- var(f)*4/(se/3)#F値
  t <- a[2]/sqrt(se/3)*sqrt(10)#t値
  Tf = c(Tf,f)
  Tt = c(Tt,t)
  #sr <- var(f)*3/(0.5)^2
  #se <- sum(e^2)/(0.5)^2#残差平方和
  Se <-c(Se,se)
  Sr <- c(Sr,sr)
}
hist(Tf,prob=T)
curve(df(x,1,3),add=T)

hist(Tt,prob=T)
curve(dt(x,3),add=T)

回帰が有意でない、即ち、説明変数は目的変数に影響を与えないとしたとき、
　　　F値は \(F(1,3\)) 　自由度(1,3)のF分布
　　　t 値は \(t(3)\) 　 自由度３のt分布
に従うことが図から確かめられます。

また、偏差平方和を母分散で割った値がカイ二乗分布に従うことも確かめられます。
　　　\(\dfrac{S_e}{(0.5)^2}\sim \chi^2(3)\)　　　（期待値は 3 ）
　　　\(\dfrac{S_R}{(0.5)^2}\sim \chi^2(1)\)　　　（期待値は 1 ）
シミュレーションの結果は、 　

mean(Se/(0.5)^2)
## [1] 2.979196
mean(Sr/(0.5)^2)
## [1] 1.013191

　　　\(\dfrac{S_e}{(0.5)^2}\) の平均は 3 　　　　\(\dfrac{S_R}{(0.5)^2}\) の平均は 1
　　　
理論上の平均と一致しています。

2 重回帰分析

3D機能を重回帰分析に利用しましょう。
重回帰分析でも単回帰と同様、説明変数は定数として固定して考えます。

2.1 ３次元プロット

Rには3Dプロットの関数があります。Rに装備されているデータセット iris を3D表示してみます
。 irisはアヤメという花で、setosa、versicolor、virginica　の3つの種類があります。 3種類のアヤメ50個体(合計150）について、がく片(sepal)の長さ、幅、および花弁(petal)の長さ、幅を測定したデータがirisです。
　Sepal.Length がく片の長さ [cm]
　Sepal.Width がく片の幅 [cm]
　Petal.Length 花弁の長さ [cm]
3つの種類によってマーカーを変えて表示したものが下図になります。

Rとウェブの融合（６）ーrglによる3次元グラフィックスー

# scatterplot3d による3次元プロット
library(scatterplot3d)
scatterplot3d(iris[, 1:3], pch = as.numeric(iris$Species))

#pch というオプションで点の形や塗りつぶしの有無

2.1.1 俯瞰図

# latticeを使った3次元プロット
library(lattice)
cloud(Petal.Length ~ Sepal.Length * Sepal.Width, data = iris, groups = Species)

2.2 重回帰モデル

真の回帰式として、
　　　\(z=1+0.2 x-0.3 y\)　を考え、誤差変動として \(\varepsilon_i \sim N(0,~0.2^2) ~:~ i.i.d.\) を加えた
　　　\(z_i=1+0.2 x_i-0.3 y_i+\varepsilon_i\)　　を回帰モデルとします。

xy平面の点として、次の9つの点を固定します。

library(tidyr)
data <- crossing(x = 1:3, y = 1:3);data
x y
1 1
1 2
1 3
2 1
2 2
2 3
3 1
3 2
3 3

\((x_i,y_i)\) に対し \(z_i=1+0.2 x_i-0.3 y_i+\varepsilon_i\)　により \(z_i\) を出力させ、3Dプロットに表示します。
次に Rの線形回帰モデル lm（linear　model）によって 回帰係数を求め回帰平面を表示します。

library(tidyr)
library(gapminder)
options(rgl.printRglwidget = TRUE) 
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
x <- data$x
y <- data$y
z <- 0.2*x - 0.3*y+ rnorm(9, sd = 0.2)
#plot3d(x = x, y = y, z = z,xlim=c(0,4),ylim=c(0,4),zlim=c(-3,3))
open3d()
## glX 
##   1
fit <- lm(z ~ x + y)
summary(fit)
## 
## Call:
## lm(formula = z ~ x + y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.19068 -0.09289 -0.04173  0.10183  0.28496 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.26777    0.20655   1.296  0.24246   
## x            0.21406    0.07016   3.051  0.02248 * 
## y           -0.39089    0.07016  -5.571  0.00142 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1719 on 6 degrees of freedom
## Multiple R-squared:  0.8705, Adjusted R-squared:  0.8274 
## F-statistic: 20.17 on 2 and 6 DF,  p-value: 0.002169
plot3d(x, y, z, type = "s", col = "blue", size = 1)

coefs <- coef(fit)
a <- coefs["x"]
b <- coefs["y"]
c <- -1
d <- coefs["(Intercept)"]
planes3d(a, b, c, d, alpha = 0.5)

重回帰モデルから作成したデータなので、当てはまりは非常によく、
t値、F値とも有意であり、決定係数も１に近い値となっています。

2.3 重回帰モデルと分散分析

２元配置分散分析を ２変量の重回帰分析によって分析し、3D表示することができます。

library(dplyr)
## 
##  次のパッケージを付け加えます: 'dplyr'
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union
library(tidyr)
library(gapminder)
options(rgl.printRglwidget = TRUE) 
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
d=c()
k1=c(0,1,2)
k2=c(0,1,2)
s1=scale(k1)
s2=scale(k2)
for (i in 1:3) {
  X_i = s1[i]
  Y_i = c()
  for (j in 1:3) {
    X_ij = 0.2*i - 0.3*j+rnorm(2,0,0.5)
    Y_i = c(Y_i,X_ij)
  }
  d=c(d,Y_i)
}
data=data.frame(x=c(rep(1,6),rep(2,6),rep(3,6)),y=c(rep(c(rep(1,2), rep(2,2), rep(3,2)),3)),z=d)
data
x y z
1 1 -0.6004066
1 1 1.4263017
1 2 -1.2181893
1 2 -1.5630003
1 3 -0.5864433
1 3 -1.3299555
2 1 0.8838414
2 1 1.1113425
2 2 0.0794811
2 2 -0.2206114
2 3 -1.1427435
2 3 -0.7042238
3 1 0.4611516
3 1 -1.1026303
3 2 0.5393277
3 2 0.3816845
3 3 -0.3873852
3 3 0.0566418 
x <- data$x
y <- data$y
z <- data$z

open3d()
## glX 
##   3
fit <- lm(z ~ x + y)
summary(fit)
## 
## Call:
## lm(formula = z ~ x + y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.72627 -0.36960  0.03536  0.46384  1.43941 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   0.1913     0.6575   0.291   0.7750  
## x             0.3184     0.2233   1.426   0.1745  
## y            -0.5228     0.2233  -2.341   0.0335 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7736 on 15 degrees of freedom
## Multiple R-squared:  0.3337, Adjusted R-squared:  0.2449 
## F-statistic: 3.756 on 2 and 15 DF,  p-value: 0.04759
plot3d(x, y, z, type = "s", col = "blue", size = 1)

coefs <- coef(fit)
a <- coefs["x"]
b <- coefs["y"]
c <- -1
d <- coefs["(Intercept)"]
planes3d(a, b, c, d, alpha = 0.5)

2.4 体重〜身長＋ウェスト

同志社大学文化情報学部 金明哲 先生の “Jin’s page” フリーソフトによるデータ解析・マイニング　第15回 「Rと重回帰分析」に、体重を身長とウェストから回帰する例題があります。
この例題の重回帰分析を3Dで空間図示します。
(現在、Jin’s pageは閲覧できません)

g<- c(50,60,65,65,70,75,80,85,90,95) 
t<- c(165,170,172,175,170,172,183,187,180,185) 
w<- c(65,68,70,65,80,85,78,79,95,97) 
taikei2<- data.frame(体重=g,身長=t,ウエスト=w) ;taikei2
体重 身長 ウエスト
50 165 65
60 170 68
65 172 70
65 175 65
70 170 80
75 172 85
80 183 78
85 187 79
90 180 95
95 185 97 
res <- lm(体重~身長+ウエスト,data=taikei2)
res
## 
## Call:
## lm(formula = 体重 ~ 身長 + ウエスト, data = taikei2)
## 
## Coefficients:
## (Intercept)         身長     ウエスト  
##   -161.6698       1.0217       0.7091
library(rgl)
options(rgl.printRglwidget = TRUE) 
plot3d(t, w, g, type = "s", col = "blue", size = 1)
coefs <- coef(res)
a <- coefs["身長"];a
##     身長 
## 1.021723
b <- coefs["ウエスト"];b
##  ウエスト 
## 0.7090631
c <- -1;c
## [1] -1
d <- coefs["(Intercept)"];d
## (Intercept) 
##   -161.6698
planes3d(a, b, c, d, alpha = 0.5)