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

KNOU/2 실험 계획과 응용 2016. 10. 17. 19:32
반응형

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
Posted by 마르띤
,

4장 이원배치법, 이원분산분석, two way anova

KNOU/2 실험 계획과 응용 2016. 10. 17. 19:19
반응형

이원산분석은 그룹변수가 하나인 일원산분석의 확장으로 두 개의 그룹변수를 가진다. 이 분석의 특징은 두 그룹변수들의 효과뿐만 아니라 두 그룹 변수들이 서로 어떤 영향을 미치는지 교호작용을 볼 수 있다는 것이다

 

1) 교호 작용이 있는 경우: lm(y~변수a*변수b, data=)

2) 교호 작용이 없는 경우: summary(lm(y~변수a+변수b, data = )

 

예시) 각기 다른 사료 a1~a4와 돼지품종 b1~3이 체중 증가에 미치는 영향을 조사

　

b1

b2

b3

a1

64

66

70

72

81

64

74

51

65

a2

65

63

58

57

43

52

47

58

67

a3

59

68

65

66

71

59

58

39

42

a4

58

41

46

57

61

53

53

59

38

 

> setwd('c:/Rwork')

> pigs<-read.csv('pigs.csv')

> head(pigs)

result breed feed

1     64    b1   a1

2     66    b1   a1

3     70    b1   a1

4     65    b1   a2

5     63    b1   a2

6     58    b1   a2

> with(pigs,tapply(result,breed,mean))

b1       b2       b3

60.25000 61.33333 54.25000

> with(pigs,tapply(result,feed,mean))

a1       a2       a3       a4

67.44444 56.66667 58.55556 51.77778

> with(pigs,tapply(result,list(breed,feed),mean))

a1       a2       a3       a4

b1 66.66667 62.00000 64.00000 48.33333

b2 72.33333 50.66667 65.33333 57.00000

b3 63.33333 57.33333 46.33333 50.00000

> boxplot(result~breed+feed,col='blue',data=pigs)



 

기술통계를 위해 tapply 함수를 활용하여 품종, 사료, 품종x사료 조합의 평균을 구하였고, boxplot 함수를 통해 각 변수 조합의 result 분포 차이를 나타냈다.


 

교호 작용 검사

> with(pigs,interaction.plot(x.factor=feed,trace.factor=breed,response=result,fun=mean,type='b',main='Interaction Plot',pch=c(1,2,16)))

 


 

#교호 작용이 있는 경우

> out=lm(result~breed*feed,data=pigs)

> out

 

 

 

> anova(out)

Analysis of Variance Table

 

Response: result

Df  Sum Sq Mean Sq F value   Pr(>F)  

breed       2  349.39  174.69  2.7926 0.081214 .

feed        3 1156.56  385.52  6.1628 0.002939 **

breed:feed   6  771.28  128.55  2.0549 0.097118 .

Residuals  24 1501.33   62.56                   

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

품종(breed)는 유의하지 않지만 p-value가 0.05에 근접하고, 사료(feed)는 유의하다. 품종에 따라 사료가 주는 효과의 크기가 조금씩 다를 수 있지만, Interaction인 breed:feed의 p-value는 0.09로 유의하지 않으므로 교호작용이 있다는 결론을 내릴 수 없다. 유의성 검정 결과를 보면 사료의 종류에 따라 체중증가에 유의한 차이가 있으나, 돼지 품종에 따라서는 체중증가에 차이가 있다는 결론은 내릴 수 없다. 품종에 대한 유의확률은 0.081로 약간의 영향력을 행사하는 듯하다.



#교호 작용이 없는 경우

 to be updated..



출처: 실험 계획과 응용, R로 하는 통계 분석


반응형

'KNOU > 2 실험 계획과 응용' 카테고리의 다른 글

6장 회귀분석과 공분산분석  (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
Posted by 마르띤
,

3장 일원배치법, 일원분산분석, one way anova

KNOU/2 실험 계획과 응용 2016. 10. 12. 10:11
반응형

제3장 일원배치 완전확률화 계획법,

 

예제) 직물의 긁힘에 대한 저항력을 측정하기 위해 원단 A1,2,3,4의 제품 마모도를 비교하조가 한다.

제품

반복

A1

1.93

2.38

2.2

2.25

A2

2.55

2.72

2.75

2.7

A3

2.4

2.68

2.32

2.28

A4

2.33

2.38

2.28

2.25

> setwd('c:/Rwork')

> fabric = read.csv('fabric.csv')

> head(fabric)

result group

1   1.93    a1

2   2.38    a1

3   2.20    a1

4   2.25    a1

5   2.55    a2

6   2.72    a2

> boxplot(result~group,data=fabric)

 





> with(fabric,tapply(result,group,mean))

a1   a2   a3   a4

2.19 2.68 2.42 2.31

> round(with(fabric,tapply(result,group,mean)))

a1 a2 a3 a4

2  3  2  2

> with(fabric,tapply(result,group,sd))

a1         a2         a3         a4

0.18920888 0.08906926 0.18036999 0.05715476

> out=lm(result~group,data=fabric)

> out

 

Call:

  lm(formula = result ~ group, data = fabric)

 

Coefficients:

  (Intercept)      groupa2      groupa3      groupa4 

2.19         0.49         0.23         0.12 

Lm의 결과물을 보면 intercept 추정치가 a1의 평균과 같으므로 a1이 reference로 사용되었음을알 수 있다. A1와 a2 평균이 각각 2.19, 2.68으로 그 차이가 groupa2 임을 알 수 있다.

 

> summary(out)

 

Call:

  lm(formula = result ~ group, data = fabric)

 

Residuals:

  Min      1Q  Median      3Q     Max

-0.2600 -0.0700  0.0150  0.0625  0.2600

 

Coefficients:

  Estimate Std. Error t value Pr(>|t|)   

(Intercept)  2.19000    0.07050  31.062 7.79e-13 ***

  groupa2      0.49000    0.09971   4.914 0.000357 ***

  groupa3      0.23000    0.09971   2.307 0.039710 * 

  groupa4      0.12000    0.09971   1.204 0.251982   

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 0.141 on 12 degrees of freedom

Multiple R-squared:  0.6871,           Adjusted R-squared:  0.6089

F-statistic: 8.785 on 3 and 12 DF,  p-value: 0.002353

 

> anova(out)

Analysis of Variance Table

 

Response: result

Df Sum Sq  Mean Sq F value   Pr(>F)  

group     3 0.5240 0.174667  8.7846 0.002353 **

Residuals 12 0.2386 0.019883                   

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

분산분석표에서 group의 F(df1=3,df2=12)=8.7846 (p-value=0.002353)으로 세 그룹의 평균들이 모두 같다는 귀무가설을 기각한다.

귀무가설 H0: μa1 = μa2 = μa3 = μa4

대립가설 H1: not H0

검정통계량 F(df1=3,df2=12)=8.7846

p-value = 0.002353 < 0.05

결정 : 모든 그룹의 평균이 같다는 귀무가설을 기각한다. 4개의 그룹 중 어느 쌍에 차이가 있는지 다중비교로 좀 더 진단할 수 있다.

 

> par(mfrow=c(2,2))

> plot(out)

> abline(out,col='blue')



Normal Q-Q Plot이 직선에 가깝고 다른 그림들도 별다른 추세를 보이지 않으므로 문제가 없다는 결론을 내린다.

 

> shapiro.test(resid(out))

 

Shapiro-Wilk normality test

 

data:  resid(out)

W = 0.97318, p-value = 0.8872

 

Shapiro.test()를 이용하여 잔차의 정규성검정 결과는 p-value=0.8872로 정규분포 가정을 만족시킨다.

 

모든 평균이 같다는 귀무가설이 기각되었다는 말은 그룹 중 최소한 하나는 0이 아니다라는 말이다. 어느 쌍의 차이로 귀무가설이 기각되었는지 조사하기 위해 다중비교를 한다. 분산분석에서 많이 쓰이는 다중비교 방법은 Dunnett와 Tukey이다. Tukey는 가능한 모든 조합의 쌍을, Dunnett는 하나의 대조군(reference)을 나머지 비교군(treatment)들과 비교한다.

 

> library(multcomp)

> dunnett=glht(out,linfct=mcp(group='Dunnett'))

> dunnett

 

General Linear Hypotheses

 

Multiple Comparisons of Means: Dunnett Contrasts

 

 

Linear Hypotheses:

  Estimate

a2 - a1 == 0     0.49

a3 - a1 == 0     0.23

a4 - a1 == 0     0.12

 

> summary(dunnett)

 

Simultaneous Tests for General Linear Hypotheses

 

Multiple Comparisons of Means: Dunnett Contrasts

 

 

Fit: lm(formula = result ~ group, data = fabric)

 

Linear Hypotheses:

  Estimate Std. Error t value Pr(>|t|)   

a2 - a1 == 0  0.49000    0.09971   4.914   <0.001 ***

a3 - a1 == 0  0.23000    0.09971   2.307   0.0962 . 

a4 - a1 == 0  0.12000    0.09971   1.204   0.5085   

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- single-step method)

> plot(dunnett)


a2-a1의 p-value만 0.001로 유의

 

> tukey=glht(out,linfct=mcp(group='Tukey'))

> tukey

 

General Linear Hypotheses

 

Multiple Comparisons of Means: Tukey Contrasts

 

 

Linear Hypotheses:

  Estimate

a2 - a1 == 0     0.49

a3 - a1 == 0     0.23

a4 - a1 == 0     0.12

a3 - a2 == 0    -0.26

a4 - a2 == 0    -0.37

a4 - a3 == 0    -0.11

 

> summary(tukey)

 

Simultaneous Tests for General Linear Hypotheses

 

Multiple Comparisons of Means: Tukey Contrasts

 

 

Fit: lm(formula = result ~ group, data = fabric)

 

Linear Hypotheses:

  Estimate Std. Error t value Pr(>|t|)  

a2 - a1 == 0  0.49000    0.09971   4.914  0.00187 **

a3 - a1 == 0  0.23000    0.09971   2.307  0.15096  

a4 - a1 == 0  0.12000    0.09971   1.204  0.63632  

a3 - a2 == 0 -0.26000    0.09971  -2.608  0.09215 .

a4 - a2 == 0 -0.37000    0.09971  -3.711  0.01375 *

a4 - a3 == 0 -0.11000    0.09971  -1.103  0.69444  

---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Adjusted p values reported -- single-step method)

> plot(tukey)


 

a2-a1, a4-a2간 95% 신뢰구간이 0을 포함하지 않으므로 유의한 차이가 있음을 알 수 있다.

 

분산분석은 모든 그룹의 평균이 같은지를 유의수준 0.05로 한 번에 검정한다. 그러나 여러 쌍의 평균이 같은지 다시 검정하면 비교하는 쌍이 하나 이상으로 늘어나 실제로 유의하지 않은데 유의하게 나올 가능성이 높아지므로 유의수준을 0.05보다 낮추거나, p-value를 높게 조종해야 한다. 아래 표처럼 실제 Tukey의 p-value가 Dunnett의 p-value보다 높게 조정되었음을 알 수 있다.

 

　

Dunnett

Tukey

a2 vs a1

0.001

0.00187

a3 vs a1

0.0962

0.15096

a4 vs a1

0.5085

0.63632

a3 vs a2

　

0.09215

a4 vs a2

　

0.01375

a4 vs a3

　

0.69444

 

출처: 실험 계획과 응용, R로 하는 통계 분석

반응형

'KNOU > 2 실험 계획과 응용' 카테고리의 다른 글

6장 회귀분석과 공분산분석  (0) 2016.10.17
4장 이원배치법, 이원분산분석, two way anova  (0) 2016.10.17
2장 짝지어진 비교 Paired T-Test  (1) 2016.10.12
2장 독립표본 두 모평균 차이 추론 two sample test  (0) 2016.10.12
Posted by 마르띤
,

2장 짝지어진 비교 Paired T-Test

KNOU/2 실험 계획과 응용 2016. 10. 12. 10:03
반응형

2.3 짝지어진 비교 Paired T-Test

Paired T-Test는 쌍을 이룬 두 변수의 차이를 보는 검정이다. 가장 흔한 예는 한 집단을 대상으로 어떤 개입의 효과를 보기 위해 개입 전-후 값을 비교하여 개입의 효과를 측정하는 것이다.

Paired t-test는 다음 순서를 따른다.


 1) 정규분포를 따르는지 검정 with(데이터, shapiro.test(사후-사전)

 2) 정규분포를 따르면 사후-사전 차이의 평균이 0인지 검정하는 one sample t-test를 이용

   with(데이터, t.test(사후-사전))

 3) 정규분포를 따르지 않으면 비모수적 방법을 적용 with(데이터, wilcox.test(사후-사전)

 

예) 운동화 밑창에 사용되는 B 재질과 A 재질 간 마모도에 차이가 있는지 알아보고자 한다. 이를 위하여 10명의 아이들을 대상으로 임의로 한 쪽 발에는 재질 A의 밑창을 단 운동화를 한 쪽 발에는 재질 B의 밑창을 단 운동화를 신겨 일정기간 사용하도록 한 다음 밑창이 얼마나 마모되었는지를 측정한 결과 다음 자료를 얻었다.

아이

재질A

재질B

차이B-A

1

13.2

14.0

0.8

2

8.2

