# 상관계수
# 공분산 : 두 변수가 같은 방향으로 움직이는 정도. 측정단위에 영향을 받는 단점.
# 상관계수 : 공분산을 표준편차로 나누 값

# cov(x,y) = SUM(x - x_bar)(y - y_bar) / n - 1
# corr(x,y) = cov(x,y) / sd(x)sd(y)

# Pearson 상관계수
# --- 직선관계의 정도를 나타냄
# --- -1 ~ 1 사이의 값. 1에 가까울수록 강한 상관관계 / 0에 가까울수록 관계 없음

# Spearman 상관계수
# --- 서열척도일 경우 사용하는 비모수적 방법
# --- 직선관계가 아니어도 상관관계가 있으면 1에 가까운 값을 갖는다.

x <- 1:10
y <- x^2
plot(x, y)

# 상관계수
cor(x, y)
## [1] 0.9745586
# 자료가 순서형/순위자료일 경우.
cor(x, y, method = "kendall")
## [1] 1
cor(x, y, method = "spearman")
## [1] 1
# 상관계수가 유의한지 검정
# H0 : cor = 0 (상관관계가 없다) : p-value < 0.05 이면 H0 기각 = 유의한 상관관계가 있다
cor.test(x, y)
## 
##  Pearson's product-moment correlation
## 
## data:  x and y
## t = 12.298, df = 8, p-value = 1.778e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8926999 0.9941604
## sample estimates:
##       cor 
## 0.9745586
# 30개 부서에서 부서당 35명의 직원 설문조사
# 데이터 숫자는 해당 질문에 긍정한 직원의 비율

attitude
##    rating complaints privileges learning raises critical advance
## 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
## 7      58         67         42       56     66       68      35
## 8      71         75         50       55     70       66      41
## 9      72         82         72       67     71       83      31
## 10     67         61         45       47     62       80      41
## 11     64         53         53       58     58       67      34
## 12     67         60         47       39     59       74      41
## 13     69         62         57       42     55       63      25
## 14     68         83         83       45     59       77      35
## 15     77         77         54       72     79       77      46
## 16     81         90         50       72     60       54      36
## 17     74         85         64       69     79       79      63
## 18     65         60         65       75     55       80      60
## 19     65         70         46       57     75       85      46
## 20     50         58         68       54     64       78      52
## 21     50         40         33       34     43       64      33
## 22     64         61         52       62     66       80      41
## 23     53         66         52       50     63       80      37
## 24     40         37         42       58     50       57      49
## 25     63         54         42       48     66       75      33
## 26     66         77         66       63     88       76      72
## 27     78         75         58       74     80       78      49
## 28     48         57         44       45     51       83      38
## 29     85         85         71       71     77       74      55
## 30     82         82         39       59     64       78      39
cov(attitude)
##               rating complaints privileges  learning    raises critical
## rating     148.17126  133.77931   63.46437  89.10460  74.68851 18.84253
## complaints 133.77931  177.28276   90.95172  93.25517  92.64138 24.73103
## privileges  63.46437   90.95172  149.70575  70.84598  56.67126 17.82529
## learning    89.10460   93.25517   70.84598 137.75747  78.13908 13.46782
## raises      74.68851   92.64138   56.67126  78.13908 108.10230 38.77356
## critical    18.84253   24.73103   17.82529  13.46782  38.77356 97.90920
## advance     19.42299   30.76552   43.21609  64.19770  61.42299 28.84598
##              advance
## rating      19.42299
## complaints  30.76552
## privileges  43.21609
## learning    64.19770
## raises      61.42299
## critical    28.84598
## advance    105.85747
round( cor(attitude), 3)
##            rating complaints privileges learning raises critical advance
## rating      1.000      0.825      0.426    0.624  0.590    0.156   0.155
## complaints  0.825      1.000      0.558    0.597  0.669    0.188   0.225
## privileges  0.426      0.558      1.000    0.493  0.445    0.147   0.343
## learning    0.624      0.597      0.493    1.000  0.640    0.116   0.532
## raises      0.590      0.669      0.445    0.640  1.000    0.377   0.574
## critical    0.156      0.188      0.147    0.116  0.377    1.000   0.283
## advance     0.155      0.225      0.343    0.532  0.574    0.283   1.000
attach(attitude)

library(psych)
pairs.panels(attitude)

# complaints & rating 상관계수가 가장 높음
plot(rating, complaints)

cor.test(rating, complaints)
## 
##  Pearson's product-moment correlation
## 
## data:  rating and complaints
## t = 7.737, df = 28, p-value = 1.988e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6620128 0.9139139
## sample estimates:
##       cor 
## 0.8254176
# 단순회귀분석 (Simple Linear Regression)
# 하나의 종속변수(y)와 하나의 설명변수(x) 간의 관계를 직선으로 표현하는 방법
# 인과관계가 아닌 상관관계!

