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