重回帰分析と分散分析の結果が違う

次のようなRのコードを書きました。
x1~x7は因子で、1がなし2がありだと思ってください。
最後のxがデータです。これが大きい程よいと言う意味だと思ってください。

data <- data.frame(
x1 = c("1","1","1","1","2","2","2","2"),
x2 = c("1","1","2","2","1","1","2","2"),
x3 = c("1","1","2","2","2","2","1","1"),
x4 = c("1","2","1","2","1","2","1","2"),
x5 = c("1","2","1","2","2","1","2","1"),
x6 = c("1","2","2","1","1","2","2","1"),
x7 = c("1","2","2","1","2","1","1","2"),
x = c(35,48,21,38,50,43,31,22)
)

data
colnames(data) <-c("B","C","B:C","D","B:D","A","7","x")
data
summary(aov(x ~　A+B+C+D+B:C+B:D,data=data))
summary(lm(x ~　A+B+C+D+B:C+B:D,data=data))


結果としましては、
分散分析では因子C以外は有意差がないという結論を出しており、
重回帰分析では因子Bと因子B:D以外は有意差がないという結論を出しています。


重回帰分析と分散分析は基本的に数学的にやっていることは同じと言われたのですが、結果が異なり戸惑っています。
なぜこのようなことになったのか心当たりがある方ご教授ください。
よろしくお願いします。

No.4 ベストアンサー 回答者： kamiyasiro 回答日時： 2014/10/18 06:32

企業でSQC推進を担当する統計家です。

まず、ご質問の趣旨から外れますが、
定石と言われるL8の割りつけは、主効果を列[1][2][4][7]の順で割り付けます。
今、A主効果が第6列に割りつけられていますので、CD交互作用と交絡しています。
マズイ実験ですよ。
CD交互作用によって、Aの効果が打ち消されている可能性があります。

では、本題です。
列名称が、「B:Cじゃダメ」という回答が、もう出ています。
B:Cは積項を作る演算記号だからです。
関数形指定の際に使用します。

さて、ご質問者はRのスクリプトもさることながら、
L8直交実験だから、列を使っても、積項を生成しても一緒でしょ。
と言いたくはありませんか。
たぶん、何かの本かサイトに、一致すると書いてあったのではないでしょうか。
実は、そうなんです。やってみえることは、何ら間違いではありません。
見方が違うだけです。
lmの結果はanovaで見なければいけません。

では、第1に、列名称をBCと変更して、
BCとB:Cの違いを、aovと、lmの結果をanovaで見たときを比較してみましょう。

第2に、L8直交表を変形して、全ての結果が同じになることを見ていきましょう。
直交表を変形する前は、第1列と第2列の積は、次のようになります。

[1][2]　　[3]とは一致しない
1　1　→　1
1　1　→　1
1　2　→　2
1　2　→　2
2　1　→　2
2　1　→　2
2　2　→　4
2　2　→　4

では、水準を-1，+1のように、コード化してみましょう。
すると、第1列と第2列の積は、符号が違いますが第3列と一致します。

[1][2]　　　[3]と一致
-1　-1　→　+1
-1　-1　→　+1
-1　+1　→　-1
-1　+1　→　-1
+1　-1　→　-1
+1　-1　→　-1
+1　+1　→　+1
+1　+1　→　+1

こうすれは、ご質問者のやり方で何ら間違いはありません。

スクリプトは、以下に示します。
コピペして実行してみて下さい。

