6장 회귀분석과 공분산분석

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로 유의.

결정계수R²는 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로 하는 통계 분석