데이터마이너를 꿈꾸며

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