# y = b0 + b1x + e

# b0 : y 절편
# b1 : 기울기
# e (epsilon) : 오차. 확률변수. 평균 0, 표준편차 sigma인 정규분포를 따른다.

# 최소자승법 (Least Squares Method)
# 선(y 추정값)으로부터 각 점(y 관찰값)이 얼마나 떨어져 있는지 수치화
# 오차제곱합을 최소로 하는 추정방법
# 각 점이 선에 가까이 붙어있을수록 추정된 회귀식이 유의미하다.
# 회귀분석 in R

# model <- lm(y ~ x, data)  : 회귀분석
# plot(model)               : 회귀분석 관련 그래프 출력
#                           (plot, residuals vs fitted, Normal QQ, scale-location, residuals vs leverage)
# summary(model)            : 회귀분석 결과 출력
# abline(model)             : 그래프에 직선 추가
# abline(intercept, slope)
# 차의 속도와 브레이크 제동거리
head(cars)
##    speed dist
## 1      4    2
## 2      4   10
## 3      7    4
## 4      7   22
## 5      8   16
summary(cars)
##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00
plot(cars)

out <- lm(dist ~ speed, data = cars)    # 회귀분석(종속변수~설명변수)
plot(out)

summary(out)
##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Residuals:
##     Min      1Q  Median      3Q     Max
## -29.069  -9.525  -2.272   9.215  43.201
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12
library(ggplot2)
ggplot(cars, aes(speed, dist)) + geom_point() + geom_smooth(method = "lm")

# 회귀계수 추정과 해석
# 
# (1) F-검정
#
# F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12
# --- 회귀모형에 대한 유의성 검정
# --- H0 : 회귀모형이 유의하지 않다.
# --- p-value < 0.05 이므로 H0 기가. 회귀모형이 유의하다.
# --- 단순회귀분석에서는 b1 = 0 (회귀계수)를 검정하는 t-test와 동일
#
# (2) t-검정
#
# Coefficients: Est Std.    Error       t value     Pr(>|t|)    
# (Intercept)   -17.5791    6.7584      -2.601      0.0123 *        --> b0
#     speed     3.9324      0.4155      9.464       1.49e-12 ***    --> b1
# 
# --- 회귀계수에 대한 t-test. H0 : b1 = 0. p-value < 0.05 이면 유의미
# --- y = -17.5791 + 3.9324 * x
# --- 속력이 0 일때 제동거리가 
#
# (3) 결정계수 (R-squared)
#
# Multiple R-squared:  0.6511,  Adjusted R-squared:  0.6438 
# --- 결정계수 = SSR / SST
# --- 1에 가까울수록 모형이 변동량에 대해 설명력을 갖는다. 
# ANOVA table (Residuals:잔차. 평균 = 0)

# Regression    SSR (회귀식에 의해 설명되는 변동량)
# Residual      SSE (회귀식에 의해 설명되지 않는 변동량)
# Total         SST (총변동량)

anova(out)
## Analysis of Variance Table
## 
## Response: dist
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## speed      1  21186 21185.5  89.567 1.49e-12 ***
## Residuals 48  11354   236.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(dist ~ speed, data = cars, col = "blue")
abline(out, col = "red")

# 회귀진단
# No Intercept Model
# 속도가 0이면 제동거리가 0 인 것이 당연하다. b0 = 0 으로 고정.

out <- lm(dist ~ speed + 0, data = cars)
out
## 
## Call:
## lm(formula = dist ~ speed + 0, data = cars)
## 
## Coefficients:
## speed  
## 2.909
summary(out)    # y = 2.909 * x
## 
## Call:
## lm(formula = dist ~ speed + 0, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -26.183 -12.637  -5.455   4.590  50.181 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## speed   2.9091     0.1414   20.58   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.26 on 49 degrees of freedom
## Multiple R-squared:  0.8963, Adjusted R-squared:  0.8942 
## F-statistic: 423.5 on 1 and 49 DF,  p-value: < 2.2e-16
# 오차항(e, 잔차)이 추세를 보인다면 무언가 중요 정보가 모형에 포함되지 않았다는 의미
# e는 평균 0, 표준편차 sigma인 정규분포를 따르는 확률변수이다.
# --> 잔차도(residual plot)의 패턴 확인 필요!
# --> 잔차가 0과 가까울수록 좋다.

par(mfcol=c(2,2))
plot(out)

par(mfcol=c(1,1))
# Residuals vs Fitted plot (잔차도) / Normal Q-Q plot (정규성 검정)
# --> 선형 패턴이 아니고, 분산이 증가하는 경향 --> 종속변수의 log 또는 sqrt 변환 시도