8.8

0.6

3

10.9

11.2

0.3

4

14.3

14.2

-0.1

5

10.7

11.8

1.1

6

6.6

6.4

-0.2

7

9.5

9.8

0.3

8

10.8

11.3

0.5

9

8.8

9.3

0.5

10

13.3

13.6

0.3

　

　

　

평균차이=0.41

 

> shoes<-read.csv('shoes.csv',header=T)

> head(shoes)

A    B

1 13.2 14.0

2  8.2  8.8

3 10.9 11.2

4 14.3 14.2

5 10.7 11.8

6  6.6  6.4

#paired t-test 1. 정규분포를 따르는지 검정

> with(shoes,shapiro.test(B-A))

 

Shapiro-Wilk normality test

 

data:  B - A

W = 0.96132, p-value = 0.8009

 

귀무가설 H0: B-A의 차이가 정규분포를 따른다.

대립가설 H1: B-A의 차이가 정규분포를 따르지 않는다.

검정통계량 W = 0.96132

p-value = 0.8009

결정 : 귀무가설 기각할 수 없다. B-A의 차이가 정규 분포를 따른다.

 

#paired t-test 2. 정규분포를 따르면 one sample t-test를 이용

> with(shoes,t.test(B-A))

 

One Sample t-test

 

data:  B - A

t = 3.3489, df = 9, p-value = 0.008539

alternative hypothesis: true mean is not equal to 0

95 percent confidence interval:

  0.1330461 0.6869539

sample estimates:

  mean of x

0.41

귀무가설 H0: μ = 0

대립가설 H1: μ ≠ 0

검정통계량 t(df=9) = 3.3489

p-value = 0.008539

결정 : 귀무가설 기각한다. 두 신발 재질의 마모도에 유의한 차이가 있다.

 

#paired t-test 3. 정규분포를 따르지 않는 경우 wilcox.test를 사용한다.

> with(shoes,wilcox.test(B-A)

귀무가설 H0: 전후 차이의 median은 0이다

대립가설 H1: 전후 차이의 median은 0이 아니다

검정통계량 V =

p-value =

결정: if p-value > 0.05 then 귀무가설 기각할 수 없다, 전후 차이가 없다.


출처: 실험계획과 응용, R로 하는 통계 분석

반응형

'KNOU > 2 실험 계획과 응용' 카테고리의 다른 글

6장 회귀분석과 공분산분석  (0) 2016.10.17
4장 이원배치법, 이원분산분석, two way anova  (0) 2016.10.17
3장 일원배치법, 일원분산분석, one way anova  (0) 2016.10.12
2장 독립표본 두 모평균 차이 추론 two sample test  (0) 2016.10.12
Posted by 마르띤
,

2장 독립표본 두 모평균 차이 추론 two sample test

KNOU/2 실험 계획과 응용 2016. 10. 12. 09:59
반응형

2.2 독립표본을 이용한 두 모평균 차이에 대한 추론. Two-sample T-test

독립표본을 바탕으로 두 개의 모집단의 평균을 비교. 가장 흔한 실험 연구는 실험군과 대조군에 서로 다른 개입(intervention)을 적용시킨 후 두 집단의 평균이 같은지를 비교하여 개입 효과의 차이를 평가하는 것이다. 이 경우 two-sample t-test를 사용하는데, 서로 독립적인 두 변수 간에 차이의 평균이 0인지를 검정한다.

Two-sample t-test는 다음 순서를 따른다.

 1) 두 집단의 분산이 같은지 검정한다. var.test(y~그룹변수)

 2) 분산이 다르면 Welch의 t-test를 적용한다. t.test(y~그룹변수)

 3) 분산이 같으면 pooled variance를 이용한 t-test를 적용한다. t.test(y~그룹변수, var.equal=TRUE)

 

예제) 제약회사에서 어떤 약을 오래 보관해도 약효가 지속되는지를 검사하려고 한다. 표본1과 2를 랜덤추출한 결과가 아래와 같다.

표본1

10.2

10.5

10.3

10.8

9.8

10.6

10.7

10.2

10.0

10.1

표본2

9.8

9.6

10.1

10.2

10.1

9.7

9.5

9.6

9.8

9.9

 

> medical<-read.csv('medical.csv',header=T)

> head(medical,3)

sample result

1 sample1   10.2

2 sample1   10.5

3 sample1   10.3

> tail(medical,3)

sample result

18 sample2    9.6

19 sample2    9.8

20 sample2    9.9

> boxplot(result~sample,data=medical)


-> 해석: Sample1의 분산이 sample2보다 더 큼을 알 수 있다.

 

#two sample test 1. 등분산 검정

> var.test(result~sample,data=medical)

F test to compare two variances

 

data:  result by sample

F = 1.7965, num df = 9, denom df = 9, p-value = 0.3959

alternative hypothesis: true ratio of variances is not equal to 1

95 percent confidence interval:

  0.4462364 7.2328801

sample estimates:

  ratio of variances

1.796545  

귀무가설 H0: σ21 = σ22

대립가설 H1: σ21 ≠ σ22

검정통계량 F(df1=9, df2=9) = 1.7965

p-value = 0.3959

결정 : 귀무가설 기각할 수 없다. 등분산을 가정한다.

 

> 1/1.7965  #F값 1.7965의 역수는 0.55로 sample2의 분산이 sample1대비 0.55배임을 알 수 있다.

[1] 0.5566379

 

#two sample test 2. 분산이 같은 경우, pooled variance사용

> t.test(result~sample,var.equal=TRUE,data=medical)

 

Two Sample t-test

 

data:  result by sample

t = 3.8511, df = 18, p-value = 0.00117

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

  0.222688 0.757312

sample estimates:

  mean in group sample1 mean in group sample2

10.32                  9.83

 

귀무가설 H0: μ1 = μ2

대립가설 H1: μ1 ≠ μ2

검정통계량 t(df=18) = 3.8511

p-value = 0.00117

결정 : 귀무가설 기각한다. 두 모집단의 평균이 다르다.

 

#two sample test 3. 분산이 다른 경우, Welch의 t-test를 한다.

> t.test(result~sample,data=medical)

 

출처: 실험계획과 응용, R로 하는 통계 분석

반응형

'KNOU > 2 실험 계획과 응용' 카테고리의 다른 글

6장 회귀분석과 공분산분석  (0) 2016.10.17
4장 이원배치법, 이원분산분석, two way anova  (0) 2016.10.17
3장 일원배치법, 일원분산분석, one way anova  (0) 2016.10.12
2장 짝지어진 비교 Paired T-Test  (1) 2016.10.12
Posted by 마르띤
,

제2장 중회귀모형 - 표준화된 중회귀분석, 추정과 검정, 변수 추가

KNOU/2 회귀모형 2016. 9. 18. 19:22
반응형

2.4 표준화된 중회귀분석

표준화된 중회귀모형

Yi* = α1Zi1 + α2Zi2 + …+ αkZik + εi’

 

 

가 되며, 표준화된 중회귀모형에서 추정된 회귀계수 αi의 절대값이 크면 클수록 설명변수 Xi가 반응변수 Yi에 주는 영향이 크게 됨.

> library(lm.beta)

> market2.lm<-lm(Y~X1+X2,data=market2)

> market2.beta<-lm.beta(market2.lm)

> market2.beta

 

Call:

lm(formula = Y ~ X1 + X2, data = market2)

 

Standardized Coefficients::

(Intercept)          X1          X2

  0.0000000   0.7015566   0.3376137

 

> summary(market2.beta)

 

Call:

lm(formula = Y ~ X1 + X2, data = market2)

 

Residuals:

     Min       1Q   Median       3Q      Max

-1.30465 -0.69561 -0.01755  0.69003  1.53127

 

Coefficients:

            Estimate Standardized Std. Error t value Pr(>|t|)   

(Intercept)  0.85041      0.00000    0.84624   1.005 0.334770   

X1         1.55811      0.70156    0.14793  10.532 2.04e-07 ***

X2         0.42736      0.33761    0.08431   5.069 0.000276 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 0.9318 on 12 degrees of freedom

Multiple R-squared:  0.9799,    Adjusted R-squared:  0.9765

F-statistic: 292.5 on 2 and 12 DF,  p-value: 6.597e-11

-> R 결과에서 변수 X1, X2에 대한 표준화 계수는 0.7016, 0.3376이 됨을 알 수 있다. 따라서 적합된 표준화된계수 모형은

 

Ŷ* = 0.7016Z1 + 0.3376Z2

 

가 된다.여기서 X1의 계수가 X2의 계수보다 크므로 상대적으로 X1의 영향이 더 큼을 알 수 있다.

 

 

2.5 추정과 검정

 

Market2 데이터에 대하여

(1) X1=10, X2=10에 E(Y)=95% , 99%를 신뢰구간으로 추정하고

(2) H0 : β­1 = 0, H0 : β­2 = 0에 대하여 유의수준 α =0.05로 가설검정을 해 보자

 

#95%의 신뢰구간

> pred.x<-data.frame(X1=10,X2=10)

> pred.x

  X1 X2

1 10 10

> pc=predict(market2.lm,int='c',newdata=pred.x)

> pc

       fit      lwr      upr

1 20.70503 19.95796 21.45209

> round(pc,3)

     fit    lwr    upr

1 20.705 19.958 21.452

-> X1=10, X2=10 에서의 추정값은 20.705이고 95%의 신뢰구간은 (19.958,21.452)가 된다.

 

#99%의 신뢰구간

> pc99<-predict(market2.lm,int='c',level=0.99,newdata=pred.x)

> pc99

       fit      lwr      upr

1 20.70503 19.65769 21.75236

> round(pc99,3)

     fit    lwr    upr

1 20.705 19.658 21.752

-> X1=10, X2=10 에서의 추정값은 20.705이고 99%의 신뢰구간은 (19.658,21.752)로 95%의 신뢰구간은 (19.958,21.452)보다 더 넓게 형성된다.

 

 

H0 : β­1 = 0, H0 : β­2 = 0에 대한 가설 검정은 회귀적합 결과를 이용하면 된다.

> summary(market2.lm)

 

Call:

lm(formula = Y ~ X1 + X2, data = market2)

 

Residuals:

     Min       1Q   Median       3Q      Max

-1.30465 -0.69561 -0.01755  0.69003  1.53127

 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)   

(Intercept)  0.85041    0.84624   1.005 0.334770   

X1         1.55811    0.14793  10.532 2.04e-07 ***

X2         0.42736    0.08431   5.069 0.000276 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 0.9318 on 12 degrees of freedom

Multiple R-squared:  0.9799,    Adjusted R-squared:  0.9765

F-statistic: 292.5 on 2 and 12 DF,  p-value: 6.597e-11

Beta hat 1 값은 1.55811 이고, 표준오차는 0.14793이므로

 

t-값은 1.55811 / 0.14793 = 10.53275

 

유의확률 p-값 2.04 X 107 이 되므로 귀무가설 H0 : β­1 = 0은 기각한다. R에서 ***는 p-값이 0.001보다 작은 경우로 표시되었음을 알 수 있다. H0 : β­2 = 0 의 가설 역시 p-값이 0.000276이므로 귀각한다.


 

2.6 변수 추가

중회귀모형을 적합할 때 어떤 특정한 변수를 회귀모형에 포함시키는 것이 바람직한가를 결정하고 싶은 경우가 있다. 이러한 경우 이 변수를 포함시키지 않고 구한 회귀제곱함(SSR)에서 이 변수를 포함시키고 구한 회귀제곱함(SSR)이 추가적으로 어느 정도 커졌는가를 검토하는것이 좋을 것이다. 이와 같은 경우에 추가적으로 증가된 제곱함을 추가제곱합(extra sum of squares)라고 부른다.

 

 

Health data 불러오기

X1: 몸무게(파운드),

X2: 분당 맥박수,

X3: 근력(들어올릴 수 있는 무게: 파운드),

X4: 1/4 마일 달리는 시간(초),

Y: 1마일 달리는 시간(초)

> health<-read.table('health.csv',header=T)

> head(health,2)

   번호.X1.X2.X3.X4.Y

1 1,217,67,260,91,481

2 2,141,52,190,66,292

> health<-read.table('health.csv',header=T,sep=',')

> head(health,2)

  번호  X1 X2  X3 X4   Y

1    1 217 67 260 91 481

2    2 141 52 190 66 292

> colnames(health[1])

[1] "번호"

> colnames(health)[1]<-'ID'

> colnames(health)[1]

[1] "ID"

> head(health)

  ID  X1 X2  X3 X4   Y

1  1 217 67 260 91 481

2  2 141 52 190 66 292

 

 

회귀분석 실시

> h4.lm=lm(Y~X1+X2+X3+X4,data=health)

> summary(h4.lm)

 

Call:

lm(formula = Y ~ X1 + X2 + X3 + X4, data = health)

 

Residuals:

   Min     1Q Median     3Q    Max

-54.49 -21.91  -5.10  18.66  45.69

 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)   

(Intercept)  -6.8606    59.0070  -0.116   0.9084   

X1            1.3826     0.2933   4.713 7.84e-05 ***

X2           -0.3745     0.8955  -0.418   0.6794   

X3           -0.5302     0.2571  -2.062   0.0497 * 

X4            3.6202     0.7573   4.781 6.58e-05 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 29.96 on 25 degrees of freedom

