ADP 실기 준비 - 머신러닝 (1) 회귀분석 패키지/함수 정리

1. 단순회귀분석, 다항회귀분석 (polynomial analysis) , 다중회귀분석 

# lm은 선형 모형을 적합시키는 데 사용됩니다. 회귀 분석, 분산의 단일 층 분석 및 
# 공분산 분석을 수행하는 데 사용할 수 있습니다
# (aov가 이러한 인터페이스를 더 편리하게 제공할 수도 있음).

result <- lm(formula, # 식 작성 - 종속변수 ~ 사용하려는 독립변수(들) 
   		data, # R에서 사용하려는 데이터 명시 
  		subset, # 필요시 조건을 부여하여 나눔 ex. subset = 컬럼 > 컬럼의 평균값 
		weights, na.action,
  		method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
 		singular.ok = TRUE, contrasts = NULL, offset, …)
   
 # 고려사항 formula에 넣을 때 심볼의 의미는 R을 이용한 통계데이터분석 책 p.152 참고 
 # formula에 제곱, 세제곱을 넣으며, 분석해볼 수 있음 
 
 # 베타회귀계수(종속변수에 대한 독립변수의 영향력을 비교할 때 사용)
library(QuantPsyc)
lm.beta(result) 
   
 # 결과 추출 함수 
 # 1) anova() : 분산분석표 
 # 2) coef() : 회귀계수
 # 3) confint() : 회귀계수에 대한 신뢰구간, 기본적으로 95% 신뢰구간이며, 필요시 level = 로 변환 가능
 # 4) fitted() : 회귀식에 대한 예측값 
 # 5) resid() : 잔차 
 # 6) summary(): 주요 분석정보 
 
 # 예측을 위한 함수 
 predict(result, # lm함수의 결과로 지정된 것 
 		newdata = newdata # 예측을 위해 넣는 데이터(일종의 테스트 데이터) 
        interval = "confidence" # 신뢰구간까지 함께 예측이 필요한 경우 
        )
        
 # 결과정리 : 논문형식으로 통계분석표를 정리하는 것 
 library(stargazer)
 stargazer(result, type= "text", no.space = T)
 
 # 베타 회귀계수를 구하는 패키지
 # 베타 회귀계수 : 종속변수에 영향을 주는 다양한 독립변수가 있을 때, 각각의 독립변수의 
 # 측정단위가 다르기 때문에 회귀계수로 각 변수의 설명력을 파악하기는 어려움
 # 이때 사용할 수 있는 방법이 표준화 계수(standardized coefficient)인 베타임 
 
 library(QuantPsyc) 
 lm.beta(result) 
# 해석: 절대값의 크기가 크면 클수록 영향력이 큰 것으로 판단함.

2. 회귀분석의 가정과 진단, 다중공산성 진단, 모델의 수정, 모델 선택, 더미변수 활용 

 

회귀분석 가정과 진단

1. 선형성(linearity) : 종속변수와 독립변수 간의 관계는 선형

2. 정규성(normality) : 독립변수값에 대해 대응되는 종속변수값의 분포는 정규분포

3. 등분산성(homoscedasticity, equality of variance) : 독립변수값에 대해 대응되는 종속변수값들의 분포는 모두 동일한 분산을 갖음 

4. 독립성(independence) : 모든 관측값은 서로 독립임. 하나의 관측값은 다른 관측값에 영향을 미치지 않음 

 

# 회귀분석의 가정 확인
plot(result)

# 다중공산성 확인
library(car)
vif(result) 
# 정확한 내용은 R을 이용한 통계데이터분석 p.165 ~ 167

# 회귀모델 수정 
# 종속변수의 변환 
library(car)
summary(powerTransform(data$종속변수 컬럼명))

# 독립변수의 변환
# 결과의 해석 : 통계적으로 변환의 결과가 유의미한 경우 변환 진행 
boxTidwell(종속변수 ~ 독립변수 1 + 독립변수 2 + ... + 독립변수 n , data = 데이터)

# 회귀모델의 선택 
# 1) ANOVA 함수
# 회귀모델 내 추가적인 변수가 예측력의 향상에 추가적인 영향이 있는지를 확인함 
result_1 <- lm(y~x1 + x2 , data = 데이터) 
result_2 <- lm(y~x1 + x2 + x3 + x4 , data = 데이터) 
anova(result_1, result_2) 
# anova의 결과 F값이 유의수준 0.05에서 통계적으로 유의미하지 않은 경우 x3, x4를 추가하지 않음 
# anova의 단점 : 중첩된 모델간의 비교만 가능함 

# 2) AIC
# anova 단점의 해소 : 두 모델 간의 AIC 값을 비교하여 AIC값이 낮을수록 좋은 모델로 판단함 
AIC(result_1, result_2) 

