重回帰分析から分散分析へ
回帰分析・分散分析、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.1897915s <- sqrt(se/(27-3));s#残差標準偏差## [1] 0.08892682ff<- (var(f)*26/2)/(se/24);ff#F値## [1] 131.9669summary(fit)##
## Call:
## lm(formula = z ~ x + y)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.234917 -0.050644 0.001932 0.047469 0.156779
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.10583 0.06171 -1.715 0.0992 .
## x 0.21904 0.02096 10.450 2.07e-10 ***
## y -0.26072 0.02096 -12.439 5.92e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08893 on 24 degrees of freedom
## Multiple R-squared: 0.9166, Adjusted R-squared: 0.9097
## F-statistic: 132 on 2 and 24 DF, p-value: 1.125e-13plot3d(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.2331851 | 0.0215345 |
| x | 0.1757752 | 0.2622991 |
| y | -0.3039865 | -0.2174626 |
95%信頼区間は、真の回帰係数 0.2 , -0.3 を含んでいます。
回帰がないという帰無仮説に対する t値、F値 は有意であり、決定係数も１に近い値となっています。
1.1.1 １つの変数を固定
coefs## (Intercept) x y
## -0.1058253 0.2190372 -0.2607246a;b;d## x
## 0.2190372## y
## -0.2607246## (Intercept)
## -0.1058253curve(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.1.2 ２変量分散分析
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.8636 0.8636 109.2 2.07e-10 ***
## y 1 1.2236 1.2236 154.7 5.92e-12 ***
## Residuals 24 0.1898 0.0079
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(x,y,z, type="b", pch=c(1,2))
interaction.plotは交差しません。
1.2 交互作用モデル
モデルに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 90.71 90.71 1354.7 < 2e-16 ***
## y 1 48.70 48.70 727.4 < 2e-16 ***
## x:y 1 13.19 13.19 197.0 9.12e-13 ***
## Residuals 23 1.54 0.07
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1fit <- lm(z ~ x * y)
f <- fit$ fitted.values;#予測値
a <- fit$ coefficients;a#回帰係数## (Intercept) x y x:y
## 0.2591776 0.1478323 -0.4521285 1.0485334e <- fit$ residuals;#e#残差
se <- sum(e^2);se#残差平方和## [1] 1.540116s <- sqrt(se/23);s#残差標準偏差## [1] 0.2587693ff<- (var(f)*26/3)/(se/23);ff#F値## [1] 759.6903summary(fit)##
## Call:
## lm(formula = z ~ x * y)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59242 -0.17367 -0.01881 0.15534 0.46226
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2592 0.3486 0.743 0.4647
## x 0.1478 0.1614 0.916 0.3691
## y -0.4521 0.1614 -2.802 0.0101 *
## x:y 1.0485 0.0747 14.037 9.12e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2588 on 23 degrees of freedom
## Multiple R-squared: 0.99, Adjusted R-squared: 0.9887
## F-statistic: 759.7 on 3 and 23 DF, p-value: < 2.2e-16plot3d(x, y, z, type = "s", col = "blue", size = 1)coefs <- coef(fit)
a <- coefs["x"];a## x
## 0.1478323b <- coefs["y"];b## y
## -0.4521285c <- coefs["x:y"];c## x:y
## 1.048533d <- coefs["(Intercept)"];d## (Intercept)
## 0.2591776df <- 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 90.71 90.71 1354.7 < 2e-16 ***
## y 1 48.70 48.70 727.4 < 2e-16 ***
## x:y 1 13.19 13.19 197.0 9.12e-13 ***
## Residuals 23 1.54 0.07
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction.plot(x,y,z, type="b", pch=c(1,2))
1.2.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.03910976hist(Vf,prob=T,xlim=c(0,0.4),breaks = "scott")
mean(Vf)## [1] 0.03968163cor(S2,Vf)## [1] 0.01335994hist(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
## 7fit <- lm(z ~ factor(x) + factor(y))
f <- fit$ fitted.values;#予測値
a <- fit$ coefficients;#回帰係数
e <- fit$ residuals;#e#残差
se <- sum(e^2);se#残差平方和## [1] 1.114323s <- sqrt(se/22);s#残差標準偏差## [1] 0.2250579ff<- (var(f)*26/4)/(se/22);ff#F値## [1] 14.54981summary(fit)##
## Call:
## lm(formula = z ~ factor(x) + factor(y))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.59336 -0.09831 0.00124 0.17634 0.31025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.09682 0.09685 -1.000 0.328304
## factor(x)2 0.25514 0.10609 2.405 0.025034 *
## factor(x)3 0.37894 0.10609 3.572 0.001704 **
## factor(y)2 -0.49515 0.10609 -4.667 0.000118 ***
## factor(y)3 -0.68964 0.10609 -6.500 1.54e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2251 on 22 degrees of freedom
## Multiple R-squared: 0.7257, Adjusted R-squared: 0.6758
## F-statistic: 14.55 on 4 and 22 DF, p-value: 5.945e-06plot3d(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)##
## 次のパッケージを付け加えます: 'dplyr'## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## filter, lag## 以下のオブジェクトは 'package:base' からマスクされています:
##
## intersect, setdiff, setequal, unionlibrary(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.0788109 |
| 1 | 1 | 0.0711652 |
| 1 | 2 | -1.0681484 |
| 1 | 2 | -0.6521538 |
| 1 | 3 | -1.3425252 |
| 1 | 3 | -0.1617509 |
| 2 | 1 | -0.3077812 |
| 2 | 1 | 0.2619091 |
| 2 | 2 | 0.0914813 |
| 2 | 2 | 0.3770143 |
| 2 | 3 | 0.2870105 |
| 2 | 3 | -1.0625904 |
| 3 | 1 | 0.1065304 |
| 3 | 1 | 0.1350579 |
| 3 | 2 | 0.0828609 |
| 3 | 2 | -0.1190660 |
| 3 | 3 | 0.2163023 |
| 3 | 3 | 0.8407757 |
x <- data$x
y <- data$y
z <- data$z
open3d()## glX
## 10fit <- lm(z ~ x + y)
summary(fit)##
## Call:
## lm(formula = z ~ x + y)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8160 -0.2864 -0.0301 0.4197 0.7128
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6430 0.4156 -1.547 0.1426
## x 0.3746 0.1412 2.653 0.0181 *
## y -0.1176 0.1412 -0.833 0.4180
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.489 on 15 degrees of freedom
## Multiple R-squared: 0.3402, Adjusted R-squared: 0.2522
## F-statistic: 3.867 on 2 and 15 DF, p-value: 0.04421plot3d(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.1065304 |
| 1 | 1 | 0.1350579 |
| 1 | 2 | 0.0828609 |
| 1 | 2 | -0.1190660 |
| 1 | 3 | 0.2163023 |
| 1 | 3 | 0.8407757 |
| 2 | 1 | 0.1065304 |
| 2 | 1 | 0.1350579 |
| 2 | 2 | 0.0828609 |
| 2 | 2 | -0.1190660 |
| 2 | 3 | 0.2163023 |
| 2 | 3 | 0.8407757 |
| 3 | 1 | 0.1065304 |
| 3 | 1 | 0.1350579 |
| 3 | 2 | 0.0828609 |
| 3 | 2 | -0.1190660 |
| 3 | 3 | 0.2163023 |
| 3 | 3 | 0.8407757 |
x <- data$x
y <- data$y
z <- data$z
open3d()## glX
## 12fit <- lm(z ~ x + y)
summary(fit)##
## Call:
## lm(formula = z ~ x + y)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32948 -0.19798 -0.01378 0.12852 0.42649
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.973e-01 2.319e-01 -0.851 0.4083
## x 2.676e-17 7.879e-02 0.000 1.0000
## y 2.039e-01 7.879e-02 2.588 0.0206 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2729 on 15 degrees of freedom
## Multiple R-squared: 0.3086, Adjusted R-squared: 0.2164
## F-statistic: 3.348 on 2 and 15 DF, p-value: 0.06278plot3d(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)