Multiple R-squared:  0.8396,    Adjusted R-squared:  0.814

F-statistic: 32.72 on 4 and 25 DF,  p-value: 1.332e-09

 

> anova(h4.lm)

Analysis of Variance Table

 

Response: Y

          Df Sum Sq Mean Sq F value    Pr(>F)   

X1         1  89117   89117 99.2926 3.444e-10 ***

X2         1   4680    4680  5.2142   0.03117 * 

X3         1   3165    3165  3.5260   0.07213 . 

X4         1  20513   20513 22.8548 6.578e-05 ***

Residuals 25  22438     898                     

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘

 

 

모형 적합 후 anova 함수 사용

> h1.lm=lm(Y~X1,data=health)

> h2.lm=lm(Y~X1+X4,data=health)

> h3.lm=lm(Y~X1+X3+X4,data=health)

> h4.lm=lm(Y~X1+X2+X3+X4,data=health)

> anova(h1.lm,h2.lm)

Analysis of Variance Table

 

Model 1: Y ~ X1

Model 2: Y ~ X1 + X4

  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    

1     28 50795                                 

2     27 26419  1     24376 24.912 3.119e-05 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1

-> 모형1의 잔차제곱합 50795, 모형2의 잔차제곱함26419, 추가제곱합 50795-26419 = 24376 이고, 검정통계랑 F0는 24.912가 된다. 이에 대한 유의확률 p-값은 3.119 X 105이므로 변수4가 유의한 변수임을 알 수 있다.

F0 =

[SSE(R) – SSE(F)] / (dfR – dfF)

=

50795-26419 / 28-27

=24.91207

SSE(F)/dfF

26419/27

 

 

> anova(h2.lm,h3.lm)

Analysis of Variance Table

 

Model 1: Y ~ X1 + X4

Model 2: Y ~ X1 + X3 + X4

  Res.Df   RSS Df Sum of Sq      F  Pr(>F) 

1     27 26419                             

2     26 22595  1    3824.4 4.4007 0.04579 *

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

> anova(h3.lm,h4.lm)

Analysis of Variance Table

 

Model 1: Y ~ X1 + X3 + X4

Model 2: Y ~ X1 + X2 + X3 + X4

  Res.Df   RSS Df Sum of Sq      F Pr(>F)

1     26 22595                          

2     25 22438  1    156.93 0.1748 0.6794

-> (X1,X4)인 모형에 X3을 추가하는 경우의 추가제곱합 3842.4 p-값은 0.04579로 X3은 유의한 변수임을 알 수 있고, (X1,X3,X4)인 모형에 X2을 추가하는 경우는 추가제곱합 156.93, p-값은 0.6794로 X2은 유의미한 변수가 아님을 알 수 있다.

 


추가변수그림

새로운 변수의 효과를 그래프로 표현할 수 있는데, 이러한 그래 중의 하나가 추가변수그림(added variable plot)이다.


> library(car)

> h4.lm=lm(Y~X1+X2+X3+X4,data=health)

> summary(h4.lm)

 

Call:

lm(formula = Y ~ X1 + X2 + X3 + X4, data = health)

 

Residuals:

   Min     1Q Median     3Q    Max

-54.49 -21.91  -5.10  18.66  45.69

 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)   

(Intercept)  -6.8606    59.0070  -0.116   0.9084   

X1            1.3826     0.2933   4.713 7.84e-05 ***

X2           -0.3745     0.8955  -0.418   0.6794   

X3           -0.5302     0.2571  -2.062   0.0497 * 

X4            3.6202     0.7573   4.781 6.58e-05 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 29.96 on 25 degrees of freedom

Multiple R-squared:  0.8396,    Adjusted R-squared:  0.814

F-statistic: 32.72 on 4 and 25 DF,  p-value: 1.332e-09

> avPlots(h4.lm)


 


-> 해석: 중회귀모형에서 어떤 특정한 회귀모형에 포함시키조가 할 떄, 변수 선택은 기존의 모형이 설명하지 못하는 부분을 새로운 변수가 들어옴으로써 추가설명력이 얼마나 유의한가에 따라 결정된다. 변수 X1의 추가변수그림에서 회귀계수는 1.3826, X4는 3.6202, 변수X1와 X에 대하 추가변수 그림이 선형성이 강한 것을 볼 수 이씅며, 따라서 이 두 변수가 회귀모형에 매우 유의하다는 것을 알 수있다. 

반응형

'KNOU > 2 회귀모형' 카테고리의 다른 글

제2장 중회귀모형 - 중회귀모형 추정, 신뢰성  (0) 2016.09.14
제1장 단순회귀모형 - 연습문제  (0) 2016.09.14
제1장 단순회귀모형 - 단손회귀의 추정과 검정  (0) 2016.09.12
제1장 단순회귀모형 - 회귀모형, 회귀선의 추정, 회귀모형의 정도  (0) 2016.08.26
Posted by 마르띤
,

제2장 회귀모형 - 로지스틱 회귀모형 연습

KNOU/2 데이터마이닝 2016. 9. 14. 10:24
반응형

목표변수가 연속형인 경우 -> 선형 회귀모델, ex) 광고비 투입 대비 매출액

목표변수가 두 개의 범주를 가진 이항형인 경우 -> 로지스틱 회귀모형, ex) 좋다1, 나쁘다0


독일신용평가 데이터 셋

독일신용평가 데이터(German Credit Data)는 머신러닝 저장소에 탑재되어 있는 데이터로 분류의 예제로 많이 활용된다.


변수명

속성

변수 설명

check

범주형

자유예금형태
Status of existing checking account

duration

수치형

기간
Duration in month

history

범주형

과거신용정보
Credit history

purpose

범주형

목적
Purpose

credit

수치형

신용대출금액
Credit amount

savings

범주형

저축예금/채권
Savings account/bonds

employment

범주형

현직장 재직기간
Present employment since

installment

수치형

가처분소득 대비 적금비율
Installment rate in percentage of disposable income

personal

범주형

결혼상황 및 성별
Personal status and sex

debtors

범주형

여타 채무/채권
Other debtors / guarantors

residence

수치형

현 거주기간
Present residence since

property

범주형

재산
Property

age

수치형

나이
Age in years

others

범주형

여타적금
Other installment plans

housing

범주형

주거형태
Housing

numcredits

수치형

해당 은행 신용계좌 수
Number of existing credits at this bank

job

범주형

직업
Job

residpeople

수치형

부양가족수
Number of people being liable to provide maintenance for

telephone

범주형

전화소유
Telephone

foreign

범주형

외국인 노동자 여부
foreign worker

y

범주형

신용등급 양호 또는 불량
credit:Good or Bad

 

 



1. 데이터 불러오기

> setwd('c:/Rwork')

> german<-read.table('germandata.txt')

> head(german,2) #열 값들의 변수명이 없음.


> names<-c("check","duration","history","purpose","credit","savings","employment","installment",         "personal",      "debtors",       "residence",     "property",       "age",   "others",         "housing",       "numcredits",    "job",   "residpeople",    "telephone",     "foreign"         ,"y")

> colnames(german)<-names

> head(german,2)


> german$y<-factor(german$y,levels=c(1,2),labels=c('good','bad'))

> head(german,2)


> summary(german)


#  residence,numcredits,residpeople는 실제 범주형이지만 수치형으로 인식. 범주형으로 변환 필요 

> class(german$residence) #integer 수치형

[1] "integer"

> class(german$check) #factor 범주형

[1] "factor"

> german$residence = factor(german$residence)

> german$numcredits = factor(german$numcredits)

> german$residpeople = factor(german$residpeople)

> class(german$residence) #factor로 변환

[1] "factor"

> class(german$numcredits) #factor로 변환

[1] "factor"

> class(german$residpeople) #factor로 변환

[1] "factor"

> table(german$residence)

1   2   3   4

130 308 149 413

> german$y<-ifelse(german$y=='good',1,0) #반응 값 good은 1로, bad는 2로 변환

 


2. 로지스틱 회귀 분석 시작

> fit.all = glm(y~.,family = binomial,data=german) #로지스틱 회귀 분석

 

또는 아래와 같은 방법도 가능하다.

> gmn<-names(german)
> f<-as.formula(paste('y~',paste(gmn[!gmn%in%y],collapse='+')))
> fit.all.1<-glm(f,family = binomial, data=german)

> fit.step = step(fit.all, direction='both') #단계적 선택방법

Start:  AIC=993.44

y ~ check + duration + history + purpose + credit + savings +

  employment + installment + personal + debtors + residence +

  property + age + others + housing + numcredits + job + residpeople +

  telephone + foreign

 

Df Deviance     AIC

- job          3   888.00  988.00

- numcredits   3   890.25  990.25

- property     3   890.70  990.70

- residpeople  1   888.52  992.52

- age          1   889.37  993.37

- telephone    1   889.40  993.40

<none>             887.44  993.44

- employment   4   895.48  993.48

- housing      2   891.63  993.63

- residence    3   894.74  994.74

- debtors      2   894.80  996.80

- others       2   895.71  997.71

- personal     3   897.80  997.80

- foreign      1   894.16  998.16

- credit       1   895.07  999.07

- duration     1   896.25 1000.25

- installment  1   900.81 1004.81

- savings      4   908.55 1006.55

- history      4   911.01 1009.01

- purpose      9   922.07 1010.07

- check        3   957.33 1057.33

 

Step:  AIC=988

y ~ check + duration + history + purpose + credit + savings +

  employment + installment + personal + debtors + residence +

  property + age + others + housing + numcredits + residpeople +

  telephone + foreign

 

Df Deviance     AIC

- numcredits   3   890.85  984.85

- property     3   891.21  985.21

- residpeople  1   889.08  987.08

- employment   4   895.67  987.67

<none>             888.00  988.00

- housing      2   892.01  988.01

- age          1   890.05  988.05

- telephone    1   890.34  988.34

- residence    3   895.32  989.32

- debtors      2   895.25  991.25

- personal     3   898.31  992.31

- others       2   896.49  992.49

- foreign      1   894.77  992.77

+ job          3   887.44  993.44

- credit       1   895.72  993.72

- duration     1   897.14  995.14

- installment  1   901.56  999.56

- savings      4   909.71 1001.71

- history      4   911.44 1003.44

- purpose      9   922.89 1004.89

- check        3   957.60 1051.60

 

Step:  AIC=984.85

y ~ check + duration + history + purpose + credit + savings +

  employment + installment + personal + debtors + residence +

  property + age + others + housing + residpeople + telephone +

  foreign

 

Df Deviance     AIC

- property     3   894.03  982.03

- employment   4   898.02  984.02

- residpeople  1   892.07  984.07

- age          1   892.85  984.85

<none>             890.85  984.85

- housing      2   895.09  985.09

- telephone    1   893.29  985.29

- residence    3   898.52  986.52

+ numcredits   3   888.00  988.00

- debtors      2   898.27  988.27

- personal     3   901.17  989.17

- others       2   899.85  989.85

- foreign      1   898.00  990.00

+ job          3   890.25  990.25

- credit       1   898.64  990.64

- duration     1   899.76  991.76

- installment  1   904.66  996.66

- history      4   911.95  997.95

- savings      4   912.53  998.53

- purpose      9   926.15 1002.15

- check        3   959.38 1047.38

 

Step:  AIC=982.03

y ~ check + duration + history + purpose + credit + savings +

  employment + installment + personal + debtors + residence +

  age + others + housing + residpeople + telephone + foreign

 

Df Deviance     AIC

- residpeople  1   895.11  981.11

- employment   4   901.94  981.94

- telephone    1   895.95  981.95

<none>             894.03  982.03

- age          1   896.10  982.10

- housing      2   898.15  982.15

- residence    3   901.53  983.53

+ property     3   890.85  984.85

+ numcredits   3   891.21  985.21

- personal     3   903.97  985.97

- debtors      2   902.35  986.35

- foreign      1   901.07  987.07

+ job          3   893.45  987.45

- others       2   903.55  987.55

- credit       1   902.94  988.94

- duration     1   903.85  989.85

- installment  1   908.62  994.62

- savings      4   915.22  995.22

- history      4   915.59  995.59

- purpose      9   930.66 1000.66

- check        3   964.51 1046.51

 

Step:  AIC=981.11

y ~ check + duration + history + purpose + credit + savings +

  employment + installment + personal + debtors + residence +

  age + others + housing + telephone + foreign

 

Df Deviance     AIC

- employment   4   903.04  981.04

- age          1   897.04  981.04

<none>             895.11  981.11

- telephone    1   897.12  981.12

- housing      2   899.31  981.31

+ residpeople  1   894.03  982.03

- residence    3   902.80  982.80

- personal     3   904.04  984.04

+ property     3   892.07  984.07

+ numcredits   3   892.19  984.19

- debtors      2   903.15  985.15

- foreign      1   902.06  986.06

+ job          3   894.59  986.59

- others       2   904.70  986.70

- credit       1   903.73  987.73

- duration     1   904.80  988.80

- installment  1   909.03  993.03

- savings      4   916.06  994.06

- history      4   916.94  994.94

- purpose      9   932.01 1000.01

- check        3   965.87 1045.87

 

Step:  AIC=981.04

y ~ check + duration + history + purpose + credit + savings +

  installment + personal + debtors + residence + age + others +

  housing + telephone + foreign

 

