1 分散分析と回帰分析

1.1 回帰分析

回帰分析では、固定されたXに対し誤差(確率変数)を付加したYという確率変数が出力されました。
関数では Xに対し出力されるYは１つですが、回帰分析ではXに対し出力されるYは変動します。

X=0 に対するYの出力値は、誤差と同じ変動を持ちます。
シミュレーションしてみましょう。
Rの lm(linear model)関数により、回帰分析を行って見ます。

X <- c(rep(0,5),rep(1,5))
e <- rnorm(10,0,0.5)
Y <- X+e
data <- data.frame(X=X,Y=Y)
data
X Y
0 0.0132615
0 -0.0957967
0 1.0445395
0 0.4936407
0 0.3634923
1 1.9895420
1 1.4047265
1 0.9724222
1 1.5007117
1 1.7184171 
plot(data)

fit = lm( formula = Y ~ X, data = data )

a <- fit$ coefficients#回帰係数
e <- fit$ residuals#残差
se <- sum(e^2)#残差平方和
#se/2#残差分散
sqrt(se/8)#残差の標準誤差
## [1] 0.416484
a[2]#回帰係数の推定値
##        X 
## 1.153336
summary(fit)
## 
## Call:
## lm(formula = Y ~ X, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.54474 -0.29103 -0.00839  0.18339  0.68071 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   0.3638     0.1863   1.953  0.08654 . 
## X             1.1533     0.2634   4.379  0.00235 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4165 on 8 degrees of freedom
## Multiple R-squared:  0.7056, Adjusted R-squared:  0.6688 
## F-statistic: 19.17 on 1 and 8 DF,  p-value: 0.002353

1.2 分散分析

X=0 に対する5つの出力値Cと、X=1 に対する5つの出力値D に対し分散分析を行います。
データは、上のデータから取り出して使用します。

C <- data[1:5,2]
D <- data[6:10,2]
d <- data.frame(C = C, D = D)
library(reshape)
x <- melt(d, variable_name = "group")
## Using  as id variables
x
group value
C 0.0132615
C -0.0957967
C 1.0445395
C 0.4936407
C 0.3634923
D 1.9895420
D 1.4047265
D 0.9724222
D 1.5007117
D 1.7184171 
boxplot(value ~ factor(group), data = x)

res <- aov(value ~ factor(group), data = x)
summary(res)
##               Df Sum Sq Mean Sq F value  Pr(>F)   
## factor(group)  1  3.325   3.325   19.17 0.00235 **
## Residuals      8  1.388   0.173                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(C,D,paired=FALSE,alternative ="two.sided",var.equal=T)
## 
##  Two Sample t-test
## 
## data:  C and D
## t = -4.3785, df = 8, p-value = 0.002353
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.7607555 -0.5459173
## sample estimates:
## mean of x mean of y 
## 0.3638274 1.5171639

分散分析のsummary のF値は、回帰分析のF値と同じになっています。

1.3 2標本t検定と分散分析

2標本t検定の値は、t=-2.8269 であり、
\(t^2=(-2.8269)^2=7.991=F\) であり、t値の平方がF値に一致しています。
つまりt検定を一般化したものがF検定であることがわかります。

C <- rnorm(5, mean = 0, sd = 0.5)
D <- rnorm(5, mean = 1, sd = 0.5)
d <- data.frame(C = C, D = D)
library(reshape)
x <- melt(d, variable_name = "group")
## Using  as id variables
x
group value
C -0.2041560
C 0.3630882
C 0.3735975
C -0.3531959
C 0.3720951
D 0.4321007
D 0.6628048
D 1.2581517
D 0.7645144
D 0.2523534 
boxplot(value ~ factor(group), data = x)

res <- aov(value ~ factor(group), data = x)
summary(res)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## factor(group)  1 0.7944  0.7944    5.77  0.043 *
## Residuals      8 1.1013  0.1377                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(C,D,paired=FALSE,alternative ="two.sided",var.equal=T)
## 
##  Two Sample t-test
## 
## data:  C and D
## t = -2.4022, df = 8, p-value = 0.04303
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.10483142 -0.02256707
## sample estimates:
## mean of x mean of y 
## 0.1102858 0.6739850

2 確率モデル

