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

回帰分析・分散分析、Rはすぐに結果を返してくれます。でも少しも分かった感じがしませんでした。

今、手元にあるデータは、無限回の調査を実施したとして、そのうちの１回に過ぎません。 無限回の結果がイメージできていないのに、今あるデータが異常かどうか判断できるはずもありません。

• 帰無仮説モデルで1000回程度のシミュレーションを行い、結果がどのように変動するか把握する.
  
• 対立仮説モデルで1000回程度のシミュレーションを行い、結果がどのように変動するか把握する.
  
• 標本数や誤差の大きさを変えると、出力にどの程度影響するか把握する.

このように、モデルからデータを生成し、aov や linear model などの関数によりモデルが正確に復元できるかみておく必要があります。
つまり、Rという機械が正しく働くか、自分の手で確かめておくということです。

1 重回帰モデル

1.1 ２変量回帰分析.

真の回帰式として、
　　　\(z=1+0.2 x-0.3 y\)　を考え、誤差変動として \(\varepsilon_i \sim N(0,~0.1^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
data <- rbind(data,data,data)#３回繰り返す
data[1:9,]

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(27, sd = 0.1)#誤差分散 0.1^2
plot3d(x = x, y = y, z = z,type = "s", col = "blue", size = 1)

#open3d()
fit <- lm(z ~ x + y)
f <- fit$ fitted.values;#予測値
a <- fit$ coefficients;#回帰係数
e <- fit$ residuals;#e#残差
se <- sum(e^2);se#残差平方和

## [1] 0.1764835

s <- sqrt(se/(27-3));s#残差標準偏差

## [1] 0.08575244

ff<- (var(f)*26/2)/(se/24);ff#F値

## [1] 140.7933

summary(fit)

## 
## Call:
## lm(formula = z ~ x + y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14400 -0.04838  0.01370  0.03735  0.18784 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02485    0.05950  -0.418     0.68    
## x            0.19061    0.02021   9.430 1.53e-09 ***
## y           -0.28054    0.02021 -13.880 5.81e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08575 on 24 degrees of freedom
## Multiple R-squared:  0.9215, Adjusted R-squared:  0.9149 
## F-statistic: 140.8 on 2 and 24 DF,  p-value: 5.507e-14

plot3d(x, y, z, type = "s", col = "blue", size = 1)

coefs <- coef(fit)
a <- coefs["x"]
b <- coefs["y"]
c <- -1
d <- coefs["(Intercept)"]
#ax+by＋cz+d=0で表される平面
planes3d(a, b, c, d, alpha = 0.5)

重回帰モデルから作成したデータですが、誤差変動を含むため、
　　　回帰係数　a=0.2　b=-0.3　とはなっていません。
回帰係数の推定値 \(\widehat{\beta_i}\) と真値 \(\beta_i\) の差を標準誤差（Std. Error）で割った値は自由度 \(n-3=27-3=24\) のt分布に従います。即ち、
　　　\(t=\dfrac{\widehat{\beta}_i-\beta_i}{std.error} \sim t(24)\)

summary のt値は、帰無仮説 \(\beta_i=0\) としたときの t値なので、
　　　\(t=\dfrac{\widehat{\beta}_i}{std.error} \sim t(24)\)　（推定値を標準誤差で割った値が記載されています）

回帰係数を区間推定する場合は、帰無仮説でない式を用いて
　　　\({\beta_i}=\widehat{\beta_i}\pm 2.064\cdot (Std.error)\)

よって、推定区間は.

summaryReg <- summary(fit)
testResults <- summaryReg$coefficients 
lower <- testResults[, "Estimate"] -2.064*testResults[, "Std. Error"]
upper <- testResults[, "Estimate"] +2.064*testResults[, "Std. Error"]
kukan <- data.frame(
  lower = testResults[, "Estimate"]-2.064*testResults[, "Std. Error"],
  upper = testResults[, "Estimate"] +2.064*testResults[, "Std. Error"])
kukan

	lower	upper
(Intercept)	-0.1476678	0.0979592
x	0.1488894	0.2323247
y	-0.3222611	-0.2388258

95%信頼区間は、真の回帰係数 0.2 , -0.3 を含んでいます。

回帰がないという帰無仮説に対する t値、F値 は有意であり、決定係数も１に近い値となっています。

1.1.1 １つの変数を固定

coefs

## (Intercept)           x           y 
## -0.02485431  0.19060709 -0.28054346

a;b;d

##         x 
## 0.1906071

##          y 
## -0.2805435

## (Intercept) 
## -0.02485431

curve(a*x+b*1+d,ylim=c(-1.1,0.1),col="blue")
curve(a*x+b*2+d,add=T,col="pink")
curve(a*x+b*3+d,add=T,col="green")

回帰平面を、\(y=1,y=2,y=3\) で切った切り口は、平行な直線になります。

1.2 パス解析

data <- data.frame(x=data$x,y=data$y,z=z)
str(data)

## 'data.frame':    27 obs. of  3 variables:
##  $ x: int  1 1 1 2 2 2 3 3 3 1 ...
##  $ y: int  1 2 3 1 2 3 1 2 3 1 ...
##  $ z: num  -0.137 -0.499 -0.72 0.123 -0.313 ...

attach(data)

## The following objects are masked _by_ .GlobalEnv:
## 
##     x, y, z

model01 <-  '
  z ~ x + y '
par(family = "HiraKakuProN-W3")  
sem01 <- sem(model01, data = data, estimator = "ML")
semPaths(sem01, what = "est", style = "lisrel", layout = "tree", rotation = 2, nCharNodes = 0, nCharEdges = 0, fade = FALSE, optimizeLatRes = TRUE, edge.width = 0.2, label.scale = FALSE, label.cex = 0.6, theme = 'gray', asize = 6.0, node.width = 1.5)

1.2.1 ２変量分散分析

aov は分散分析ですが、説明変数が数値型の場合は回帰分析を実行します。 xという要因（変数）が有意であるかの判定は、
- 回帰分析では回帰係数のt値、分散分析では回帰変動のF値が出力されます。
\(t^2=F\) となっており、判定は同じです。

res <- aov(z ~ x + y)
summary(res)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## x            1 0.6540  0.6540   88.93 1.53e-09 ***
## y            1 1.4167  1.4167  192.66 5.81e-13 ***
## Residuals   24 0.1765  0.0074                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

interaction.plot(x,y,z, type="b",  pch=c(1,2))

interaction.plotは交差しません。

1.3 交互作用モデル

モデルにx,yの積がある場合、

#data,x,yは継承
z <- 0.2*x - 0.3*y+ x*y+rnorm(27, sd = 0.2)#誤差分散 0.2^2
plot3d(x = x, y = y, z = z,type = "s", col = "blue", size = 1)

#open3d()
res <- aov(z ~ x * y)
summary(res)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## x            1  85.15   85.15    2446  < 2e-16 ***
## y            1  47.43   47.43    1363  < 2e-16 ***
## x:y          1  12.67   12.67     364 1.35e-15 ***
## Residuals   23   0.80    0.03                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fit <- lm(z ~ x * y)
f <- fit$ fitted.values;#予測値
a <- fit$ coefficients;a#回帰係数

## (Intercept)           x           y         x:y 
##   0.2938350   0.1196235  -0.4320214   1.0276775

e <- fit$ residuals;#e#残差
se <- sum(e^2);se#残差平方和

## [1] 0.8006946

s <- sqrt(se/23);s#残差標準偏差

## [1] 0.1865819

ff<- (var(f)*26/3)/(se/23);ff#F値

## [1] 1390.837

summary(fit)

## 
## Call:
## lm(formula = z ~ x * y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29842 -0.09923 -0.00170  0.11728  0.37950 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.29384    0.25135   1.169  0.25437    
## x            0.11962    0.11635   1.028  0.31460    
## y           -0.43202    0.11635  -3.713  0.00114 ** 
## x:y          1.02768    0.05386  19.080 1.35e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1866 on 23 degrees of freedom
## Multiple R-squared:  0.9945, Adjusted R-squared:  0.9938 
## F-statistic:  1391 on 3 and 23 DF,  p-value: < 2.2e-16

plot3d(x, y, z, type = "s", col = "blue", size = 1)

coefs <- coef(fit)
a <- coefs["x"];a

##         x 
## 0.1196235

b <- coefs["y"];b

##          y 
## -0.4320214

c <- coefs["x:y"];c

##      x:y 
## 1.027677

d <- coefs["(Intercept)"];d

## (Intercept) 
##    0.293835

df <- data.frame(x=data$x,y=data$y)
X <- seq(min(df$x),max(df$x),len=100)#定義域方眼
Y <- seq(min(df$y),max(df$y),len=100)#定義域方眼
plot.df <- expand.grid(x=X,y=Y)#定義域方眼
#定義域方眼の座標を(A,B,C)とする
A <- plot.df$x
B <- plot.df$y
C <- a*A+b*B+c*A*B#z座標の列を追加
plot3d(A,B,C,col="lightblue",
       alpha = 0.2,xlab="x",ylab="y",zlab="z")
plot3d(x, y, z, type = "s", col = "blue", size = 1,add=T)

curve(a*x+b*1+c*x,xlim=c(-0.2,3.2),ylim=c(-0.2,2),col="blue")
curve(a*x+b*2+2*c*x,add=T,col="pink")
curve(a*x+b*3+3*c*x,add=T,col="green")

回帰平面を、\(y=1,y=2,y=3\) で切った切り口の直線は平行にならない。 このとき、多重共線性があるという。

res <- aov(z ~ x * y)
summary(res)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## x            1  85.15   85.15    2446  < 2e-16 ***
## y            1  47.43   47.43    1363  < 2e-16 ***
## x:y          1  12.67   12.67     364 1.35e-15 ***
## Residuals   23   0.80    0.03                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

interaction.plot(x,y,z, type="b",  pch=c(1,2))

1.3.1 帰無仮説のもとで

S2=c();Vf=c();FF=c()
for (i in 1:1000){
  x <- data$x
  y <- data$y
  z <- 0*x - 0*y+ rnorm(27, sd = 0.2)#誤差分散 0.2^2
  #plot3d(x = x, y = y, z = z,type = "s", col = "blue", size = 1)
  #open3d()
  fit <- lm(z ~ x + y)
  f <- fit$ fitted.values;#f#予測値
  e <- fit$ residuals;#e#残差
  se <- sum(e^2);se#残差平方和
  s2 <- se/24;#残差分散
  S2=c(S2,s2)
  vf<- var(f)*26/2
  Vf=c(Vf,vf)
  FF=c(FF,vf/s2)
}
hist(S2,prob=T,xlim=c(0,0.1),breaks = "scott")

mean(S2)

## [1] 0.03970672

hist(Vf,prob=T,xlim=c(0,0.4),breaks = "scott")

mean(Vf)

## [1] 0.03882001

cor(S2,Vf)

## [1] 0.006371104

hist(FF,prob=T,xlim=c(0,10),breaks = "scott")
curve(df(x, 2, 24),add=T)

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

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(27, sd = 0.2)
plot3d(x = x, y = y, z = z,type = "s", col = "blue", size = 1)

open3d()

## glX 
##   7

fit <- lm(z ~ factor(x) + factor(y))
f <- fit$ fitted.values;#予測値
a <- fit$ coefficients;#回帰係数
e <- fit$ residuals;#e#残差
se <- sum(e^2);se#残差平方和

## [1] 1.484355

s <- sqrt(se/22);s#残差標準偏差

## [1] 0.2597512

ff<- (var(f)*26/4)/(se/22);ff#F値

## [1] 13.67173

summary(fit)

## 
## Call:
## lm(formula = z ~ factor(x) + factor(y))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.46419 -0.13181  0.00621  0.13896  0.39377 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.02778    0.11178  -0.248  0.80607    
## factor(x)2   0.18067    0.12245   1.475  0.15426    
## factor(x)3   0.41335    0.12245   3.376  0.00272 ** 
## factor(y)2  -0.21590    0.12245  -1.763  0.09176 .  
## factor(y)3  -0.77965    0.12245  -6.367 2.09e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2598 on 22 degrees of freedom
## Multiple R-squared:  0.7131, Adjusted R-squared:  0.661 
## F-statistic: 13.67 on 4 and 22 DF,  p-value: 9.58e-06

plot3d(x, y, z, type = "s", col = "blue", size = 1)

coefs <- coef(fit)
a <- coefs["x"]
b <- coefs["y"]
c <- -1
d <- coefs["(Intercept)"]
#ax+by＋cz+d=0で表される平面
planes3d(a, b, c, d, alpha = 0.5)

z <- 0.2x - 0.3y+ rnorm(9, sd = 0.2)

library(dplyr)
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)
    #print(X_ij)
  }
  #print(Y_i)
  d=c(d,Y_i)
}
#d
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.2148792
1	1	-0.3084408
1	2	0.4633932
1	2	-1.0774330
1	3	-1.1873358
1	3	-0.8105574
2	1	0.5482433
2	1	0.2941053
2	2	0.1156236
2	2	0.7362176
2	3	-0.2483747
2	3	-0.8225333
3	1	-0.5366544
3	1	0.1606842
3	2	0.7743664
3	2	0.2161830
3	3	-0.1616249
3	3	-0.3364377