Df Deviance     AIC

- age          1   904.91  980.91

<none>             903.04  981.04

+ employment   4   895.11  981.11

- telephone    1   905.28  981.28

- housing      2   907.58  981.58

+ residpeople  1   901.94  981.94

- residence    3   910.50  982.50

+ property     3   899.28  983.28

+ numcredits   3   900.64  984.64

- foreign      1   909.67  985.67

- debtors      2   912.24  986.24

+ job          3   902.89  986.89

- personal     3   915.04  987.04

- others       2   913.21  987.21

- duration     1   911.34  987.34

- credit       1   911.50  987.50

- installment  1   917.92  993.92

- savings      4   925.25  995.25

- history      4   925.74  995.74

- purpose      9   939.70  999.70

- check        3   975.57 1047.57

 

Step:  AIC=980.91

y ~ check + duration + history + purpose + credit + savings +

  installment + personal + debtors + residence + others + housing +

  telephone + foreign

 

Df Deviance     AIC

<none>             904.91  980.91

+ age          1   903.04  981.04

+ employment   4   897.04  981.04

- telephone    1   907.69  981.69

+ residpeople  1   903.95  981.95

- housing      2   910.11  982.11

- residence    3   912.96  982.96

+ property     3   901.18  983.18

+ numcredits   3   902.60  984.60

- foreign      1   911.56  985.56

- debtors      2   914.35  986.35

- others       2   914.61  986.61

+ job          3   904.63  986.63

- credit       1   913.18  987.18

- personal     3   917.50  987.50

- duration     1   914.06  988.06

- installment  1   919.35  993.35

- savings      4   927.70  995.70

- history      4   928.79  996.79

- purpose      9   940.82  998.82

- check        3   978.40 1048.40

 

> fit.step$anova #제거된 변수 보기

Step Df  Deviance Resid. Df Resid. Dev      AIC

1               NA        NA       947   887.4372 993.4372

2         - job  3 0.5588674       950   887.9960 987.9960

3  - numcredits  3 2.8582392       953   890.8543 984.8543

4    - property  3 3.1777611       956   894.0320 982.0320

5 - residpeople  1 1.0747973       957   895.1068 981.1068

6  - employment  4 7.9298736       961   903.0367 981.0367

7         - age  1 1.8704615       962   904.9072 980.9072

> summary(fit.step) #최종모델

 

Call:

  glm(formula = y ~ check + duration + history + purpose + credit +

      savings + installment + personal + debtors + residence +

      others + housing + telephone + foreign, family = binomial,

      data = german)

 

Deviance Residuals:

  Min       1Q   Median       3Q      Max 

-2.7904  -0.7290   0.3885   0.6911   2.1780 

 

Coefficients:

              Estimate Std. Error z value Pr(>|z|)   

(Intercept)   -9.736e-01  7.032e-01  -1.385 0.166204   

checkA12       3.863e-01  2.136e-01   1.809 0.070468 . 

checkA13       1.055e+00  3.636e-01   2.902 0.003714 **

checkA14       1.782e+00  2.308e-01   7.721 1.15e-14 ***

duration      -2.726e-02  9.034e-03  -3.018 0.002546 **

historyA31     1.290e-01  5.297e-01   0.244 0.807596   

historyA32     8.608e-01  4.104e-01   2.097 0.035956 * 

historyA33     9.975e-01  4.675e-01   2.133 0.032889 * 

historyA34     1.564e+00  4.329e-01   3.612 0.000303 ***

purposeA41     1.591e+00  3.684e-01   4.320 1.56e-05 ***

purposeA410    1.397e+00  7.732e-01   1.806 0.070849 . 

purposeA42     6.766e-01  2.529e-01   2.675 0.007467 **

purposeA43     8.867e-01  2.443e-01   3.629 0.000284 ***

purposeA44     5.231e-01  7.546e-01   0.693 0.488206   

purposeA45     1.335e-01  5.388e-01   0.248 0.804301   

purposeA46    -2.006e-01  3.883e-01  -0.517 0.605426   

purposeA48     2.060e+00  1.202e+00   1.714 0.086523 . 

purposeA49     7.396e-01  3.318e-01   2.229 0.025820 * 

credit        -1.230e-04  4.314e-05  -2.852 0.004351 **

savingsA62     3.126e-01  2.805e-01   1.115 0.264984   

savingsA63     4.303e-01  3.887e-01   1.107 0.268291   

savingsA64     1.396e+00  5.184e-01   2.692 0.007106 **

savingsA65     1.004e+00  2.606e-01   3.852 0.000117 ***

installment   -3.218e-01  8.621e-02  -3.733 0.000189 ***

personalA92    2.159e-01  3.754e-01   0.575 0.565268   

personalA93    8.302e-01  3.672e-01   2.261 0.023766 * 

personalA94    3.551e-01  4.434e-01   0.801 0.423122   

debtorsA102   -4.978e-01  4.005e-01  -1.243 0.213967   

debtorsA103    1.074e+00  4.205e-01   2.555 0.010628 * 

residence2    -7.181e-01  2.796e-01  -2.568 0.010223 * 

residence3    -3.929e-01  3.246e-01  -1.210 0.226104   

residence4    -2.893e-01  2.806e-01  -1.031 0.302546   

othersA142     5.959e-02  4.061e-01   0.147 0.883344   

othersA143     6.787e-01  2.355e-01   2.882 0.003955 **

housingA152    5.098e-01  2.271e-01   2.245 0.024799 * 

housingA153    2.464e-01  3.288e-01   0.749 0.453710   

telephoneA192  3.051e-01  1.838e-01   1.660 0.096958 . 

foreignA202    1.439e+00  6.253e-01   2.301 0.021383 * 

  ---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

(Dispersion parameter for binomial family taken to be 1)

 

Null deviance: 1221.73  on 999  degrees of freedom

Residual deviance:  904.91  on 962  degrees of freedom

AIC: 980.91

 

Number of Fisher Scoring iterations: 5



->  <!--[endif]-->해석: fit.step = step(fit.all, direction='both'를 통해 AIC가 가장 작은 모형을 찾는다.

check는 4개의 범주(checkA11 계좌 없음 / A12 잔액 없음 / A13 잔액 200 이하 / A14 잔액 200 이상)를 가지므로 3개의 가변 수 생성. 추정된 회귀계수는 모두 양수이므로, A12~A14 즉 계좌가 있는 경우 계좌 없음(A11)대비 신용이 좋을 확률(Y=1)이 더 높다. 대출기간인 duration은 마이너스의 값을 지니므로 대출 기간이 오래 될 수록 신용도는 낮아진다. 모델의 AIC는 980.91로, AIC가 클 경우 그 모형은 적합하지 않기 때문에, 여러 후보 모형 중에서 AIC가 가장 작은 모형을 선택한다.

 

 

단계적선택법의 AIC는 980.91

[참고] 후진소거법의 AIC는 980.91

> fit.step.back = step(fit.all,direction='backward')

Step:  AIC=980.91

y ~ check + duration + history + purpose + credit + savings +

    installment + personal + debtors + residence + others + housing +

    telephone + foreign

 

              Df Deviance     AIC

<none>             904.91  980.91

- telephone    1   907.69  981.69

- housing      2   910.11  982.11

- residence    3   912.96  982.96

- foreign      1   911.56  985.56

- debtors      2   914.35  986.35

- others       2   914.61  986.61

- credit       1   913.18  987.18

- personal     3   917.50  987.50

- duration     1   914.06  988.06

- installment  1   919.35  993.35

- savings      4   927.70  995.70

- history      4   928.79  996.79

- purpose      9   940.82  998.82

- check        3   978.40 1048.40

 

> fit.step.back$anova #제거된 변수 보기

           Step Df  Deviance Resid. Df Resid. Dev      AIC

1               NA        NA       947   887.4372 993.4372

2         - job  3 0.5588674       950   887.9960 987.9960

3  - numcredits  3 2.8582392       953   890.8543 984.8543

4    - property  3 3.1777611       956   894.0320 982.0320

5 - residpeople  1 1.0747973       957   895.1068 981.1068

6  - employment  4 7.9298736       961   903.0367 981.0367

7         - age  1 1.8704615       962   904.9072 980.9072

 

 

[참고] 전진선택법 AIC : 993.44

> fit.step.forward = step(fit.all, direction = 'forward')

Start:  AIC=993.44

y ~ check + duration + history + purpose + credit + savings +

    employment + installment + personal + debtors + residence +

    property + age + others + housing + numcredits + job + residpeople +

telephone + foreign

 

> fit.step.forward$anova #제거된 변수 보기

  Step Df Deviance Resid. Df Resid. Dev      AIC

1      NA       NA       947   887.4372 993.4372

 


 

3. 예측함수 및 정오분류표 작성

> p = predict(fit.step, newdata=german,type='response')

> threshold = 0.5 #cutoff기준 0.5로 정함

> yhat = ifelse(p>threshold,1,0)

> head(yhat)

  1 2 3 4 5 6

  1 0 1 1 0 1

> class.tab = table(german$y,yhat,dnn=c("Actual","Predicted"))#실값과 예측값 배열

> class.tab

       Predicted

Actual   0   1

     0 158 142

     1  82 618

->  해석: 1로 예측할 확률이 임계치(threshold) 0.5보다 클 경우에는 1로, 0.5이하일 경우에는 0으로 예측. 실제로는 0인데 0으로 예측한 경우가 158개, 1인데 1로 분류한 경우가 618개이다.반면에 0인데 1로 오분류한 경우가 142개, 1인데 0으로 오분류한 경우가 82개이다.


 

4. 예측력 측도

> sum(german$y==yhat)/length(german$y) #Prediction Accuracy 예측정확도

[1] 0.776

> sum(german$y!=yhat)/length(german$y) #Misclassification Rate 오분류율

[1] 0.224

> class.tab[1,1]/apply(class.tab,1,sum)[1] #Specificity 특이도

0

0.5266667

> class.tab[2,2]/apply(class.tab,1,sum)[2] #Sensitivity 민감도

1

0.8828571

-> 해석: 민감도는 실제 양성(Y=1)일 때 양성으로 예측할 확률, 특이도는 실제 음성(Y=0)일 때 음성으로 예측할 확률이다. 예측정확도(prediction accuracy)는 실제 양서일 때 양성으로, 음성일 때 음성으로 제대로 예측할 확률로 민감도와 특이도의 가중평균이다. 오분류율(misclassification rate)는 양성일 때 음성으로, 음성일 때 양성으로 잘못 예측할 확률이다.


 

5. ROC 곡선 및 AUC 생성

 

> library(ROCR)

> pred<-prediction(p,german$y)

> perf<-performance(pred,'tpr','fpr') #민감도와 1-특이도 계산 과정

> plot(perf,lty=1,col=2,xlim=c(0,1),ylim=c(0,1),xlab='1-Specificity',ylab='Sensitivity',main='ROC Curve')

> lines(x=c(0,1),y=c(0,1),col='grey')


> performance(pred,'auc')@y.values #면적 계산

[[1]]

[1] 0.8312286


->  민감도와 특이도는 임계치에 다라 달라지고 임계치는 상황에 따라 다르게 결정할 수 이다. 여러 가능한 임계치에 대해 ‘1-특이도(Specificity)’를 가로축에, 민감도를 세로축에 놓고 그린 그래프를 ROC(Receiver operating characteristic) 곡선이라 한다. 민감도와 특이도가 높을수록 예측력이 좋다고 할 수 있기 때문에 ROC 곡선이 좌상단에 가까울수록 ROC 곡선 아래 면적인 AUC(area under the ROC curve)가 커지고, 예측력이 좋다고 할 수 있다.이 독일신용평가 데이터에 적합한 로지스틱 회귀모형에 대한 예측력의 측도인 AUC는 최대치 1보다 다소 작은 0.831로 상당히 높음을 알 수 있다.

반응형

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제3장 나무모형 - 회귀나무모형  (0) 2016.10.26
제3장 나무모형 - 분류나무모형  (4) 2016.10.18
제2장 회귀모형 - 선형회귀 연습  (0) 2016.09.14
제2장 회귀모형 - 선형회귀, 로지스틱회귀  (0) 2016.09.14
1장 데이터과학과 데이터마이닝  (0) 2016.08.26
Posted by 마르띤
,

제2장 회귀모형 - 선형회귀 연습

KNOU/2 데이터마이닝 2016. 9. 14. 10:10
반응형

목표변수가 연속형인 경우 -> 선형 회귀모델, ex) 광고비 투입 대비 매출액

목표변수가 두 개의 범주를 가진 이항형인 경우 -> 로지스틱 회귀모형, ex) 좋다1, 나쁘다0



보스턴 하우징 데이터 Housing Values in Suburbs of Boston

(출처: http://127.0.0.1:31866/library/MASS/html/Boston.html)


변수명

속성

변수 설명

crim

수치형(numeric)

per capita crime rate by town
타운별 1인당 범죄율

zn

수치형(numeric)

proportion of residential land zoned for lots over 25,000 sq.ft.
25,000평방피트를 초과하는 거주지역 비율

indus

수치형(numeric)

proportion of non-retail business acres per town.
비소매 사업지역의 토지 비율

chas

범주형(integer)

Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).
찰스강 더비 변수 (강의 경계에 위치 = 1, 아니면 = 0)

nox

수치형(numeric)

