Linear Regression

######################################################################################################
# Linear Regression & GLM                                                                            #
# ref : machine learning with R, CH 6                                                                #
# ref : https://onlinecourses.science.psu.edu/stat504/node/216                                       #
# ref : http://statmath.wu-wien.ac.at/courses/heather_turner/glmCourse_001.pdf                       #
# ref : https://www.sagepub.com/sites/default/files/upm-binaries/21121_Chapter_15.pdf                #
# ref : https://rstudio-pubs-static.s3.amazonaws.com/41074_62aa52bdc9ff48a2ba3fb0f468e19118.html     #
# ref : https://rpubs.com/FelipeRego/MultipleLinearRegressionInRFirstSteps                           #
# ref : http://blog.naver.com/PostView.nhn?blogId=jindog2929&logNo=10161293384                       #
# ref : https://stats.idre.ucla.edu/r/modules/coding-for-categorical-variables-in-regression-models/ #
# ref : http://rstudio-pubs-static.s3.amazonaws.com/190997_40fa09db8e344b19b14a687ea5de914b.html     #
######################################################################################################
################################################################################
# terms                                                                        #
# dependent var / independent var                                              #
# covariance / correlation                                                     #
# simple linear regression / multiple linear regression                        #
# OLS ( Ordinary Least Squares ) Estimation                                    #
# dummy coding / dummy variable                                                #
# GLM / link function                                                          #
# Transformation ( to binary )                                                 #
# parameter estimate                                                           #
# Term : Interaction / non-linear                                              #
# Evaluation : Residual , RMSE , p-value , significance level, MAE             #
# residuals                                                                    #
# SDR                                                                          #
# dummy variable / dummy coding                                                #
################################################################################
################################################################################
# package & methods                                                            #
# glm : https://stat.ethz.ch/R-manual/R-devel/library/stats/html/glm.html      #
# corrplot : https://cran.r-project.org/web/packages/corrplot/corrplot.pdf     #
################################################################################

# 빠르게 Linear Regression 을 훑고. GLM 을 중심으로 살펴보면 좋겠습니다;
# Linear Regression 의 경우, 몇가지 통계적 가정을 바탕으로 하고 있기때문에, 
# 서비스단에서는 생각보다 적용하기가 쉽지않고, GLM 이 좋은 선택안이 되는 경우가 있다. 
# 기본개념과 evaluation 과정보다는, 가정이 잘 만족하는지 검증하는 부분에 방점을 찍었으면 한다.

0. Basic Terms

# dependent variable  : the value to be predicted , 종속변수
# indepndent variable : the predictors , 독립변수 , 예측 변수 , 설명변수 , feature

1. Simple Linear Regression

# a single independent variable
# OLS
## In order to determine the optimal estimates of α and β, \
## an estimation method known as Ordinary Least Squares (OLS) was used
## error : the difference between the actual y value and the predicted y value

# 다소 쓸데없지만, 직접 simple linear regression 을 수행해보고, 이를 library method 결과와 비교해본다.
# a. 데이터 분포 
plot(faithful$waiting, faithful$eruptions)

# b. 직접 계산 
b <- stats::cov(faithful$waiting, faithful$eruptions) / stats::var(faithful$waiting)
a <- mean(faithful$eruptions) - b * mean(faithful$waiting)
plot(faithful$waiting, faithful$eruptions); abline(a = a, b = b, col = "blue") 

# c. 비교 
model <- lm(eruptions ~ waiting, data = faithful)
intercept <- paste("intercept  manual : library = " , round(a,5) , " : ", round(model$coefficients[[1]],5))
slope     <- paste("slope manual : library = " , round(b,5) , " : ", round(model$coefficients[[2]],5))
plot(faithful$waiting, faithful$eruptions, main = paste(intercept, slope, sep = "\n")); abline(model, col = "red")

2. Multiple Linear Regression

# two or more independent variables 
# The goal is now to solve for β, 
# the vector of regression coef cients that minimizes the sum of the squared errors 
# between the predicted and actual Y values.
# psuedo inverse : http://b.mytears.org/tag/pseudo-inverse