確率モデルが分かればデータを生成することでシミュレーションが可能となる。

2.1 一元配置分散分析

　　\(\displaystyle x_{ij}=\mu_0+\delta_j+\varepsilon_{ij},~~~~~~\varepsilon_{ij} \sim N(0,\sigma^2)~i.i.d.~~~~~~~\sum_j\delta_j=0\)
　　 効果は合計が０となるよう基準化しておく必要がある。

2.1.1 基準化.

a<-　seq(1:5);scale(a)
##            [,1]
## [1,] -1.2649111
## [2,] -0.6324555
## [3,]  0.0000000
## [4,]  0.6324555
## [5,]  1.2649111
## attr(,"scaled:center")
## [1] 3
## attr(,"scaled:scale")
## [1] 1.581139

2.1.2 1要因3水準の分散分析

効果を d=(0,1,2)
基準化した効果は、scale(d) で計算される。
簡単のため、全平均は０とする。

k=c(0,1,2)
s=scale(k)
X_1 <- s[1]+rnorm(4,0,0.5)
X_2 <- s[3]+rnorm(4,0,0.5)
X_3 <- s[2]+rnorm(4,0,0.5)

d<-data.frame(X_1=X_1,X_2=X_2,X_3=X_3)
d
X_1 X_2 X_3
-0.4382725 1.449227 -0.7255023
-0.5999813 1.779837 -1.0234081
-0.6494128 0.136971 0.5362841
-1.8001031 0.825855 -0.1620181 
library(reshape)
x <- melt(d, variable_name = "group")
## Using  as id variables
boxplot(value ~ factor(group), data = x)

res <- aov(value ~ factor(group), data = x)
summary(res)
##               Df Sum Sq Mean Sq F value  Pr(>F)   
## factor(group)  2  7.869   3.935   8.504 0.00844 **
## Residuals      9  4.164   0.463                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

効果が有意であることが検出される。
一般に、誤差変動が効果に比べて小さければ効果は検出されやすく、
誤差変動が効果に比べて大きければ効果は検出されにくい。

誤差を大きくしてシミュレーションしてみよう。

k=c(0,1,2)
s=scale(k)
X_1 <- s[1]+rnorm(4,0,2)
X_2 <- s[3]+rnorm(4,0,2)
X_3 <- s[2]+rnorm(4,0,2)
d<-data.frame(X_1=X_1,X_2=X_2,X_3=X_3)
library(reshape)
x <- melt(d, variable_name = "group")
## Using  as id variables
boxplot(value ~ factor(group), data = x)

res <- aov(value ~ factor(group), data = x)
summary(res)
##               Df Sum Sq Mean Sq F value Pr(>F)
## factor(group)  2   2.11   1.057   0.215  0.811
## Residuals      9  44.31   4.923

誤差変動を sd=2 とすると、効果は検出されなくなる。

2.2 二元配置分散分析

2.2.1 繰り返しがないとき

