[R 데이터분석] 15장. 회귀분석

# Chap15_1_LinearRegression & Chap15_2_LogisticRegression 

# 선형회귀분석(LinearRegration)
# 1. 단순회귀분석 : 독립변수(x) : 1 -> 종속변수(y)
# 2. 다중회귀분석 : 독립변수(x) : n -> 종속변수(y)

# 3. 변수선택법

# 4. 다중공선성과 기계학습 

#  5. 선형회귀분석 잔차검정과 모형진단 

 

# 로지스틱 회귀분석(Logistic Regression) 

# 1. 오즈비(odds ratio) 

# 2. 로짓변환

# 3. sigmoid function

#########################
## 1. 단순회귀분석
#########################
# 독립변수와 종속변수가 1개인 경우 

 

### 사용된 파일은 https://joyfuls.tistory.com/4 에서 다운 받으실 수 있습니다.
# 1) 데이터셋 가져오기 
setwd("c:/2_Rwork/Part-IV")

product <- read.csv("product.csv")
str(product)
# 'data.frame': 264 obs. of  3 variables:
#$ 제품_친밀도: int  3 3 4 2 2 3 4 2 3 4 ...
#$ 제품_적절성: int  4 3 4 2 2 3 4 2 2 2 ... : x
#$ 제품_만족도: int  3 2 4 2 2 3 4 2 3 3 ... : y

# 2) 변수 선택 
x <- product$제품_적절성 # 독립변수 : 영향을 주는 변수 
y <- product$제품_만족도 # 종속변수 : 영향을 받는 변수 
df <- data.frame(x, y)

# 3) 회귀모델 생성 
product_lm <- lm(y ~ x, data = df)
product_lm
# Coefficients: 회귀계수 
#(Intercept) : 절편    x : 기울기   
#      0.7789(y절편)   0.7393(x 기울기)


# 회귀방정식 : Y = aX + b(a:기울기, b:절편)
head(df)
dim(df) # 264   2

x <- 4 # 입력 
y <- 3 # 출력(정답:관측치)
a <- 0.7393
b <- 0.7789

# 모델 예측치 
y_pred <- a*x + b # 회귀방정식 
y_pred # 예측치 : 3.7361

# 오차(error)
err <- y - y_pred # 관측치 - 예측치 
err # -0.7361
abs(err) # 0.7361

names(product_lm)
# "coefficients" : 계수 
# "residuals" : 오차(잔차)
# "fitted.values" : 예측치 

product_lm$coefficients
product_lm$fitted.values[1] # 3.735963 
product_lm$residuals[1] # -0.735963

# 4) 회귀모델 결과 분석 
summary(product_lm)

# <분석 절차>
# 1. F검정 통계량 : 모델의 통계적 유의검정(p-value < 0.05) 
# 2. 모델의 설명력 : Adjusted R-squared: 0.59(예측력 60%)
# 3. x변수의 통계적 유의성 검정 : p-value < 0.05

# R^2 -> R
sqrt(0.5865) # 0.7658329

0.7658329 ** 2 # 0.5865

# 5) 회귀선 시각화 
plot(y ~ x, data = df)
abline(product_lm, col='red')


#########################
## 2. 다중회귀분석
#########################
# 독립변수(n)와 종속변수 1개인 경우 

# 1) dataset 가져오기 
install.packages("car")
library(car)

Prestige
str(Prestige)
# 'data.frame': 102 obs. of  6 variables:
# 102개 직업군에 대한 교육수준, 수입, 여성비율, 평판,직원수  
# $ education : 다소 양의 상관성 
# $ income : y
# $ women : 다소 음의 상관성 
# $ prestige : 양의 상관성 
# $ census

head(Prestige)

# 2) 상관분석 
newdata <- Prestige[c(1,3:4,2)]
head(newdata)

cor(newdata)
# income    0.57758023 -0.44105927  0.7149057  1.0000000


# 3) 회귀모델 
#lm(formula = y ~ x1 + x2 + xn , data = )
#pre_lm <- lm(formula = income ~ education + women + prestige, data = newdata)
pre_lm <- lm(formula = income ~ ., data = newdata)
pre_lm # Coefficients:
# Coefficients:
#(Intercept)          education      women     prestige  
#   -253.8(절편)        177.2        -50.9        141.4

