ウォッチ
次のような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ベストアンサー
- 回答日時：
企業で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)))
#
#
#
-
-
- 0
- 件
-
この回答へのお礼
お礼日時：2014/10/18 15:52
非常に丁寧で分かりやすい回答をしていただきありがとうございました！
コードの方も非常に参考になりました。
#1,#2#3様もありがとうございました。
No.3
- 回答日時：
いろいろ修正していたらおかしなまま回答してしまったので、#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でないかをを検定しているので、検定の結果も当然異なります。
-
-
- 0
- 件
-
No.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による変動、εは偶然誤差です。
すると、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{・}は期待値を意味します。
従って、違うものの推定値なので検定の結果も当然異なります。
-
-
- 0
- 件
-
No.1
- 回答日時：
「:」（コロン）は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
-
-
- 0
- 件
-
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう！
似たような質問が見つかりました
- 統計学 加重最小二乗法＝①「変数を自然対数変換」＝②「誤差項の分散の逆数を重み付け」？ 8 2022/11/26 11:15
- 統計学 回帰分析検証不一致 5 2022/05/04 20:48
- Excel（エクセル） SUMIFSと日付変換 10 2023/04/16 15:38
- 統計学 Rを用いた「繰り返しがある直交表実験計画法」の分析方法 8 2023/08/01 17:58
- その他（プログラミング・Web制作） データ解析ソフトRでのファイル入力read.csvがエラーになります 7 2022/03/27 22:11
- 統計学 ダミー変数の因子分析について質問です 1 2022/03/27 14:31
- 統計学 卒論でエクセルの単回帰分析を使用して分析を行いました。結果は出たのですが、どの数値をどのように論文に 2 2023/01/10 03:17
- 心理学 卒業論文の分析について 1 2022/09/05 14:19
- 大人・中高年 苦しい時って落ち込んだり暗くならないといけないんですか？ 1 2022/04/22 19:13
- その他（Microsoft Office） Excel 2列分のDATAの並べ替え 1 2023/01/01 17:12
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
因数分解の問題です a^2(b+c)+b...
-
a(b二乗−c二乗)＋b(c二乗−a二乗...
-
行列の二次方程式
-
数学の問題
-
高１の因数分解の問題ですが、...
-
判別式D/4の公式は b²ｰac です...
-
よろしくおねがいします
-
この展開、因数分解ができませ...
-
文字を使った行列式の計算
-
(a+b)(b+c)(c+a)+abcの因数分解
-
因数分解を教えてください
-
ab(b+c)+bc(b+c)+ca...3abc
-
数学の問題
-
数学です。 P=a^2+b^2+c^2-ab-b...
-
数学のあまり難しくない問題です
-
a^3(b-c)+b^3(c-a)+c^3(a-b)
-
複素関数論、１次変換
-
難解な等式の証明：α^6+β^6+γ^6=?
-
因数分解の問題
-
mm3とμｌ
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
因数分解を教えてください
-
(a+b)(b+c)(c+a)+abcの因数分解
-
数学の問題
-
a(b二乗−c二乗)＋b(c二乗−a二乗...
-
因数分解の問題です a^2(b+c)+b...
-
数学のあまり難しくない問題です
-
（a+b+c）³を展開して a³+b³+c³...
-
一つの文字について整理とは？...
-
二次関数の難問です。 P=x^2+y^...
-
高１の因数分解の問題ですが、...
-
a^2 + b^2 + c^2 - 2ab - 2bc -...
-
a三乗(bーc)＋b三乗(cーa)＋c三...
-
因数分解の問題
-
数学です。 P=a^2+b^2+c^2-ab-b...
-
aとbが互いに素 aとcが互いに素...
-
この因数分解がわかりません＞＜
-
a>0 b<0 c<0 のとき√(a^2 b c^3...
-
（a-b)^2+(b-c)^2+(c-a)^2を簡...
-
山梨大学の数Iの入試問題が回答...
-
1つの頂点が原点にある三角形の...
おすすめ情報