x <- data$x
y <- data$y
z <- data$z

open3d()

## glX 
##  10

fit <- lm(z ~ x + y)
summary(fit)

## 
## Call:
## lm(formula = z ~ x + y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.96705 -0.25150 -0.05517  0.23083  0.86930 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept) -0.09006    0.44632  -0.202   0.8428  
## x            0.27098    0.15161   1.787   0.0941 .
## y           -0.29249    0.15161  -1.929   0.0728 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5252 on 15 degrees of freedom
## Multiple R-squared:  0.3156, Adjusted R-squared:  0.2243 
## F-statistic: 3.458 on 2 and 15 DF,  p-value: 0.05819

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)

#res <- aov(value~A*B,data=data)
#summary(res)

#attach(data) 
#interaction.plot(A,B,value)

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

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

library(dplyr)
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)
Y= c()
e <- rnorm(18,0,0.2)
for (i in 1:3) {
  #X_i = s1[i]

  for (j in 1:3) {
    X_ij = 0.2*i - 0.3*j
    Y = c(Y,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.5366544
1	1	0.1606842
1	2	0.7743664
1	2	0.2161830
1	3	-0.1616249
1	3	-0.3364377
2	1	-0.5366544
2	1	0.1606842
2	2	0.7743664
2	2	0.2161830
2	3	-0.1616249
2	3	-0.3364377
3	1	-0.5366544
3	1	0.1606842
3	2	0.7743664
3	2	0.2161830
3	3	-0.1616249
3	3	-0.3364377

x <- data$x
y <- data$y
z <- data$z

open3d()

## glX 
##  12

fit <- lm(z ~ x + y)
summary(fit)

## 
## Call:
## lm(formula = z ~ x + y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.58660 -0.32533 -0.01989  0.19676  0.75495 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  8.047e-02  3.974e-01   0.202    0.842
## x           -7.259e-17  1.350e-01   0.000    1.000
## y           -3.052e-02  1.350e-01  -0.226    0.824
## 
## Residual standard error: 0.4677 on 15 degrees of freedom
## Multiple R-squared:  0.003396,   Adjusted R-squared:  -0.1295 
## F-statistic: 0.02556 on 2 and 15 DF,  p-value: 0.9748

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)