This report is a summary of the lesson by Richie Cotton, Datacamp

1. Linear Regression

Basic

Response variable(y) = dependent variable, 응답변수, 종속변수
Explanatory variable(x) = indepedent variable, 설명변수, 독립변수

Linear regression

  • The response variable is numeric

Logistic regression

  • The response variable is logical

numeric vs. numeric -> scatter plot

head(taiwan_real_estate)
##   dist_to_mrt_m n_convenience house_age_years price_twd_msq
## 1      84.87882            10        30 to 45     11.467474
## 2     306.59470             9        15 to 30     12.768533
## 3     561.98450             5         0 to 15     14.311649
## 4     561.98450             5         0 to 15     16.580938
## 5     390.56840             5         0 to 15     13.040847
## 6    2175.03000             3         0 to 15      9.712557
ggplot(taiwan_real_estate, aes(x = n_convenience, y = price_twd_msq)) +
  geom_jitter(alpha = 0.5) +
  geom_smooth(method = "lm", se = F)

mdl_price_vs_conv <- lm(price_twd_msq ~ n_convenience, data = taiwan_real_estate)
summary(mdl_price_vs_conv)
## 
## Call:
## lm(formula = price_twd_msq ~ n_convenience, data = taiwan_real_estate)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.7132  -2.2213  -0.5409   1.8105  26.5299 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.22424    0.28500   28.86   <2e-16 ***
## n_convenience  0.79808    0.05653   14.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.384 on 412 degrees of freedom
## Multiple R-squared:  0.326,  Adjusted R-squared:  0.3244 
## F-statistic: 199.3 on 1 and 412 DF,  p-value: < 2.2e-16

numeric vs. factor -> histogram

범주형 회귀식의 경우 설명변수에 + 0을 해줘야 제대로된 회귀 계수 계산이 가능함.

## 범주형 설명변수가 1개만 있을 경우 회귀식 계수는 각 범주들의 평균

ggplot(taiwan_real_estate, aes(x = price_twd_msq)) +
  geom_histogram(bins = 10) +
  facet_wrap(vars(house_age_years))

mdl_price_vs_age <- lm(price_twd_msq ~ house_age_years + 0, taiwan_real_estate)
summary(mdl_price_vs_age)
## 
## Call:
## lm(formula = price_twd_msq ~ house_age_years + 0, data = taiwan_real_estate)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.3379  -2.9119   0.2218   2.5544  22.9147 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## house_age_years0 to 15   12.6375     0.2866   44.10   <2e-16 ***
## house_age_years15 to 30   9.8767     0.3478   28.40   <2e-16 ***
## house_age_years30 to 45  11.3933     0.4053   28.11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.95 on 411 degrees of freedom
## Multiple R-squared:  0.896,  Adjusted R-squared:  0.8953 
## F-statistic:  1180 on 3 and 411 DF,  p-value: < 2.2e-16

Prediction

predict

## 테스트 위한 데이터 임의 생성
explanatory_data <- tibble(n_convenience = c(0:10))

prediction_data <- explanatory_data %>% 
  mutate(
    response_var = predict(mdl_price_vs_conv, explanatory_data)
  )
head(prediction_data)
## # A tibble: 6 × 2
##   n_convenience response_var
##           <int>        <dbl>
## 1             0         8.22
## 2             1         9.02
## 3             2         9.82
## 4             3        10.6 
## 5             4        11.4 
## 6             5        12.2
ggplot(taiwan_real_estate, aes(n_convenience, price_twd_msq)) +
  geom_jitter() +
  geom_smooth(method = "lm", se = F) +
  geom_point(data = prediction_data, aes(n_convenience, response_var), col = "red")

Object

summary(mdl_price_vs_conv)
## 
## Call:
## lm(formula = price_twd_msq ~ n_convenience, data = taiwan_real_estate)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.7132  -2.2213  -0.5409   1.8105  26.5299 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    8.22424    0.28500   28.86   <2e-16 ***
## n_convenience  0.79808    0.05653   14.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.384 on 412 degrees of freedom
## Multiple R-squared:  0.326,  Adjusted R-squared:  0.3244 
## F-statistic: 199.3 on 1 and 412 DF,  p-value: < 2.2e-16