#
#計画行列はエクセルから貼り付けることも考慮して、こう書こう
#
dat <- data.frame(
matrix(c(
1,1,1,1,1,1,1,35,
1,1,1,2,2,2,2,48,
1,2,2,1,1,2,2,21,
1,2,2,2,2,1,1,38,
2,1,2,1,2,1,2,50,
2,1,2,2,1,2,1,43,
2,2,1,1,2,2,1,31,
2,2,1,2,1,1,2,22)
,ncol=8,nrow=8,byrow=T)
)
names(dat)=c("B","C","BC","D","BD","A","e","x")
dat
#
#正しい変数指定と間違った変数指定（B:Dとすると内部で生成）
#
summary(aov(x ~ A+B+C+D+BC +BD, data=dat)) #正
summary(aov(x ~ A+B+C+D+B:C+B:D,data=dat)) #誤だけど一致
summary(lm (x ~ A+B+C+D+BC +BD, data=dat)) #正
summary(lm (x ~ A+B+C+D+B:C+B:D,data=dat)) #誤一致しない
#
#内部生成してもANOVAを取れば結果は一致
#
result1 <- lm (x ~ A+B+C+D+B:C+B:D,data=dat)
anova(result1)
#
#回帰なら変数選択しておこう（正しい変数指定で）
#
result2 <-lm (x ~ A+B+C+D+BC+BD,data=dat)
result2.step <- step(result2)
summary(result2.step)
#
#変数選択したモデルでの分散分析結果は一致する
#
summary(aov(x ~ B+C+D+BC+BD,data=dat))
#
#コード化，中心化を行って、計画行列を正しい表現にする
#
dat.s <- dat[,-8]
dat.s <- (dat.s-1.5)*2
dat[,1:7] <- dat.s
dat
#
#そうすれば結果は全て一致する
#
summary(aov(x ~ A+B+C+D+BC +BD, data=dat)) #正
summary(aov(x ~ A+B+C+D+B:C+B:D,data=dat)) #誤だけど一致
summary(lm (x ~ A+B+C+D+BC +BD, data=dat)) #正
summary(lm (x ~ A+B+C+D+B:C+B:D,data=dat)) #正
#
#変数選択しておこう
#
result3 <-(lm (x ~ A+B+C+D+BC+BD,data=dat))
result3.step <- step(result3)
summary(result3.step)
#
#そうすれば、予測も可能
#
x <- data.frame(1,-1,1,1,1,1,1) #予測したいセット，ここでは－1～1の値を取る
index <- names(dat[1:7]) %in% names(result3.step$coefficients)
x <- x[index==TRUE]
p <- length(result3.step$coefficients)
y.hat <- (result3.step$coefficients[2:p] %*% t(x)) + mean(dat$x)
y.hat
#
#もちろん、元データを代入して、重相関を見ることも可能
#
x <- dat[,1:7]
index <- names(dat[1:7]) %in% names(result3.step$coefficients)
x <- x[index==TRUE]
p <- length(result3.step$coefficients)
y.hat <- (result3.step$coefficients[2:p] %*% t(x)) + mean(dat$x)
plot(y.hat,dat$x,pch=16)
r <- cor(t(y.hat),dat$x)
lm.obj <- lm(dat$x~t(y.hat))
abline(lm.obj,col=2)
mtext(substitute(paste(R,"=",r1),list(r1=r)))
#
#
#

この回答へのお礼

非常に丁寧で分かりやすい回答をしていただきありがとうございました！
コードの方も非常に参考になりました。
#1,#2#3様もありがとうございました。

お礼日時：2014/10/18 15:52

No.3 回答者： quaestio 回答日時： 2014/10/18 00:00

いろいろ修正していたらおかしなまま回答してしまったので、#2は無視してください。

#1さんも指摘していますが、
> summary(aov(x ~　A+B+C+D+B:C+B:D,data=data))
のB:CとB:Dはdataの3列目と5列目を意味するものではありませんので注意してください。

まずは、Rで次のコードを実行してみてください。

data <- data.frame(
x1 = c("1","1","1","1","2","2","2","2"),
x2 = c("1","1","2","2","1","1","2","2"),
x3 = c("1","1","2","2","2","2","1","1"),
x4 = c("1","2","1","2","1","2","1","2"),
x5 = c("1","2","1","2","2","1","2","1"),
x6 = c("1","2","2","1","1","2","2","1"),
x7 = c("1","2","2","1","2","1","1","2"),
x = c(35,48,21,38,50,43,31,22)
)
colnames(data) <-c("B","C","BC","D","BD","A","7","x")
data
summary(aov(x ~ A+B+C+D+B:C+B:D,data=data)) # summary(aov(x ~　A+B+C+D+B:C+B:D,data=data))と同じ
summary(aov(x ~ A+B+C+D+BC+BD,data=data))
summary(lm(x ~ A+B+C+D+B:C+B:D,data=data)) # summary(lm(x ~　A+B+C+D+B:C+B:D,data=data))と同じ
summary(lm(x ~ A+B+C+D+BC+BD,data=data))

