앤디 필드의 유쾌한 R 통계학 챕터 7 - 회귀분석

 

회귀 분석은 lm함수를 활용해서 간단히 실행할 수 있다.

다만 변수 선택과 이상치를 확인하는 작업을 신경 써서 진행해야 한다.

 

 

library(boot);library(car);library(QuantPsyc)

# 데이터를 불러와 lm(linear model) 함수 사용하기

album1 <- read.delim('datafiles/Album Sales 1.dat', header = TRUE)
albumSales.1 <- lm(sales~adverts, album1, na_action = na.omit) # na_action = na.omit은 결측치를 제외하라는 뜻이다.

summary(albumSales.1) # 함수의 결과 확인하기

# 결정계수 R^2는 상관계수의 제곱이므로 결정계수의 제곱근으로 상관계수를 찾을 수 있다.
sqrt(0.3346)

# Y 절편(intercept)과 각 변수가 유의한지는 t 검정을 통해 판단한다.
# 이 때의 귀무가설은 회귀계수(b0, b1 등)가 0이다. 이며 이를 검증하기 위해 (추정지-0)/표준편차 식을 계산한다.
# 값은 lm 함수를 사용하면 R이 알아서 계산해준다.

# F 통계량은 모형의 평균제곱/잔차의 평균제곱으로 모형이 얼마나 좋은지 평가할 때 사용한다.
# 평균제곱은 제곱합을 자유도로 나눈 것이다.

# 다중회귀 분석하기

album2 <- read.delim('datafiles/Album Sales 2.dat', header = TRUE)
albumSales.2 <- lm(sales~adverts, album2)
albumSales.3 <- lm(sales~adverts + airplay + attract, album2)
# 혹은 albumSales3 <- update(albumSales.2, .~. + airplay + attract)

summary(albumSales.2);summary(albumSales.3)

# 각 변수들의 회귀계수를 표준화하여 직접적으로 효과를 비교할 수 있다.
# 이때 사용하는 lm.beta 함수는 QuantPsyc 패키지에 포함되어 있는데 이 패키지는 현재는 cran에서 제외되어 있다.
# lm.beta 함수의 결과는 각 변수가 1 표준편차 변화했을때 종속변수가 몇 표준편차 변화하는지를 나타낸다.

lm.beta(albumSales.3)

# 각 변수들의 신뢰구간 확인하기

confint(albumSales.3)

# 위계적인 순서의 모델(변수를 추가해가는 등)의 비교에는 anova함수를 사용한다.
# 변화량과 변화의 유의성을 확인할 수 있다.

anova(albumSales.2, albumSales.3)


# 이상치 혹은 영향력이 큰 사례(관측치) 진단하기

album2$residuals <- resid(albumSales.3) # 잔차
album2$standardized.residuals <- rstandard(albumSales.3) # 표준화 잔차 : 차이가 3 초과시 이상치
album2$studentized.residuals <- rstudent(albumSales.3) # 스튜던트화 잔차
album2$cooks.distance <- cooks.distance(albumSales.3) # 쿡 거리 : 1 초과시 영향을 미칠 수 있다
album2$dfbeta <- dfbeta(albumSales.3) # 베타 차이
album2$dffit <- dffits(albumSales.3) # 적합도 차이
album2$leverage <- hatvalues(albumSales.3) # 모자 값 : 평균의 2배 또는 3배 이상인 경우 이상치
album2$covariance.ratios <- covratio(albumSales.3) # 공분산비 : 1에서 모자값의 3배 이상 차이나는 경우 문제

album2 <- round(album2, digits = 3)
head(album2)

# 표준화 잔차 기준으로 이상치 확인하기

album2$large.residual <- album2$standardized.residuals > 2 | album2$standardized.residuals < -2
sum(album2$large.residual) # 개수 확인하기
album2[album2$large.residual, c('sales','airplay','attract','adverts','standardized.residuals')] # 사례 확인하기
# 표준화 잔차에서 이상치였던 사례들의 다른 진단량 확인하기
album2[album2$large.residual, c('cooks.distance','leverage','covariance.ratios','standardized.residuals')]

# dwt 함수로 오차의 독립성 검증하기
# 더빈왓슨 통계량은 2에 가까울수록 좋다
# 영가설: 잔차 사이에 관계가 없다

dwt(albumSales.3)

# vif 함수로 다중공선성 확인하기
# vif가 10 이상이면 다중공선성이 있는 것으로 판단한다.
# 평균 vif가 1과 차이가 크면 과적합한 것으로 판단한다.

vif(albumSales.3)
mean(vif(albumSales.3))

# 시각적으로 잔차 확인하기

plot(albumSales.3)
hist(album2$standardized.residuals)

# 부트스트랩 기법 활용하기

bootReg <- function(formula, data, indices){
  d <- data [indices,]
  fit <- lm(formula, data = d)
  return(coef(fit))
}

bootResults <- boot(statistic = bootReg, formula = sales~adverts+airplay+attract, data = album2, R = 2000)
boot.ci(bootResults, type = 'bca', index = 1)
boot.ci(bootResults, type = 'bca', index = 2)
boot.ci(bootResults, type = 'bca', index = 3)
boot.ci(bootResults, type = 'bca', index = 4)

# 더미 변수 만들기

gfr <- read.delim('datafiles/GlastonburyFestivalRegression.dat', header = TRUE)
head(gfr)
gfr$music <- as.factor(gfr$music) # contrasts는 요인(factor)에만 사용할 수 있다.
contrasts(gfr$music) <- contr.treatment(4, base = 4) # contr.treatment(그룹 수, base = 베이스라인으로 삼을 그룹 번호)
# 혹은 이름 지정을 위해 직접 만든다
# crusty_v_NMA <- c(1, 0, 0, 0)
# indie_v_NMA <- c(0, 1, 0, 0)
# metal_v_NMA <- c(0, 0, 1, 0)
# contrasts(gfr$music) <- cbind(crusty_v_NMA, indie_v_NMA, metal_v_NMA)

gfr$music

# 더미 변수를 포함한 lm 함수

gfModel <- lm(change~music, gfr, na.action = na.omit)
summary(gfModel)

# 각 변수의 평균 변화량 확인하기(회귀계수로도 확인할 수 있다)
round(tapply(gfr$change, gfr$music, mean, na.rm = TRUE), 3)