데이터마이너를 꿈꾸며

2.4 표준화된 중회귀분석

표준화된 중회귀모형

Y_i^* = α₁Z_i1 + α₂Z_i2 + …+ α_kZ_ik + ε_i’

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

> 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.7016Z₁ + 0.3376Z₂

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

2.5 추정과 검정

Market2 데이터에 대하여

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

(2) H₀ : β­₁ = 0, H₀ : β­₂ = 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)보다 더 넓게 형성된다.

H₀ : β­₁ = 0, H₀ : β­₂ = 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 10⁷ 이 되므로 귀무가설 H₀ : β­₁ = 0은 기각한다. R에서 ***는 p-값이 0.001보다 작은 경우로 표시되었음을 알 수 있다. H₀ : β­₂ = 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 이고, 검정통계랑 F₀는 24.912가 된다. 이에 대한 유의확률 p-값은 3.119 X 10⁵이므로 변수4가 유의한 변수임을 알 수 있다.

> 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에 대하 추가변수 그림이 선형성이 강한 것을 볼 수 이씅며, 따라서 이 두 변수가 회귀모형에 매우 유의하다는 것을 알 수있다. 

종속 변수의 변화를 설명하기 위하여 두 개 이상의 독립 변수가 사용되는 선형회귀모형을 중선형회귀(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)이다.

H² = HH = H

H’ = H

잔차의 성질

(1)잔차의 합은 0

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

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

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

2.3 회귀방정식의 신뢰성

<중회귀의 분산분석표>

귀무가설 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값이 유의미함