# 회귀방정식 : Y = a1*x1 + a2*x2 + a3*x3 + b
head(newdata)
education <- 13
women <- 11.6
prestige <- 69
Y = 177.2*education + -50.9*women + 141.4*prestige -253.8
Y # 11215.96


# 4) 회귀모델 분석 : 분석절차 
summary(pre_lm)
# F-statistic: p-value: < 2.2e-16
# Adjusted R-squared:  0.6323 
# education : 영향 없음 
# women, prestige : 강한 영향 있음 


#########################
## 3. 변수 선택법
#########################
# - 최적의 모델을 위한 x변수 선택법(3가지) 

library(MASS)

newdata2 <- Prestige[1:5]
head(newdata2)
# x(4) -> y
model2 <- lm(income ~ ., data = newdata2)

step <- stepAIC(model2, direction = "both")
# Step:  AIC=1604.96
# income ~ women + prestige

# 수정된 model 
model3 <- lm(income ~ women + prestige, data = newdata2)

# x(2) -> y
mean(model3$residuals) # 오차의 평균 : -7.53319e-15
# x(3) -> y
mean(pre_lm$residuals) # 오차의 평균 : 1.704083e-14


##############################
## 4. 다중공선성과 기계학습 
##############################

# - 다중공선성 : 독립변수 간의 강한 상관관계로 인해서 발생하는 문제
# - 기계학습 : 훈련셋과 검정셋을 이용하여 모델을 학습하고, 평가하는 방법 

# 1) model vs 다중공선성 문제 확인 
iris
str(iris)

# model 생성 
iris_model <- lm(Sepal.Length ~ Sepal.Width+Petal.Length+Petal.Width,
                  data = iris)

# 다중공선성 문제 확인 
library(car)
vif(iris_model)
#Sepal.Width Petal.Length  Petal.Width 
#  1.270815    15.097572    14.234335

sqrt(vif(iris_model)) > 2
# Sepal.Width Petal.Length  Petal.Width 
#      FALSE         TRUE         TRUE

cor(iris[-5])


# 2) 훈련셋(7)과 검정셋(3) 분류 : 홀드아웃 방식 
row_num <- sample(x=nrow(iris), size = nrow(iris)*0.7, replace = FALSE)
row_num

trainset <- iris[row_num, ] # 훈련셋 
testset <- iris[-row_num, ] # 검정셋 
dim(trainset) # 105   5
dim(testset) # 45  5

# 3) model 학습 : trainset
iris_lm <- lm(Sepal.Length ~ Sepal.Width+Petal.Length, data = trainset)
iris_lm

# 4) model 결과 분석 
summary(iris_lm)

# 5) model 평가 : 상관계수, MSE
y_pred <- predict(iris_lm, testset) # 예측치 
y_pred

# 관측치(정답) vs 예측치
testset$Sepal.Length[1] # 5.1
y_pred[1] # 5.02295

real_value <- testset$Sepal.Length

# (1) 상관계수 
cor(real_value, y_pred) # 0.905183

# (2) MSE(Mean Square Error)
mean((real_value - y_pred)^2) # 0.1186105


##########################################
##  5. 선형회귀분석 잔차검정과 모형진단
##########################################

# 1. 변수 모델링  
# 2. 회귀모델 생성 
# 3. 모형의 잔차(error)검정 
#   1) 잔차의 등분산성 검정
#   2) 잔차의 정규성 검정 
#   3) 잔차의 독립성(자기상관) 검정 
# 4. 다중공선성 검사 
# 5. 회귀모델 생성/ 평가 


names(iris)

# 1. 변수 모델링 : y:Sepal.Length <- x:Sepal.Width,Petal.Length,Petal.Width
formula = Sepal.Length ~ Sepal.Width + Petal.Length + Petal.Width

# 2. 회귀모델 생성 
model <- lm(formula = formula,  data=iris)
model
names(model)

