데이터분석응용 in R_회귀분석 pt2 실습260401
# 선형회귀
Supervisor <- read.table("./datasets/Supervisor.txt", sep = "\t", header = TRUE)
head(Supervisor, 5)
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
reg <- lm(Y ~ ., data = Supervisor)
reg <- lm(Y ~ X1+X2+X3+X4+X5+X6, data = Supervisor)
# 둘 다 같은 코드이다.
# 컬럼이 Y 제외 X1~X6인데,
# 전부 불러 Y와의 관계를 볼거면 '.'
# 일부를 볼거면 보고 싶은 컬럼들만 +로 연결해서 적으면 된다.
?lm
# lm 사용법에 대한 도움말 부르기 기능
# Value 섹션에 가면 아래와 같은 열들이 있다.
# 이들은 lm을 실행하는 순간 자동으로 생성되어 저장되는 값들이다.
| coefficients | a named vector of coefficients |
| residuals | the residuals, that is response minus fitted values. |
| fitted.values | the fitted mean values. |
| rank | the numeric rank of the fitted linear model. |
| weights | (only for weighted fits) the specified weights. |
| df.residual | the residual degrees of freedom. |
| call | the matched call. |
| terms | the terms object used. |
| contrasts | (only where relevant) the contrasts used. |
| xlevels | (only where relevant) a record of the levels of the factors used in fitting. |
| offset | the offset used (missing if none were used). |
| y | if requested, the response used. |
| x | if requested, the model matrix used. |
| model | if requested (the default), the model frame used. |
| na.action | (where relevant) information returned by model.frame on the special handling of NAs. |
reg$coef # hat 추정 계수
(Intercept) X1 X2 X3
10.78707639 0.61318761 -0.07305014 0.32033212
X4 X5 X6
0.08173213 0.03838145 -0.21705668 # b hat의 계수를 추정한 것이다.
# 추가적인 해석을 가미해보자면,
# 하나의 변수를 나머지 변수들을 통제한 상황에서 1 증가시키면
# 그 변수의 계수만큼 값이 증가한다고도 말 할 수 있겠다.
reg$fitted # hat계수들로 추정했을때의 값
1 2 3 4 5
51.11030 61.35277 69.93944 61.22684 74.45380
6 7 8 9 10
53.94185 67.14841 70.09701 79.53099 59.19846
11 12 13 14 15
57.92572 55.40103 59.58168 70.21401 76.54933
16 17 18 19 20
84.54785 76.15013 61.39736 68.01656 55.62014
21 22 23 24 25
42.60324 63.81902 63.66400 44.62475 57.31710
26 27 28 29 30
67.84347 75.14036 56.04535 77.66053 76.87850
reg$res # y-(hat_y) : 잔차 계산하기
1 2 3 4
-8.1102953 1.6472337 1.0605589 -0.2268416
5 6 7 8
6.5462010 -10.9418499 -9.1484140 0.9029929
9 10 11 12
-7.5309862 7.8015424 6.0742817 11.5989723
13 14 15 16
9.4183197 -2.2140147 0.4506705 -3.5478519
17 18 19 20
-2.1501319 3.6026355 -3.0165587 -5.6201442
21 22 23 24
7.3967582 0.1809831 -10.6639999 -4.6247464
25 26 27 28
5.6828983 -1.8434727 2.8596385 -8.0453540
29 30
7.3394730 5.1215016
summary(reg)
Call:
lm(formula = Y ~ X1 + X2 + X3 + X4 + X5 + X6, data = Supervisor)
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
(Intercept)
X1 ***
X2
X3 .
X4
X5
X6
---
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
# 여기서 유심히 봐야하는 지표들은 아래와 같다.
# 1. coefficients 오른쪽 열의 Pr값 & 해당 부분 밑의 ***이 체크되는 부분
# 2. Adjusted R-squared값
# 3. F-statistic의 p-value값
# 1. coefficients 오른쪽 열의 Pr값 & 해당 부분 밑의 ***이 체크되는 부분
# 첫 행이 ***이다. 즉 X1의 회귀계수 유의수준이 0.001이상이면 기존 귀무가설 기각한다는 뜻이다.
# 즉, X1의 추정된 회귀계수 0.61319는 유의미하다.
# 세번째 X3의 회귀계수는 그럼 유의수준이 0.1이상일때 유의미
# 아참. 이들의 귀무가설은 각각의 회귀계수가 0인것이다.
# 따라서 유의수준을 0.05라고한다면, X1 제외하고는 모든 회귀계수들을 0이다가 결론이며,
# X1의 회귀계수는 0이 아니다. 즉, X1은 종속변수에 영향을 준다. 나머지 변수들은 영향을 안준다.
# 2. Adjusted R-squared값
# 0.6628이니 분명한 해석 기준선이 있지 않지만, 이정도면~ 설명력이 있다.
# 3. F-statistic의 p-value값
# 1.24e-05는 0.0000124. 즉 대부분의 유의수준이 0.05이니 이정도면~ 유의하다.
# 귀무가설은 H0는 모형이 적합하지 않다. H1은 모형이 적합하다.
# 위에서 봤던 p-value는 한개 변수의 회귀계수의 유의미한 정도를 본다면,
# f-statistic의 p-value는 전체 모형의 유의미한 정도를 본다.
# 이 p-value값의 귀무가설은 설명변수 전부 의미없다이며,
# 이 귀무가설의 대립가설은 최소 하나의 변수는 유의미하다.(그게 뭔지는 모르지만!)
confint(reg, level=0.95) #각각의 hat들의 신뢰구간 첫과 끝
2.5 % 97.5 %
(Intercept) -13.18712881 34.7612816
X1 0.28016866 0.9462066
X2 -0.35381806 0.2077178
X3 -0.02827872 0.6689430
X4 -0.37642935 0.5398936
X5 -0.26570179 0.3424647
X6 -0.58571106 0.1515977# 신뢰구간을 봐도 유의미한 변수를 예측할수 있다.
# 예를 들어, X1의 신뢰구간안에 0이 없어서 h1은 0이 아니다.
# X2의 신뢰구간안에는 0이 있어서 가능성이 있다.
# ols_regress
# lm은 모델을 만드는 함수라면,
# ols_regress는 결과를 자세히 분석하는 함수이다.
# lm 결과를 기반으로 ANOVA, F-test, R^2 등 상세 분석을 출력한다.
# summary와 동일한 기능을 수행하지만,
# summary와는 다르게(?) 조금 더 친절하다. 더 보기 좋게 정리했다.
if(!require("olsrr")) install.packages("olsrr") # 필요한 패키지를 설치안했다면 설치하기
library(olsrr)
ols <- ols_regress(Y ~ ., data = Supervisor)
print(ols)
Model Summary
---------------------------------------------------------------
R 0.856 RMSE 6.189
R-Squared 0.733 MSE 38.300
Adj. R-Squared 0.663 Coef. Var 10.936
Pred R-Squared 0.547 AIC 210.500
MAE 5.179 SBC 221.709
---------------------------------------------------------------
RMSE: Root Mean Square Error
MSE: Mean Square Error
MAE: Mean Absolute Error
AIC: Akaike Information Criteria
SBC: Schwarz Bayesian Criteria
ANOVA
--------------------------------------------------------------------
Sum of
Squares DF Mean Square F Sig.
--------------------------------------------------------------------
Regression 3147.966 6 524.661 10.502 0.0000
Residual 1149.000 23 49.957
Total 4296.967 29
--------------------------------------------------------------------
Parameter Estimates
-----------------------------------------------------------------------------------------
model Beta Std. Error Std. Beta t Sig lower upper
-----------------------------------------------------------------------------------------
(Intercept) 10.787 11.589 0.931 0.362 -13.187 34.761
X1 0.613 0.161 0.671 3.809 0.001 0.280 0.946
X2 -0.073 0.136 -0.073 -0.538 0.596 -0.354 0.208
X3 0.320 0.169 0.309 1.901 0.070 -0.028 0.669
X4 0.082 0.221 0.070 0.369 0.715 -0.376 0.540
X5 0.038 0.147 0.031 0.261 0.796 -0.266 0.342
X6 -0.217 0.178 -0.183 -1.218 0.236 -0.586 0.152
-----------------------------------------------------------------------------------------# 현재 우리가 반드시 챙겨가야할 값은 회귀계수 / 표준오차 / 결정계수이다.
# 화면상으로는 Beta, Std. Error, R-Squared가 되겠다.
# 이들만 따로 보는 방법도 있다. 각각 회귀계수 / 표준오차 / 결정계수이다.
ols$betas
(Intercept) X1 X2 X3 X4
10.78707639 0.61318761 -0.07305014 0.32033212 0.08173213
X5 X6
0.03838145 -0.21705668
ols$std_errors
(Intercept) X1 X2 X3 X4
11.5892572 0.1609831 0.1357247 0.1685203 0.2214777
X5 X6
0.1469954 0.1782095
ols$rsq
[1] 0.732602
# anova_ 결과 분석 함수
# 위의 결과에서 유의수준 0.1이상에서는 1,3 계수만 유의미하다는 결론이 나왔다.
# 그럼 전체 변수를 사용한 모형과 두개의 변수만을 사용한 모형의 적절성을 평가해보자.
# anova 모형 비교는 분산분석과는 다르다.
# aov()는 분산분석 모델을 "만드는 함수"라면
# anova()는 모델을 "분석(ANOVA 테이블 출력)"하는 함수이다.
reg1 <- lm(Y ~ X1+X2+X3+X4+X5+X6, data = Supervisor)
reg0 <- lm(Y ~ X1+X3, data = Supervisor)
anova(reg0, reg1)
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.7158# 봐야할 지표는 F Pr값이다. 슬슬 예상이 가지 않는가?
# pr, 즉 p-value값은 0.7158로 일반 유의수준 0.05보다 한참 크다.
# 귀무가설 채택. 귀무가설은 뭐지? b2=b4=b5=b6=0이다.
# 대립가설은 b2, b4, b5, b6 중 어느 하나라도 0이 아니다.
# 결론은 귀무가설 채택이므로 b2=b4=b5=b6=0이다.
# 즉 b2,b4,b5,b6 변수들은 종속변수에 영향을 주지 않는다.
# 잠깐 복습
Hamilton <- read.table("./datasets/Hamilton.Data.txt", sep = "\t", header = TRUE)
head(Hamilton, 5)
Y X1 X2
1 12.37 2.23 9.66
2 12.66 2.57 8.94
3 12.00 3.87 4.40
4 11.93 3.10 6.64
5 11.06 3.39 4.91
# 복습1겸
# regression 모델을 lm을 이용해 만들고
# 모델의 y 추정값, 잔차, 개별 데이터의 영향력, 표준화잔차, 보정된 표준화 잔차를
# 기존 Hamilton 데이터셋에 합쳐서 만들어보자.
# 참고
# 회귀계수(β)와 "개별 데이터의 영향력" (hatvalues) 차이
| 대상 | 변수 (X1, X2) | 데이터 한 개 (i번째 관측치) |
| 의미 | 영향 크기 | 영향력 |
| 개수 | 변수 개수 | 데이터 개수 |
reg <- lm(Y ~ X1+X2, data = Hamilton)
Y.hat <- fitted(reg)
e <- residuals(reg)
h <- hatvalues(reg)
r <- rstandard(reg) # 표준화된 잔차
rs <- rstudent(reg) # (my)표준화+더 정확하게 보정된,,,
resid <- data.frame(Hamilton, Y.hat, e, h, r, rs)
# 복습2겸
# Y, X1, X2를 이용한 plot도 만들어보자.
par(mar = c(3, 3, 2, 1)) # 마진 줄이기
plot(Hamilton[c("Y", "X1", "X2")], pch = 19, cex = 1)
# pch = 19 : 점모양 설정. 19는 채워진 동그라미로 설정하라는 의미
# cex = 1 : 점크기 설정. 1은 기본, 2는 두배, 0.5는 절반을 의미
# 그림에서 봐야할것은 크게 3가지다.
# 1. Y와 X1,X2간의 관계. Y와 선형관계가 있는지. 글쎄다?
# 2. X1과 X2간의 관계도 체크. 둘은 선형관계가 있어보인다. 다중공선성 의심!!
# 3. 이상치 확인. 이건 겸사겸사~
# 복습3겸
# 상관계수도 구해보자. Y, X1, X2는 모두 수치형이므로 cor()로 구할 수 있다.
cor <- cor(Hamilton[c("Y", "X1", "X2")])
# 언젠가 했던 말인것 같은데 cor<-cor~에 전체 괄호를 쳐두면 해당 cor값 출력도 가능하다.
# 출력하면 다음과 같다.
Y X1 X2
Y 1.000000000 0.002497966 0.4340688
X1 0.002497966 1.000000000 -0.8997765
X2 0.434068758 -0.899776481 1.0000000
# 복습 끝
# qq plot / 잔차 plot
# 참고로 복습을 한 이유가 있었다.
# 위에서 쓴 변수를 여기서 사용할 계획..
# 실행 안한 코드가 있다면 실행하고 오자.
# 이제 위의 Hamilton 데이터로 qqplot과 잔차plot을 구해보자.
# 그나저나 qqplot은 왜 구하는걸까?
# QQplot으로 우리는 정규성 검사를 할 수 있기 때문이다.
# 잔차plot으로는 등분산, 독립, 선형성을 검사할 수 있다.
# 참고로 이러한 특징을 검사하는 이유는
# 이들이 회귀모형의 가정사항이기 때문이다.
# 1. 잘 만들어진 모형은 "표준화된 잔차와 정규점수를 함께 순서대로 y,x축에 배치한" qqplot을 y=x에 가깝게 그려지게 한다.
# 2. 잘 만들어진 모형은 또한 잔차에 각 변수별 아직 설명 못한 패턴이 남겨두지 않는다.
# 패턴이 남아있는지는 잔차plot을 그려보고 특정한 패턴이 보이는지를 직접 눈으로 확인하면 된다.
# 1. qqplot
qqnorm(resid$r, ylim = c(-3,3))
qqline(resid$r, lwd = 1, col = "red")
# qqplot을 그리고 해당 그래프에 y=x 선을 그려보자.
# 2. 잔차플롯(resid vs X)
plot(r ~ X1, data = resid, pch = 19)
# 랜덤 패턴인게 좋다.
abline(h = c(-2, 0, 2), lty=2)
# 수평선 하나 긋기. -2,0,2에 대해 선긋기
# -2와 +2선이 안보인다. 안그린게 아니라 안보이는것이다!..!!
# 이 plot은 크게 벗어나는게 없는거임
plot(r ~ X2, data = resid, pch = 19)
abline(h = c(-2,0,2), lty=2)
# 이것도 마찬가지,,
# 3. resid vs fitted
#y와 잔차간의 잔차플롯도 그려보자,,
plot(r ~ Y.hat, data = resid, pch = 19)
abline(h = c(-2,0,2), lty=2)
# 상관관계 시각화_corrplot
# 또 다른 데이터셋을 하나 더 구경하자.
Salary <- read.table("./datasets/Salary.Survey.txt", sep = "\t", header = TRUE)
head(Salary, 5)
S X E M
1 13876 1 1 1
2 11608 1 3 0
3 18701 1 3 1
4 11283 1 2 0
5 11767 1 3 0
summary(Salary)
S X
Min. :10535 Min. : 1.0
1st Qu.:13321 1st Qu.: 3.0
Median :16436 Median : 6.0
Mean :17270 Mean : 7.5
3rd Qu.:20720 3rd Qu.:11.0
Max. :27837 Max. :20.0
E M
Min. :1.000 Min. :0.0000
1st Qu.:1.000 1st Qu.:0.0000
Median :2.000 Median :0.0000
Mean :1.978 Mean :0.4348
3rd Qu.:3.000 3rd Qu.:1.0000
Max. :3.000 Max. :1.0000
(cor <- cor(Salary))
# 언젠가 얘기했다는게 먼저 이 내용을 복기할적의 머릿속에서였나 보다.. 이제 알았네..
# 무튼 ()를 씌워주면 결과 출력도 함께 가능하다.
S X E M
S 1.0000000 0.53888588 0.2792746 0.72952541
X 0.5388859 1.00000000 -0.1914649 -0.05144053
E 0.2792746 -0.19146486 1.0000000 0.19668422
M 0.7295254 -0.05144053 0.1966842 1.00000000
# 이 상관관계를 시각화해보자.
# install.packages("corrplot")
# 안깔았으면 깔기
library(corrplot)
corrplot(cor, method = "number", number.digits = 3)
# 소수점 셋째짜리까지 노출하기. plot에 숫자를 제시하기
# 분석 주객전도(?)하기
# 위에서 썼던 데이터로 lm의 주객을 바꿔보자.
# 주객이라 쓰는게 적절한 표현인지는 모르겠다.
str(Salary)
'data.frame': 46 obs. of 4 variables:
$ S: num 13876 11608 18701 11283 11767 ...
$ X: num 1 1 1 1 1 2 2 2 2 3 ...
$ E: num 1 3 3 2 3 2 2 1 3 2 ...
$ M: num 1 0 1 0 0 1 0 0 0 0 ...# 구조+타입보기
# E와 M은 factor로 설정해볼까?
# 이들을 임시로 as.factor로 변환해 lm에 넣으면 된다.
# 실제로는 안바뀌게!
reg <- lm(S ~ X+as.factor(E)+as.factor(M), data = Salary)
# 모형을 돌릴때만, E와 M을 카테고리로 받아들여줘
# 기존 data를 바꾸진 않는다.
Call:
lm(formula = S ~ X + as.factor(E) + as.factor(M), data = Salary)
Residuals:
Min 1Q Median 3Q Max
-1884.60 -653.60 22.23 844.85 1716.47
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8035.60 386.69 20.781 < 2e-16 ***
X 546.18 30.52 17.896 < 2e-16 ***
as.factor(E)2 3144.04 361.97 8.686 7.73e-11 ***
as.factor(E)3 2996.21 411.75 7.277 6.72e-09 ***
as.factor(M)1 6883.53 313.92 21.928 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1027 on 41 degrees of freedom
Multiple R-squared: 0.9568, Adjusted R-squared: 0.9525
F-statistic: 226.8 on 4 and 41 DF, p-value: < 2.2e-16# 봐야될건 3가지. 늘 똑같다.
# 1. 각 변수들 오른쪽의 Pr, 2. Adjusted R-squared값, 3. F-statistic의 p-value값
# 1. 모든 변수 p-value값 유의미. 즉 종속변수에 모두 영향 준다. 즉, 회귀계수 모두 그대로 사용하자.
# 2. 0.95정도면 설명력이 높다고 확신에 차서 말할 수 있을것 같다.
# 3. F의 p-value는 엄청 작다. 2.2*10^-16이면.. 따라서 유의미하다. 모형 유의미하다.
# 근데 이제 살펴봐야할것.
# 첫째.
# 변수가 3개였던것 같은데 왜 더 많아졌지?
# as.factor를 하면, 각 열들이 가지는 (값 종류)-1만큼 해당 열의 변수가 생긴다.
# 둘째.
# 해석 방법도 수치형때와 달라진다.
# as.factor(E)2와 as.factor(E)3가 있는데,
# 이들의 Estimate값은 여기에 없는 as.factor(E)1을 기준으로 봐야한다.
# 즉 이 상황은 as.factor(E)1을 baseline으로 잡은거다.
# as.factor(M)1에 대해서도 마찬가지로 여기에 없는 as.factor(M)0을 기준으로 봐야한다.
# 참고로 눈치챘을 수도 있지만
# 회귀식의 범주형 변수들은 as.factor(컬럼)(값)
# <- 이런식으로 구성되어 있다.
# 쨋든 as.factor(E)1을 기준으로 as.factor(E)2의 변수가 1증가하면,
# y값이 3144.04가 증가한다.
# 마찬가지로 as.factor(E)1을 기준으로 as.factor(E)3의 변수가 1증가하면,
# y값이 2996.21가 증가한다.
# 만약 주객전도하고싶다면??
# 다시말해 baseline을 바꾸고 싶다면??
E1 <- as.numeric(Salary$E ==1)
E2 <- as.numeric(Salary$E ==2)
E3 <- as.numeric(Salary$E ==3)
Salary <- cbind(Salary, E1, E2, E3)
head(Salary, 5)
# 우선 새로운 열을 만들자.
# 예를 들어, X3을 중심으로 보고 싶다면 다음과 같이 쓰자.
reg2 <- lm(S ~ X+E1+E2+as.factor(M), data = Salary)
summary(reg2)
Call:
lm(formula = S ~ X + E1 + E2 + as.factor(M), data = Salary)
Residuals:
Min 1Q Median 3Q Max
-1884.60 -653.60 22.23 844.85 1716.47
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11031.81 383.22 28.787 < 2e-16 ***
X 546.18 30.52 17.896 < 2e-16 ***
E1 -2996.21 411.75 -7.277 6.72e-09 ***
E2 147.82 387.66 0.381 0.705
as.factor(M)1 6883.53 313.92 21.928 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1027 on 41 degrees of freedom
Multiple R-squared: 0.9568, Adjusted R-squared: 0.9525
F-statistic: 226.8 on 4 and 41 DF, p-value: < 2.2e-16
# 해석은,, 이제 여기쯤 왔으면 알아서~~~~~~~~~~~~~~~~~