- 중상관계수(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

분석사례

어떤 슈퍼마켓에서 고객이 구입하는 상품의 금액과 카운터에서 값을 치르는 데 걸리는 시간 사이에 회귀함수 관계가 있는가를 알아보기 위해 10명의 고객을 임의로 추출하여 데이터를 얻었다. R을 이용하여 회귀 모형을 적합해 보자.

1. 자료를 읽어 산점도를 그려보자.

> setwd('c:/Rwork')

> super=read.table('supermarket.txt',header=T)

> head(super)

  price time

1   6.4  1.7

2  16.1  2.7

3  42.1  4.9

4   2.1  0.3

5  30.7  3.9

6  32.1  4.1

> length(price)

[1] 10

> a<-1:10

> attach(super)

> plot(price,time,pch=a)

2. 회귀모형과 분석분석표 구하기

> super.lm=lm(time~price,data=super)

> plot(super$price,super$time)

> abline(super.lm,lty=2,col='BLUE')

> summary(super.lm)

Call:

lm(formula = time ~ price, data = super)

Residuals:

     Min       1Q   Median       3Q      Max

-0.37928 -0.32771 -0.04431  0.32231  0.56126

Coefficients:

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

(Intercept)  0.396460   0.191488    2.07   0.0722 . 

price             0.115982   0.008979   12.92 1.22e-06 ***

---

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

Residual standard error: 0.3925 on 8 degrees of freedom

Multiple R-squared:  0.9542,    Adjusted R-squared:  0.9485

F-statistic: 166.9 on 1 and 8 DF,  p-value: 1.221e-06

> anova(super.lm)

Analysis of Variance Table

Response: time

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

price        1 25.7036 25.7036  166.85 1.221e-06 ***

Residuals  8  1.2324  0.1541                     

---

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

> qf(0.95,1,8)

[1] 5.317655

> 1-pf(166.85,1,8)

[1] 1.220851e-06

데이터 해석을 해 보자

① 회귀계수 추정값은 time hat = 0.396 + 0.116 X price 

② 기울기에 대한 t-값은 12.92이고, p-값은 1.221e-06 = 1.221 X 110⁶ 로 그 값이 매우 적기 때문에 유의하다. 

③ F-값은 166.85 이고, 이에 대한 p-값은 1.221e-06으로서 적합된 회귀직선이 유의하다는 것을 알 수 있다.

F₀ = 166.85 > F 기각역 F(1,8;0.05)의 값 qf(0.95,1,8) = 5.317655 

이므로 귀무가설 H0 : β1 = 0을 기각한다. 따라서 구해진 회귀직선이 유의하다. 

유의 확률 값은 1-pf(166.85,1,8) = 1.22081 X 106 

로서, 유의확율 p-값을 이용한 검정은 유의확률 p-값 < 유의수준 α(=0.05)이면 회귀방정식이 유의하다는 것을 나타낸다.

 ④ 추정값의 표준오차(standard error of estimate) = 0.3925 = MSE의 제곱근 sqrt(0.1541)

 ⑤ 결정 계수는 0.9542 = SSR / SST = 25.7036 / (25.7036 + 1.2324), 총 변동 중에서 회귀 방정식으로 설명되는 부분이 95.42%라는 것을 나타낸다.

3. 잔차, 추정값 보기 및 잔차 그림 그리기

> names(super.lm)

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

 [5] "fitted.values" "assign"        "qr"            "df.residual" 

 [9] "xlevels"       "call"          "terms"         "model"       

> cbind(super,super.lm$resi,super.lm$fit)

   price time super.lm$resi super.lm$fit

1    6.4  1.7    0.56125840     1.138742

2   16.1  2.7    0.43623742     2.263763

3   42.1  4.9   -0.37928275     5.279283

4    2.1  0.3   -0.34002095     0.640021

5   30.7  3.9   -0.05709314     3.957093

6   32.1  4.1   -0.01946730     4.119467

7    7.2  1.2   -0.03152683     1.231527

8    3.4  0.5   -0.29079696     0.790797

9   20.8  3.3    0.49112416     2.808876

10   1.5  0.2   -0.37043203     0.570432

> plot(price,super.lm$resi,pch=19)

> abline(h=0,lty=2)

-> time은 실제 소요 시간, super.lm$fit은 추정값, super.lm$resi는 잔차. 잔차는 실제값과 추정값 차이

   abline(h=0,lty=2)은 잔차가 0인 선을 그리되, 라인 타입을 점선(lty=2)으로 그리라는 명령이다. 그림에서 보면 잔차는 0을 중심으로 일정한 범위 내에 있으므로 회귀에 대한 기본 가정을 만족한다고 할 수 있으나, X가 증가함에 따라 곡선관계를 보여주고 있다. 따라서 2차 곡선 회귀식 Ŷ = b₀ + b₁X = b₂X²을 구해보는 것도 의미가 있으리라 생각된다.

4. 추정값의 신뢰대 그리기

> p.x = data.frame(price=c(1,45))

> pc = predict(super.lm,int='c',newdata=p.x)

> pred.x=p.x$price

> plot(super$price,super$time,ylim=range(super$time,pc))

> matlines(pred.x,pc,lty=c(1,2,2),col='BLUE')

걸음만 단계. 무한 반복과 많은 연습이 필요.

단순 회귀 직선의 기본 가정 중 하나는 주어진 X와 Y의 기대값을 E(Y)라고 하면,

E(Y)= β₀ + β₁X

가 성립되는 것이다. 이러한 의미는 관심이 되고 있는 모집단에서 변수 X와 Y사이에 선형관계가 성립된다고 하는 것이다.

1. β_0, β₁의 신뢰구간

> summary(market.lm)

Call:

lm(formula = Y ~ X, data = market)

Residuals:

   Min     1Q Median     3Q    Max

-3.600 -1.502  0.813  1.128  4.617

Coefficients:

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

(Intercept)    -2.2696     3.2123  -0.707 0.499926   

X                  2.6087     0.3878   6.726 0.000149 ***

---

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

Residual standard error: 2.631 on 8 degrees of freedom

Multiple R-squared:  0.8497,    Adjusted R-squared:  0.831

F-statistic: 45.24 on 1 and 8 DF,  p-value: 0.0001487

# β₀ 의 95%신뢰구간

> q.val=qt(0.975,8)

> -2.2696-3.2123*q.val;-2.2696+3.2123*q.val

[1] -9.677177

[1] 5.137977

# β1의 95%신뢰구간

q.val=qt(0.975,8)

 2.6084-q.val*0.3878 ; 2.6084+q.val*0.3878

[1] 1.714132

[1] 3.502668 

2. 추정값의 신뢰 구간을 그려보자

> range(market$X)

[1]  4 12

> pred.frame=data.frame(X=seq(3.5,13.5,0.2))

> pc=predict(market.lm,int='c',newdata=pred.frame) #기대값의 신뢰구간

> pp=predict(market.lm,int='p',newdata=pred.frame) # 새로운 값의 신뢰구간

> head(pc,3);head(pp,3)

       fit      lwr      upr

1 6.860870 2.402421 11.31932

2 7.382609 3.084944 11.68027

3 7.904348 3.765982 12.04271

       fit         lwr      upr

1 6.860870 -0.66731914 14.38906

2 7.382609 -0.05148732 14.81670

3 7.904348  0.56119219 15.24750

> pred.X=pred.frame$X

> plot(market$X, market$Y,ylim=range(market$Y,pp))

> matlines(pred.X,pc,lty=c(1,2,2),col='BLUE')

> matlines(pred.X,pp,lty=c(1,3,3),col='RED')

-> 안쪽 파란색 점섬은 pc의 신뢰구간, 바깥쪽 빨간색 점선은 pp의 신뢰구간. Pp의 신뢰대가 pc의 신뢰대보다 더 넓은 영역을 차지하는 것을 알 수 있다.

> identify(market$X,market$Y)

> market[c(1,2,6),]

  NUMBER  X  Y

1      1  4  9

2      2  8 20

6      6 12 30

> market[c(4,7,9),]

  NUMBER X  Y

4      4 8 15

7      7 6 18

9      9 6 10

-> 1,2,6번 데이터는 pc 신뢰구간 내 있는 데이터, 4,7,9번은 pc 신뢰구간 바깥에 있는 데이터

3. β1의 검정

귀무가설 H₀ : β₁ = 0

대립가설 H1 : β₁ ≠ 0

위 가설에 대한 검정통계량은, 양측검정이므로 t(n-2; α/2)인 기각값을 구한 후, 만약 |t₀| > t(n-2; α/2)이면 귀무가설을 기각하게 된다.

> summary(market.lm)

Call:

lm(formula = Y ~ X, data = market)

Residuals:

   Min     1Q Median     3Q    Max

-3.600 -1.502  0.813  1.128  4.617

Coefficients:

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

(Intercept)    -2.2696     3.2123  -0.707 0.499926   

X                  2.6087     0.3878   6.726 0.000149 ***

---

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

Residual standard error: 2.631 on 8 degrees of freedom

Multiple R-squared:  0.8497,    Adjusted R-squared:  0.831

F-statistic: 45.24 on 1 and 8 DF,  p-value: 0.0001487

이 결과에서 기울기 β1의 추정값 b1 = 2.6087이고, t-값

t₀= 2.6087 / 0.3878 = 6.726921

임을 알 수 있다. P – 값은 다음과 같이 구해진다.

p-값 = Pr(T > |6.726|) = 2 X (1-Pr (T<6.726)) 이므로

> 2*(1-pt(6.726921,8))

[1] 0.0001485319

-> 결론: p-값 = 0.0001485319로서 매우 작으므로 귀무가설 H₀ : β₁ = 0을 기각하게 된다.

4. 추정값의 신뢰구간

Ŷ의 분산은 X의 함수로서 X=X bar일 경우 최소. X=X bar를 대칭으로 X의 값이 X bar에서 멀어질수록 커짐. 또한 n이 커져도 Var(Ŷ)이 작아짐을 알 수 있음.

문제: 광고비를 투입할수록 매출액이 커진다고 한다. 이 자료를 바탕올 회귀식을 구했더니 매출액(Ŷ) = -2.2696 + 2.6087 * 광고비가 나왔고, 이 회귀식을 이용하여 광고비 7을 투입할 경우 매출액 15.993을 구하였다. 이 때 매출액 15.993에 대하여 95%의 신뢰구간을 설정하시오.

t(n-2 ; α/2) 은 신뢰수준이 95%이므로 α/2 = 0.025이고, 자유도는 10-2 = 8이다. 해당 값을 t분포표 값에서 찾으면 2.306이다. 이를 R에서 계산하면 아래와 같다.

> anova(market.lm)

Analysis of Variance Table

Response: Y

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

X             1 313.043  313.04   45.24 0.0001487 ***

Residuals  8  55.357    6.92                     

---

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

> qt(0.975,8)

[1] 2.306004

-> MSE = 6.92, t(10-2 ; 5%/2) = 2.306

= 15.993 ± 2.306 x SQRT ( MSE x (1 + 1/(10-2) + (7-8)^2 / 46))

= 15.993 ± 2.306 x SQRT ( 55.357/8 x (1 + 1/(10-2) + (7-8)^2 / 46))

= 15.993 ± 6.49

따라서 매출액 15.993에 대하여 95%의 신뢰구간은 9.503에서 22.483로 예측할 수 있다.

회귀분석(regression analysis): 독립변수와 종속변수 간의 함수 관계를 규명하는 통계적인 분석방법

Ŷ=f(X)+ε

 - 독립변수(independent variable) 또는 설명변수(explanatory variable): 다른 변수에 영향을 주는 변수, 흔히 Y= β₀ + β₁X 공식에서 X

 - 종속변수(dependent variable) 또는 반응변수(response variable): 독립 변수에 의해 영향을 받는다는 변수, 흔히 Y= β₀ + β₁X 공식에서 Y

회귀(回歸)라는 말은 ‘다시 본디의 자리로 돌아온다’라는 뜻으로 통계 분석에 처음 사용한 사람은 영국의 우생학자 Galton. 완두콩 실험을 통해 부모콩의 무게를 X축, 자식콩의 무게를 Y축으로 산점도를 그리자, 이들의 관계식은 양이 관계이나 1보다 작아서 자식의 무게는 평균 무게로 회귀하려는 경향이 있다는 사실을 발견하고 이를 회귀(regression)으로 표현. 당시에 Galton의 연구실에서 일하던 동료 연구원 Karl Pearson이 이를 계량적으로 처음으로 분석하여 발표.

1. 데이터를 불러와서 산점도 그래프를 그리기

> market=read.table('market-1.txt',header=T)

> head(market,3)

NUMBER X  Y

1      1 4  9

2      2 8 20

3      3 9 22

> plot(market$X,market$Y,xlab='광고료',ylab='총판매액',pch=19)

> title('광고료와 판매액의 산점도')

2. 단순 회귀 분석 실시

> market.lm=lm(Y~X,data=market)

> summary(market.lm)

해석

추정값은 Coefficients: Estimate에서 확인. 추정된 회귀식은 Ŷ=-2.27 + 2.6 X

3. 산점도 위에 회귀직선을 그리자 - 회귀선의 추정

> abline(market.lm)

> identify(market$X,market$Y)

[1]  4  5 10

# Identify는 재미있는 함수인데, 본 함수를 입력하고 마우스로 점을 클릭하면 그림처럼 값을 알 수 있다.

> xbar = mean(market$X)

> ybar = mean(market$Y)

> xbar

[1] 8

> ybar

[1] 18.6

> points(xbar,ybar,pch=17,cex=2.0,col='RED')

> text(xbar,ybar,"(8,18.6)")

> fx <- "Y-hat = -2.27 + 2.6*X "

> text(locator(1),fx) #locator(1)은 마우스로 클릭하면서 지정

4. 회귀식 특징

> names(market.lm)

[1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values"

[6] "assign"        "qr"            "df.residual"   "xlevels"       "call"        

[11] "terms"         "model"       

> market.lm$resid

1          2          3          4          5          6          7          8

0.8347826  1.4000000  0.7913043 -3.6000000 -1.6000000  0.9652174  4.6173913  1.1826087

9         10

-3.3826087 -1.2086957

> resid=market.lm$residual

> sum(resid) #특징1. 잔차의합은 0이다

[1] 0

> sum(market$X*resid) #특징2. 잔차들의 Xi에 의한 가중합은 0이다

[1] 2.220446e-15

> sum(market.lm$fitted*resid) #특징3. 잔차들의 Yi에 의한 가중합은 0이다

[1] -1.24345e-14

> names(market.lm)

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

 [5] "fitted.values" "assign"        "qr"            "df.residual"  

 [9] "xlevels"       "call"          "terms"         "model"        

> sum(market.lm$fitted.values)

[1] 186

> sum(market$Y)

[1] 186

#특징4. 추정값Yhat의 값과 관찰값Yi의 값은 같다.

5. 회귀모형의 정도

 - 산점도 위에 회귀직선을 그려 회귀선의 정도를 대략 짐작할 수 있으나, 이러한 경우는 독립변수가 하나인 경우에만 유용하게 쓰일 수 있다. 추정된 회귀선의 정도를 측정하는 여러 가지 측도(measure)들 중에서 널리 이용되는 세가지를 알아보자

 ① 분산분석표에 의한 F-검정

 ② 결정계수

 ③ 추정값의 표준오차 

 ④ 상관계수와 결정계수

① 분산분석표에 의한 F-검정

- SST(Total sum of squares): 총 제곱합

- SSR(Sum of squares due to regression): 회귀제곱합, 설명되는 편차

- SSE(Sum of squares due to residual errors): 잔차제곱합, 설명되지 않는 편차

귀무가설 H₀ : β₁ = 0

대립가설 H1 : β₁ ≠ 0

F₀ > F(1, n-2 ; α )이면 귀무가설 H₀ : β₁ = 0을 기각하고, 회귀직선이 유의하다고 말한다. R분석 결과에서는 검정통계량 F0에 대한 유의확률 p값이 제공된다. 즉 p값 < 유의확률 α이면 귀무가설 H₀ : β₁ = 0을 기각한다. 

분산분석표

> anova(market.lm) #분산분석표

p값은0.0001487로 p값 < 유의확률 α 이므로 귀무가설  H₀ : β₁ = 0을 기각한다, 따라서 본 회귀식은 유의하다.

> qf(0.95,1,8)

[1] 5.317655

→ 유의수준 α =0.05에서 F-기각역 F(1,8;0.05)의 값은 5.32, " F₀ = 45.24 > 5.32"이므로 귀무가설 기각, 회귀선은 유희하다

> 1-pf(45.25,1,8)

[1] 0.0001485485

→ 유의확률 p값을 이용한 검정은 0.0001485485 이 값이 주어진 유의수준 α =0.05 보다 작을수록 귀무가설을 기각한다.

② 결정계수

R²=SSR/SST =1-SSE/SST

R² 을 결정계수(coefficient of determination)라고 부른다.

R²=SSR/SST = 313.04/368.4= 84.97%

이는 총변동 주에서 회귀직선에 의하여 설명되는 부분이 84.97%라는 의미로서, 추정된 회귀선의 정도가 높다는 것을 알 수 있다.

R을 통해서도 알 수 있는데, Multiple R-squared가 0.8497임을 알 수 있다.

> summary(market.lm)  

③ 추정값의 표준오차

선형회귀모형 Y= β₀ + β₁X + ε 을 표본의 자료로부터 적합시킬 때,

Y의 기댓값은 E(Y) = β₀ + β₁X, 분산은 σ²로 가정,

하였다. 따라서 Y의 측정값들이 회귀선 주위에 가깝게 있다면 σ의 추정값은 작아질 것이다. 

분산분석표에서 잔차평균제곱 MSE는 σ2의 불편추정량이 된다. 따라서 MSE의 제곱근을 추정값의 표준오차(standard error of estimate)라고 부르며, 다음과 같이 표현한다.

S_Y•X = SQRT(MSE) = SQRT(SSE/n-2) = 2.63

R을 통해서도 알 수 있는데, Residual standard error가 2.631임을 알 수 있다.

> summary(market.lm) 

④ 상관계수와 결정계수

상관계수는 연속인 두 변수 간의 선형관계(linear relationship)가 어느 정도인가를 재는 측도로서, 단순 회귀 분석에서는 상관계수 r을 다음과 같이 구할 수 있다.

r = ±SQRT(R² )

즉 상관계수는 결정계수 R²의 제곱근이며, 만약 추정된 회귀선의 기울기 b1가 양이면 양의 상관계수를 갖고, 기울기b1이 음이면 음의 상관계수를 갖는다. 회귀식 Ŷ=-2.27 + 2.6 X에서 b1가 2.6으로 양이므로 상관계수는 0.8497의 sqrt 값, 0.9217임을 알 수 있다.