library(dplyr)
## 
##  次のパッケージを付け加えます: 'dplyr'
##  以下のオブジェクトは 'package:reshape' からマスクされています: 
## 
##      rename
##  以下のオブジェクトは 'package:stats' からマスクされています: 
## 
##      filter, lag
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      intersect, setdiff, setequal, union
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 = X_i+s2[j]+rnorm(1,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(A=factor(c(rep("a1",3),rep("a2",3),rep("a3",3))),
                B=factor(rep(c(rep("b1",1), rep("b2",1), rep("b3",1)),3)),value=d)
data
A B value
a1 b1 -2.0768658
a1 b2 -1.4889651
a1 b3 0.2718914
a2 b1 -0.8065674
a2 b2 0.3066276
a2 b3 0.9277435
a3 b1 -0.1586784
a3 b2 0.5819345
a3 b3 1.9121813 
#boxplot(value~A+B,data=data,col="lightblue")
res <- aov(value~A+B,data=data)
summary(res)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## A            2  5.464   2.732   32.86 0.00329 **
## B            2  6.401   3.201   38.49 0.00244 **
## Residuals    4  0.333   0.083                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2.2.2 繰り返しがあるとき

library(dplyr)
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 = X_i+s2[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(A=c(rep(1,6),rep(2,6),rep(3,6)),B=c(rep(c(rep(1,2), rep(2,2), rep(3,2)),3)),value=d)
data
A B value
1 1 -1.7687349
1 1 -3.1209229
1 2 -1.0957711
1 2 -0.8533220
1 3 -1.2033367
1 3 -0.1725810
2 1 -0.8894118
2 1 -0.9672258
2 2 0.8752919
2 2 0.9919933
2 3 0.4846233
2 3 0.6598565
3 1 0.7940276
3 1 0.1614898
3 2 0.8906812
3 2 1.0615645
3 3 1.2800124
3 3 2.2141413 
#boxplot(value~A+B,data=data,col="lightblue")
res <- aov(value~A*B,data=data)
summary(res)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            1 17.804  17.804  52.332 4.33e-06 ***
## B            1  6.830   6.830  20.078 0.000518 ***
## A:B          1  0.119   0.119   0.349 0.563898    
## Residuals   14  4.763   0.340                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
attach(data) 
interaction.plot(A,B,value)

#3D表示
options(rgl.printRglwidget = TRUE) 
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
plot3d(A, B, value, type = "s", col = "blue", size = 1)

coefs <- coef(res)
a <- coefs["A"]
b <- coefs["B"]
c <- -1
d <- coefs["(Intercept)"]
planes3d(a, b, c, d, alpha = 0.5)
library(dplyr)
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 = X_i+s2[j]+rnorm(2,0,0.5)
    Y_i = c(Y_i,X_ij)
  }
  d=c(d,Y_i)
}
data=data.frame(A=c(rep(1,6),rep(2,6),rep(3,6)),B=c(rep(c(rep(1,2), rep(2,2), rep(3,2)),3)),value=d)
data
A B value
1 1 -2.6175749
1 1 -2.5422282
1 2 -2.1628725
1 2 -0.6392673
1 3 -0.2171366
1 3 1.4540408
2 1 -1.3380676
2 1 -0.9909216
2 2 -0.2699353
2 2 0.4284128
2 3 0.8523626
2 3 1.0262649
3 1 0.1321705
3 1 0.5251828
3 2 0.9837612
3 2 0.4439254
3 3 2.2023613
3 3 1.3890992 
res <- aov(value~A+B,data=data)
summary(res)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            1 12.817  12.817   35.48 2.63e-05 ***
## B            1 15.274  15.274   42.28 9.99e-06 ***
## Residuals   15  5.419   0.361                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#3D表示
options(rgl.printRglwidget = TRUE) 
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
plot3d(A, B, value, type = "s", col = "blue", size = 1)

coefs <- coef(res)
a <- coefs["A"]
b <- coefs["B"]
c <- -1
d <- coefs["(Intercept)"]
planes3d(a, b, c, d, alpha = 0.5)

2.3 交互作用モデル

2.3.1 交互作用なし

library(rgl)
df <- trees
fit <- lm(Volume ~ Girth + Height, df)
df$pred <- predict(fit)
Girth <- seq(min(df$Girth),max(df$Girth),len=100)
Height <- seq(min(df$Height),max(df$Height),len=100)
plot.df <- expand.grid(Girth=Girth,Height=Height)
plot.df$z <- predict(fit,newdata=plot.df)
library(reshape2)
## 
##  次のパッケージを付け加えます: 'reshape2'
##  以下のオブジェクトは 'package:reshape' からマスクされています: 
## 
##      colsplit, melt, recast
z <- dcast(plot.df,Height~Girth,value.var="z")[-1]
open3d(scale=c(1,1,0.2))
## glX 
##   4
points3d(df$Girth,df$Height,df$Volume,col="blue")
surface3d(Girth,Height,as.matrix(z),col="gray",alpha=.2)
apply(df,1,function(row)lines3d(rep(row[1],2),rep(row[2],2),c(log(row[3]),row[4]),col="lightblue"))
##  [1] 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
## [26] 70 71 72 73 74 75
axes3d()
title3d(xlab="Girth",ylab="Height",zlab="Volume")

2.3.2 交互作用あり