residual이 정규분포에 수렴하는지 확인해야 한다.(median이 0에 수렴하는지, 1Q와 3Q의 절대값이 비슷한 수준인지)

coefficients(mdl_price_vs_conv)
fitted(mdl_price_vs_conv)
residuals(mdl_price_vs_conv)

Response value = fitted value + residual

coefficients : 회귀식 계수를 출력
fitted : the stuff you explained
residuals : the stuff you couldn’t explain

잔차가 발생하는 이유는 보통 모델이 좋지 않기 때문이지만 일반적으로 완벽한 모델을 갖는 것은 불가능하다.

보통 randomness에 기인한 잔차들로 이것들은 장기적으로는 지속하지 못하기에 회귀식은 결국 평균에 수렴한다.

summary 출력값은 분석하기에 용이한 데이터 타입이 아니므로 다음 함수를 사용하여 분석한다.

tidy(mdl_price_vs_conv)
## # A tibble: 2 × 5
##   term          estimate std.error statistic   p.value
##   <chr>            <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)      8.22     0.285       28.9 5.81e-101
## 2 n_convenience    0.798    0.0565      14.1 3.41e- 37
augment(mdl_price_vs_conv)
## # A tibble: 414 × 8
##    price_twd_msq n_convenience .fitted .resid    .hat .sigma  .cooksd .std.resid
##            <dbl>         <dbl>   <dbl>  <dbl>   <dbl>  <dbl>    <dbl>      <dbl>
##  1         11.5             10   16.2  -4.74  0.0121    3.38  1.22e-2     -1.41 
##  2         12.8              9   15.4  -2.64  0.00913   3.39  2.83e-3     -0.783
##  3         14.3              5   12.2   2.10  0.00264   3.39  5.10e-4      0.621
##  4         16.6              5   12.2   4.37  0.00264   3.38  2.21e-3      1.29 
##  5         13.0              5   12.2   0.826 0.00264   3.39  7.92e-5      0.244
##  6          9.71             3   10.6  -0.906 0.00275   3.39  9.91e-5     -0.268
##  7         12.2              7   13.8  -1.62  0.00477   3.39  5.50e-4     -0.479
##  8         14.1              6   13.0   1.12  0.00343   3.39  1.88e-4      0.331
##  9          5.69             1    9.02 -3.33  0.00509   3.38  2.49e-3     -0.988
## 10          6.69             3   10.6  -3.93  0.00275   3.38  1.87e-3     -1.16 
## # ℹ 404 more rows
glance(mdl_price_vs_conv)
## # A tibble: 1 × 12
##   r.squared adj.r.squared sigma statistic  p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>    <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1     0.326         0.324  3.38      199. 3.41e-37     1 -1091. 2188. 2200.
## # ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

RMS(residual standard error) : glance sigma, 자유도를 포함해서 계산 RMSE(root mean square error) : RMS와 비슷하지만 자유도 포함하지 않음 대부분 RMS를 사용한다.

Mean reversion

일반적으로 대부분의 데이터들, 특히 극단에 위치한 데이터들은 평균에 수렴하는 모습을 보인다.

head(sp500_yearly_returns)
##   symbol return_2018 return_2019
## 1   MSFT  0.20795279   0.5755810
## 2   AAPL -0.05390179   0.8895780
## 3   AMZN  0.28431683   0.2302776
## 4     FB -0.25711215   0.5657183
## 5   GOOG -0.01031158   0.2910459
## 6    JNJ -0.05131476   0.1622003
ggplot(sp500_yearly_returns, aes(return_2018, return_2019)) +
  geom_point() +
  geom_abline(col = "green", size = 1) +
  geom_smooth(method = "lm", se = F) +
  coord_fixed()

