分散分析を行う
下記のような表があったとき（はてな記法でcolspanのやつとかできないのね！めんどくさい！）
要因１，要因２の影響について調べる。
| 要因１ | ||||
|---|---|---|---|---|
| 要因２ | A | B | C | D |
| a | 34 | 41 | 52 | 56 |
| b | 35 | 46 | 55 | 62 |
| c | 40 | 53 | 56 | 63 |
このデータ、その形のままだとRに取り込みにくいのでこんな形の単純な表（データフレーム）にする。
factor1 factor2 data 1 a A 34 2 a B 41 3 a C 52 4 a D 56 5 b A 35 6 b B 46 7 b C 55 8 b D 62 9 c A 40 10 c B 53 11 c C 56 12 c D 63
そのためには下記のように。関数repは同じ回数分繰り返したデータを作ってくれるので、aを連続したければ rep("a", n)のようにする。
> factor1 = c(rep("a", 4), rep("b", 4), rep("c", 4))
> factor2 = rep(c("A", "B", "C", "D"),3)
> data = c(34, 41,52,56,35,46,55,62,40,53,56,63)
> data
> f = data.frame(f1=factor1, f2=factor2, data)これに対してaovで分散分析を行う
> summary(aov(data~f1+f2, data=f))
Df Sum Sq Mean Sq F value Pr(>F)
f1 2 105.2 52.6 12.37 0.00743 **
f2 3 966.3 322.1 75.78 3.68e-05 ***
Residuals 6 25.5 4.2
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1すると、P値が5%より低く、f1もf2もこの要因として適している、と分析できる。多分この***の多いf2の要因のほうが強いのだろう。
試しに、明らかに要因１の影響しか受けないこんなデータでやると、
> data = c(1, 1,1,1,2,2,2,2,3,3,3,3)
> f = data.frame(f1=factor1, f2=factor2, data)
> f
f1 f2 data
1 a A 1
2 a B 1
3 a C 1
4 a D 1
5 b A 2
6 b B 2
7 b C 2
8 b D 2
9 c A 3
10 c B 3
11 c C 3
12 c D 3
> summary(aov(data~f1+f2, data=f))
Df Sum Sq Mean Sq F value Pr(>F)
f1 2 8 4 3.865e+31 <2e-16 ***
f2 3 0 0 1.000e+00 0.455
Residuals 6 0 0
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘f1のP値のみが小さく、f2は45%とかなり大きいのでf2の影響はないと言える
このデータを回帰分析してみよう
要するにa,b, A,B,C のみでこのデータは決まるわけなので、数量化理論一類にて、下記のような表と見なせる
| a | b | A | B | C | data |
| 1 | 0 | 1 | 0 | 0 | 34 |
| 1 | 0 | 0 | 1 | 0 | 41 |
| 1 | 0 | 0 | 0 | 1 | 52 |
| 1 | 0 | 0 | 0 | 0 | 56 |
| 0 | 1 | 1 | 0 | 0 | 35 |
| 0 | 1 | 0 | 1 | 0 | 46 |
| 0 | 1 | 0 | 0 | 1 | 55 |
| 0 | 1 | 0 | 0 | 0 | 62 |
| 0 | 0 | 1 | 0 | 0 | 40 |
| 0 | 0 | 0 | 1 | 0 | 53 |
| 0 | 0 | 0 | 0 | 1 | 56 |
| 0 | 0 | 0 | 0 | 0 | 63 |
これをデータフレームにして
> a = c(1,1,1,1,0,0,0,0,0,0,0,0) > b = c(0,0,0,0,1,1,1,1,0,0,0,0) > A = c(1,0,0,0,1,0,0,0,1,0,0,0) > B = rep(c(0,1,0,0), 3) # 早速使ってみた > C = rep(c(0,0,1,0), 3) > data.frame(a=a,b=b,A=A,B=B,C=C, data=data) a b A B C data 1 1 0 1 0 0 34 2 1 0 0 1 0 41 3 1 0 0 0 1 52 4 1 0 0 0 0 56 5 0 1 1 0 0 35 6 0 1 0 1 0 46 7 0 1 0 0 1 55 8 0 1 0 0 0 62 9 0 0 1 0 0 40 10 0 0 0 1 0 53 11 0 0 0 0 1 56 12 0 0 0 0 0 63
線形回帰分析する
> summary(lm(data~a+b+A+B+C))
Call:
lm(formula = data ~ a + b + A + B + C)
Residuals:
Min 1Q Median 3Q Max
-2.0000 -1.0417 -0.2917 1.3333 2.7500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 63.917 1.458 43.846 9.42e-09 ***
a -7.250 1.458 -4.973 0.002518 **
b -3.500 1.458 -2.401 0.053224 .
A -24.000 1.683 -14.258 7.44e-06 ***
B -13.667 1.683 -8.119 0.000187 ***
C -6.000 1.683 -3.565 0.011862 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.062 on 6 degrees of freedom
Multiple R-squared: 0.9768, Adjusted R-squared: 0.9574
F-statistic: 50.42 on 5 and 6 DF, p-value: 8.03e-05 b,Cを抜いて
> summary(lm(data~a+A+B))
Call:
lm(formula = data ~ a + A + B)
Residuals:
Min 1Q Median 3Q Max
-4.1667 -2.6667 -0.1667 2.4583 4.5000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 59.167 1.646 35.952 3.92e-10 ***
a -5.500 2.208 -2.491 0.03746 *
A -21.000 2.550 -8.237 3.54e-05 ***
B -10.667 2.550 -4.184 0.00306 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.606 on 8 degrees of freedom
Multiple R-squared: 0.9052, Adjusted R-squared: 0.8696
F-statistic: 25.46 on 3 and 8 DF, p-value: 0.0001912 data = -5.5 * a - 21 * A - 10 * B + 59.167 の回帰式を得ることができた。
たとえば aとCの要因があるデータは、 data = -5.5 + 59.167 = 53.667 と求めることができるし、（実データは52）、cとDの場合は59.167（実データは63）となる。
