본문 바로가기

Deep Learning + AI

데이터분석응용 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)

hamilton data


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)

resid

 

 

# 복습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는 절반을 의미

plot

 

# 그림에서 봐야할것은 크게 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 선을 그려보자.

qqplot

 

 

# 2. 잔차플롯(resid vs X)

 

plot(r ~ X1, data = resid, pch = 19)

# 랜덤 패턴인게 좋다.

 

abline(h = c(-2, 0, 2), lty=2) 

# 수평선 하나 긋기. -2,0,2에 대해 선긋기

잔차플롯1

 

# -2와 +2선이 안보인다. 안그린게 아니라 안보이는것이다!..!!

# 이 plot은 크게 벗어나는게 없는거임

 

plot(r ~ X2, data = resid, pch = 19)
abline(h = c(-2,0,2), lty=2)

잔차플롯2

 

# 이것도 마찬가지,,

 

 

 

# 3. resid vs fitted

 

#y와 잔차간의 잔차플롯도 그려보자,,


plot(r ~ Y.hat, data = resid, pch = 19)
abline(h = c(-2,0,2), lty=2)

잔차플롯3

 

 

 

 

# 상관관계 시각화_corrplot

 

# 또 다른 데이터셋을 하나 더 구경하자.

 

Salary <- read.table("./datasets/Salary.Survey.txt", sep = "\t", header = TRUE)

salary

 

 

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에 숫자를 제시하기

corplot

 

 

 

 # 분석 주객전도(?)하기

# 위에서 썼던 데이터로 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

 

# 해석은,, 이제 여기쯤 왔으면 알아서~~~~~~~~~~~~~~~~~

반응형