Transforming Variables

만약 회귀식과 데이터들간의 관계가 선형이 아닌경우 변수를 변형해서 선형으로 만들어 모델 적합성을 높일수 있다.

제곱근, 지수화가 대표적으로 사용된다.

Square Root

ggplot(taiwan_real_estate, aes(dist_to_mrt_m, price_twd_msq)) +
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  labs(title = "NOT Applying Square root")

일반 회귀식으로 데이터들이 회귀식과 linear하지 않다. 그럼 여기에 dist_to_mrt_m에 제곱근(sqrt)을 해보자.

ggplot(taiwan_real_estate, aes(sqrt(dist_to_mrt_m), price_twd_msq)) +
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  labs(title = "Applying Square root")

lm(price_twd_msq ~ sqrt(dist_to_mrt_m), taiwan_real_estate)
## 
## Call:
## lm(formula = price_twd_msq ~ sqrt(dist_to_mrt_m), data = taiwan_real_estate)
## 
## Coefficients:
##         (Intercept)  sqrt(dist_to_mrt_m)  
##             16.7098              -0.1828

데이터들이 훨씬 회귀식과 선형 관계를 보이고 있다. 더 적합한 모델이라고 말할 수 있다.

Exponentialization

head(ad_conversion)
##   spent_usd n_impressions n_clicks
## 1      1.43          7350        1
## 2      1.82         17861        2
## 3      1.25          4259        1
## 4      1.29          4133        1
## 5      4.77         15615        3
## 6      1.27         10951        1
ggplot(ad_conversion, aes(x = n_impressions, y = n_clicks)) +
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  labs(title = "NOT Applying Exponentialization")

이번에는 데이터들이 좌측 하단으로 몰려있어 선형관계를 보이지 않고 있다. 두 변수 모두에게 지수화를 적용해보자.

주의해야할 점은 lm함수에서 지수화를 적용할 때 I함수를 사용해야 한다는 것이다!!!

ggplot(ad_conversion, aes(x = n_impressions ^ 0.25, y = n_clicks ^ 0.25)) +
  geom_point() +
  geom_smooth(method = "lm", se = F) +
  labs(title = "Applying Exponentialization")

lm(data = ad_conversion, I(n_clicks ^ 0.25) ~ I(n_impressions ^ 0.25))
## 
## Call:
## lm(formula = I(n_clicks^0.25) ~ I(n_impressions^0.25), data = ad_conversion)
## 
## Coefficients:
##           (Intercept)  I(n_impressions^0.25)  
##               0.07175                0.11153

Assessing model fit 1

mdl_click_vs_impression_orig %>% 
  glance() %>% 
  pull(r.squared)
## [1] 0.8916135
mdl_click_vs_impression_trans %>% 
  glance() %>% 
  pull(r.squared)
## [1] 0.9445273

mdl_click_vs_impression_orig R-squared 89% vs. mdl_click_vs_impression_trans R-squared 94% 로 지수화를 적용한 모델이 better fit model이다.

mdl_click_vs_impression_orig %>% 
  glance() %>% 
  pull(sigma)
## [1] 19.90584
mdl_click_vs_impression_trans %>% 
  glance() %>% 
  pull(sigma)
## [1] 0.1969064

mdl_click_vs_impression_orig sigma 19.9 vs. mdl_click_vs_impression_trans sigma 0.19 로 지수화를 적용한 모델이 better fit model이다.

Assessing model fit 2 - Outliers, leverage and influence

Outliers

이상치(outliers)는 보통 2가지로 분류가능

  1. 설명변수가 극단적인 경우
  2. 응답변수가 회귀식으로부터 멀리 떨어진 경우

Leverage

leverage is a measure of how extreme the explanatory variable values are.

설명변수가 얼마나 극단에 위치한 정도에 대한 척도로 augment의 .hat 칼럼이 leverage 이다.

Influence

Influence measures how much the model would change if you left the observation out of the dataset when modeling.