# 다소 쓸데없지만, 직접 multiple linear regression 을 수행해보고, 이를 library method 결과와 비교해본다.
# formula : mpg ~ cyl + wt 
# a. 직접 계산 값과 lm function 결과 비교 
reg <- function(y, x) {
  x <- as.matrix(x)
  x <- cbind(Intercept = 1, x)
  b <- solve(t(x) %*% x) %*% t(x) %*% y
  colnames(b) <- "estimate"
  print(b)
}
model <- lm(mpg ~ cyl + wt, data = mtcars)
reg(mtcars$mpg, mtcars[c("cyl", "wt")]); print(model)

           estimate
Intercept 39.686261
cyl       -1.507795
wt        -3.190972

Call:
lm(formula = mpg ~ cyl + wt, data = mtcars)

Coefficients:
(Intercept)          cyl           wt  
     39.686       -1.508       -3.191  

# 이를 역시, 다소 쓸데없지만 시각화를 해보자.
library(rgl)
library(rglwidget)
library(knitr)
knit_hooks$set(webgl = hook_webgl)
test_data <- expand.grid(cyl=seq(4,8,by=1),wt=seq(1.5,5.4,by=.1))
test_data$mpg <- predict(model, newdata = test_data, type = "response")
with(mtcars, plot3d(cyl,wt,mpg, col="blue", size=1, type="s", main="3D Linear Model Fit"))
with(test_data, surface3d(unique(cyl),unique(wt), mpg, alpha=0.3,front="line", back="line"))
rglwidget()

3. Testing of Assumptions

# ref : http://r-statistics.co/Assumptions-of-Linear-Regression.html
# 위의 링크는 조금 더 많은 가정을 체크하지만, 일단 아래 네개정도만 봐도 헉헉이라 
# Assumptions
## 정규성(Normality) 
## 독립성(Independence) 
## 선형성(Linearity) 
## 등분산성(Homoscedasticity) 
# a. 일단 종합 검진 gvlma
# ref : https://cran.r-project.org/web/packages/gvlma/gvlma.pdf
# Global Stat 값을 보면되는데, 이 값이 통과과 되지 않는다고 하면, 뭐가 문제인지 하나하나 봐야한다. 
# 예는 아깝게 탈락 ;; ( 0.042 > 0.05 )
library(gvlma)
gvmodel<-gvlma(model)
summary(gvmodel)

Call:
lm(formula = mpg ~ cyl + wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2893 -1.5512 -0.4684  1.5743  6.1004 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)  39.6863     1.7150  23.141 < 0.0000000000000002 ***
cyl          -1.5078     0.4147  -3.636             0.001064 ** 
wt           -3.1910     0.7569  -4.216             0.000222 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.568 on 29 degrees of freedom
Multiple R-squared:  0.8302,    Adjusted R-squared:  0.8185 
F-statistic: 70.91 on 2 and 29 DF,  p-value: 0.000000000006809


ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = model) 

                    Value p-value                   Decision
Global Stat        9.8891 0.04234 Assumptions NOT satisfied!
Skewness           3.3162 0.06860    Assumptions acceptable.
Kurtosis           0.1101 0.74000    Assumptions acceptable.
Link Function      5.5472 0.01851 Assumptions NOT satisfied!
Heteroscedasticity 0.9156 0.33864    Assumptions acceptable.

# ref : http://data.library.virginia.edu/diagnostic-plots/
# 기본 plot 으로 체크해보자 
par(mfrow = c(2,2))
plot(model)

## 정규성(Normality) : 실제값(종속변수)가 정규분포이면, Residuals 도 정규분포이고 , 평균은 0
## 아래 그림에서 CASE1 과 같이 선위에 올라와 있으면, 만족. ( CASE1 : Good , CASE2 : Bad )
## 보통 Normal Q-Q plot 으로 확인한다. 
## car::qqPlot 은 신뢰구간 95% 선을 추가로 보여준다. 
par(mfrow = c(2,2))
plot(model, which = c(2))
plot(out, which = c(2))

car::qqPlot(model, id.method = "identify", simulate = T) # only lm model , not working for glm model

## 선형성(Linearity) : 종속변수와 독립변수가 선형관계에 있다면 잔차와 예측치 사이에 어떤 체계적인 관계가 있으면 안 된다. random noise 만 보여야 한다. 
## 아래 그림에서 CASE1 과 같이 어떤 패턴이 안보이면, 만족. ( CASE1 : Good , CASE2 : Bad )
## 보통 Residuals vs Fitted 으로 확인한다. 
par(mfrow = c(1,2))
plot(model, which = c(2))
plot(out, which = c(2))