# 3. 모형의 잔차검정
plot(model)
#Hit  to see next plot: 잔차 vs 적합값 -> 패턴없이 무작위 분포(포물선 분포 좋지않은 적합) 
#Hit  to see next plot: Normal Q-Q -> 정규분포 : 대각선이면 잔차의 정규성 
#Hit  to see next plot: 척도 vs 위치 -> 중심을 기준으로 고루 분포 
#Hit  to see next plot: 잔차 vs 지렛대값 -> 중심을 기준으로 고루 분포 

# (1) 등분산성 검정 
plot(model, which =  1) 
methods('plot') # plot()에서 제공되는 객체 보기 

# (2) 잔차 정규성 검정
attributes(model) # coefficients(계수), residuals(잔차), fitted.values(적합값)
res <- residuals(model) # 잔차 추출 
shapiro.test(res) # 정규성 검정 - p-value = 0.9349 >= 0.05
# 귀무가설 : 정규성과 차이가 없다.

# 정규성 시각화  
hist(res, freq = F) 
qqnorm(res)

# (3) 잔차의 독립성(자기상관 검정 : Durbin-Watson) 
install.packages('lmtest')
library(lmtest) # 자기상관 진단 패키지 설치 
dwtest(model) # 더빈 왓슨 값
# DW = 2.0604, p-value = 0.6013 >= 0.05

# 4. 다중공선성 검사 
library(car)
sqrt(vif(model)) > 2 # TRUE 

# 5. 모델 생성/평가 
formula = Sepal.Length ~ Sepal.Width + Petal.Length 
model <- lm(formula = formula,  data=iris)
summary(model) # 모델 평가


-------------------------------------------------------------------------------------------------------------

# Chap15_2_LogisticRegression 

 

# 로지스틱 회귀분석

# 1. 오즈비(odds ratio) 

# 2. 로짓변환

# 3. sigmoid function

###############################################
# 15_2. 로지스틱 회귀분석(Logistic Regression) 
###############################################

# 목적 : 일반 회귀분석과 동일하게 종속변수와 독립변수 간의 관계를 나타내어 
# 향후 예측 모델을 생성하는데 있다.

# 차이점 : 종속변수가 범주형 데이터를 대상으로 하며 입력 데이터가 주어졌을 때
# 해당 데이터의결과가 특정 분류로 나눠지기 때문에 분류분석 방법으로 분류된다.
# 유형 : 이항형(종속변수가 2개 범주-Yes/No), 다항형(종속변수가 3개 이상 범주-iris 꽃 종류)
# 다항형 로지스틱 회귀분석 : nnet, rpart 패키지 이용 
# a : 0.6,  b:0.3,  c:0.1 -> a 분류 

# 분야 : 의료, 통신, 기타 데이터마이닝

# 선형회귀분석 vs 로지스틱 회귀분석 
# 1. 로지스틱 회귀분석 결과는 0과 1로 나타난다.(이항형)
# 2. 정규분포 대신에 이항분포를 따른다.
# 3. 로직스틱 모형 적용 : 변수[-무한대, +무한대] -> 변수[0,1]사이에 있도록 하는 모형 
#    -> 로짓변환 : 출력범위를 [0,1]로 조정
# 4. 종속변수가 2개 이상인 경우 더미변수(dummy variable)로 변환하여 0과 1를 갖도록한다.
#    예) 혈액형 AB인 경우 -> [1,0,0,0] AB(1) -> A,B,O(0)


# 1. 오즈비(odds ratio) : [실패(0)에 대한] 성공(1) 비율
# p : 성공확률
# 1-p : 실패확률(~p)

p = 0.5 # success 
1-p # fail

odds_ratio = p / (1-p) # 실패에 대한 성공비율 
odds_ratio # 1
p = 0.6 # success
odds_ratio = p / (1-p) # 실패에 대한 성공비율 
odds_ratio # 1.5

p = 0.4 # success
odds_ratio = p / (1-p) # 실패에 대한 성공비율 
odds_ratio # 0.6666667
# [정리] 성공확률이 큰 경우 오즈비 증가 