最後から二行目の結果だけ違いますよね。
ですが、次のコードで分散分析表とfitted valueを確認してください。

data.lm <- lm(x ~ A+B+C+D+B:C+B:D,data=data)
data.aov <- aov(x ~ A+B+C+D+B:C+B:D,data=data)

anova(data.lm)
anova(data.aov)
fitted(data.lm)
fitted(data.aov)

おや、全く同じですね。
ということで、質問文中にもあるように、重回帰分析と分散分析は基本的に数学的にやっていることは同じです。

では、何が違うのかというと、検定しているものが違うからです。

データの構造を
x[ijkl] = m + a[i] + b[j] + c[k] + d[l] + (bc)[jk] + (bd)[jl] + ε[ijkl]
としましょう。
ただし、mは総平均、a,b,c,dは要因A,B,C,Dによる変動、(bc),(bd)は交互作用B×C,B×Dによる変動、εは偶然誤差です。

すると、lmの係数はそれぞれ次の推定値となります。
ただし、E{・}は期待値を意味します。

A2 a[2]-a[1]　(= E{x[2111]-x[1111]})
B2 b[2]-b[1] + 2(bc)[12] + 2(bd)[12]　(= E{x[1211]-x[1111]})
C2 c[2]-c[1] + 2(bc)[12]　(= E{x[1121]-x[1111]})
D2 d[2]-d[1] + 2(bd)[12]　(= E{x[1112]-x[1111]})
B2:C2 -4(bc)[12]　(= E{x[1221]-x[1111]}-E{x[1211]-x[1111]}-E{x[1121]-x[1111]})
B2:D2 -4(bd)[12]　(= E{x[1212]-x[1111]}-E{x[1211]-x[1111]}-E{x[11121-x[1111]})
(Intercept)の値はi=j=k=l=1の推定値です。

そして、検定はこれらが0でないかを検定しています。

一方、aovの方は

A2 a[2]-a[1]
B2 b[2]-b[1]
C2 c[2]-c[1]
D2 d[2]-d[1]
B2:C2 (bc)[12] = (bc)[21] = -(bc)[11] = -(bc)[22]
B2:D2 (bd)[12] = (bd)[21] = -(bd)[11] = -(bd)[22]

が0でないかをを検定しているので、検定の結果も当然異なります。

No.2 回答者： quaestio 回答日時： 2014/10/17 23:36

#1さんも指摘していますが、

> summary(aov(x ~　A+B+C+D+B:C+B:D,data=data))
のB:CとB:Dはdataの3列目と5列目を意味するものではありませんので注意してください。

まずは、Rで次のコードを実行してみてください。

data <- data.frame(
x1 = c("1","1","1","1","2","2","2","2"),
x2 = c("1","1","2","2","1","1","2","2"),
x3 = c("1","1","2","2","2","2","1","1"),
x4 = c("1","2","1","2","1","2","1","2"),
x5 = c("1","2","1","2","2","1","2","1"),
x6 = c("1","2","2","1","1","2","2","1"),
x7 = c("1","2","2","1","2","1","1","2"),
x = c(35,48,21,38,50,43,31,22)
)
colnames(data) <-c("B","C","BC","D","BD","A","7","x")
data
summary(aov(x ~ A+B+C+D+B:C+B:D,data=data)) # summary(aov(x ~　A+B+C+D+B:C+B:D,data=data))と同じ
summary(aov(x ~ A+B+C+D+BC+BD,data=data))
summary(lm(x ~ A+B+C+D+B:C+B:D,data=data)) # summary(lm(x ~　A+B+C+D+B:C+B:D,data=data))と同じ
summary(lm(x ~ A+B+C+D+BC+BD,data=data))

最後から二行目の結果だけ違いますよね。
ですが、次のコードで分散分析表とfitted valueを確認してください。

data.lm <- lm(x ~ A+B+C+D+B:C+B:D,data=data)
data.aov <- aov(x ~ A+B+C+D+B:C+B:D,data=data)

anova(data.lm)
anova(data.aov)
fitted(data.lm)
fitted(data.aov)

おや、全く同じですね。
ということで、質問文中にもあるように、重回帰分析と分散分析は基本的に数学的にやっていることは同じです。

では、何が違うのかというと、
それぞれの係数の意味が違うからです。

データの構造を
x[ijkl] = m + a[i] + b[j] + c[k] + d[l] + (bc)[jk] + (bd)[jl] + ε[ijkl]
としましょう。
ただし、mは総平均、a,b,c,dは要因A,B,C,Dによる変動、(bc),(bd)は交互作用B×C,B×Dによる変動、εは偶然誤差です。

すると、aovの方の(Intercept)以外の係数はそれぞれ次の推定値になります。

A2 a[2]-a[1]
B2 b[2]-b[1]
C2 c[2]-c[1]
D2 d[2]-d[1]
BC2 2(bc)[12] = 2(bc)[21] = -2(bc)[11] = -2(bc)[22]
BD2 2(bd)[12] = 2(bd)[21] = -2(bd)[11] = -2(bd)[22]
(Intercept)の値はi=j=k=l=1の推定値です。

一方、lmはどうかというと、

A2 a[2]-a[1]　(= E{x[2111]-x[1111]})
B2 b[2]-b[1] + 2(bc)[12] + 2(bd)[12]　(= E{x[1211]-x[1111]})
C2 c[2]-c[1] + 2(bc)[12]　(= E{x[1121]-x[1111]})
D2 d[2]-d[1] + 2(bd)[12]　(= E{x[1112]-x[1111]})
B2:C2 -4(bc)[12]　(= E{x[1221]-x[1111]}-E{x[1211]-x[1111]}-E{x[1121]-x[1111]})
B2:D2 -4(bd)[12]　(= E{x[1212]-x[1111]}-E{x[1211]-x[1111]}-E{x[11121-x[1111]})

となるのです。
ただし、E{・}は期待値を意味します。
従って、違うものの推定値なので検定の結果も当然異なります。

No.1 回答者： toukeier 回答日時： 2014/10/17 15:57

「:」（コロン）はformula中で意味がある記号ですので、変数名に使うのは適切でない気がします。

交互作用項を投入するときに使います。それを織り込み済みなのかどうかちょっとわかりかねました。変数名をx1からx6のままで実施するのはいかがでしょうか。ちなみに「*」（アスタリスク）は、主効果と交互作用項をともに投入するときに使います。

> summary(lm(x ~ B+C, data=data))

Call:
lm(formula = x ~ B + C, data = data)

Residuals:
1 2 3 4 5 6 7 8
-8.5 4.5 -6.5 10.5 5.5 -1.5 2.5 -6.5

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 43.500 4.945 8.797 0.000315 ***
B2 1.000 5.710 0.175 0.867839
C2 -16.000 5.710 -2.802 0.037892 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.075 on 5 degrees of freedom
Multiple R-squared: 0.6119, Adjusted R-squared: 0.4567
F-statistic: 3.942 on 2 and 5 DF, p-value: 0.09383

> summary(lm(x ~ B:C, data=data))

Call:
lm(formula = x ~ B:C, data = data)

Residuals:
1 2 3 4 5 6 7 8
-6.5 6.5 -8.5 8.5 3.5 -3.5 4.5 -4.5

Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.500 6.062 4.371 0.012 *
B1:C1 15.000 8.573 1.750 0.155
B2:C1 20.000 8.573 2.333 0.080 .
B1:C2 3.000 8.573 0.350 0.744
B2:C2 NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.573 on 4 degrees of freedom
Multiple R-squared: 0.65, Adjusted R-squared: 0.3875
F-statistic: 2.476 on 3 and 4 DF, p-value: 0.2008

> summary(lm(x ~ B*C, data=data))

Call:
lm(formula = x ~ B * C, data = data)

Residuals:
1 2 3 4 5 6 7 8
-6.5 6.5 -8.5 8.5 3.5 -3.5 4.5 -4.5

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 41.500 6.062 6.846 0.00238 **
B2 5.000 8.573 0.583 0.59105
C2 -12.000 8.573 -1.400 0.23418
B2:C2 -8.000 12.124 -0.660 0.54542
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.573 on 4 degrees of freedom
Multiple R-squared: 0.65, Adjusted R-squared: 0.3875
F-statistic: 2.476 on 3 and 4 DF, p-value: 0.2008