nitrogen oxides concentration (parts per 10 million).
10ppm당 농축 일산화질소

rm

수치형(numeric)

average number of rooms per dwelling.
주택 1가구등 방의 평균 개수

age

수치형(numeric)

proportion of owner-occupied units built prior to 1940.
1940년 이전에 건축된 소유자 주택 비율

dis

수치형(numeric)

weighted mean of distances to five Boston employment centres.
5개의 보스턴 고용센터까지의 접근성 지수

rad

범주형(integer)

index of accessibility to radial highways.
방사형 도로까지의 접근성 지수

tax

수치형(numeric)

full-value property-tax rate per \$10,000.
10,000달러당 재산세율

ptratio

수치형(numeric)

pupil-teacher ratio by town.
타운별 학생/교사 비율

black

수치형(numeric)

1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.
타운별 흑인의 비율

lstat

수치형(numeric)

lower status of the population (percent).
모집단의 하위계층의 비율

medv
(목표변수)

수치형(numeric)

median value of owner-occupied homes in \$1000s.
본인 소유의 주택가격(중앙값)

 

 

 

 

1. 데이터 불러오기

> library(MASS)

> range(Boston$medv)

[1]  5 50

> stem(Boston$medv)

The decimal point is at the | 

4 | 006

6 | 30022245

8 | 1334455788567

10 | 2224455899035778899

12 | 013567778011112333444455668888899

14 | 0111233445556689990001222344666667

16 | 01112234556677880111222344455567888889

18 | 01222334445555667778899990011112233333444444555566666778889999

20 | 0000011111223333444455566666677888990001122222444445566777777788999

22 | 00000001222223344555666667788889999000011111112222333344566777788889

24 | 001112333444455566777888800000000123

26 | 24456667011555599

28 | 01244567770011466889

30 | 111357801255667

32 | 0024579011223448

34 | 679991244

36 | 01224502369

38 | 78

40 | 37

42 | 38158

44 | 084

46 | 07

48 | 358

50 | 0000000000000000

 

> i=which(Boston$medv==50)#본인 소유의 주택가격(중앙값)

> Boston[i,]

        crim zn indus chas    nox    rm   age    dis rad tax ptratio  black lstat medv

162 1.46336  0 19.58    0 0.6050 7.489  90.8 1.9709   5 403    14.7 374.43  1.73   50

163 1.83377  0 19.58    1 0.6050 7.802  98.2 2.0407   5 403    14.7 389.61  1.92   50

164 1.51902  0 19.58    1 0.6050 8.375  93.9 2.1620   5 403    14.7 388.45  3.32   50

167 2.01019  0 19.58    0 0.6050 7.929  96.2 2.0459   5 403    14.7 369.30  3.70   50

187 0.05602  0  2.46    0 0.4880 7.831  53.6 3.1992   3 193    17.8 392.63  4.45   50

196 0.01381 80  0.46    0 0.4220 7.875  32.0 5.6484   4 255    14.4 394.23  2.97   50

205 0.02009 95  2.68    0 0.4161 8.034  31.9 5.1180   4 224    14.7 390.55  2.88   50

226 0.52693  0  6.20    0 0.5040 8.725  83.0 2.8944   8 307    17.4 382.00  4.63   50

258 0.61154 20  3.97    0 0.6470 8.704  86.9 1.8010   5 264    13.0 389.70  5.12   50

268 0.57834 20  3.97    0 0.5750 8.297  67.0 2.4216   5 264    13.0 384.54  7.44   50

284 0.01501 90  1.21    1 0.4010 7.923  24.8 5.8850   1 198    13.6 395.52  3.16   50

369 4.89822  0 18.10    0 0.6310 4.970 100.0 1.3325  24 666    20.2 375.52  3.26   50

370 5.66998  0 18.10    1 0.6310 6.683  96.8 1.3567  24 666    20.2 375.33  3.73   50

371 6.53876  0 18.10    1 0.6310 7.016  97.5 1.2024  24 666    20.2 392.05  2.96   50

372 9.23230  0 18.10    0 0.6310 6.216 100.0 1.1691  24 666    20.2 366.15  9.53   50

373 8.26725  0 18.10    1 0.6680 5.875  89.6 1.1296  24 666    20.2 347.88  8.88   50

 

> boston=Boston[-i,] #최대값 50인 관측치 16개를 찾아 제거

> boston$chas = factor(boston$chas) #범주형으로 변경

> boston$rad = factor(boston$rad) #범주형으로 변경

> table(boston$rad)

1   2   3   4   5   6   7   8  24

19  24  37 108 109  26  17  23 127

 

> boston$chas <- as.factor(boston$chas)

> boston$rad <- as.factor(boston$rad)

> class(boston$rad);class(boston$chas)
[1] "factor"
[1] "factor"

 

[참고] 아래와 같은 방법으로 이용하면 모든 변수를 수치로 변경할 수 있다.

> for(i in 1:ncol(boston))if(!is.numeric(boston[,i])) boston[,i]=as.numeric(boston[,i])
> str(boston)
'data.frame':   490 obs. of  14 variables:
 $ crim   : num  0.00632 0.02731 0.02729 0.03237 0.06905 ...
 $ zn     : num  18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
 $ indus  : num  2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
 $ chas   : num  1 1 1 1 1 1 1 1 1 1 ...
 $ nox    : num  0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
 $ rm     : num  6.58 6.42 7.18 7 7.15 ...
 $ age    : num  65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
 $ dis    : num  4.09 4.97 4.97 6.06 6.06 ...
 $ rad    : num  1 2 2 3 3 3 5 5 5 5 ...
 $ tax    : num  296 242 242 222 222 222 311 311 311 311 ...
 $ ptratio: num  15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
 $ black  : num  397 397 393 395 397 ...
 $ lstat  : num  4.98 9.14 4.03 2.94 5.33 ...
 $ medv   : num  24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...

 

 

 

 

 

 

2. 선형 회귀 모형 만들기

#선형회귀모형 만들기

> fit1 = lm(medv~.,data=boston) #목표변수 = medv, 선형회귀모형 함수, ~.는 목표 변수를 제외한 모든 변수를 입력변수로 사용

> summary(fit1)

 

  Call:

    lm(formula = medv ~ ., data = boston)

 

  Residuals:

    Min      1Q  Median      3Q     Max

  -9.5220 -2.2592 -0.4275  1.6778 15.2894

 

  Coefficients:

    Estimate Std. Error t value Pr(>|t|)   

  (Intercept)  30.120918   4.338656   6.942 1.29e-11 ***

    crim         -0.105648   0.025640  -4.120 4.47e-05 ***

    zn            0.044104   0.011352   3.885 0.000117 ***

    indus        -0.046743   0.051044  -0.916 0.360274   

    chas1         0.158802   0.736742   0.216 0.829435   

    nox         -11.576589   3.084187  -3.754 0.000196 ***

    rm            3.543733   0.356605   9.937  < 2e-16 ***

    age          -0.026082   0.010531  -2.477 0.013613 * 

    dis          -1.282095   0.160452  -7.991 1.05e-14 ***

    rad2          2.548109   1.175012   2.169 0.030616 * 

    rad3          4.605849   1.064492   4.327 1.85e-05 ***

    rad4          2.663393   0.950747   2.801 0.005299 **

    rad5          3.077800   0.962725   3.197 0.001483 **

    rad6          1.314892   1.157689   1.136 0.256624   

    rad7          4.864208   1.241760   3.917 0.000103 ***

    rad8          5.772296   1.194221   4.834 1.82e-06 ***

    rad24         6.195415   1.417826   4.370 1.53e-05 ***

    tax          -0.009396   0.003070  -3.061 0.002333 **

    ptratio      -0.828498   0.114436  -7.240 1.85e-12 ***

    black         0.007875   0.002084   3.779 0.000178 ***

    lstat        -0.354606   0.041901  -8.463 3.36e-16 ***

    ---

    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

  Residual standard error: 3.671 on 469 degrees of freedom

  Multiple R-squared:  0.7911,     Adjusted R-squared:  0.7821

  F-statistic: 88.78 on 20 and 469 DF,  p-value: < 2.2e-16

 

또는 아래와 같은 방법도 가능하다

> names(boston)
 [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"   
 [8] "dis"     "rad"     "tax"     "ptratio" "black"   "lstat"   "medv"  
> bn <- names(boston)

> f <- as.formula(paste('medv~',paste(bn[!bn %in% 'medv'],collapse='+')))
> f
medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +
    tax + ptratio + black + lstat
> fit2 <- lm(f,data=boston)
> summary(fit2)

Call:
lm(formula = f, data = boston)

Residuals:
    Min      1Q  Median      3Q     Max
-9.5220 -2.2592 -0.4275  1.6778 15.2894 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept)  30.120918   4.338656   6.942 1.29e-11 ***
crim         -0.105648   0.025640  -4.120 4.47e-05 ***
zn            0.044104   0.011352   3.885 0.000117 ***
indus        -0.046743   0.051044  -0.916 0.360274   
chas2         0.158802   0.736742   0.216 0.829435   
nox         -11.576589   3.084187  -3.754 0.000196 ***
rm            3.543733   0.356605   9.937  < 2e-16 ***
age          -0.026082   0.010531  -2.477 0.013613 * 
dis          -1.282095   0.160452  -7.991 1.05e-14 ***
rad2          2.548109   1.175012   2.169 0.030616 * 
rad3          4.605849   1.064492   4.327 1.85e-05 ***
rad4          2.663393   0.950747   2.801 0.005299 **
rad5          3.077800   0.962725   3.197 0.001483 **
rad6          1.314892   1.157689   1.136 0.256624   
rad7          4.864208   1.241760   3.917 0.000103 ***
rad8          5.772296   1.194221   4.834 1.82e-06 ***
rad9          6.195415   1.417826   4.370 1.53e-05 ***
tax          -0.009396   0.003070  -3.061 0.002333 **
ptratio      -0.828498   0.114436  -7.240 1.85e-12 ***
black         0.007875   0.002084   3.779 0.000178 ***
lstat        -0.354606   0.041901  -8.463 3.36e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.671 on 469 degrees of freedom
Multiple R-squared:  0.7911,    Adjusted R-squared:  0.7821
F-statistic: 88.78 on 20 and 469 DF,  p-value: < 2.2e-16 

  

 #가장 적절한 모형 선택 위한 변수 선택

> fit.step = step(fit1,direction='both') #both는 단계적 선택법 적용

  Start:  AIC=1295.03

  medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +

    tax + ptratio + black + lstat

 

  Df Sum of Sq    RSS    AIC

  - chas     1      0.63 6321.5 1293.1

  - indus    1     11.30 6332.2 1293.9

  <none>                 6320.9 1295.0

  - age      1     82.67 6403.5 1299.4

  - tax      1    126.28 6447.1 1302.7

  - nox      1    189.88 6510.7 1307.5

  - black    1    192.42 6513.3 1307.7

  - zn       1    203.44 6524.3 1308.5

  - crim     1    228.82 6549.7 1310.5

  - rad      8    721.85 7042.7 1332.0

  - ptratio  1    706.41 7027.3 1344.9

  - dis      1    860.51 7181.4 1355.6

  - lstat    1    965.26 7286.1 1362.7

  - rm       1   1330.92 7651.8 1386.7

 

  Step:  AIC=1293.08

  medv ~ crim + zn + indus + nox + rm + age + dis + rad + tax +

    ptratio + black + lstat

 

  Df Sum of Sq    RSS    AIC

  - indus    1     11.00 6332.5 1291.9

  <none>                 6321.5 1293.1

  + chas     1      0.63 6320.9 1295.0

  - age      1     82.48 6404.0 1297.4

  - tax      1    130.45 6451.9 1301.1

  - nox      1    189.27 6510.8 1305.5

  - black    1    193.59 6515.1 1305.9

  - zn       1    203.76 6525.2 1306.6

  - crim     1    230.58 6552.1 1308.6

  - rad      8    738.26 7059.8 1331.2

  - ptratio  1    719.40 7040.9 1343.9

  - dis      1    861.64 7183.1 1353.7

  - lstat    1    965.11 7286.6 1360.7

  - rm       1   1333.37 7654.9 1384.9

 

  Step:  AIC=1291.93

  medv ~ crim + zn + nox + rm + age + dis + rad + tax + ptratio +

    black + lstat

 

  Df Sum of Sq    RSS    AIC

  <none>                 6332.5 1291.9

  + indus    1     11.00 6321.5 1293.1

  + chas     1      0.32 6332.2 1293.9

  - age      1     81.09 6413.6 1296.2

  - tax      1    192.78 6525.3 1304.6

  - black    1    196.55 6529.0 1304.9

  - zn       1    220.63 6553.1 1306.7

  - crim     1    225.50 6558.0 1307.1

  - nox      1    239.09 6571.6 1308.1

  - rad      8    791.09 7123.6 1333.6

  - ptratio  1    732.81 7065.3 1343.6

  - dis      1    857.27 7189.8 1352.1

  - lstat    1    987.73 7320.2 1361.0

  - rm       1   1380.21 7712.7 1386.5