# 2. 로짓변환 : 오즈비에 log(자연로그) 함수 적용 
p = 0.5
odds_ratio = p / (1-p)
odds_ratio # 1
logit = log(odds_ratio)
logit # 0

p = 1
odds_ratio = p / (1-p)
odds_ratio # Inf
logit = log(odds_ratio)
logit # Inf

p = 0
odds_ratio = p / (1-p)
odds_ratio # Inf
logit = log(odds_ratio)
logit # -Inf
# [정리] 성공확률 0.5 이상 -> Inf, 이하면 -Inf

# 3. sigmoid function : 로짓값에 지수(exp)함수 적용 
# sigmoid function = 1 / (1 + exp(-로짓값))
logit = 0
sig = 1 / (1 + exp(-logit)) # 0
sig # 0.5

logit = Inf
sig = 1 / (1 + exp(-logit)) # 0
sig # 1

logit = -Inf
sig = 1 / (1 + exp(-logit)) # 0
sig # 0
# [정리] logit값이 무한이 큰 값(Inf)은 1로, 무한이 작은 값(-Inf) 0


# 단계1. 데이터 가져오기
weather = read.csv("C:/2_Rwork/Part-IV/weather.csv", 
                   stringsAsFactors = F) 
dim(weather)  # 366  15
head(weather)
str(weather)

# chr 칼럼, Date, RainToday 칼럼 제거 
weather_df <- weather[, c(-1, -6, -8, -14)]
str(weather_df)

# RainTomorrow 칼럼 -> 로지스틱 회귀분석 결과(0,1)에 맞게 더미변수 생성      
weather_df$RainTomorrow[weather_df$RainTomorrow=='Yes'] <- 1
weather_df$RainTomorrow[weather_df$RainTomorrow=='No'] <- 0
weather_df$RainTomorrow <- as.numeric(weather_df$RainTomorrow)
head(weather_df)

#  단계2.  데이터 셈플링 : 70 vs 30
set.seed(123)
idx <- sample(1:nrow(weather_df), nrow(weather_df)*0.7)
train <- weather_df[idx, ]
test <- weather_df[-idx, ]
dim(train) # 256  11
dim(test) # 110  11

#  단계3. 로지스틱 회귀모델 생성 : 학습데이터 
weater_model <- glm(RainTomorrow ~ ., data = train, family = 'binomial')
weater_model 
summary(weater_model) 


# 단계4. 로지스틱  회귀모델 예측치 생성 : 검정데이터 
# newdata=test : 새로운 데이터 셋, type="response" : 0~1 확률값으로 예측 
pred <- predict(weater_model, newdata=test, type="response")  
pred # 110
length(pred)
range(pred, na.rm = T) # 0.0004382011 ~ 0.9953997957
summary(pred)
str(pred)

# cutoff = 0.5
cpred <- ifelse(pred >= 0.5, 1, 0)
table(cpred)
#cpred
# 0  1 
#96 13

real_value <- test$RainTomorrow

# model 평가 : 교차분할표(confusion matrix)
table(real_value, cpred)
#           cpred
#real_value  0  1
#         0 79  3 = 82 
#         1 17 10 = 27

nrow(test) # 100
table(test$RainTomorrow) # 83(no) 27(yes)

# 분류정확도 
acc <- (79 + 10) / (79+17+3+10)
acc # 0.8165138

# 특이도 : NO -> NO 예측 비율 
79 / (79  + 3) # 0.9634146

# 재현율 : YES -> YES 예측 비율 
Recall <- 10 / (17 + 10) # 0.3703704

# 정확률 : 예측치(yes) -> yes 예측 비율 
Precision <- 10 / (3 + 10) # 0.7692308

# F1 score 
F1_score <- 2 * ((Recall * Precision) / (Recall + Precision))
F1_score # 0.5

### ROC Curve를 이용한 모형평가(분류정확도)  ####
# Receiver Operating Characteristic

install.packages("ROCR")
library(ROCR)

# ROCR 패키지 제공 함수 : prediction() -> performance
pr <- prediction(pred, test$RainTomorrow)
prf <- performance(pr, measure = "tpr", x.measure = "fpr")
plot(prf)