## 등분산성(Homoscedasticity) : 분산이 일정하다는 가정. 
## 아래 그림에서 CASE1 과 같이 수평선 주위로 random band 를 나타내면, 만족. ( CASE1 : Good , CASE2 : Bad )
## 보통 Spread-Location plot 으로 확인한다. 
par(mfrow = c(1,2))
plot(model, which = c(3))
plot(out, which = c(3))

## 독립성(Independence) : 종속 변수의 독립성 검정 
##                  ( 시계열 데이터에서 짧은 간격으로 수집된 데이터등에서는, 높은 자기상관을 가질 수 있음 )
# 그래프로는 확인이 안되고 별도의 함수를 사용하자. 
# H1 이 자기상관 ( Alternative hypothesis: rho != 0 )
# 높은 자기상관을 갖는 경우, Cochrane - Orcutt 방법과 Prais - winsten 방법을 사용한다. 
durbinWatsonTest(model)

 lag Autocorrelation D-W Statistic p-value
   1       0.1302185      1.671096   0.296
 Alternative hypothesis: rho != 0

## 기본 그래프의 마지막인 

4. Model Evaluation

# precision / recall 등 model evaluation 에서 수행한 모델 평가외에 
# summary function 을 통한 리포트를 해석해보자 ( python summary 는 한큐에 다 보여준다;; )
summary(model)

Call:
lm(formula = mpg ~ cyl + wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2893 -1.5512 -0.4684  1.5743  6.1004 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)  39.6863     1.7150  23.141 < 0.0000000000000002 ***
cyl          -1.5078     0.4147  -3.636             0.001064 ** 
wt           -3.1910     0.7569  -4.216             0.000222 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.568 on 29 degrees of freedom
Multiple R-squared:  0.8302,    Adjusted R-squared:  0.8185 
F-statistic: 70.91 on 2 and 29 DF,  p-value: 0.000000000006809

# Residuals 
## Real Value - Predicted Value ( y - ŷ )
plot(mtcars$mpg, mtcars$pred, ylim = c(10, 35)); abline(a = 0, b = 1, col = "red")

# ref : http://blog.yhat.com/posts/r-lm-summary.html
# Coefficients
## p-value : probability for a null hypothesis that the coefficients have values of zero
## Signif. codes : p-value 에 의한 유의성 
## Estimate : slope var.
## Std. Error of the Coefficient Estimate : the average amount that the coefficient estimates vary from the actual average value of predicted value ( 적을수록 좋음 )
## t value of the Coefficient Estimate : standard deviations our coefficient estimate is far away from 0 ( 클 수록 좋은 값 / p-value 와 반대 )
# Residual standard error 
# R-squared
## Metric for evaluating the goodness of fit of your model. Higher is better with 1 being the best

5. (etc) Covariance & Correlation

# a. covariance
# Cov(X, Y) > 0    X가 증가 할 때 Y도 증가한다.
# Cov(X, Y) < 0    X가 증가 할 때 Y는 감소한다.
# Cov(X, Y) = 0    공분산이 0이라면 두 변수간에는 아무런 선형관계가 없으며 두 변수는 서로 독립적인 관계에 있음을 알 수 있다.
#                  그러나 두 변수가 독립적이라면 공분산은 0이 되지만, 공분산이 0이라고 해서 항상 독립적이라고 할 수 없다.
                
# Cov(X, X) = Var(X)
# Cov(X, Y) = Cov(Y, X)
# Cov(aX, bY) = ab * Cov(X, Y)
# b. Correlation
# indicates how closely their relationship follows a straight line
# covariance 는 unboundary value , 편의를 위해 -1 ~ 1 구간으로 한정되는 Correlation 을 사용 ( pearson )
# weak 0.1 ~ 0.3 / moderate 0.3 ~ 0.5 / string 0.5 ~ 

# correlation 은 feature selection 이 redundant 한 feature 를 제거하는데 사용되며, 
# 시각화해서 확인을 많이 하는데, 몇가지 방법은 아래와 같다. 

cor(mtcars[c("mpg", "cyl", "hp", "disp", "drat")])

            mpg        cyl         hp       disp       drat