> summary(fit.step) #최종모형, rad는 범주형 변수를 가변수로 변환한 것.#AIC가 가장 작은 변수가 단계적 선택법에 의해 변수들이 정의 됨

  Call:

    lm(formula = medv ~ crim + zn + nox + rm + age + dis + rad +

         tax + ptratio + black + lstat, data = boston)

 

  Residuals:

    Min      1Q  Median      3Q     Max

  -9.5200 -2.2850 -0.4688  1.7535 15.3972

 

  Coefficients:

    Estimate Std. Error t value Pr(>|t|)   

  (Intercept)    30.252522   4.329907   6.987 9.64e-12 ***

    crim         -0.104568   0.025533  -4.095 4.96e-05 ***

    zn            0.045510   0.011235   4.051 5.97e-05 ***

    nox         -12.366882   2.932651  -4.217 2.97e-05 ***

    rm            3.583130   0.353644  10.132  < 2e-16 ***

    age          -0.025822   0.010514  -2.456 0.014412 * 

    dis          -1.253903   0.157029  -7.985 1.08e-14 ***

    rad2          2.387130   1.160735   2.057 0.040278 * 

    rad3          4.644091   1.062157   4.372 1.51e-05 ***

    rad4          2.608777   0.944668   2.762 0.005977 **

    rad5          3.116933   0.960550   3.245 0.001258 **

    rad6          1.422890   1.150280   1.237 0.216705   

    rad7          4.868388   1.240114   3.926 9.94e-05 ***

    rad8          5.872144   1.180865   4.973 9.26e-07 ***

    rad24         6.420553   1.393304   4.608 5.24e-06 ***

    tax          -0.010571   0.002792  -3.787 0.000172 ***

    ptratio      -0.837356   0.113420  -7.383 7.08e-13 ***

    black         0.007949   0.002079   3.823 0.000149 ***

    lstat        -0.357576   0.041718  -8.571  < 2e-16 ***

    ---

    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

  Residual standard error: 3.667 on 471 degrees of freedom

  Multiple R-squared:  0.7907,     Adjusted R-squared:  0.7827

F-statistic: 98.83 on 18 and 471 DF,  p-value: < 2.2e-16

-> 결과 해석: 최초 만든 회귀함수 fit1 = lm(medv~.,data=boston)에서, 가장 적절한 모형 선택 위한 변수 선택을 위해 step 함수를 사용한다. fit.step = step(fit1,direction='both') 이 함수에서 both는 단계적 선택법 적용을 의미한다. 결과적으로 lm(formula = medv ~ crim + zn + nox + rm + age + dis + rad + tax + ptratio + black + lstat, data = boston)라는 모형이 만들어졌고, 최초 만든 모형 대비 indus, chas1 변수가 사라졌음을 알 수 있다. 또한 최종 모형에서 범주 rad2,3,4 등은 범주형 범수 중 특정 변수를 의미한다.


    입력변수 crim의 회귀계수 추정치는 음수이므로 crim이 증가함에 따라 목표변수medv는 감소한다. nox 변수의 회귀곗수는 -12인데, nox 변수가 올라갈 때 마다 medv 값은 내려간다. nox 변수는 10ppm당 농축 일산화질소를 뜻한다.

 

      rad변수는 9개 범주로 구성되어 있기 때문에 8개의 가변수가 생성되었다. 각 입력 변수의 t값의 절대값으 커서 대응하는 p-값은 0.05보다 작아서 유의하다고 할 수 있다. 단, rad6는 유의하지 않지만 다른 가변수가 유의하므로 제거되지 않고 여전히 모형에 포함된다. 

     

     R2은 79.07%로 적합한 선형 회귀모형으로 데이터를 설명할 수 있는 부분이 약 80%로 높고, F-검정의 p-value도 2.2e-16로 아주 작은 것도 모형이 적합하다는 것을 지지하다.  

 

 

 

[참고] 단계적선택법(stepwise selection)의 AIC는 1291.93이다. 후진소거법과 전친선택법은??

후진소거법(backward elimination)의 AIC는 1291.93

> fit.step.back = step(fit1,direction='backward')

Start:  AIC=1295.03

medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +

    tax + ptratio + black + lstat

 

          Df Sum of Sq    RSS    AIC

- chas     1      0.63 6321.5 1293.1

- indus    1     11.30 6332.2 1293.9

<none>                 6320.9 1295.0

- age      1     82.67 6403.5 1299.4

- tax      1    126.28 6447.1 1302.7

- nox      1    189.88 6510.7 1307.5

- black    1    192.42 6513.3 1307.7

- zn       1    203.44 6524.3 1308.5

- crim     1    228.82 6549.7 1310.5

- rad      8    721.85 7042.7 1332.0

- ptratio  1    706.41 7027.3 1344.9

- dis      1    860.51 7181.4 1355.6

- lstat    1    965.26 7286.1 1362.7

- rm       1   1330.92 7651.8 1386.7

 

Step:  AIC=1293.08

medv ~ crim + zn + indus + nox + rm + age + dis + rad + tax +

    ptratio + black + lstat

 

          Df Sum of Sq    RSS    AIC

- indus    1     11.00 6332.5 1291.9

<none>                 6321.5 1293.1

- age      1     82.48 6404.0 1297.4

- tax      1    130.45 6451.9 1301.1

- nox      1    189.27 6510.8 1305.5

- black    1    193.59 6515.1 1305.9

- zn       1    203.76 6525.2 1306.6

- crim     1    230.58 6552.1 1308.6

- rad      8    738.26 7059.8 1331.2

- ptratio  1    719.40 7040.9 1343.9

- dis      1    861.64 7183.1 1353.7

- lstat    1    965.11 7286.6 1360.7

- rm       1   1333.37 7654.9 1384.9

 

Step:  AIC=1291.93

medv ~ crim + zn + nox + rm + age + dis + rad + tax + ptratio +

    black + lstat

 

          Df Sum of Sq    RSS    AIC

<none>                 6332.5 1291.9

- age      1     81.09 6413.6 1296.2

- tax      1    192.78 6525.3 1304.6

- black    1    196.55 6529.0 1304.9

- zn       1    220.63 6553.1 1306.7

- crim     1    225.50 6558.0 1307.1

- nox      1    239.09 6571.6 1308.1

- rad      8    791.09 7123.6 1333.6

- ptratio  1    732.81 7065.3 1343.6

- dis      1    857.27 7189.8 1352.1

- lstat    1    987.73 7320.2 1361.0

- rm       1   1380.21 7712.7 1386.5

 

> summary(fit.step.back )

 

Call:

lm(formula = medv ~ crim + zn + nox + rm + age + dis + rad +

    tax + ptratio + black + lstat, data = boston)

 

Residuals:

    Min      1Q  Median      3Q     Max

-9.5200 -2.2850 -0.4688  1.7535 15.3972

 

Coefficients:

              Estimate Std. Error t value Pr(>|t|)   

(Intercept)  30.252522   4.329907   6.987 9.64e-12 ***

crim         -0.104568   0.025533  -4.095 4.96e-05 ***

zn            0.045510   0.011235   4.051 5.97e-05 ***

nox         -12.366882   2.932651  -4.217 2.97e-05 ***

rm            3.583130   0.353644  10.132  < 2e-16 ***

age          -0.025822   0.010514  -2.456 0.014412 * 

dis          -1.253903   0.157029  -7.985 1.08e-14 ***

rad2          2.387130   1.160735   2.057 0.040278 * 

rad3          4.644091   1.062157   4.372 1.51e-05 ***

rad4          2.608777   0.944668   2.762 0.005977 **

rad5          3.116933   0.960550   3.245 0.001258 **

rad6          1.422890   1.150280   1.237 0.216705   

rad7          4.868388   1.240114   3.926 9.94e-05 ***

rad8          5.872144   1.180865   4.973 9.26e-07 ***

rad9          6.420553   1.393304   4.608 5.24e-06 ***

tax          -0.010571   0.002792  -3.787 0.000172 ***

ptratio      -0.837356   0.113420  -7.383 7.08e-13 ***

black         0.007949   0.002079   3.823 0.000149 ***

lstat        -0.357576   0.041718  -8.571  < 2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 3.667 on 471 degrees of freedom

Multiple R-squared:  0.7907,    Adjusted R-squared:  0.7827

F-statistic: 98.83 on 18 and 471 DF,  p-value: < 2.2e-16

 

 

[참고] 전진선택법(forward selection)의 AIC는 1295.03

> fit.step.forward = step(fit1,direction='forward')

Start:  AIC=1295.03

medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +

tax + ptratio + black + lstat

 

> summary(fit.step.forward)

 

Call:

lm(formula = medv ~ crim + zn + indus + chas + nox + rm + age +

    dis + rad + tax + ptratio + black + lstat, data = boston)

 

Residuals:

    Min      1Q  Median      3Q     Max

-9.5220 -2.2592 -0.4275  1.6778 15.2894

 

Coefficients:

              Estimate Std. Error t value Pr(>|t|)   

(Intercept)  30.120918   4.338656   6.942 1.29e-11 ***

crim         -0.105648   0.025640  -4.120 4.47e-05 ***

zn            0.044104   0.011352   3.885 0.000117 ***

indus        -0.046743   0.051044  -0.916 0.360274   

chas2         0.158802   0.736742   0.216 0.829435   

nox         -11.576589   3.084187  -3.754 0.000196 ***

rm            3.543733   0.356605   9.937  < 2e-16 ***

age          -0.026082   0.010531  -2.477 0.013613 * 

dis          -1.282095   0.160452  -7.991 1.05e-14 ***

rad2          2.548109   1.175012   2.169 0.030616 * 

rad3          4.605849   1.064492   4.327 1.85e-05 ***

rad4          2.663393   0.950747   2.801 0.005299 **

rad5          3.077800   0.962725   3.197 0.001483 **

rad6          1.314892   1.157689   1.136 0.256624   

rad7          4.864208   1.241760   3.917 0.000103 ***

rad8          5.772296   1.194221   4.834 1.82e-06 ***

rad9          6.195415   1.417826   4.370 1.53e-05 ***

tax          -0.009396   0.003070  -3.061 0.002333 **

ptratio      -0.828498   0.114436  -7.240 1.85e-12 ***

black         0.007875   0.002084   3.779 0.000178 ***

lstat        -0.354606   0.041901  -8.463 3.36e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 3.671 on 469 degrees of freedom

Multiple R-squared:  0.7911,    Adjusted R-squared:  0.7821

F-statistic: 88.78 on 20 and 469 DF,  p-value: < 2.2e-16

 

 

 

 

 

3. 어떤 변수들이 제거 되었을까?

> fit.all = lm(medv~.,data=boston)

> fit.step = step(fit.all, direction = "both")

  Start:  AIC=1295.03

  medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad +

    tax + ptratio + black + lstat

 

  Df Sum of Sq    RSS    AIC

  - chas     1      0.63 6321.5 1293.1

  - indus    1     11.30 6332.2 1293.9

  <none>                 6320.9 1295.0

  - age      1     82.67 6403.5 1299.4

  - tax      1    126.28 6447.1 1302.7

  - nox      1    189.88 6510.7 1307.5

  - black    1    192.42 6513.3 1307.7

  - zn       1    203.44 6524.3 1308.5

  - crim     1    228.82 6549.7 1310.5

  - rad      8    721.85 7042.7 1332.0

  - ptratio  1    706.41 7027.3 1344.9

  - dis      1    860.51 7181.4 1355.6

  - lstat    1    965.26 7286.1 1362.7

  - rm       1   1330.92 7651.8 1386.7

 

  Step:  AIC=1293.08

  medv ~ crim + zn + indus + nox + rm + age + dis + rad + tax +

    ptratio + black + lstat

 

  Df Sum of Sq    RSS    AIC

  - indus    1     11.00 6332.5 1291.9

  <none>                 6321.5 1293.1

  + chas     1      0.63 6320.9 1295.0

  - age      1     82.48 6404.0 1297.4

  - tax      1    130.45 6451.9 1301.1

  - nox      1    189.27 6510.8 1305.5

  - black    1    193.59 6515.1 1305.9

  - zn       1    203.76 6525.2 1306.6

  - crim     1    230.58 6552.1 1308.6

  - rad      8    738.26 7059.8 1331.2

  - ptratio  1    719.40 7040.9 1343.9

  - dis      1    861.64 7183.1 1353.7

  - lstat    1    965.11 7286.6 1360.7

  - rm       1   1333.37 7654.9 1384.9

 

  Step:  AIC=1291.93

  medv ~ crim + zn + nox + rm + age + dis + rad + tax + ptratio +

    black + lstat

 

  Df Sum of Sq    RSS    AIC

  <none>                 6332.5 1291.9

  + indus    1     11.00 6321.5 1293.1

  + chas     1      0.32 6332.2 1293.9

  - age      1     81.09 6413.6 1296.2

  - tax      1    192.78 6525.3 1304.6

  - black    1    196.55 6529.0 1304.9

  - zn       1    220.63 6553.1 1306.7

  - crim     1    225.50 6558.0 1307.1

  - nox      1    239.09 6571.6 1308.1

  - rad      8    791.09 7123.6 1333.6

  - ptratio  1    732.81 7065.3 1343.6

  - dis      1    857.27 7189.8 1352.1

  - lstat    1    987.73 7320.2 1361.0

  - rm       1   1380.21 7712.7 1386.5

 

> names(fit.step)
 [1] "coefficients"  "residuals"     "effects"       "rank"        
 [5] "fitted.values" "assign"        "qr"            "df.residual" 
 [9] "contrasts"     "xlevels"       "call"          "terms"       