leverage가 높은 관측치들을 제외하고 모델을 다시 실행했을 때 얼마나 모델이 변하는지에 대한 척도

이는 residual과 leverage의 함수로 이들의 단순한 곱이 아니라 cook's distance라는 척도를 사용하며 잔차와 레버리지가 클수록 영향도 커진다.

augment의 .cooksd 에서 확인이 가능하다.

또한 autoplot함수를 통해 fit 척도들을 그래프로 한눈에 확인할 수 있다.

## which arg value
## 1. residuals vs. fitted values : residuals이 0에 근접하는지
## 2. Q-Q plot
## 3. scale-location
## 4. Cook's distance
## 5. residuals vs. leverage
## 6. Cook's distance vs. leverage

autoplot(mdl_click_vs_impression_trans, which = 1:6)

2. Logistic Regression

응답변수가 logical일 때 사용하며 응답변수는 logistic curve(S shaped)를 따른다.

Basic

head(churn)
##   has_churned time_since_first_purchase time_since_last_purchase
## 1           0               -1.08922097               -0.7213215
## 2           0                1.18298297                3.6344354
## 3           0               -0.84615637               -0.4275823
## 4           0                0.08694165               -0.5356717
## 5           0               -1.16664155               -0.6726400
## 6           0                0.49339968               -0.7700030
mdl_recency_glm <- glm(has_churned ~ time_since_last_purchase, data = churn, family = binomial)

summary(mdl_recency_glm)
## 
## Call:
## glm(formula = has_churned ~ time_since_last_purchase, family = binomial, 
##     data = churn)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              -0.03502    0.10152  -0.345  0.73013   
## time_since_last_purchase  0.26921    0.09812   2.744  0.00607 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 554.52  on 399  degrees of freedom
## Residual deviance: 546.40  on 398  degrees of freedom
## AIC: 550.4
## 
## Number of Fisher Scoring iterations: 4

logistic regression이 일반 regression과 다른 점은 다음과 같다.

  1. glm 사용
  2. famaily = binomial 변수 설정
ggplot(churn, aes(x = time_since_last_purchase, y = has_churned)) +
  geom_point() +
  geom_smooth(method = "glm", se = F, method.args = list(family = binomial))

ggplot2에서 사용할 때에도 역시 method.args = list(family = binomial) 을 개별적으로 설정해줘야 한다.

Odds ratio

  • odds ration = prob / (1 -prob)
  • 특정 사건이 발생할 확률이 0.75 일때, odds ratio = 0.75 / 0.25 = 3
  • 즉, 사건이 발생할 확률이 발생하지 않을 확률보다 3배 높음을 의미

comparing scales

  1. most likely outcome : yes, no 만 있기에 이해하기 쉽지만 정확도가 떨어짐
  2. odds ratio & probability : 이해하기 쉽지만 비선형 관계에서 설명력이 떨어짐
  3. log odds ratio : 이해하기는 어렵지만 설명변수와 선형관계로 설명 가능
mdl_churn_vs_relationship <- glm(has_churned ~ time_since_first_purchase, data = churn, family = binomial)
summary(mdl_churn_vs_relationship)
## 
## Call:
## glm(formula = has_churned ~ time_since_first_purchase, family = binomial, 
##     data = churn)
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)   
## (Intercept)               -0.01518    0.10156  -0.150  0.88115   
## time_since_first_purchase -0.35479    0.11095  -3.198  0.00139 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 554.52  on 399  degrees of freedom
## Residual deviance: 543.73  on 398  degrees of freedom
## AIC: 547.73
## 
## Number of Fisher Scoring iterations: 4
## 임의 데이터 생성
explanatory_data <- tibble(
  time_since_first_purchase = seq(-1.5, 0.75, 0.25)
)