par(mfcol=c(1,2))
plot(log(dist) ~ speed + 0, data = cars)
plot(sqrt(dist) ~ speed + 0, data = cars)  # sqrt (root) 변환이 적절한 패턴을 보임

par(mfcol=c(1,1))

out2 <- lm(sqrt(dist) ~ speed + 0, data = cars)
str(out2)
## List of 12
##  $ coefficients : Named num 0.397
##   ..- attr(*, "names")= chr "speed"
##  $ residuals    : Named num [1:50] -0.173 1.575 -0.777 1.913 0.826 ...
##   ..- attr(*, "names")= chr [1:50] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:50] -45.631 1.581 -0.767 1.923 0.838 ...
##   ..- attr(*, "names")= chr [1:50] "speed" "" "" "" ...
##  $ rank         : int 1
##  $ fitted.values: Named num [1:50] 1.59 1.59 2.78 2.78 3.17 ...
##   ..- attr(*, "names")= chr [1:50] "1" "2" "3" "4" ...
##  $ assign       : int 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:50, 1] -115.013 0.0348 0.0609 0.0609 0.0696 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:50] "1" "2" "3" "4" ...
##   .. .. ..$ : chr "speed"
##   .. ..- attr(*, "assign")= int 1
##   ..$ qraux: num 1.03
##   ..$ pivot: int 1
##   ..$ tol  : num 1e-07
##   ..$ rank : int 1
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 49
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = sqrt(dist) ~ speed + 0, data = cars)
##  $ terms        :Classes 'terms', 'formula'  language sqrt(dist) ~ speed + 0
##   .. ..- attr(*, "variables")= language list(sqrt(dist), speed)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "sqrt(dist)" "speed"
##   .. .. .. ..$ : chr "speed"
##   .. ..- attr(*, "term.labels")= chr "speed"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 0
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(sqrt(dist), speed)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "sqrt(dist)" "speed"
##  $ model        :'data.frame':   50 obs. of  2 variables:
##   ..$ sqrt(dist): num [1:50] 1.41 3.16 2 4.69 4 ...
##   ..$ speed     : num [1:50] 4 4 7 7 8 9 10 10 10 11 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language sqrt(dist) ~ speed + 0
##   .. .. ..- attr(*, "variables")= language list(sqrt(dist), speed)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "sqrt(dist)" "speed"
##   .. .. .. .. ..$ : chr "speed"
##   .. .. ..- attr(*, "term.labels")= chr "speed"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 0
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(sqrt(dist), speed)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "sqrt(dist)" "speed"
##  - attr(*, "class")= chr "lm"
summary(out2)   # sqrt 변환후 회귀분석 --> R-squared 값 높아짐
## 
## Call:
## lm(formula = sqrt(dist) ~ speed + 0, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2781 -0.6972  0.0208  0.7965  3.3898 
## 
## Coefficients:
##       Estimate Std. Error t value Pr(>|t|)    
## speed  0.39675    0.01015   39.09   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.167 on 49 degrees of freedom
## Multiple R-squared:  0.9689, Adjusted R-squared:  0.9683 
## F-statistic:  1528 on 1 and 49 DF,  p-value: < 2.2e-16
par(mfcol=c(2,2))
plot(out2)      # Residuals vs Fitted plot / Normal Q-Q plot 좋은 패턴을 보임

par(mfcol=c(1,1))

shapiro.test(out2$residuals)    # 잔차가 정규분포를 따른다 (p-value > 0.05)
## 
##  Shapiro-Wilk normality test
## 
## data:  out2$residuals
## W = 0.97919, p-value = 0.5185
# 최종 모형으로 추정된 회귀식
# sqrt(dist) =  0.397 * speed
# dist = (0.397 * speed)^2
# 추정과 예측
# 속도가 10 또는 30일때의 제동거리 예측
nspeed <- data.frame("speed" = c(10,30))
nspeed
##   speed
## 1    10
## 2    30
predict(out2, nspeed)
##         1         2 
##  3.967494 11.902483
# predict(out2, nspeed) ^ 2   # 결과값은 sqrt(dist) 이므로 실제값은 변환 필요