# 3) Step 방법

# 4) regsubsets 방법
# 가지고 있는 독립변수를 기준으로 수행할 수 있는 모든 경우의 수를 해봄
library(leaps) 
regsub_result <- regsubsets(x = 종속변수 ~ 독립변수 1 + 독립변수 2 + ... , data= 데이터, 
                           nbest = 독립변수의 각 서브셋별 크기) 
# 서브셋별 몇 개의 최적의 모델을 산출할 것인지 
# 결과 산출
library(RColorBrewer) 
plot(regsub_result, scale = "adjr2", col = brewer.pal(9, "Pastel1"), main = "제목")

# 더미변수가 있는 경우
# 단, str()을 통해 확인 간 독립변수가 factor화가 되어 있지 않은 경우 
# data$더미변수 <- as.factor(data$더미변수)을 통해 factor화를 먼저 진행 
dummy_result <- lm(종속변수 ~ 독립변수(더미변수), data = 데이터) 
summary(dummy_result)
# 결과값은 자동으로 각각의 더미별 결과가 산출됨
# 이에 대한 해석은 더미화되지 않은 변수의 평균값이 절편값으로 표현됨 

# 참고 : 일원분산분석을 통해서도 각 더미변수간의 차이가 유의미한지 확인할 수 있음 
ano_result <- aov(종속변수 ~ 더미형 독립변수, data = 데이터) 
summary(ano_result)
TukeyHSD(ano_result) 
# 결과값을 통해 각 더미변수값의 차이가 유의미한지 확인할 수 있음. 
# 통계적으로 유의미한 것들만 유의미한 차이가 있는 것으로 판단함

 

3. 매개효과분석과 조절효과 분석 

매개효과 분석(mediation effect analysis) : X 변수가 K 변수에 영향을 미치고, K 변수가 다시 Y 변수에 영향을 미치는 경우(연쇄적인 영향관계를 검증하는 것)에 대한 분석

ex. 소득이 높은 사람이 기대 수명이 높다 → 소득이 높으면 더 나은 의료환경을 기대할 수 있고, 더 나은 의료환경이 기대 수명을 높이게 됨 

이를 수행하는 방법은 바론 & 케니 분석을 통해 진행됨. 바론 & 케니 분석은 통계적 분석 결과를 제공하지 않으므로, 이를 확인하기 위해서는 소벨검정(Sobel test)나 부트스트랩핑(bootstrapping)을 통해 통계적 유의성 검정을 진행함.

# 바론 & 케니 방법론
total <- lm(종속변수 ~ 독립변수1, data = 데이터) 
summary(total) 
# 하나의 독립변수가 종속변수에 미치는 총효과(total effect)를 확인함 

model_1 <- lm(독립변수2 ~ 독립변수1, data = 데이터) 
summary(model_1) 
# 독립변수 2와 독립변수 1이 유의미한 상관관계를 갖는지 확인함
# 만약 유의미한 상관관계가 없는 경우 더 이상 X, 유의미한 상관관계 존재 시 아래 절차 수행

model_2 <- lm(종속변수 ~ 독립변수 1 + 독립변수 2, data = 데이터) 
summary(model_2) 
# 결과를 볼 때 독립변수 1이 더 이상 유의미하지 않다면, 독립변수 2가 독립변수 1을 
# 완전 매개한다고 볼 수 있음. 독립변수 1과 독립변수 2의 간접효과는 
# model_1의 estimate(독립변수 1) * model_2의 estimate(독립변수 2) 

# 통계적 검정 
# 1) 소벨검정
# 통계적 유의성 파악 가능 
library(multilevel)
model_sob <- sobel(pred = 데이터$독립변수1, med = 데이터$독립변수2, out = 데이터$종속변수)
# 결과값 확인
# Mod 1~3을 통한 검증 
model_sob

#간접 효과에 대한 p값 계산 
pnorm(abs(model_sob$z.value), lower.tail =F )* 2

# sobel 검정의 단점 : 간접효과가 정규분포를 따라야 한다는 가정을 충족해야하고, 표본크기가
# 충분히 커야한다는 단점이 있음.
# 이에 따라 최근에는 해당 가정 없이도 사용할 수 있는 bootstraping 기법을 선호함 
library(mediation)
set.seed(1234)
model_1 <- lm(독립변수2~ 독립변수1, data= 데이터) 
model_2 <- lm(종속변수 ~ 독립변수1 + 독립변수2, data = 데이터) 
model_3 <- mediate(model = model.m = model_1, model.y = model_2, treat = "독립변수1", 
mediator = "독립변수2", boot = T, sims = 500)
summary(model_3)

