4장 다중회귀분석
유병주
2023 3 28
준비할 사항
필요한 패키지 설치 : install.packages(c(‘car’, ‘psych’,
‘leaps’))
leaps(Regression Subset Selection package) 패키지 참고 : https://cran.r-project.org/web/packages/leaps/index.html
데이터파일 : supervisorperformance.rdata
참고교재 : Samprit Chatterjee, Ali S. Hadi(2012), Regression Analysis By
Example Fifth Edition, John Wiley & Sons.
http://www1.aucegypt.edu/faculty/hadi/RABE5/ 3장
1. 다중선형회귀모형
자료의 형태
\[\begin{equation} \left[ \begin{array}{c} y_1 \\ \vdots\\ y_n \end{array} \right] = \left[ \begin{array}{cccc} 1 & x_{11} & \cdots & x_{1p} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & x_{n1} & \cdots & x_{np} \end{array} \right] \left[ \begin{array}{c} \beta_0 \\ \vdots \\ \beta_p \end{array} \right] + \left[ \begin{array}{c} \varepsilon_1 \\ \vdots \\ \varepsilon_n \end{array} \right] \end{equation}\]
일반선형모형
\[\begin{equation} y = X\beta + \varepsilon ,\quad \varepsilon \sim N(0,\sigma^2 I) \end{equation}\] \[\begin{equation} → y|x \sim N(X\beta, \sigma^2 I) \end{equation}\]
다중선형회귀모형
\[\begin{equation} y = \beta_0 + \beta_1 x_1 + … + \beta_p x_p + \varepsilon \end{equation}\]
\[\begin{equation} y_i = \beta_0 + \beta_1 x_{i1} + … + \beta_p x_{ip} + \varepsilon_i ,\quad i =1, 2, …, n \end{equation}\]
\[\begin{equation} \varepsilon_i \sim N(0,\sigma^2),\quad cov(\varepsilon_i , \varepsilon_j )=0,\quad i≠j \end{equation}\]
2. 사례분석 자료
아래의 자료는 감독자의 직무수행 능력 자료이며 금융기관에서 근무하는 사무직원의 감독자에 대한 만족도 설문조사 결과이다.
supervisor = read.table('https://raw.githubusercontent.com/byungjooyoo/Dataset/main/supervisorperformance.csv',sep='', header = TRUE)
data=supervisor
str(data)## 'data.frame': 30 obs. of 7 variables:
## $ Y : int 43 63 71 61 81 43 58 71 72 67 ...
## $ X1: int 51 64 70 63 78 55 67 75 82 61 ...
## $ X2: int 30 51 68 45 56 49 42 50 72 45 ...
## $ X3: int 39 54 69 47 66 44 56 55 67 47 ...
## $ X4: int 61 63 76 54 71 54 66 70 71 62 ...
## $ X5: int 92 73 86 84 83 49 68 66 83 80 ...
## $ X6: int 45 47 48 35 47 34 35 41 31 41 ...head(data)## Y X1 X2 X3 X4 X5 X6
## 1 43 51 30 39 61 92 45
## 2 63 64 51 54 63 73 47
## 3 71 70 68 69 76 86 48
## 4 61 63 45 47 54 84 35
## 5 81 78 56 66 71 83 47
## 6 43 55 49 44 54 49 34총 30개 부서에 대한 만족도 조사 결과이고, 하나의 관찰치는 5점척도로 35명의 직원과 1명의 감독자로 구성된 부서별 설문조사 결과를 합한 접수이며 만족도 1-2는 만족, 4-5는 불만족을 의미한다. 각 항목별 조사내용은 아래와 같다.
- y : 상사의 업무수행에 대한 전반적인 평가
- x1 : 피고용인의 불만 처리
- x2 : 특권을 허용하지 않음
- x3 : 새로운 것을 배울 기회
- x4 : 업무성과에 따른 승진
- x5 : 과실에 대한 지나친 비판
- x6 : 더 나은 일로의 진급
3. 자료의 탐색
선형모형을 분석하기 위한 자료를 탐색하는 방법으로 종속변수와 독립변수간의 관계를 대략적으로 알아보기 위해 psych::pairs.panels()방법이 있다.
#install.packages("psych")
psych::pairs.panels(data)
자료를 대략적으로 살펴보면 종속변수와 x1변수가 가장 상관관계가 크고,
독립변수간에는 x1과 x2, x3, x4의 관계와, x3와 x4의 관계 등이 크다.
단순회귀모형의 경우 일대 일의 인과관계를 입증하면 되지만 다중회귀모형의
경우 다대 일의 인과관계를 설명해야하므로 설명변수간의 독립이 전제되어야
한다. 그래서 서로 교락되어 있는 효과가 있는 경우 모형에 어떻게 반영하는
것이 좋은지 탐색하고 처리(주효과와 교호작용등으로 분리)하는 절차가
필요하다.
4. 초기모형 선택
초기 모형은 모든 독립변수를 고려한 완전모형을 고려하기로 한다.
y=data[,1]
xdata=data[,-1]
n=nrow(data);p=ncol(xdata)
fit=lm(y ~ . , data=xdata)
summary(fit)##
## Call:
## lm(formula = y ~ ., data = xdata)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.9418 -4.3555 0.3158 5.5425 11.5990
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.78708 11.58926 0.931 0.361634
## X1 0.61319 0.16098 3.809 0.000903 ***
## X2 -0.07305 0.13572 -0.538 0.595594
## X3 0.32033 0.16852 1.901 0.069925 .
## X4 0.08173 0.22148 0.369 0.715480
## X5 0.03838 0.14700 0.261 0.796334
## X6 -0.21706 0.17821 -1.218 0.235577
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.068 on 23 degrees of freedom
## Multiple R-squared: 0.7326, Adjusted R-squared: 0.6628
## F-statistic: 10.5 on 6 and 23 DF, p-value: 1.24e-05분석결과 전반적인 모형은 만족스럽지만 유의한 변수는 x1만 만족하고 있으며 각 변수별로 \(\beta_j = 0\)라는 귀무가설에 대한 95% 신뢰구간 추정결과를 보면 알수 있다.
confint(fit,"X1"); confint(fit,"X2")## 2.5 % 97.5 %
## X1 0.2801687 0.9462066## 2.5 % 97.5 %
## X2 -0.3538181 0.2077178독립변수를 모형에 포함시킬 것인지 여부도 영향이 있지만 어떤 변수를 먼저 모형에 포함시킬 것인가도 다음 변수선택에 영향을 미친다.
fit6=lm(y ~ X6+X5+X4+X3+X2+X1, data=data)
anova(fit);anova(fit6)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 2927.58 2927.58 58.6026 9.056e-08 ***
## X2 1 7.52 7.52 0.1505 0.7016
## X3 1 137.25 137.25 2.7473 0.1110
## X4 1 0.94 0.94 0.0189 0.8920
## X5 1 0.56 0.56 0.0113 0.9163
## X6 1 74.11 74.11 1.4835 0.2356
## Residuals 23 1149.00 49.96
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## X6 1 103.35 103.35 2.0688 0.1638099
## X5 1 59.13 59.13 1.1836 0.2878999
## X4 1 1561.92 1561.92 31.2656 1.089e-05 ***
## X3 1 653.08 653.08 13.0729 0.0014532 **
## X2 1 45.69 45.69 0.9145 0.3488557
## X1 1 724.80 724.80 14.5086 0.0009029 ***
## Residuals 23 1149.00 49.96
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 15.모형 선택(변수 선택)
가. Testing Based Procedures
모형을 선택하는 방법은 Testing Based Procedures와 Criterion Based
Procedures가 있다.
Testing Based Procedures는 아래와 같은 방법이 있으며 이는 step()함수의
direction 옵션을 forward, backward, both로 선택하여 진행할 수
있다.
* Backward Elimination(후진제거법) : Full model에서 변수를 단계적으로
하나씩 제거하는데 이때 유의수준을 5%보다 큰 15~20%를 고려
* Forward Selection(전진선택법) : 어떠한 설명변수가 추가되지 않을 때까지
진행
* Stepwise Regression(단계적 방법) : 설명변수가 작은 모형을 선호하는
경향
#step(fit, direction = 'backward')
#step(fit, direction = 'forward')
step(fit, direction = 'both')## Start: AIC=123.36
## y ~ X1 + X2 + X3 + X4 + X5 + X6
##
## Df Sum of Sq RSS AIC
## - X5 1 3.41 1152.4 121.45
## - X4 1 6.80 1155.8 121.54
## - X2 1 14.47 1163.5 121.74
## - X6 1 74.11 1223.1 123.24
## <none> 1149.0 123.36
## - X3 1 180.50 1329.5 125.74
## - X1 1 724.80 1873.8 136.04
##
## Step: AIC=121.45
## y ~ X1 + X2 + X3 + X4 + X6
##
## Df Sum of Sq RSS AIC
## - X4 1 10.61 1163.0 119.73
## - X2 1 14.16 1166.6 119.82
## - X6 1 71.27 1223.7 121.25
## <none> 1152.4 121.45
## + X5 1 3.41 1149.0 123.36
## - X3 1 177.74 1330.1 123.75
## - X1 1 724.70 1877.1 134.09
##
## Step: AIC=119.73
## y ~ X1 + X2 + X3 + X6
##
## Df Sum of Sq RSS AIC
## - X2 1 16.10 1179.1 118.14
## - X6 1 61.60 1224.6 119.28
## <none> 1163.0 119.73
## + X4 1 10.61 1152.4 121.45
## + X5 1 7.21 1155.8 121.54
## - X3 1 197.03 1360.0 122.42
## - X1 1 1165.94 2328.9 138.56
##
## Step: AIC=118.14
## y ~ X1 + X3 + X6
##
## Df Sum of Sq RSS AIC
## - X6 1 75.54 1254.7 118.00
## <none> 1179.1 118.14
## + X2 1 16.10 1163.0 119.73
## + X4 1 12.54 1166.6 119.82
## + X5 1 7.18 1171.9 119.96
## - X3 1 186.12 1365.2 120.54
## - X1 1 1259.91 2439.0 137.94
##
## Step: AIC=118
## y ~ X1 + X3
##
## Df Sum of Sq RSS AIC
## <none> 1254.7 118.00
## + X6 1 75.54 1179.1 118.14
## - X3 1 114.73 1369.4 118.63
## + X2 1 30.03 1224.6 119.28
## + X4 1 1.19 1253.5 119.97
## + X5 1 0.00 1254.7 120.00
## - X1 1 1370.91 2625.6 138.16##
## Call:
## lm(formula = y ~ X1 + X3, data = xdata)
##
## Coefficients:
## (Intercept) X1 X3
## 9.8709 0.6435 0.2112나. Criterion Based Procedures
- 잔차의 MSE가 작거나, 또는 \(R^2\)나
\(R_a^2\)이 큰 모형
- AIC(Akaike’s Information Criterion) = \(n*log(SSE/n)+2p\)이 작은 모형
- BIC(Bayesian Information Criterion) = \(n*log(SSE/n)+p*log(n)\)이 작은 모형
- Mallow’s \(C_p =\frac{SSE}{\hat\sigma^2}+(2p+2-n)\)의 Full model 기대값 p+1 선과 근접한 것
require(leaps) ## 필요한 패키지를 로딩중입니다: leaps## Warning: 패키지 'leaps'는 R 버전 4.2.2에서 작성되었습니다modelset=regsubsets(y ~ .,xdata)
x.set=summary(modelset)
x.set$which## (Intercept) X1 X2 X3 X4 X5 X6
## 1 TRUE TRUE FALSE FALSE FALSE FALSE FALSE
## 2 TRUE TRUE FALSE TRUE FALSE FALSE FALSE
## 3 TRUE TRUE FALSE TRUE FALSE FALSE TRUE
## 4 TRUE TRUE TRUE TRUE FALSE FALSE TRUE
## 5 TRUE TRUE TRUE TRUE TRUE FALSE TRUE
## 6 TRUE TRUE TRUE TRUE TRUE TRUE TRUEx.set$rss## [1] 1369.382 1254.649 1179.109 1163.012 1152.406 1149.000plot(1:p, x.set$adjr2,ylab="Ajusted R-Square", xlab="Number of Predictors")
which.max(x.set$adjr2)## [1] 3plot(modelset, scale = "r2")
plot(modelset, scale = "adjr2")
plot(modelset, scale = "bic")
plot(modelset, scale = "Cp")
다. 완전모형(Full Model)과 축소모형(Reduced Model)의 비교
전반적인 모형에 대한 검정통계량은 \(F(p,n-p-1)=MSR/MSE\)으로 확인하며 이때의
귀무가설은 \(H_0 : \beta_1 =\beta_2 = \cdots
=\beta_p =0\)이며 대립가설은 \(H_1
:\) 최소한 하나의 \(j\)에 대해
\(\beta_j ≠ 0\)이다.
여기서 모형에 대한 SSR은 변수가 진입하면서 순차적으로 계산이 되기때문에
변수선택 순서에 따라 결과가 다를 수 있다. 기존 full model fit의 경우를
계산하면 아래와 같다.
\[\begin{equation}
SSR = SSR(x_1 )+SSR(x_2 |x_1 )+ \cdots + SSR(x_p | x_{p-1}, \cdots,x_1 )
\end{equation}\]
anova(fit)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 2927.58 2927.58 58.6026 9.056e-08 ***
## X2 1 7.52 7.52 0.1505 0.7016
## X3 1 137.25 137.25 2.7473 0.1110
## X4 1 0.94 0.94 0.0189 0.8920
## X5 1 0.56 0.56 0.0113 0.9163
## X6 1 74.11 74.11 1.4835 0.2356
## Residuals 23 1149.00 49.96
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1축소모형(RM)과 완전모형(FM)을 비교할 경우 예를 들어 \(x_1\)과 \(x_3\)만 포함하는 축소모형과 완전모형을 비교 검증하고자 한다면 축소모형의 변수 개수를 \(k\)라고 하면 그 모형과 검정통계량은 아래와 같다. \[\begin{equation} RM : y = \beta_0 + \beta_1 x_1 + \beta_3 x_3 + \varepsilon \\ H_0 : \beta_2 = \beta_4 = \beta_5 =\beta_6 =0 \\ F(p-k,n-p-1) = \frac{(SSE_{RM} -SSE_{FM})/(p-k)}{SSE_{FM}/(n-p-1)} \\ \end{equation}\]
fit13=lm(y ~ X1+X3, data=data)
anova(fit13,fit)## Analysis of Variance Table
##
## Model 1: y ~ X1 + X3
## Model 2: y ~ X1 + X2 + X3 + X4 + X5 + X6
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 27 1254.7
## 2 23 1149.0 4 105.65 0.5287 0.7158fit2=lm(y ~ X2, data=data)
anova(fit2,fit)## Analysis of Variance Table
##
## Model 1: y ~ X2
## Model 2: y ~ X1 + X2 + X3 + X4 + X5 + X6
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 3516.7
## 2 23 1149.0 5 2367.7 9.4792 5.185e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1이때 귀무가설이 기각되는 경우 다른 변수를 선택하거나 더 많은 변수를 선택해야 한다는 의미이다.
라. 최선의 모형선택 방법은?
| 변수선택방법 | 상수항 | X1 | X2 | X3 | X4 | X5 | X6 |
|---|---|---|---|---|---|---|---|
| Testing Based | O | O | O | ||||
| \(R_a^2\) | O | O | O | O | |||
| AIC | O | O | O | ||||
| BIC | O | O | |||||
| \(C_p\) | O | O | O |
이 경우 어떤 모델을 선택할 것인가?
* 너무 수치적인 것(검정통계량, p-value 등)에 집착하지 말 것
* 연구하고자 하는 방향과 목적이 더 중요
* 해당 분야의 전문가 의견이 중요
fit1=lm(y ~ X1, data=data)
anova(fit1,fit13)## Analysis of Variance Table
##
## Model 1: y ~ X1
## Model 2: y ~ X1 + X3
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 1369.4
## 2 27 1254.7 1 114.73 2.4691 0.1278- y : 상사의 업무수행에 대한 전반적인 평가
- x1 : 피고용인의 불만 처리
- x3 : 새로운 것을 배울 기회
설문조사의 경우 전반적인 부서 만족도가 여러가지 질문에 대해 비슷하게 답하였을 경우가 있을 수 있어 문항별 차별성있는 질문지 작성이 중요하게 생각된다.
// knitr::purl(“~/lmr_2023/4장_다중회귀분석.Rmd”)