6.1 회귀분석
예) 페인트의 불순도는 페인트를 얼마나 빨리 저어주느냐에 따라 달라진다. 아래표는 휘젓는 장치의 회전율과 불순도를 측정한 데이터이다.
x(회전율) |
20 |
22 |
24 |
26 |
28 |
30 |
32 |
34 |
36 |
38 |
40 |
42 |
y(불순도) |
8.4 |
9.5 |
11.8 |
10.4 |
13.3 |
14.8 |
13.2 |
14.7 |
16.4 |
16.5 |
18.9 |
18.5 |
> setwd('c:/Rwork')
> paint=read.csv('paint.csv')
> paint
x y
1 20 8.4
2 22 9.5
3 24 11.8
4 26 10.4
5 28 13.3
6 30 14.8
7 32 13.2
8 34 14.7
9 36 16.4
10 38 16.5
11 40 18.9
12 42 18.5
> out=lm(y~x,data=paint)
> plot(y~x,data=paint)
> abline(out)
> out
Call:
lm(formula = y ~ x, data = paint)
Coefficients:
(Intercept) x
-0.2893 0.4566
> summary(out)
Call:
lm(formula = y ~ x, data = paint)
Residuals:
Min 1Q Median 3Q Max
-1.1834 -0.5432 -0.3233 0.8333 1.3900
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.28928 1.22079 -0.237 0.817
x 0.45664 0.03844 11.880 3.21e-07 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.9193 on 10 degrees of freedom
Multiple R-squared: 0.9338, Adjusted R-squared: 0.9272
F-statistic: 141.1 on 1 and 10 DF, p-value: 3.211e-07
회귀식은 Y = -0.28928 + 0.45664X
기울기의 0.45664의 p-value는 3.21e-07로 유의.
결정계수R2는 0.9338로, 총 변동 중 93.38%가 회귀모형에 의해 설명되고 있다. 이는 두 변수 사이의 피어슨 상관계수의 제곱이다.
> with(paint,cor(x,y))^2
[1] 0.933832
F값(df1=1,df2=10)은 141.1로 p-value는 3.21e-07로 유의.
cor.test()로 구한 Pearson 상관계수의 t(df=10)=11.88, p-value 값 3.211e-07은 위의 p-value값과 일치한다. T 값인 11.88의 제곱 11.88^2 값은 141.1344로 모형의 적합성을 나타내는 F값 141.1과 일치하고 당연히 p-value도 3.21e-07로 동일하다.
> with(paint,cor.test(x,y))
Pearson's product-moment correlation
data: x and y
t = 11.88, df = 10, p-value = 3.211e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.8810937 0.9907768
sample estimates:
cor
0.9663498
out으로 저장된 결과를 plot()을 이용하여 회귀진단에 필요한 그래프를 얻을 수 있다.
> par(mfrow=c(2,2))
> plot(out)
해석:표준화된 잔차의 Normal Q-Q plot이 직선에 가깝고 다른 그래프들도 별다른 추세를 보이지 않는다.
resid()을 이용하면 회귀분석 결과물로부터 잔차를 어을 수 있다. 이렇게 받은 잔차에 아래의 명령어를 이용하여 normal Q-Q Plot을 그릴 수 있다.
> qqnorm(resid(out))
> qqline(resid(out))
정확한 p-value를 알고 싶다면 shapiro.test()를 이용하여 잔차가 정규분포를 따르는지 검정한다.
> shapiro.test(resid(out))
Shapiro-Wilk normality test
data: resid(out)
W = 0.9196, p-value = 0.2826
H0 정규분포를 따른다
H1 정규분포를 따르지 않는다
p-value는 0.2826로 잔차는 정규분포 가정을 만족시킨다. 분산분석을 통하여도 회귀의 유의성을 검정할 수 있다.
> anova(out)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 119.275 119.275 141.13 3.211e-07 ***
Residuals 10 8.451 0.845
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
6.2 공분산분석
실험이 잘 통제되어 종속변수Y의 변동을 설명하는 데 그룹변수 이외에 다른 변인이 없다면 Two-sample t-test나 일원분산분석을 하면 된다. 그러나 현실적으로 통제는 쉽지 않기 때문에 이 경우에는 여러가지 다른 변수들을 통제해주어 조사하고자 하는 변수만의 효과를 조사해야 한다. 공분산분석은 분산분석에 연속형 변수를 추가한 것이다. 궁극적인 목적인 각 그룹 간 평균들의 차이가 있는 검정하는 것으로 분산분석과 동일하나, 통제가 안 되는 연속형 변수(covariate)를 추가하여 오차를 줄이고 검정력을 높이는 것이 차이점이다.
기계1 |
기계2 |
기계3 |
|||
y 강도 |
x 섬유 두께 |
y 강도 |
x 섬유 두께 |
y 강도 |
x 섬유 두께 |
36 |
20 |
40 |
22 |
35 |
21 |
41 |
25 |
48 |
28 |
37 |
23 |
39 |
24 |
39 |
22 |
42 |
26 |
42 |
25 |
45 |
30 |
34 |
21 |
49 |
32 |
44 |
28 |
32 |
15 |
lm(y~공변량변수 + 그룹변수)
> machine<-read.csv('machine.csv')
> head(machine)
machine y x
1 m1 36 20
2 m1 41 25
3 m1 39 24
4 m1 42 25
5 m1 49 32
6 m2 40 22
> levels(machine$machine)
[1] "m1" "m2" "m3"
해석: m1기계1이 대조군이고 m2와 m3은 비교군이다.
> out=lm(y~x+machine,data=machine)
> anova(out)
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
x 1 305.130 305.130 119.9330 2.96e-07 ***
machine 2 13.284 6.642 2.6106 0.1181
Residuals 11 27.986 2.544
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
해석: p-value 값 0.1181로 기계 간의 차이가 있다고 결론을 내릴 수 없다.
summary()로 자세한 결과물을 보자.
> summary(out)
Call:
lm(formula = y ~ x + machine, data = machine)
Residuals:
Min 1Q Median 3Q Max
-2.0160 -0.9586 -0.3841 0.9518 2.8920
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.360 2.961 5.862 0.000109 ***
x 0.954 0.114 8.365 4.26e-06 ***
machinem2 1.037 1.013 1.024 0.328012
machinem3 -1.584 1.107 -1.431 0.180292
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.595 on 11 degrees of freedom
Multiple R-squared: 0.9192, Adjusted R-squared: 0.8972
F-statistic: 41.72 on 3 and 11 DF, p-value: 2.665e-06
해석: m2를 reference인 m1와 비교한 machinem2의 p-value값이 0.328012이고 m3을 m1과 비교한 machinem3의 p-value는 0.180292이다. 그러나 이 p-value들을 조정없이 사용하여 두 machine 모두 m1와 유의한 차이가 없다고 결론을 내리지는 않는다.
1개의 대조군(m1)을 2개의 비교군과 비교하므로 Dunnett의 방법으로 p-value를 조정한다.
> library(multcomp)
> dunnett=glht(out,linfct=mcp(machine='Dunnett'))
> summary(dunnett)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Dunnett Contrasts
Fit: lm(formula = y ~ x + machine, data = machine)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
m2 - m1 == 0 1.037 1.013 1.024 0.518
m3 - m1 == 0 -1.584 1.107 -1.431 0.304
(Adjusted p values reported -- single-step method)
> plot(dunnett)
해석: 두 p-value 모두가 원래 p-value보다 커졌으며, 여전히 유의한 차이가 있다고 결론을 지을 수 없다. Plot 그래프로 확인할 수도 있다. 신뢰구간이 0을 포함하고 있으므로 유의하다고 볼 수 없다.
|
summary(out) |
dunnett |
m2 |
0.328012 |
0.518 |
m3 |
0.180292 |
0.304 |
> with(data=machine,tapply(y,machine,mean))
m1 m2 m3
41.4 43.2 36.0
> with(data=machine,tapply(y,machine,sd))
m1 m2 m3
4.827007 3.701351 3.807887
해석: 각 기계에서 생산된 섬유 제품의 강도 평균은 41.4, 43.2, 36.0으로 얻어져 혹시 차이가 있다면 기계3의 강도가 떨어지는 듯 하다.
출처: 실험 계획과 응용, R로 하는 통계 분석
'KNOU > 2 실험 계획과 응용' 카테고리의 다른 글
4장 이원배치법, 이원분산분석, two way anova (0) | 2016.10.17 |
---|---|
3장 일원배치법, 일원분산분석, one way anova (0) | 2016.10.12 |
2장 짝지어진 비교 Paired T-Test (1) | 2016.10.12 |
2장 독립표본 두 모평균 차이 추론 two sample test (0) | 2016.10.12 |