#이에 대한 결과 해석
plot.mediate(model_3, cex = 1.2, col = "royalblue", lwd = w, main = "제목")

조절효과 분석(moderation effect analysis) = 상호작용 효과 (interaction effect) 

종속변수에 대해 독립변수1과 2 각각이 아닌 이들 간의 상호작용 하 어떤 식으로 모델이 생성되는지 확인 

model_1 <- lm(종속변수 ~ 독립변수1+ 독립변수2 + 독립변수1:독립변수2, data = 데이터) 
#결과 확인 
summary(model_1) 
#결과에 대해 그래프로 확인해보기
# effect 함수활용 
library(effect)
m <- round(mean(데이터$독립변수1), 1)
s <- round(sd(데이터$독립변수1) ,1)
plot(effect(term = "독립변수2:독립변수1", mod = model_1, xlevels = list(독립변수1 = c(m-s, m, m+s)),
lines = list(multiline = T, lwd = 2, lty = c(3,2,1), col = c("royalblue", "violet", "maroon"),
main = "제목") 

#plotSlopes 함수 활용
library(rockchalk)
plotSlopes(model = model_1, plotx = "독립변수1", modx = "독립변수2", # 조절효과를 분석하고자 하는 변수
modxVals = "std.dev.",pch = 21, col = rainbow(3), cex = 1 , bg = "dimgray", main = "제목")

조절매개효과분석(moderated mediation effect analysis)는 매개변수에 의해 매개된 두 변수(독립변수, 종속변수) 간 직접적인 혹은 간접적인 관계가 조절변수에 미치는 영향을 검정함 

# 앞서 독립변수2 → 독립변수1 → 종속변수로 이어지는 매개효과가 존재하는지 여부를 확인함 
# 이때 새로운 독립변수 3(1,0으로만 이루어진 변수로 가정)에 따라 매개효과가 달라지는지 확인하고자 함

model_1 <- lm(독립변수1 ~ 독립변수2*독립변수3, data = 데이터) # 이때 독립변수3은 factor형 변수 
model_2 <- lm(종속변수 ~ 독립변수2*독립변수3 + 독립변수1*독립변수3, data= 데이터) 

library(mediation)
set.seed(12) 
model_med1 <- mediate(model.m = model_1, model.y = model_2, covariates = list(독립변수3 = 0),
treat = "독립변수2", mediator = "독립변수1", boot = T, sims = 500)
summary(model_med1)

model_med2 <- mediate(model.m = model_1, model.y = model_2, covariates = list(독립변수3 = 1),
treat = "독립변수2", mediator = "독립변수1", boot = T, sims = 500)
summary(model_med2)

# summary에서 나오는 ACME 결과를 비교함 
# 통계적 유의성도 함께 확인

# 간접효과의 차이가 통계적으로 유의미한지 검정하는 코드
set.seed(123)
model_med3 <- mediate(model.m = model_1, model.y = model_2,
treat = "독립변수2", mediator = "독립변수1", sims = 500)

test.modmed(object = model_3, covariates.1 = list(am = 0), covariates.2 = list(am=1), 
sims = 500)
# 만약 조절변수가 연속형 변수인 경우(본 예시에서는 범주형 변수) 
# covariates.1과 covariates2를 임의로 설정하기

 

4. 패널티 회귀분석 - 릿지(Ridge), 랏소(Lasso), 일래스틱넷(Elasticnet) 회귀분석

기본 목적 : 모델의 단순화를 통해 과적합 피하기, 다중공산성 방지 등을 달성하기 위함 

# 활용하는 패키지와 함수
library(glmnet)
glmnet(x, y, family, alpha = 1, lambda = NULL)
# alpha = 0 이면 릿지, alpha = 1이면 랏쏘 
# 주의할 점 
# 수식 형태로 모델을 지원하지 않아서 별도 지정 필요
# 범주형 변수는 더미변수화하여 넣어야 함 
x <- model.matrix(종속변수 ~ . 훈련 데이터셋)[,-1] # model.matrix 첫 열은 의미가 없어서 제외
y <- 훈련 데이터셋$종속변수

 

1) 릿지 회귀분석 

모델의 설명력에 기여하지 못하는 독립변수의 회귀계수를 0에 가깝게 만듦(L2-norm이라는 패널티항을 통해 회귀모델에 패널티 부여하며, 람다를 통해 패널티의 정도를 조절함/ 아예 0으로 만드는 것은 아니기 때문에 모델 단순화의 효과는 없음) 

주의 사항 : 독립변수의 척도에 영향을 받기 때문에 이에 대해 표준화를 선행하는 것이 좋음 

library(glmnet)
set.seed(1234)
#  최적의 람다를 찾기 위해 cross-validation (최적의 람다 : 어느 수준으로 패널티를 줘야 가장 좋은 모델이 되는지) 
ridge_model_cv <- cv.glment(x = x, y= y, family = "gaussian", alpha = 0) 
plot(ridge_model_cv)

# 최적의 람다 확인 
ridge_model_cv$lambda.min 

# 실제 모델링 
ridge_model <- glmnet(x,y,family = "gaussian", alpha = 0, lambda = ridge_model_cv$lambda.min)

# 모델링에 대한 평가 
# ridge_test_data <- model.matrix를 통한 변환 
ridge_predict <- predict(ridge_model, newx = ridge_test_data) 
library(caret) 
postResample(pred = ridge_predict, obs = ridge_test_data$종속변수)

 

2) 랏소 회귀분석