mpg   1.0000000 -0.8521620 -0.7761684 -0.8475514  0.6811719
cyl  -0.8521620  1.0000000  0.8324475  0.9020329 -0.6999381
hp   -0.7761684  0.8324475  1.0000000  0.7909486 -0.4487591
disp -0.8475514  0.9020329  0.7909486  1.0000000 -0.7102139
drat  0.6811719 -0.6999381 -0.4487591 -0.7102139  1.0000000

pairs(mtcars[c("mpg", "cyl", "hp", "disp", "drat")])

library(psych)
pairs.panels(mtcars[c("mpg", "cyl", "hp", "disp", "drat")])

library(corrplot)
mcor <- cor(mtcars[c("mpg", "cyl", "hp", "disp", "drat")])
round(mcor, 2)

       mpg   cyl    hp  disp  drat
mpg   1.00 -0.85 -0.78 -0.85  0.68
cyl  -0.85  1.00  0.83  0.90 -0.70
hp   -0.78  0.83  1.00  0.79 -0.45
disp -0.85  0.90  0.79  1.00 -0.71
drat  0.68 -0.70 -0.45 -0.71  1.00

corrplot(mcor, order="AOE",type="upper",tl.pos="d")
corrplot(mcor, add=TRUE, type="lower", method="number",order="AOE",diag=FALSE,tl.pos="n", cl.pos="n")

sf <- as.data.frame(scale(faithful, center = F))
sm <- lm(sf$eruptions ~ sf$waiting)
plot(sf$waiting, sf$eruptions, xlim = c(0, 1.5), ylim = c(0,1.5))
abline(sm, col = "red")

abline(a = 0, b = method_cor, col = "blue")

6. (etc) Dummy Coding & Variable

# 막할때는 안보이다가, 조금 살펴보면 보이는 것이 categorical type feature 이다. 
# 그냥 numeric 으로 처리할때와 factor 로 처리할 때의 차이가 있는데, dummy coding 기법이 그것이다.
# dummy coding 은 factor type 의 categorical feature ( independent var ) 를 
# binary type 의 dummy var 로 쪼개어 분류하는 기법을 말하며,
# 이때 생성된 binary type 의 variable 을 dummy variable 이라 한다.
# 아래 예제의 data.frame 은 race.f 라는 categorical feature 를 factor 타입으로 가지고 있다.
hsb2 <- read.csv("https://stats.idre.ucla.edu/stat/data/hsb2.csv")
hsb2$race.f <- factor(hsb2$race)
table(hsb2$race.f); levels(hsb2$race.f)

  1   2   3   4 
 24  11  20 145 
[1] "1" "2" "3" "4"

# 이 경우, lm 을 돌려보면, 
# independent var 로 추가하지 않은, race.f2 / race.f3 / race.f4 가 
# dummy coding 에 의해 dummy variable 로 추가됨을 알 수 있다. 
summary(lm(write ~ race.f, data = hsb2))

Call:
lm(formula = write ~ race.f, data = hsb2)

Residuals:
     Min       1Q   Median       3Q      Max 
-23.0552  -5.4583   0.9724   7.0000  18.8000 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)   46.458      1.842  25.218 < 0.0000000000000002 ***
race.f2       11.542      3.286   3.512             0.000552 ***
race.f3        1.742      2.732   0.637             0.524613    
race.f4        7.597      1.989   3.820             0.000179 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.025 on 196 degrees of freedom
Multiple R-squared:  0.1071,    Adjusted R-squared:  0.0934 
F-statistic: 7.833 on 3 and 196 DF,  p-value: 0.00005785

# race.f 의 타입이 factor 가 아니면, dummy coding 은 적용되지 않는다. 
# dummy coding 이 발생한 경우, race.f 의 category 별 significance 가 각각 확인되지만, 
# 일반 numeric feature 로 취급되는 경우 퉁쳐진다. 이것이 차이다. 
hsb2 <- transform(hsb2, race.f = as.integer(as.character(race.f)))
summary(lm(write ~ race.f, data = hsb2))

Call:
lm(formula = write ~ race.f, data = hsb2)

Residuals:
    Min      1Q  Median      3Q     Max 
-22.919  -5.912   1.091   8.082  17.100 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)  45.8941     2.2652  20.260 < 0.0000000000000002 ***
race.f        2.0061     0.6322   3.173              0.00175 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.27 on 198 degrees of freedom
Multiple R-squared:  0.0484,    Adjusted R-squared:  0.04359 
F-statistic: 10.07 on 1 and 198 DF,  p-value: 0.001747