# "평균" 제동거리의 95% 신뢰구간
predict(out2, nspeed, interval = "confidence")
##         fit       lwr      upr
## 1  3.967494  3.763518  4.17147
## 2 11.902483 11.290554 12.51441
# 새로운 한 차량에 대한 95% 예측구간
predict(out2, nspeed, interval = "prediction")
##         fit      lwr       upr
## 1  3.967494 1.612650  6.322338
## 2 11.902483 9.477995 14.326970
# 모든 관측치에 대한 추정치
pred_dist <- fitted(out2)^2
cbind(cars, pred_dist)
##    speed dist pred_dist
## 1      4    2  2.518562
## 2      4   10  2.518562
## 3      7    4  7.713095
## 4      7   22  7.713095
## 5      8   16 10.074246
## 6      9   10 12.750218
## 7     10   18 15.741010
## 8     10   26 15.741010
## 9     10   34 15.741010
## 10    11   17 19.046622
## 11    11   28 19.046622
## 12    12   14 22.667055
## 13    12   20 22.667055
## 14    12   24 22.667055
## 15    12   28 22.667055
## 16    13   26 26.602307
## 17    13   34 26.602307
## 18    13   34 26.602307
## 19    13   46 26.602307
## 20    14   26 30.852380
## 21    14   36 30.852380
## 22    14   60 30.852380
## 23    14   80 30.852380
## 24    15   20 35.417273
## 25    15   26 35.417273
## 26    15   54 35.417273
## 27    16   32 40.296986
## 28    16   40 40.296986
## 29    17   32 45.491519
## 30    17   40 45.491519
## 31    17   50 45.491519
## 32    18   42 51.000873
## 33    18   56 51.000873
## 34    18   76 51.000873
## 35    18   84 51.000873
## 36    19   36 56.825047
## 37    19   46 56.825047
## 38    19   68 56.825047
## 39    20   32 62.964041
## 40    20   48 62.964041
## 41    20   52 62.964041
## 42    20   56 62.964041
## 43    20   64 62.964041
## 44    22   66 76.186489
## 45    23   54 83.269944
## 46    24   70 90.668218
## 47    24   92 90.668218
## 48    24   93 90.668218
## 49    24  120 90.668218
## 50    25   85 98.381313
plot(cars$speed, pred_dist)

# 관측치 속도의 최대값 25 --> 데이터 범위 밖에서 예측하는 것은 주의 해야한다!
newspeed2 <- data.frame("speed"= c(50, 70))
predict(out2, newspeed2)^2
##        1        2 
## 393.5253 771.3095
# Outlier (이상점) / Influential Points (영향점)
# Outlier
# - 측정상/실험상의 과오로 모집단에 속하지 않는다고 의심이 될 정도로 정상범위 밖에 떨어진 점
# - 대개 큰 잔차를 가짐

# Influential Points
# 소수의 관측치이지만 통계량에 큰 영향을 미침
# Leverage plot 에서 점선 영역 밖에 위치 (Cook' distance)
# practice 4
# 중고차. Odometer (주행거리 / 100 mile)에 따른 Price (가격 / $1000)

ucar <- read.csv("data/sonata.csv")
head(ucar)
##   Price Odometer Color
## 1  14.6     37.4     1
## 2  14.1     44.8     1
## 3  14.0     45.8     3
## 4  15.6     30.9     3
## 5  15.6     31.7     2
## 6  14.7     34.0     2
summary(ucar)
##      Price          Odometer         Color     
##  Min.   :13.60   Min.   :19.10   Min.   :1.00  
##  1st Qu.:14.47   1st Qu.:32.17   1st Qu.:1.00  
##  Median :14.70   Median :36.20   Median :2.00  
##  Mean   :14.84   Mean   :36.01   Mean   :1.96  
##  3rd Qu.:15.20   3rd Qu.:40.27   3rd Qu.:3.00  
##  Max.   :16.40   Max.   :49.20   Max.   :3.00
ucar <- ucar[-3]
attach(ucar)

# 산점도와 상관계수를 통해 선형관계를 판단
plot(Price ~ Odometer, ucar)

cor(Price, Odometer)
## [1] -0.805168
pairs.panels(ucar)    # cor = -0.81

# 회귀식 추정
model <- lm(Price ~ Odometer, data = ucar)
summary(model)
## 
## Call:
## lm(formula = Price ~ Odometer, data = ucar)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68679 -0.27263  0.00521  0.23210  0.70071 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.248727   0.182093   94.72   <2e-16 ***
## Odometer    -0.066861   0.004975  -13.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3265 on 98 degrees of freedom
## Multiple R-squared:  0.6483, Adjusted R-squared:  0.6447 
## F-statistic: 180.6 on 1 and 98 DF,  p-value: < 2.2e-16
# price = 17.249 - 0.067 * Odometer
# F-test : 회귀모형 유의함
# t-test : 회귀계수 유의함
# R-squared : 0.648

plot(model)     # 잔차도 적절

# 주행거리 3600 마일인 소나타
nOdometer <- data.frame("Odometer" = c(36, 46))
predict(model, nOdometer, interval = "confidence") # 평균가격
##        fit      lwr      upr
## 1 14.84174 14.77694 14.90653
## 2 14.17313 14.05513 14.29112
predict(model, nOdometer, interval = "prediction") # 특정한 차를 팔때 예측가격
##        fit      lwr      upr
## 1 14.84174 14.19060 15.49287
## 2 14.17313 13.51456 14.83169