모델의 설명력에 기여하지 못하는 독립변수의 회귀계수를 0으로 만듦(모델에서 해당 독립변수를 제거하여 모델의 단순화 가능) 

set.seed(1234)
lasso_model_cv <- cv.glmnet(x = x, y= y, family = "gaussian", alpha = 1) 
# 예측 오차를 최소화시켜주는 람다값 확인 
lasso_model_cv$lamnda.min
plot(lasso_model_cv)
# 릿지 모형과의 차이점은 예측변수의 갯수가 상단에서 람다 값에 따라 변경됨  
# cv.glmnet에서 찾고자 하는 최적의 람다의 경우 최소한의 예측변수로 적정 수준의 정확도를 제공하는 모델임
# 이에 따라 예측 오차를 최소화하는 최적의 람다 및 최소 예측 오차 내 1개 sd 내 있으면서, 
# 예측변수의 갯수를 최소화 시키는(즉, 과적합의 위험을 낮추는) 람다를 찾음 
lasso_model_cv$lambda.1se

# 두 가지 버전으로 Ridge와 같이 모델을 평가해보면, lambda.min의 예측도는 더 좋음
# 단, 독립변수의 갯수를 줄여서 간명하면서도 적당한 예측도를 제공해줄 수 있는 건 lambda.1se

3) 엘라스틱넷 회귀분석

L1-norm과 L2-norm을 모두 사용하여 회귀모델에 패널티를 부여함(릿지와 랏소의 혼합) 

엘라스틱넷에서는 위에서 조절한 람다와 alpha값도 조절해야함(L1-norm과 L2-norm의 패널티 비율 결정을 위해) 

library(caret) 
set.seed(1234)
# alpha와 lambda의 조합 각각 10개씩을 생성하여 검증에 사용
elastic_cv <- train(form = 종속변수 ~ . , data = 훈련데이터셋, method = "glmnet",
trControl = trainControl(method = "cv", number = 10), tuneLength = 10) 
# 가장 나은 결과확인
elastic_cv$bestTune

elastic_model <- glmnet(x, y, family= "gaussian", alpha = elastic_cv$bestTune$alpha,
lambda = elastic_cv$bestTune$lambda) 

# 모델 평가 방식은 위와 동일함

모델 비교해보기 

# Ridge 
set.seed(1234)
ridge <- train(종속변수 ~ . , data = 훈련데이터셋, method = "glmnet", 
trControl = trainControl(method = "cv", number= 10), 
tuneGrid = expand.grid(alpha = 0, lambda = lambda))
coef(ridge$finalModel, ridge$bestTune$lambda)
ridge_pred <- predict(ridge, 테스트데이터셋) 
postResample(pred = ridge_pred, obs = 테스트데이터셋$종속변수) 

# Lasso 
set.seed(1234)
lasso <- train(종속변수 ~ . , data = 훈련데이터셋, method = "glmnet", 
trControl = trainControl(method = "cv", number= 10), 
tuneGrid = expand.grid(alpha = 1, lambda = lambda))
coef(lasso$finalModel, lasso$bestTune$lambda)
lasso_pred <- predict(lasso, 테스트데이터셋) 
postResample(pred = lasso_pred, obs = 테스트데이터셋$종속변수) 

# Elasticnet  
set.seed(1234)
elastic <- train(종속변수 ~ . , data = 훈련데이터셋, method = "glmnet", 
trControl = trainControl(method = "cv", number= 10), tuneLength = 10)
coef(elastic$finalModel, elastic$bestTune$lambda)
elastic_pred <- predict(elastic, 테스트데이터셋) 
postResample(pred = elastic_pred, obs = 테스트데이터셋$종속변수) 

# 최종 비교 
models <- list(ridge = ridge, lasso = lasso, elastic = elastic) 
summary(resamples(models), metrics = "RMSE")

# 통계적으로 차이가 유의미한지 확인 
summary(diff(resampels(models), metric = "RMSE")