prediction_data <- explanatory_data %>% 
  mutate(
    has_churned = predict(mdl_churn_vs_relationship, explanatory_data, type = "response"),
    most_likely_outcome = round(has_churned),
    odds_ratio = has_churned / (1 - has_churned),
    log_odds_ratio = log(odds_ratio),
    ## response 설정하지 않으면 자동으로 log 로 계산됨
    log_odds_ratio2 = predict(mdl_churn_vs_relationship, explanatory_data)
    
  )

knitr::kable(prediction_data, align = "c", label = "scales 비교")
time_since_first_purchase has_churned most_likely_outcome odds_ratio log_odds_ratio log_odds_ratio2
-1.50 0.6264479 1 1.6770025 0.5170080 0.5170080
-1.25 0.6054699 1 1.5346606 0.4283092 0.4283092
-1.00 0.5840959 1 1.4044004 0.3396105 0.3396105
-0.75 0.5624009 1 1.2851966 0.2509117 0.2509117
-0.50 0.5404646 1 1.1761107 0.1622130 0.1622130
-0.25 0.5183703 1 1.0762839 0.0735142 0.0735142
0.00 0.4962039 0 0.9849302 -0.0151845 -0.0151845
0.25 0.4740525 0 0.9013305 -0.1038833 -0.1038833
0.50 0.4520027 0 0.8248267 -0.1925820 -0.1925820
0.75 0.4301398 0 0.7548164 -0.2812808 -0.2812808

odds ratio에 log를 취해 log odds ratio 게산이 가능하다.

또한, predict 함수에서 type = "reponse"를 사용하지 않아도 log odds ratio 계산이 된다. 위 테이블에서 log odds ratio와 log odds ratio2의 값이 일치함을 확인 할 수 있다.

Confusion matrix

mdl_recency <- glm(has_churned ~ time_since_last_purchase, data = churn, family = binomial)
actual_response <- churn$has_churned
predicted_response <- round(fitted(mdl_recency))
outcomes <- table(predicted_response, actual_response)

knitr::kable(outcomes)
0 1
0 141 111
1 59 89
  • model 예측값 기준으로 positive negative 구분
  • 141명이 이탈하지 않고 89명이 이탈한다고 정확히 예측
    • true negative(141, TN) : 실제 이탈하고 모델도 이탈 예측
    • true positve(89, TP) : 실제 이탈하지 않고 모델도 이탈하지 않는다고 예측
    • false negative (111, FN) : 실제 이탈했지만 모델은 이탈하지 않는다고 예측
    • false positive(59, FP) : 실제 이탈하지 않았지만 모델은 이탈했다고 예측
library(yardstick)

confusion <- yardstick::conf_mat(outcomes)
confusion
##                   actual_response
## predicted_response   0   1
##                  0 141 111
##                  1  59  89
autoplot(confusion)

## factor 대신에 0, 1을 사용했기에 "second" 사용
summary(confusion, event_level = "second")
## # A tibble: 13 × 3
##    .metric              .estimator .estimate
##    <chr>                <chr>          <dbl>
##  1 accuracy             binary         0.575
##  2 kap                  binary         0.15 
##  3 sens                 binary         0.445
##  4 spec                 binary         0.705
##  5 ppv                  binary         0.601
##  6 npv                  binary         0.560
##  7 mcc                  binary         0.155
##  8 j_index              binary         0.150
##  9 bal_accuracy         binary         0.575
## 10 detection_prevalence binary         0.37 
## 11 precision            binary         0.601
## 12 recall               binary         0.445
## 13 f_meas               binary         0.511

Accuracy

summary(confusion) %>% 
  slice(1)
## # A tibble: 1 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.575

올바른 예측의 비율(TRUE, FALSE 모두 포함)
TN + TP / TN + FN + FP + TP
클수록 GOOD

Sensitivity

summary(confusion) %>% 
  slice(3)
## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 sens    binary         0.705

예측관측치가 TRUE인 비율
TP / TP + FN
클수록 GOOD

Specificity

summary(confusion) %>% 
  slice(4)
## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 spec    binary         0.445

예측관측치가 FALSE인 비율
TN / TN + FP
클수록 GOOD