[13] "model"         "anova" 

 

 > fit.step$anova #최종모형에서 제거된 변수를 알 수 있다.

  Step Df   Deviance Resid. Df Resid. Dev      AIC

  1         NA         NA       469   6320.865 1295.031

  2  - chas  1  0.6261633       470   6321.491 1293.079

3 - indus  1 10.9964825       471   6332.487 1291.931

->  해석: fit.step$anova라는 함수를 통해 최종 모형에서 제거된 변수를 알 수 있다. 여러 후보 모형 중에서 AIC가 가장 작은 모형을 선택하게 되는데, 여기서는 chas와 indus가 제거되었음을 일목요연하게 알 수 있다.


 

4. 목표 예측값을 알아보자

> yhat=predict(fit.step,newdata=boston,type='response') #목표값 예측 시, type='response'

> head(yhat) #예측된 값 산출

  1        2        3        4        5        6

  26.59831 24.00195 28.99396 29.60018 29.07676 26.41636

> plot(fit.step$fitted,boston$medv, xlim=c(0,50),ylim=c(0,50),xlab="Fitted",ylab="Observed")#실제값과 가까운지 평가

> abline(a=0,b=1) # or abline(0,1)

> mean((boston$medv-yhat)^2) #MSE

  [1] 12.92344

         ->   함수 predict 는 다양한 모형 적합결과로부터 예측값을 계산할 때 사용하고, 이중 type 

           옵션은 예측 형태를 입력하는 것으로, 목표값을 예측할 때 ‘response’를 사용한다.


->   목표변수가 연속형인 경우에 모형의 예측력 측도로서 MSE(mean squared error)를 주로 사용한다. 관측치(yi)와 예측치 Ŷi의 차이가 적을수록 그 모형의 예측력은 높다고 할 수 있다. 이를 시각적으로 확인하기 위해서는 이들을 가로축 및 세로축에 놓고 그린 산점도가 45도 대각선을 중심으로 모여 있으면 예측력이 좋다고 할 수 있다. 

 

출처: 데이터마이닝(장영재, 김현중, 조형준 공저)

반응형

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제3장 나무모형 - 회귀나무모형  (0) 2016.10.26
제3장 나무모형 - 분류나무모형  (4) 2016.10.18
제2장 회귀모형 - 로지스틱 회귀모형 연습  (0) 2016.09.14
제2장 회귀모형 - 선형회귀, 로지스틱회귀  (0) 2016.09.14
1장 데이터과학과 데이터마이닝  (0) 2016.08.26
Posted by 마르띤
,

제2장 회귀모형 - 선형회귀, 로지스틱회귀

KNOU/2 데이터마이닝 2016. 9. 14. 10:05
반응형

목표변수가 연속형인 경우 -> 선형 회귀모델, ex) 광고비 투입 대비 매출액

목표변수가 두 개의 범주를 가진 이항형인 경우 -> 로지스틱 회귀모형, ex) 좋다1, 나쁘다0


 

2.1 선형회귀모형(linear regression model)

① 모형의 정의


Y= β0 + β1X1i + β2X2i +…..+ βpXpi + εi,   i = 1,…, n


- β0, β1, β2, …. Βp 를 회귀 모수(regression parameters) 또는 회귀 계수(regression coefficients)로서 알려지지 않는 상수이다.

-  εi는 Yi의 근사에서 오차(error)

 

② 회귀 모수의 추정

오차 εi= Y - β0 - β1X1i - β2X2i -…..- βpXpi


의 제곱합을 최소화하는 추정값을 이용한 최소제곱회귀직선(least square regression line)은 다음과 같다.



 

③ 회귀 계수의 해석

회귀계수 βj 다른 입력 변수들이 일정할 때 j번째 입력변수가 한 단위 변동할 때 대응하는 Y의 변동 양으로 해석. 즉, βj는 다른 입력변수를 보정한 후에 Y에 대한 Xj의 기여도. 회귀계수 βj가 양수이면 Xj가 증가할 때 Y도 증가하고, 반대로 βj가 음수이면, Xj가 증가할 때 Y는 감소함을 의미한다.

 


④ 입력 변수의 중요도


Tj =

Beta hatj

SE(Beta hatj)


 


t값의 절대 값이 클수록 영향력이 크다고 할 수 있다. P-값이 유의수준 (보통 0.05)보다 작을 때


귀무가설 H0 : β1 = 0

대립가설 H1 : β1 ≠ 0

 

귀무가설을 기각하여, Xj의 영향력이 있다고 말할 수 있다.

 

⑤ 모형의 적합도 – F값

모형의 상수항 β0을 제외한 모든 회귀계수가 0인지 아닌지를 검정하는 측도를 F-값이라 한다.

F-값은 회귀직선에 의해 평균적으로 설명할 수 있는 부분(MSR: mean squared regression)을 설명할 수 없는 부분(MSE: meas squared erroe)으로 나눈 값.

F =

MSR

=

SSR/p

MSE

SSE/(n-p-1)

 

F-값이 크면 p개 입력변수 중에 최소한 하나는 유의하다(회귀계수가 0이 아니다)라는 뜻이고, F-값이 작아서 p-값(보통 0.05) 보다 크면 모든 입력변수가 유의하지 않아서 회귀선이 쓸모가 없다.

 

⑤ 모형의 적합도 – 결정계수 R2

모형의 적합도(goodness-of-fit)를 결정계수(coefficient of determination) R2으로 측정할 수 있다.

결정계수 R2는 설명할 수 있는 부분의 총합을 변동의 총합으로 나눈 값으로 0과 1사의 값을 지닌다.

R2=

SSR

= 1 -

SSE

SST

SST

 

R2이 1에 가까울수록 모형이 데이터에 더 잘 적합되었다(fitted), 또는 회귀직선이 설명력이 높다라고 말할 수 있다.

 

다만, 모형에 포함된 변수의 수(p)가 증가하면 할수록 R2은 증가하므로 변수의 수가 다른 모형을 비교할 때는 수정된(adjusted) R2를 사용

Ra2=

Adjusted R2

= 1 -

n-1

n-p-1

 

Ra2은 변수의 수가 증가한다고 항상 증가하지는 않는다.

 

⑤ 모형의 적합도 – AIC

입력변수의 수가 다른 모형을 비교 평가하는 기준으로 AIC(Akaike information criterion)을 사용한다.


AIC = nlog(SSE/n) + 2p

 

SSE는 오차제곱합으로 작을수록 모형이 적합이 잘 되었다고 할 수 있다. 입력변수의 수가 증가할수록 SSE는 감소하지만, 벌점 2p를 더한 AIC는 항상 감소하지는 않는다. 여러 후보 모형 중에서 AIC가 가장 작은 모형을 선택한다.

 

⑥ 모형을 이용한 예측

주어진 데이터에 기반하여 회귀식 Ŷ을 얻었다고 하자. 임의의 객체 i*에 대해 관측한 입력변수의 값 x1i*, x2i*, …, xpi*를 그 회귀식에 대입하여 목표변수의 예측치 Ŷix를 얻을 수 있다.


 

⑦ 예측력

목표변수가 연속형인 경우에 모형의 예측력 측도로서 MSE(mean squared error)를 주로 사용. 시각적으로 관측치(yi)와 예측치 Ŷi의 차이를 확인하기 위해서는 이들을 가로축 및 세로축에 놓고 그린 산점도가 45도 대각선을 중심으로 모여 있으면 예측력이 좋다고 할 수 있다.


 

2.2 로지스틱 회귀모형(logistic regression model)

목표변수가 두 개의 범주를 가진 이항형인 경우, 가령 목표변수의 두 범주 값 ‘신용이 좋다’는 1, ‘신용이 나쁘다’는 0인 경우.

① 모형의 정의

총 n의 객체(subject) 중에 i번째 객체에 대한 두 개의 범주(성공 또는 실패)를 가지는 이항형 목표변수 값을 Yi, 입력변수들의 값을 X1i, X2i, …, Xpi라고 하자. 이항형 목표변수는 이항분포(binomial distribution)를 따른다. 두 개의 범주값을 1과 0으로 표시하고, 목표변수가 ‘성공’ 1을 가질 확률을 πi = pr( Yi = 1)이라고 하자. 로지스틱 회귀모형을 다음과 같이 나타낸다.

 

πi =

exp(0 + β1X1i + β2X2i +…..+ βpXpi)

, i = 1,..,n

1 + exp(0 + β1X1i + β2X2i +…..+ βpXpi)

 

단, X1i, X2i, …, Xpi 입력변수의 값이고, 이항형 목표변수는 이항분포(binomial distribution)을 따른다고 가정. 위의 식에서 목표변수가 범주형에서 연속형 변수로 바뀌지만, πi가 0~1사이의 값만을 가짐. πi가 1에 가까워질수록 입력변수의 값은 ∞로 증가, 0에 가까워질수록 0으로 수렴. 따라서 로지스틱 회귀모형은 다음과 같이 변환하여 표시할 수 있음.

log(

πi

) =

exp(0 + β1X1i + β2X2i +…..+ βpXpi)

, i = 1,..,n

1- πi

 

- β0, β1, β2, …. Βp 를 회귀 모수(regression parameters) 또는 회귀 계수(regression coefficients)로서 알려지지 않는 상수이다.


성공확률 πi와 입력변수 관계는 로지스틱 반응 함수로 표현할 수 있다. 입력변수가 증가함에 따라 초기에는 천천히 증가하다가 증가속도가 점차 빨라지고 확률 1/2 이후에는 다시 증가속도가 줄어드는 성장곡선 (growth curve) 형태이다(좌측 도형). 성공 확률과 실패 확률의 비를 오즈비(odds ratio)라고 하고 오즈비에 로그(log)를 취한 것을 로짓변화(logit transformation)이라고 부른다. 입력변수와 로짓의 관계는 직선이고, πi가 0~1의 값만 취하는 반면, 로짓변화는 -4, 6 등 다양한 값을 가진다.

 

② 회귀 모수의 추정

회귀모수는 최대우도추정법(MLE: maximum likelihood estimation method)에 이해 추정. 데이터의 확률함수를 모수β의 함수로 취급한 것을 우도함수(likelihood function) L(β)라고 하고, 이 우도함수가 최대가 될 때 모수의 추정치를 최대우도추정치(MLE)라고 한다. 적합된 로지스틱 회귀식(logistic regression line)은 다음과 같다.



 

③ 회귀계수의 해석