library(rgl)
df <- trees
fit <- lm(Volume ~ Girth * Height, df)
df$pred <- predict(fit)
Girth <- seq(min(df$Girth),max(df$Girth),len=100)
Height <- seq(min(df$Height),max(df$Height),len=100)
plot.df <- expand.grid(Girth=Girth,Height=Height)
plot.df$z <- predict(fit,newdata=plot.df)
library(reshape2)
z <- dcast(plot.df,Height~Girth,value.var="z")[-1]
open3d(scale=c(1,1,0.2))
## glX 
##   6
points3d(df$Girth,df$Height,df$Volume,col="blue")
surface3d(Girth,Height,as.matrix(z),col="gray",alpha=.2)
apply(df,1,function(row)lines3d(rep(row[1],2),rep(row[2],2),c(log(row[3]),row[4]),col="gray"))
##  [1]  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108 109 110 111 112
## [20] 113 114 115 116 117 118 119 120 121 122 123 124
axes3d()
title3d(xlab="Girth",ylab="Height",zlab="Volume")

2.3.3 中心化項

library(dplyr)
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 = X_i+s2[j]+s1[i]*s2[j]+rnorm(2,0,0.3)
    Y_i = c(Y_i,X_ij)
  }
  d=c(d,Y_i)
}
data=data.frame(A=c(rep(1,6),rep(2,6),rep(3,6)),B=c(rep(c(rep(1,2), rep(2,2), rep(3,2)),3)),value=d)
data
A B value
1 1 -0.4120396
1 1 -1.2894149
1 2 -1.0288215
1 2 -0.9947029
1 3 -0.7579738
1 3 -1.2705877
2 1 -1.2980142
2 1 -0.9101919
2 2 -0.4765643
2 2 -0.2101470
2 3 0.6448207
2 3 0.9178943
3 1 -0.7190937
3 1 -1.1331140
3 2 1.1352804
3 2 0.4471658
3 3 3.4693884
3 3 3.2772862 
res <- aov(value~A+B,data=data)
summary(res)
##             Df Sum Sq Mean Sq F value  Pr(>F)   
## A            1  12.46  12.465   16.11 0.00113 **
## B            1  12.09  12.086   15.62 0.00128 **
## Residuals   15  11.61   0.774                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#3D表示
options(rgl.printRglwidget = TRUE) 
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
plot3d(A, B, value, type = "s", col = "blue", size = 1)

coefs <- coef(res)
a <- coefs["A"]
b <- coefs["B"]
c <- -1
d <- coefs["(Intercept)"]
planes3d(a, b, c, d, alpha = 0.5)
library(dplyr)
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 = X_i+s2[j]+s1[i]*s2[j]+rnorm(2,0,0.3)
    Y_i = c(Y_i,X_ij)
  }
  d=c(d,Y_i)
}
data=data.frame(A=c(rep(1,6),rep(2,6),rep(3,6)),B=c(rep(c(rep(1,2), rep(2,2), rep(3,2)),3)),value=d)
data
A B value
1 1 -0.9625852
1 1 -0.4701273
1 2 -0.6005661
1 2 -1.5332573
1 3 -1.3851606
1 3 -0.9258685
2 1 -1.5361835
2 1 -0.0888921
2 2 -0.2901109
2 2 0.1974197
2 3 0.5752382
2 3 0.8018991
3 1 -1.5119947
3 1 -1.0577640
3 2 1.3396724
3 2 0.4529975
3 3 2.9539458
3 3 3.0678735 
res <- aov(value~A*B,data=data)
summary(res)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## A            1 10.309  10.309   57.67 2.50e-06 ***
## B            1  9.568   9.568   53.52 3.81e-06 ***
## A:B          1 11.210  11.210   62.71 1.54e-06 ***
## Residuals   14  2.503   0.179                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
attach(data) 
##  以下のオブジェクトは data (pos = 5) からマスクされています: 
## 
##      A, B, value
interaction.plot(A,B,value)

#3D表示
options(rgl.printRglwidget = TRUE) 
#WebGLに変換してRStudio, VSCode, ウェブブラウザなど出力
library(rgl)
plot3d(A, B, value, type = "s", col = "blue", size = 1)

coefs <- coef(res)
a <- coefs["A"]
b <- coefs["B"]
c <- -1
d <- coefs["(Intercept)"]
planes3d(a, b, c, d, alpha = 0.5)