CH 04. 회귀와 예측
setting
library(MASS)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(ggplot2)
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:dplyr':
##
## intersect, setdiff, union## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(splines)
library(mgcv)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapse## This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.## Import datasets needed for chapter 4
PSDS_PATH <- file.path('/Users/jinameliachoi/Documents/statistics-for-data-scientists')
lung <- read.csv(file.path(PSDS_PATH, 'data', 'LungDisease.csv'))
zhvi <- read.csv(file.path(PSDS_PATH, 'data', 'County_Zhvi_AllHomes.csv'))
zhvi <- unlist(zhvi[13,-(1:5)])
dates <- parse_date_time(paste(substr(names(zhvi), start=2, stop=8), "01", sep="."), "Ymd")
zhvi <- data.frame(ym=dates, zhvi_px=zhvi, row.names = NULL) %>%
mutate(zhvi_idx=zhvi/last(zhvi))
house <- read.csv(file.path(PSDS_PATH, 'data', 'house_sales.csv'), sep='\t')단순선형회귀
단순선형회귀를 통해 X가 얼만큼 변하면 Y가 어느 정도 변하는지를 정확히 추정할 수 있음
example)
노동자들이 면진(Exposure)에 노출된 연수와 폐활량(PEFR) 사이의 어떤 관계가 있을 지 선형회귀분석을 통해 알아보자.
plot(formula=PEFR ~ Exposure, data=lung)
model <- lm(PEFR ~ Exposure, data=lung)
model##
## Call:
## lm(formula = PEFR ~ Exposure, data = lung)
##
## Coefficients:
## (Intercept) Exposure
## 424.583 -4.185par(mar=c(4,4,0,0)+.1)
plot(lung$Exposure, lung$PEFR, xlab="Exposure", ylab="PEFR", ylim=c(300,450), type="n", xaxs="i")
abline(a=model$coefficients[1], b=model$coefficients[2], col="blue", lwd=2)
text(x=.3, y=model$coefficients[1], labels=expression("b"[0]), adj=0, cex=1.5)
x <- c(7.5, 17.5)
y <- predict(model, newdata=data.frame(Exposure=x))
segments(x[1], y[2], x[2], y[2] , col="red", lwd=2, lty=2)
segments(x[1], y[1], x[1], y[2] , col="red", lwd=2, lty=2)
text(x[1], mean(y), labels=expression(Delta~Y), pos=2, cex=1.5)
text(mean(x), y[2], labels=expression(Delta~X), pos=1, cex=1.5)
text(mean(x), 400, labels=expression(b[1] == frac(Delta ~ Y, Delta ~ X)), cex=1.5)
노동자가 면진에 노출되는 연수가 1씩 증가할 때마다, 폐활량은 -4.185의 비율로 줄어든다고 해석할 수 있다.
fitted <- predict(model)
resid <- residuals(model)par(mar=c(4,4,0,0)+.1)
lung1 <- lung %>%
mutate(Fitted=fitted,
positive = PEFR>Fitted) %>%
group_by(Exposure, positive) %>%
summarize(PEFR_max = max(PEFR),
PEFR_min = min(PEFR),
Fitted = first(Fitted)) %>%
ungroup() %>%
mutate(PEFR = ifelse(positive, PEFR_max, PEFR_min)) %>%
arrange(Exposure)
plot(lung$Exposure, lung$PEFR, xlab="Exposure", ylab="PEFR")
abline(a=model$coefficients[1], b=model$coefficients[2], col="blue", lwd=2)
segments(lung1$Exposure, lung1$PEFR, lung1$Exposure, lung1$Fitted, col="red", lty=3)
다중선형회귀
예측변수가 여러 개라면 다중선형회귀 분석으로 들어간다.
example) 주택 가격 평가사는 세금을 산정할 목적으로 주택 가치를 추정해야 한다.
킹 카운티의 주택 가격 데이터에서 다중회귀분석을 통해 주택 가치를 추정해보자.
head(house[, c('AdjSalePrice', 'SqFtTotLiving', 'SqFtLot', 'Bathrooms', 'Bedrooms', 'BldgGrade')])## AdjSalePrice SqFtTotLiving SqFtLot Bathrooms Bedrooms BldgGrade
## 1 300805 2400 9373 3.00 6 7
## 2 1076162 3764 20156 3.75 4 10
## 3 761805 2060 26036 1.75 4 8
## 4 442065 3200 8618 3.75 5 7
## 5 297065 1720 8620 1.75 4 7
## 6 411781 930 1012 1.50 2 8house_lm <- lm(AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms + BldgGrade, data=house, na.action = na.omit)
house_lm##
## Call:
## lm(formula = AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms +
## Bedrooms + BldgGrade, data = house, na.action = na.omit)
##
## Coefficients:
## (Intercept) SqFtTotLiving SqFtLot Bathrooms Bedrooms
## -5.219e+05 2.288e+02 -6.051e-02 -1.944e+04 -4.778e+04
## BldgGrade
## 1.061e+05계수를 해석하는 방식은 단순선형회귀와 같다.
예를 들어, 주택에 1제곱피트를 추가하면 가격이 대략 229달러 정도 증가할 것이다(SqFtTotLiving). 1,000제곱피트를 추가하면 228,800달러 증가할 것이다.
모형 평가에서 가장 중요한 성능 지표는 제곱근 평균제곱오차 (RMSE), 예측된 y햇 값들의 평균제곱오차의 제곱근을 말한다.
또 다른 유용한 지표는 결정계수라 부르는 R 제곱 통계량이다.
모델이 데이터에 얼마나 적합하지 평가하고자 할 때, 회귀분석을 설명하기 위한 용도로 활용된다.
summary(house_lm)##
## Call:
## lm(formula = AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms +
## Bedrooms + BldgGrade, data = house, na.action = na.omit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1199508 -118879 -20982 87414 9472982
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.219e+05 1.565e+04 -33.349 < 2e-16 ***
## SqFtTotLiving 2.288e+02 3.898e+00 58.699 < 2e-16 ***
## SqFtLot -6.051e-02 6.118e-02 -0.989 0.323
## Bathrooms -1.944e+04 3.625e+03 -5.362 8.32e-08 ***
## Bedrooms -4.778e+04 2.489e+03 -19.194 < 2e-16 ***
## BldgGrade 1.061e+05 2.396e+03 44.287 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 261200 on 22683 degrees of freedom
## Multiple R-squared: 0.5407, Adjusted R-squared: 0.5406
## F-statistic: 5340 on 5 and 22683 DF, p-value: < 2.2e-16해당 모델은 R-squared 값이 0.5407로 해당 변수들이 모델을 약 54% 정도 설명할 수 있다고 해석한다.
모형 선택 및 단계적 회귀
오컴의 면도날(모든 것이 동일한 조건에서는 복잡한 모델보다는 단순한 모델을 우선 사용해야 한다) 원리를 사용하여, 변수를 지정해야한다.
모델에 항을 추가할 수록 불이익을 주는 AIC를 최소화하는 모델을 검색한다.
부분집합회귀로서 모든 가능한 모델을 검색하거나,
단계적회귀 stepwise regression (예측변수를 연속적으로 추가/삭제하여 AIC를 낮추는 모델)을 사용한다.
house_full <- lm(AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms + BldgGrade + PropertyType + NbrLivingUnits + SqFtFinBasement + YrBuilt + YrRenovated + NewConstruction, data=house, na.action = na.omit)
step_lm <- stepAIC(house_full, direction="both")## Start: AIC=563193.9
## AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms +
## BldgGrade + PropertyType + NbrLivingUnits + SqFtFinBasement +
## YrBuilt + YrRenovated + NewConstruction
##
## Df Sum of Sq RSS AIC
## - NbrLivingUnits 1 6.4936e+09 1.3662e+15 563192
## - NewConstruction 1 1.0188e+10 1.3662e+15 563192
## - YrRenovated 1 2.4839e+10 1.3662e+15 563192
## - SqFtLot 1 1.0643e+11 1.3663e+15 563194
## <none> 1.3662e+15 563194
## - SqFtFinBasement 1 1.3984e+11 1.3664e+15 563194
## - PropertyType 2 4.4129e+12 1.3706e+15 563263
## - Bathrooms 1 7.6340e+12 1.3739e+15 563318
## - Bedrooms 1 2.8246e+13 1.3945e+15 563656
## - YrBuilt 1 1.2905e+14 1.4953e+15 565240
## - SqFtTotLiving 1 1.3265e+14 1.4989e+15 565294
## - BldgGrade 1 1.9059e+14 1.5568e+15 566155
##
## Step: AIC=563192.1
## AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms +
## BldgGrade + PropertyType + SqFtFinBasement + YrBuilt + YrRenovated +
## NewConstruction
##
## Df Sum of Sq RSS AIC
## - NewConstruction 1 1.0394e+10 1.3662e+15 563190
## - YrRenovated 1 2.5399e+10 1.3663e+15 563190
## - SqFtLot 1 1.0717e+11 1.3663e+15 563192
## <none> 1.3662e+15 563192
## - SqFtFinBasement 1 1.3781e+11 1.3664e+15 563192
## + NbrLivingUnits 1 6.4936e+09 1.3662e+15 563194
## - PropertyType 2 4.4221e+12 1.3706e+15 563261
## - Bathrooms 1 7.7519e+12 1.3740e+15 563318
## - Bedrooms 1 2.8308e+13 1.3945e+15 563655
## - YrBuilt 1 1.3011e+14 1.4963e+15 565254
## - SqFtTotLiving 1 1.3288e+14 1.4991e+15 565296
## - BldgGrade 1 1.9186e+14 1.5581e+15 566172
##
## Step: AIC=563190.2
## AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms +
## BldgGrade + PropertyType + SqFtFinBasement + YrBuilt + YrRenovated
##
## Df Sum of Sq RSS AIC
## - YrRenovated 1 2.5655e+10 1.3663e+15 563189
## - SqFtLot 1 1.1468e+11 1.3664e+15 563190
## <none> 1.3662e+15 563190
## - SqFtFinBasement 1 1.4477e+11 1.3664e+15 563191
## + NewConstruction 1 1.0394e+10 1.3662e+15 563192
## + NbrLivingUnits 1 6.7000e+09 1.3662e+15 563192
## - PropertyType 2 4.5236e+12 1.3708e+15 563261
## - Bathrooms 1 7.7505e+12 1.3740e+15 563317
## - Bedrooms 1 2.8304e+13 1.3945e+15 563653
## - SqFtTotLiving 1 1.3391e+14 1.5001e+15 565310
## - YrBuilt 1 1.3757e+14 1.5038e+15 565365
## - BldgGrade 1 1.9253e+14 1.5588e+15 566179
##
## Step: AIC=563188.7
## AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms +
## BldgGrade + PropertyType + SqFtFinBasement + YrBuilt
##
## Df Sum of Sq RSS AIC
## - SqFtLot 1 1.1399e+11 1.3664e+15 563189
## <none> 1.3663e+15 563189
## - SqFtFinBasement 1 1.4938e+11 1.3664e+15 563189
## + YrRenovated 1 2.5655e+10 1.3662e+15 563190
## + NewConstruction 1 1.0651e+10 1.3663e+15 563190
## + NbrLivingUnits 1 7.2737e+09 1.3663e+15 563191
## - PropertyType 2 4.5012e+12 1.3708e+15 563259
## - Bathrooms 1 7.7815e+12 1.3740e+15 563316
## - Bedrooms 1 2.8286e+13 1.3945e+15 563652
## - SqFtTotLiving 1 1.3389e+14 1.5002e+15 565308
## - YrBuilt 1 1.5088e+14 1.5171e+15 565563
## - BldgGrade 1 1.9253e+14 1.5588e+15 566178
##
## Step: AIC=563188.5
## AdjSalePrice ~ SqFtTotLiving + Bathrooms + Bedrooms + BldgGrade +
## PropertyType + SqFtFinBasement + YrBuilt
##
## Df Sum of Sq RSS AIC
## <none> 1.3664e+15 563189
## + SqFtLot 1 1.1399e+11 1.3663e+15 563189
## - SqFtFinBasement 1 1.4058e+11 1.3665e+15 563189
## + YrRenovated 1 2.4965e+10 1.3664e+15 563190
## + NewConstruction 1 1.8206e+10 1.3664e+15 563190
## + NbrLivingUnits 1 8.1478e+09 1.3664e+15 563190
## - PropertyType 2 4.4353e+12 1.3708e+15 563258
## - Bathrooms 1 7.7135e+12 1.3741e+15 563314
## - Bedrooms 1 2.8588e+13 1.3950e+15 563656
## - SqFtTotLiving 1 1.3750e+14 1.5039e+15 565362
## - YrBuilt 1 1.5077e+14 1.5171e+15 565561
## - BldgGrade 1 1.9243e+14 1.5588e+15 566176step_lm##
## Call:
## lm(formula = AdjSalePrice ~ SqFtTotLiving + Bathrooms + Bedrooms +
## BldgGrade + PropertyType + SqFtFinBasement + YrBuilt, data = house,
## na.action = na.omit)
##
## Coefficients:
## (Intercept) SqFtTotLiving
## 6.178e+06 1.993e+02
## Bathrooms Bedrooms
## 4.240e+04 -5.197e+04
## BldgGrade PropertyTypeSingle Family
## 1.372e+05 2.285e+04
## PropertyTypeTownhouse SqFtFinBasement
## 8.438e+04 7.032e+00
## YrBuilt
## -3.565e+03lm(AdjSalePrice ~ Bedrooms, data=house)##
## Call:
## lm(formula = AdjSalePrice ~ Bedrooms, data = house)
##
## Coefficients:
## (Intercept) Bedrooms
## 117436 132968함수 실행 결과, SqFtLot, NbrLivingUnits, YrRenovated, NewConstruction라는 변수가 삭제된 모델을 선택하였음
더 단순한 선택 방법에는 (1)전진 선택 (2)후진 선택이 있다.
벌점회귀는 개념적으로는 AIC와 같지만, 여기서는 계수 크기를 감소시키거나 경우에 따라 거의 0으로 만들어 벌점을 적용한다. (제거가 아니라 수를 줄이는 것 / 능형회귀와 라소 참고)
회귀에서의 요인 변수
범주형 변수로 불리는 요인변수 factor variable은 개수가 제한된 이산값을 뜻한다.
회귀분석에는 수치 입력이 필요하기 때문에, 요인변수를 이진 가변수들의 집합(dummy variable)로 변환하는 작업이 필요하다.
head(house[, 'PropertyType'])## [1] Multiplex Single Family Single Family Single Family Single Family
## [6] Townhouse
## Levels: Multiplex Single Family Townhouse이 때 가능한 값은 Multiplex, Single Family, Townhouse
이진변수들로 변환하여 나타낸다.
prop_type_dummies <- model.matrix(~PropertyType -1, data=house)
head(prop_type_dummies)## PropertyTypeMultiplex PropertyTypeSingle Family PropertyTypeTownhouse
## 1 1 0 0
## 2 0 1 0
## 3 0 1 0
## 4 0 1 0
## 5 0 1 0
## 6 0 0 1이런 식의 표현법을 머신러닝 커뮤니티에서는 원-핫 인코딩이라 칭한다.
회귀분석에서 P개의 개별 수준을 갖는 요인변수는 보통 P-1개의 열을 갖는 행렬로 표시된다.
R의 기본 표현법은 첫번째 요인 수준을 기준으로 하고, 나머지 수준을 이 기준에 상대적인 것으로 해석한다.
lm(formula = AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms +
BldgGrade + PropertyType, data = house)##
## Call:
## lm(formula = AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms +
## Bedrooms + BldgGrade + PropertyType, data = house)
##
## Coefficients:
## (Intercept) SqFtTotLiving
## -4.469e+05 2.234e+02
## SqFtLot Bathrooms
## -7.041e-02 -1.597e+04
## Bedrooms BldgGrade
## -5.090e+04 1.094e+05
## PropertyTypeSingle Family PropertyTypeTownhouse
## -8.469e+04 -1.151e+05R 회귀분석 결과에서 두 PropertyType에 대응하는 두 계수 PropertyTypeSingle Family와 PropertyTypeTownhouse를 볼 수 있음.
PropertyTypeSingle Family == 0 & PropertyTypeTownhouse == 0 일 때, Multiplex는 암묵적으로 정의되기 때문에 표현하지 않는 것이다.
다수의 수준을 갖는 요인 변수들
table(house$ZipCode)##
## 9800 89118 98001 98002 98003 98004 98005 98006 98007 98008 98010 98011 98014
## 1 1 358 180 241 293 133 460 112 291 56 163 85
## 98019 98022 98023 98024 98027 98028 98029 98030 98031 98032 98033 98034 98038
## 242 188 455 31 366 252 475 263 308 121 517 575 788
## 98039 98040 98042 98043 98045 98047 98050 98051 98052 98053 98055 98056 98057
## 47 244 641 1 222 48 7 32 614 499 332 402 4
## 98058 98059 98065 98068 98070 98072 98074 98075 98077 98092 98102 98103 98105
## 420 513 430 1 89 245 502 388 204 289 106 671 313
## 98106 98107 98108 98109 98112 98113 98115 98116 98117 98118 98119 98122 98125
## 361 296 155 149 357 1 620 364 619 492 260 380 409
## 98126 98133 98136 98144 98146 98148 98155 98166 98168 98177 98178 98188 98198
## 473 465 310 332 287 40 358 193 332 216 266 101 225
## 98199 98224 98288 98354
## 393 3 4 9해당 변수의 경우 주택 가격에 대한 위치의 효과를 볼 수 있는 중요한 변수이다.
하지만, 모든 수준을 포함하려면 자유도 81에 해당하는 81개의 계수가 필요하여 사실상 불가능
-> 대안으로 매매가격과 같은 다른 변수에 따라 우편번호를 그룹으로 묶거나, 초기 모델의 잔차를 사용하여 우편번호 그룹을 만드는 것이 좋다.
zip_groups <- house %>%
mutate(resid = residuals(house_lm)) %>%
group_by(ZipCode) %>%
summarize(med_resid = median(resid),
cnt = n()) %>%
# sort the zip codes by the median residual
arrange(med_resid) %>%
mutate(cum_cnt = cumsum(cnt),
ZipGroup = factor(ntile(cum_cnt, 5)))house <- house %>%
left_join(select(zip_groups, ZipCode, ZipGroup), by='ZipCode')회귀방정식 해석
예측변수 간 상관
다중회귀분석에서 예측변수는 종종 서로 상관성이 있다.
step_lm$coefficients## (Intercept) SqFtTotLiving Bathrooms
## 6.177658e+06 1.992747e+02 4.240306e+04
## Bedrooms BldgGrade PropertyTypeSingle Family
## -5.197477e+04 1.371811e+05 2.285488e+04
## PropertyTypeTownhouse SqFtFinBasement YrBuilt
## 8.437893e+04 7.032179e+00 -3.564935e+03여기서 침실 개수를 뜻하는 Bedrooms 변수가 음수로 나타남.
이는 논리적으로 맞지 않는 결과. 예측변수들이 서로 연관되어 있기 때문에 위와 같은 값이 등장한다.
-> 상호 연관된 예측 변수들을 제거 후에 다시 알아보자.
update(step_lm, .~. -SqFtTotLiving - SqFtFinBasement - Bathrooms)##
## Call:
## lm(formula = AdjSalePrice ~ Bedrooms + BldgGrade + PropertyType +
## YrBuilt, data = house, na.action = na.omit)
##
## Coefficients:
## (Intercept) Bedrooms
## 4913025 27136
## BldgGrade PropertyTypeSingle Family
## 249019 -19932
## PropertyTypeTownhouse YrBuilt
## -47424 -3211update 함수를 이용하여 변수 추가, 제거가 가능하다
현재 Bedrooms 변수는 양수로 변환된 것을 확인할 수 있다.
다중공선성
교란변수는 회귀방정식에서 중요한 변수가 포함되지 못해서 생기는 누락의 문제
방정식 계수에 대한 순진한 해석은 잘못된 결론으로 이어질 수 있음
교란변수
교란변수는 회귀방정식에 중요한 변수가 포함되지 못해서 생기는 누락의 문제.
이 경우, 방정식 계수에 대한 순진한 해석은 잘못된 결론으로 이어질 수 있다.
예를 들어, house_lm 모형에서 SqFtLot, Bathrooms, Bedrooms는 모두 음수였다.
여기에 변수 ZipGroup을 포함해보자.
lm(AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms + BldgGrade + PropertyType + ZipGroup, data=house, na.action = na.omit)##
## Call:
## lm(formula = AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms +
## Bedrooms + BldgGrade + PropertyType + ZipGroup, data = house,
## na.action = na.omit)
##
## Coefficients:
## (Intercept) SqFtTotLiving
## -6.709e+05 2.112e+02
## SqFtLot Bathrooms
## 4.692e-01 5.537e+03
## Bedrooms BldgGrade
## -4.139e+04 9.893e+04
## PropertyTypeSingle Family PropertyTypeTownhouse
## 2.113e+04 -7.741e+04
## ZipGroup2 ZipGroup3
## 5.169e+04 1.142e+05
## ZipGroup4 ZipGroup5
## 1.783e+05 3.391e+05ZipGroup이 분명 중요한 변수라는 사실을 알 수 있다.
가장 비싼 우편번호 그룹의 주택 가격이 약 340,000달러나 더 높다.
이제 SqFtLot과 Bathrooms의 계수는 양수이고, 욕실을 하나 추가하면 판매 가격이 5,537달러 정도 증가한다.
Bedrooms의 경우, 살기 좋은 지역에서는 욕실 수가 같은 주택의 경우 작은 침실이 여러 개 있으면 오히려 가격이 낮아지기 때문에 어느 정도 맞는 결과로 나온다.
상호작용과 주효과
통계학자는 주효과(독립변수)와 주효과 사이의 상호작용을 구별하길 좋아한다.
주효과는 회귀방정식에서 종종 예측변수라고 불린다.
여기에는 예측변수와 응답변수 간의 관계가 다른 예측변수들에 대해 독립적이라는 암묵적인 가정이 있다.
예를 들어, 해당 모델은 ZipCode를 비롯한 여러 변수를 주효과로 포함한다.
부동산에서는 위치가 가장 중요하기 때문에,
주택 크기와 매매 가격 간의 관계가 위치에 달려 있다고 가정하는 것은 자연스러운 일이다.
연산자 *를 사용하여 모델에 변수의 상호작용을 포함시킬 수 있다.
lm(AdjSalePrice ~ SqFtTotLiving*ZipGroup + SqFtLot + Bathrooms + Bedrooms + BldgGrade + PropertyType, data = house, na.action = na.omit)##
## Call:
## lm(formula = AdjSalePrice ~ SqFtTotLiving * ZipGroup + SqFtLot +
## Bathrooms + Bedrooms + BldgGrade + PropertyType, data = house,
## na.action = na.omit)
##
## Coefficients:
## (Intercept) SqFtTotLiving
## -4.919e+05 1.176e+02
## ZipGroup2 ZipGroup3
## -1.342e+04 2.254e+04
## ZipGroup4 ZipGroup5
## 1.776e+04 -1.555e+05
## SqFtLot Bathrooms
## 7.176e-01 -5.130e+03
## Bedrooms BldgGrade
## -4.181e+04 1.053e+05
## PropertyTypeSingle Family PropertyTypeTownhouse
## 1.603e+04 -5.629e+04
## SqFtTotLiving:ZipGroup2 SqFtTotLiving:ZipGroup3
## 3.165e+01 3.893e+01
## SqFtTotLiving:ZipGroup4 SqFtTotLiving:ZipGroup5
## 7.051e+01 2.298e+02이 결과를 볼 때, 주택의 위치와 크기는 강한 상호작용이 있는 것으로 보인다.
가격대가 가장 낮은 ZipGroup에 집에 대한 기울기 (SqFtTotLiving)은 제곱당 118달러.
가장 비싼 ZipGroup에 집에 대한 기울기 (SqFtTotLiving:ZipGroup5)은 제곱당 118 + 230 = 348달러.
다시 말해, 가격대가 가장 비싼 지역에서는 주택의 크기가 1제곱 피트 늘어날 때마다
가장 낮은 지역에 비해 예상 매매가가 거의 3배가 차이난다.
가정 검정: 회귀 진단
특잇값
킹 카운티 주택 매매 데이터에서 우편번호가 98105인 데이터만 가지고 회귀모형 만들기
house_98105 <- house[house$ZipCode == 98105, ]
lm_98105 <- lm(AdjSalePrice ~ SqFtTotLiving + SqFtLot + Bathrooms + Bedrooms + BldgGrade, data=house_98105)다음과 같이 rstandard 함수 이용해 표준화잔차를 구하고, order 함수로 가장 작은 잔차 위치 구할 수 있음
sresid <- rstandard(lm_98105)
idx <- order(sresid)
sresid[idx[1]]## 20431
## -4.326732결과를 보면 표준오차의 4배 이상이나 회귀식과 차이를 보이는 데, 이에 해당하는 추정치는 757,753달러다. 이 특잇값에 해당하는 레코드는 다음과 같다.
house_98105[idx[1], c('AdjSalePrice', 'SqFtTotLiving', 'SqFtLot', 'Bathrooms', 'Bedrooms', 'BldgGrade')]## AdjSalePrice SqFtTotLiving SqFtLot Bathrooms Bedrooms BldgGrade
## 20431 119748 2900 7276 3 6 7이 경우는 레코드가 조금 이상한 것을 알 수 있다. 이 평수에는 적어도 119,748 달러보다는 가격이 비싸야 정상이다. 이와 같이 비정상적인 판매 데이터가 있을 수 있고, 이는 회귀에 포함해서는 안된다.
영향값
회귀모형에서 제외됐을 때 모델에 중요한 변화를 가져오는 값을 주영향관측값 infulential observation이라고 한다.
seed <- 11
set.seed(seed)
x <- rnorm(25)
y <- -x/5 + rnorm(25)
x[1] <- 8
y[1] <- 8
par(mar=c(3,3,0,0)+.1)
plot(x, y, xlab='', ylab='', pch=16)
model <- lm(y~x)
abline(a=model$coefficients[1], b=model$coefficients[2], col="blue", lwd=3)
model <- lm(y[-1]~x[-1])
abline(a=model$coefficients[1], b=model$coefficients[2], col="red", lwd=3, lty=2)
해당 그래프를 보자. 실선은 모든 데이터를 고려한 회귀선에 해당하며, 점선은 오른쪽 위의 한 점을 제거하였을 때의 회귀선을 보여준다.
이 데이터값은 특잇값은 아니지만, 회귀에 대해 높은 레버리지를 가진 것으로 볼 수 있다.
레버리지를 측정하는 일반적인 척도는 햇 값 hat-value이다.
2(P+1)/n 이상의 값들은 레버리지가 높은 데이터 값을 나타낸다.
또 다른 측정 지표는 쿡의 거리 cook’s distance 이다.
이것은 레버리지와 잔차의 크기를 합쳐서 영향력을 판단한다.
쿡의 거리가 4/(n-P-1)보다 크면 영향력이 높다고 보는 편이다.
영향력그림 또는 거품그림 bubble plot은 표준화잔차, 햇 값, 쿡의 거리를 모두 한 그림에 표현한다.
std_resid <- rstandard(lm_98105)
cooks_D <- cooks.distance(lm_98105)
hat_values <- hatvalues(lm_98105)
plot(hat_values, std_resid, cex=10*sqrt(cooks_D))
abline(h=c(-2.5, 2.5), lty=2)
햇 값은 x축, 전차 정보는 y축에 위치하며, 쿡의 거리에 해당하는 값은 원의 크기로 나타낸다.
데이터의 수가 클 경우, 어떤 한 값이 회귀방정식에 엄청난 변화를 가져오기란 쉽지 않다/
물론 이상 검출이 목적이라면 높은 영향력의 값들을 찾는 것이 큰 도움이 된다.
이분산성 / 비정규성 / 오차 간 상관
통계학자들은 잔차 분포에 상당한 주의를 기울인다.
보통최소제곱 추정은 다양한 분포 가정하에서 편향성도 없고 경우에 따라 ’최적’이라 할 수 있는 추정을 제공한다.
즉, 대부분의 문제에서 데이터 과학자라면 잔차 분포에 너무 많은 신경을 쓸 필요는 없다.
잔차의 분포는 주로 공식적인 통계적 추론의 유효성과 관련이 있으므로, 예측 정확도를 중요하기 생각하는 데이터 과학자들에게는 별로 중요하지 않다.
형식적 추론이 완전히 유효하려면 잔차는 1) 동일한 분산을 가지며 2) 정규분포를 따르고 3) 서로 독립이라는 가정이 필요하다.
데이터 과학자가 신경 쓰는 것 한 가지는, 잔차에 대한 가정을 기반으로 예상 값에 대한 신뢰구간을 계산하는 방법이다.
먼저 이분산성은 다양한 범위의 예측값에 따라 잔차의 분산이 일정하지 않은 것을 의미.
어떤 일부분에서 오차가 다른 데보다 훨씬 크게 나타나는 것을 말한다.
df <- data.frame(
resid = residuals(lm_98105),
pred = predict(lm_98105))
ggplot(df, aes(pred, abs(resid))) +
geom_point() +
geom_smooth()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
# histogram
par(mar=c(4,4,0,0)+.1)
hist(std_resid, main='')
주택의 가격이 높아질 수록 잔차의 분산이 높아지는 경향이 있다.
이를 통해, 회귀모형 lm_98105는 이분산성 오차를 갖고 있다고 볼 수 있다.
통계학에서는 오차가 독립적이라는 가정을 점검하기 위해,
시계열 데이터를 다루는 회귀분석에서 유의미한 자기상관이 있는지를 탐지하기 위해 더빈-왓슨 통계량을 사용한다.
하지만, 데이터 과학에서는 ’예측 정확도’가 가장 중요하기 때문에 이분산성에 대한 검토는 그다음으로 이뤄질 수 있다.
편잔차그림과 비선형성
편잔차그림은 예측 모델이 예측변수와 결과변수 간의 관계를 얼마나 잘 설명하는 지 시각화하는 방법
하나의 예측변수와 응답변수 사이의 관계를 모든 다른 예측 변수로부터 분리하는 것이 기본 개념
terms <- predict(lm_98105, type='terms')
partial_resid <- resid(lm_98105) + termsdf <- data.frame(SqFtTotLiving = house_98105[, 'SqFtTotLiving'],
Terms = terms[, 'SqFtTotLiving'],
PartialResid = partial_resid[, 'SqFtTotLiving'])
ggplot(df, aes(SqFtTotLiving, PartialResid)) +
geom_point(shape=1) + scale_shape(solid = FALSE) +
geom_smooth(linetype=2) +
geom_line(aes(SqFtTotLiving, Terms))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
편잔차 그래프는 SqFtTotLiving 변수가 주택 가겨에 얼마나 영향을 미치는지 보여준다.
회귀선에 따르면, 1000제곱피트보다 작은 평수에서는 가격을 원래보다 낮게,
2000 ~ 3000제곱피트에서는 원래보다 높게 추정하고 있다. (4000피트 이상은 데이터가 적어서 결론 내리기 힘듬)
다항회귀와 스플라인 회귀
응답변수와 예측변수 간의 관계가 반드시 선형일 필요가 없다.
비선형 효과를 회귀분석에 담기 위해 회귀모형을 확장하는 몇 가지 방법이 있다.
다항식
다항회귀란 회귀식에 다항 항을 포함한 것을 말한다.
lm(AdjSalePrice ~ poly(SqFtTotLiving, 2) + SqFtLot + BldgGrade + Bathrooms + Bedrooms, data=house_98105)##
## Call:
## lm(formula = AdjSalePrice ~ poly(SqFtTotLiving, 2) + SqFtLot +
## BldgGrade + Bathrooms + Bedrooms, data = house_98105)
##
## Coefficients:
## (Intercept) poly(SqFtTotLiving, 2)1 poly(SqFtTotLiving, 2)2
## -402530.47 3271519.49 776934.02
## SqFtLot BldgGrade Bathrooms
## 32.56 135717.06 -1435.12
## Bedrooms
## -9191.94결과를 통해 SqFtTotLiving에 대한 두 가지 계수가 있다는 것을 볼 수 있다 (일차항과 이차항).
편잔차그림을 비교해보면, 회귀 결과가 선형회귀일 때보다 평활곡선에 더 가까운 것을 확인할 수 있다.
lm_poly <- lm(AdjSalePrice ~ poly(SqFtTotLiving, 2) + SqFtLot +
BldgGrade + Bathrooms + Bedrooms,
data=house_98105)
terms <- predict(lm_poly, type='terms')
partial_resid <- resid(lm_poly) + terms
df <- data.frame(SqFtTotLiving = house_98105[, 'SqFtTotLiving'],
Terms = terms[, 1],
PartialResid = partial_resid[, 1])
ggplot(df, aes(SqFtTotLiving, PartialResid)) +
geom_point(shape=1) + scale_shape(solid = FALSE) +
geom_smooth(linetype=2) +
geom_line(aes(SqFtTotLiving, Terms))+
theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
스플라인
다항회귀는 비선형 관계에 대해 어느 정도의 곡률을 담아낼 수 있다.
하지만 3, 4차 다항식의 경우 회귀방정식에 바람직하지 않은 ’흔들림’을 초래한다.
이 때는 스플라인을 활용한다.
스플라인은 고정된 점들 사이를 부드럽게 보간하는 방법이다.
R 패키지 splines는 회귀모형에서 b-스플라인 b-spline 항을 만드는 데 사용할 수 있는 bs 함수를 포함하고 있다.
library(splines)
knots <- quantile(house_98105$SqFtTotLiving, p=c(.25, .5, .75))
lm_spline <- lm(AdjSalePrice ~ bs(SqFtTotLiving, knots=knots, degree=3) + SqFtLot +
Bathrooms + Bedrooms + BldgGrade, data=house_98105)terms <- predict(lm_spline, type='terms')
partial_resid <- resid(lm_spline) + terms
df <- data.frame(SqFtTotLiving = house_98105[, 'SqFtTotLiving'],
Terms = terms[, 1],
PartialResid = partial_resid[, 1])
ggplot(df, aes(SqFtTotLiving, PartialResid)) +
geom_point(shape=1) + scale_shape(solid = FALSE) +
geom_smooth(linetype=2) +
geom_line(aes(SqFtTotLiving, Terms))+
theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
선형 항에서는 계수가 변수에 대한 직접적인 의미를 갖는 반면,
스플라인 항의 계수는 해석하기 어렵다.
대신 적합도를 확인하기 위해 시각화 방법을 사용하는 것이 더 유용하다.
다항회귀 모델과는 달리 스플라인 모델은 좀 더 매끄럽게 매칭되어 유연성이 뛰어난 것을 볼 수 있다.
일반화가법모형
일반화가법모형(GAM)은 스플라인 회귀를 자동으로 찾는 기술이다.
lm_gam <- gam(AdjSalePrice ~ s(SqFtTotLiving) + SqFtLot + Bathrooms + Bedrooms + BldgGrade, data=house_98105)terms <- predict.gam(lm_gam, type='terms')
partial_resid <- resid(lm_gam) + terms
df <- data.frame(SqFtTotLiving = house_98105[, 'SqFtTotLiving'],
Terms = terms[, 5],
PartialResid = partial_resid[, 5])
ggplot(df, aes(SqFtTotLiving, PartialResid)) +
geom_point(shape=1) + scale_shape(solid = FALSE) +
geom_smooth(linetype=2) +
geom_line(aes(SqFtTotLiving, Terms))+
theme_bw()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'