회귀계수 βj는 다른 입력변수들을 보정한 후 성공(Y=1)의 로그오즈(log odds = log(π/(1-π))에 미치는 Xj의 효과. 다른 입력 변수가 일정할 때, exp(βj)는 j번째 입력변수 Xj가 한 단위 변동할 때 오즈에 미치는 기여도. 회귀계수 βj가 양수이면 Xj가 증가할 때 성공확률 π와 로짓 log(π/(1-π)는 증가하고, 반대로 βj가 음수이면 Xj가 증가할 때 이들은 감소한다.

 

④ 변수의 중요도

선형회귀모형에서와 유사하게 로지스틱 회귀모형에서는 변수의 중요도는 z값으로 측정할 수 있다.



⑤ 모형의 적합도

모형의 적합도의 측도로서 이탈도(deviance)를 사용할 수 있다. 이탈도란 어떤 모형 M의 최대로그 우도(maximized log-likelihood) log(LM)에서 포화모형(saturated model) S의 최대로그우도 log(LS)를 뺀 것에 -2를 곱한 값이다.


이탈도 = -2[log(LM)-log(LS)]

 

포화모형은 각 관측에 모수 하나씩을 사용하여 완벽한 모형을 의미하며, 이탈도가 클 경우에는 그 모형은 적합하지 않다고 한다. 데이터를 모형에 적합하여 얻은 이탈도에 대응하는 p-값(보통 > 0.05)이 클 때 우리는 그 모형 M이 의미있다고 한다. 입력변수의 수가 다른 모형을 비교 평가하는 기준으로 AIC(Akaike Information Criterion)를 종종 사용한다. LM은 모형 M에 대한 우도함수의 최대값, p는 모수의 수이다.

AIC = -2log(LM) + 2P

 

AIC는 입력변수 또는 모수의 수가 증가한다고 항상 작아지지는 않으므로, 여러 후보 모형들 중에서 가장 작은 AIC를 가지는 모형을 선택한다.

 

⑥ 모형을 이용한 예측

임의의 객체 i*에 대해 관측한 입력변수의 값 x1i*,x2i*,…,xpi*를 그 로지스틱 회귀식에 대입하여 성공확률 πi* = Pr(Yi* =1)의 예측치를 얻을 수 있다.


예측치가 크면 1, 작으면 0으로 분류한다. 크고 작음을 분류하는 임계치(π0)는 보통 0.5~0.7을 사용하지만, 적영 분야에 따라 달리 결정할 수 있다. 기존 고객 데이터로 회귀계수를 추정하여 로지스틱 회귀식을 얻은 후, 수집한 새로운 고객의 입력변수값을 로지스틱 회귀식에 대입하여 새로운 고객의 성공확률을 예측한다.

 

⑦ 예측력 – 정오분류표

목표변수가 이항형(신용도가 좋으면 1, 나쁘면 0)인 경우에, 정오분류표를 만들어 예측력을 평가


 

예측범주 Ŷ

합계

1

0

실제범주 Y

1

n11

n10

n1+

0

n01

n00

n0+

합계

n+1

n+0

n

 

- 민감도(Sensitivity): Pr(Ŷ = 1 | Y =1) = n11/n1+

- 특이도(Specificity): Pr(Ŷ = 0 | Y =0) = n00/n0+

- 예측정확도(Prediction Accuracy): Pr(Ŷ =1 | Y =1) + Pr(Ŷ = 0 | Y =0) = (n11+n00)/n

- 오분류율(Misclassification Rate): Pr(Ŷ≠1 | Y =1) + Pr(Ŷ≠0 | Y =0) = (n10+n01)/n

 

민감도: 실제 양성(Y=1)일 때, 양성으로 예측할 확률(Ŷ = 1)

특이도: 실제 음성(Y=0)일 때 음성으로 예측할 확률(Ŷ = 0)

예측정확도: 실제 양성인데 양성으로, 음성일 때 음성으로 제대로 예측할 확률로 민감도와 특이도의 가중평균

오분류율: 양성인데 음성으로, 음성일 때 양성으로 잘못 예측할 확률

 


⑦ 예측력 – ROC 곡선

여러 가능한 임계치에 대해 ‘1-특이도(Specificity)’를 가로축에, 민감도를 세로축에 놓고 그린 그래프를 ROC(Receiver Operating Characteristic)곡선이라고 한다. 민감도와 특이도가 높을수록 예측력이 좋다고 할 수 있기 때문에 ROC 곡선이 좌상단에 가까울수록 ROC 곡선 아래 면적인 AUC(Area Under the ROC curve)가 커지고, AUC가 커질수록 예측력이 좋다고 할 수 있다.

2.3 범주형 입력변수 처리

입력변수가 범주형일 경우에는 가변수(dummy variable)로 변환하여 처리한다. 예를 들어, 어떤 입력변수 X가 3개의 범주(a,b,c)를 가진다고 하자. 그러면 두 개의 가변수를 다음과 같이 새롭게 정의한다.

 

X’= 1, X=a  or  0, X≠a

X”= 1, X=b  or  0, X≠b

 

X가 범주 a를 가지는 경우, X’=1, X”=0, X가 범주 b를 가지는 경우 X’=0, X”=1, X가 범주 c를 가지는 경우 X’=0, X”=0. 따라서 X가 범주를 L개 가지는 경우 L-1개의 가변수를 새롭게 생성한다.

 


2.4 모형 구축을 위한 변수 선택

① 후진소거법(backward elimination): 모든 변수를 포함시킨 모형부터 시작하여 가장 유의하지 않은 변수를 하나씩 제거

② 전진선택법(forward selection): 상수항만 가진 모형부터 시작하여 가장 유의한 변수를 하나씩 포함시켜 포함되지 않고 남은 변수가 모두 유의하지 않을 때까지 추가

③ 단계적 선택법(stepwise selection): 전진선택법처럼 상수항부터 시작하여 가장 유의한 변수를 하나씩 모형에 포함시킨다. 하지만 어떤 변수가 포함된 이후에 기존에 포함된 변수 중에 유의하지 않은 변수를 제거. 즉 전진선택법 + 후진소거법

반응형

'KNOU > 2 데이터마이닝' 카테고리의 다른 글

제3장 나무모형 - 회귀나무모형  (0) 2016.10.26
제3장 나무모형 - 분류나무모형  (4) 2016.10.18
제2장 회귀모형 - 로지스틱 회귀모형 연습  (0) 2016.09.14
제2장 회귀모형 - 선형회귀 연습  (0) 2016.09.14
1장 데이터과학과 데이터마이닝  (0) 2016.08.26
Posted by 마르띤
,

제2장 중회귀모형 - 중회귀모형 추정, 신뢰성

KNOU/2 회귀모형 2016. 9. 14. 09:50
반응형

종속 변수의 변화를 설명하기 위하여 두 개 이상의 독립 변수가 사용되는 선형회귀모형을 중선형회귀(multiple linear regression model)라 부르며, 간단히 중회귀모형(multiple regression model)이라함.




 

- 중회귀모형에서독립변수가2개인경우


 

- 벡터 표현

 

중회귀모형의 행렬 표현은

 

 


2.2 중회귀모형의 추정

-최소제곱법

 

R활용한 행렬 연산

> setwd("c:/Rwork")

> market2<-read.table("market-2.txt",header=T)

> head(market2,3)

ID  X1   X2    Y

1  1 4.2  4.5  9.3

2  2 8.5 12.0 18.5

3  3 9.3 15.0 22.8

> X=market2[,c(2:3)]

> head(X)

X1   X2

1  4.2  4.5

2  8.5 12.0

3  9.3 15.0

4  7.5  8.5

5  6.3  7.4

6 12.2 18.5

> X=cbind(1,X)

> head(X)

1   X1   X2

1 1  4.2  4.5

2 1  8.5 12.0

3 1  9.3 15.0

4 1  7.5  8.5

5 1  6.3  7.4

6 1 12.2 18.5

> Y=market2[,4]

> X=as.matrix(X)

> head(Y)

[1]  9.3 18.5 22.8 17.7 14.6 27.9

> Y=as.matrix(Y)

> head(Y)

[,1]

[1,]  9.3

[2,] 18.5

[3,] 22.8

[4,] 17.7

[5,] 14.6

[6,] 27.9

> XTX=t(X)%*%X

> XTX

1      X1      X2

1   15.0  132.80  165.50

X1 132.8 1280.80 1610.68

X2 165.5 1610.68 2149.49

> XTXI=solve(XTX)

> XTXI

1          X1           X2

1   0.82483283 -0.09810472  0.010004922

X1 -0.09810472  0.02520628 -0.011334282

X2  0.01000492 -0.01133428  0.008188029

> XTY=t(X)%*%Y

> XTY

[,1]

1   290.40

X1 2796.89

X2 3568.95

> beta=XTXI%*%XTY

> beta=round(beta,3)

> beta

[,1]

1  0.850

X1 1.558

X2 0.427

-> 결론: 마지막 beta 함수를 통해 선형회귀식을 도출할 수 있다.




 

광고료가 1000만원 (X1=10) 이고 상점의 크기가 100 (X2=10) 인 상점의 평균 총판매액의 추정값은

0.850 + 1.558*10 + 0.427*10 = 20.7 즉 207백만원으로 추정

 

 

 

 

잔차의 성질

추정된회귀식의값 Ŷi과 관찰값Yi의차이를잔차εi = Yi – Ŷi라 함.

추정값과잔차벡터의행렬표현은 아래와 같다.

 

 

여기서 H는 햇행렬(hat matrix)라고 한다. 햇행렬은 다음을 만족하는 멱등행렬(idempotent matrix)이다.

 

H2 = HH = H

H’ = H

 

잔차의 성질

(1)잔차의 합은 0

(2) 잔차의독립변수에대한가중합은0

(3) 잔차의추정값에대한가중합도0

(4) 중회귀모형 Y = Xβ + ε에서, 잔차 εi간에는 상관관계가 일반적으로 존재함.



2.3 회귀방정식의 신뢰성

<중회귀의 분산분석표>

요인

자유도

제곱합

평균제곱

F0

회귀

K

SSR(회귀제곱합)

MSR(회귀 평균 제곱)=SSR/k

MSR/MSE

잔차

n-k-1

SSE(잔차제곱합)

MSE(잔차 평균제곱)

= SSE/ n-k-1

 

계

n-1

SST(총 제곱합)

 

 

 

 

F0 =

MSR

회귀방정식이 유의한가를 검정하기 위한 검정통계량

MSE

 

귀무가설 H0 : β1 β2 = ……. = βk = 0

대립가설 H1 : 최소한 하나의 βi ≠ 0

유의수준 a에서 만약 F0의 값이 F0>F(K,n-k-1;a)이면 귀무가설을 기각하여, 회귀방정식이 유의(significant)하다는 것을 의미

 

 

R을 활용한 회귀분석

> head(market2)

ID   X1   X2    Y

1  1  4.2  4.5  9.3

2  2  8.5 12.0 18.5

3  3  9.3 15.0 22.8

4  4  7.5  8.5 17.7

5  5  6.3  7.4 14.6

6  6 12.2 18.5 27.9

> market2.lm = lm(Y~X1+X2,data=market2)

> summary(market2.lm)

 

Call:

  lm(formula = Y ~ X1 + X2, data = market2)

 

Residuals:

  Min       1Q   Median       3Q      Max

-1.30465 -0.69561 -0.01755  0.69003  1.53127

 

Coefficients:

  Estimate Std. Error t value Pr(>|t|)   

(Intercept)  0.85041    0.84624   1.005 0.334770   

X1           1.55811    0.14793  10.532 2.04e-07 ***

X2           0.42736    0.08431   5.069 0.000276 ***

  ---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 0.9318 on 12 degrees of freedom

Multiple R-squared:  0.9799,           Adjusted R-squared:  0.9765

F-statistic: 292.5 on 2 and 12 DF,  p-value: 6.597e-11

->  해석

회귀식: Ŷ = 0.85041 + 1.55811 * X1 + 0.42736 * X2

결정계수: 0.9799, 수정 결정계수는 0.9765

F값은 292.5이고 유의확률 6.597e-11로서 귀무가설을 기각하고 적합된 중회귀모형이 이 데이터를 설명하는데 유의함.  

     

 

 

R 활용한 분산분석표

> anova(market2.lm)

Analysis of Variance Table

 

Response: Y

Df Sum Sq Mean Sq F value    Pr(>F)   

X1         1 485.57  485.57 559.283 1.955e-11 ***

X2         1  22.30   22.30  25.691 0.0002758 ***

  Residuals 12  10.42    0.87                     

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  -> 해석:

SS(X1) = 485.57

SS(X2|X1) = 22.30, 변수X1이 적합된 후, 변수X2가 추가되었을 대의 추가제곱합을 의미.

X1, X2 모두 P값이 유의미함

 

요인

자유도

제곱합

평균제곱

F0

Pr(>F)

회귀

K = 2

SSR(회귀제곱합)

= 507.87

MSR=SSR/k

= 253.94

MSR/MSE

= 292

6.597e-11

잔차

n-k-1 = 12

SSE(잔차제곱합)

= 10.42

MSE(잔차 평균제곱)

= SSE/ n-k-1

= 0.87

 

 

계

n-1 = 14

SST(총 제곱합)

= 518.29

 

 

 

 

- 중상관계수(multiple correlation coefficient)

 단순 회귀 모형에서 결정계수: 두 변수의 상관계수의 제곱과 같음

 중회귀모형의 결정계수: 반응변수 Yi와 추정값Ŷi의 상관계수의 제곱과 같음. 따라서 결정계수의 제곱근을 중산관계수(multiple correlation coefficient)라고 함. 

 

 

> names(market2.lm)

[1] "coefficients"  "residuals"     "effects"     

[4] "rank"          "fitted.values" "assign"      

[7] "qr"            "df.residual"   "xlevels"     

[10] "call"          "terms"         "model"       

> yhat=market2.lm$fitted

> cor(market2$Y,yhat)

[1] 0.9898983

> cor(market2$Y,yhat)^2

[1] 0.9798986

 

->  cor(market2$Y,yhat)^2의 값이 0.9798이 결정계수 0.9799와 동일

 

 

 

 

잔차평균제곱(Residual mean squares)

 

MES = SSE / n-k-1, MSE의 값이 작으면 작을수록 관찰값 Yi들이 추정값Ŷi과 차이가 거의 업다는 것을 의미하며, 추정된 회귀방정식을 믿을 수 있게 됨.

 

> summary(market2.lm)

 

Call:

  lm(formula = Y ~ X1 + X2, data = market2)

 

Residuals:

  Min       1Q   Median       3Q      Max

-1.30465 -0.69561 -0.01755  0.69003  1.53127

 

Coefficients:

  Estimate Std. Error t value Pr(>|t|)   

(Intercept)  0.85041    0.84624   1.005 0.334770   

X1           1.55811    0.14793  10.532 2.04e-07 ***

X2           0.42736    0.08431   5.069 0.000276 ***

  ---

  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 

Residual standard error: 0.9318 on 12 degrees of freedom

Multiple R-squared:  0.9799,           Adjusted R-squared:  0.9765

F-statistic: 292.5 on 2 and 12 DF,  p-value: 6.597e-11

 

> anova(market2.lm)

Analysis of Variance Table

 

Response: Y

Df Sum Sq Mean Sq F value    Pr(>F)   

X1         1 485.57  485.57 559.283 1.955e-11 ***

X2         1  22.30   22.30  25.691 0.0002758 ***

  Residuals 12  10.42    0.87                     

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

-> 잔차평균제곱근 0.9318, sqrt(10.42/12) = 0.931844

반응형

'KNOU > 2 회귀모형' 카테고리의 다른 글

제2장 중회귀모형 - 표준화된 중회귀분석, 추정과 검정, 변수 추가  (0) 2016.09.18
제1장 단순회귀모형 - 연습문제  (0) 2016.09.14
제1장 단순회귀모형 - 단손회귀의 추정과 검정  (0) 2016.09.12
제1장 단순회귀모형 - 회귀모형, 회귀선의 추정, 회귀모형의 정도  (0) 2016.08.26
Posted by 마르띤
,

티스토리툴바