filnal exam
jinbum KIM
2020 12 1
융합심화R통계분석 기말고사 레포트
iris 데이터셋을 이용하여 다음 문제를 해결하시오.
데이터 셋 불러오기
data(iris)
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa우선 iris 데이터 셋을 불러온다. iris 데이터 셋은 R 자체에 내장되어 있기 때문에 별도의 패키지를 실행시킬 필요가 없다.(head()함수는 데이터 셋이 잘 불러들여졌는지 확인하기 위해 사용한다.)
1. EDA를 수행하시오.(25점)
EDA(:탐색적 데이터 분석)에 필요한 패키지 실행하기
library(tidyverse)
library(gridExtra)
library(RColorBrewer)
library(stargazer)
library(caret)
library(knitr)
library(car)-탐색적 데이터 분석에 필요한 패키지를 실행시킨다.-
tidyverse : ggplot2, dplyr, stringr 등 총 8개의 패키지들이 내장되어 있다.
gridExtra : 그래프를 출력한 후 한번에 비교할 수 있도록 도와준다.
RColorBrewer : 데이터 시각화에 더 다양한 색감을 부여한다.
나머진 분석의 편의상 사용했다.
결측치 확인하기
table(is.na(iris))##
## FALSE
## 750전반적인 변수 및 속성 확인하기
str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...iris는 Sepal.Length, Sepal.Width, Petal.Length, Petal.Width, Species로 총 5변수로 구성되어 있다. 이 변수들 중 Species는 Factor의 속성이며, 나머지는 numeric의 속성을 보이고 있다.
요약 통계량 확인하기
summary(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100
## 1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300
## Median :5.800 Median :3.000 Median :4.350 Median :1.300
## Mean :5.843 Mean :3.057 Mean :3.758 Mean :1.199
## 3rd Qu.:6.400 3rd Qu.:3.300 3rd Qu.:5.100 3rd Qu.:1.800
## Max. :7.900 Max. :4.400 Max. :6.900 Max. :2.500
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
## 각 변수들의 전반적인 요약 통계량은 위와 같다. 다만 Species 변수의 경우 factor의 속성이기 때문에 다른 속성들과 다르게 각 value의 갯수가 출력이 된다.
종속변수의 분포 살펴보기
iris %>%
select(Sepal.Width) %>%
ggplot(mapping = aes(x = Sepal.Width))+
stat_density(alpha = 0.7,
fill = 'steelblue',
color = 'blue')+
labs(title = 'Distribution:Sepal.Width')
모형을 추정하기 위해서는 종속변수의 분포 형태를 살펴보는 것이 중요한데, Sepal.width의 경우 완벽하진 않지만 전반적으로 대칭적인 정규분포와 유사한 분포를 보이는 것으로 나타난다.
나머지 변수들의 분포 확인하기
A <- iris %>%
select(Sepal.Length) %>%
ggplot(mapping = aes(x = Sepal.Length))+
stat_density(alpha = 0.6,
fill = 'orange')
B<- iris %>%
select(Petal.Length) %>%
ggplot(mapping = aes(x = Petal.Length))+
stat_density(alpha = 0.6,
fill = 'green')
C <- iris %>%
select(Petal.Width) %>%
ggplot(mapping = aes(x = Petal.Width))+
stat_density(alpha = 0.6,
fill = 'gray')
D <- iris %>%
select(Species) %>%
ggplot(mapping = aes(x = Species))+
geom_bar(alpha = 0.6,
fill = 'firebrick')
gridExtra::grid.arrange(A, B, C, D, ncol = 2)
나머지 변수들의 분포의 형태는 위 그래프들을 참고하면 알 수 있다.
종속변수와 개별 독립변수들 간의 관계 파악하기
AA <- iris %>%
ggplot(mapping = aes(x = Sepal.Length, y = Sepal.Width))+
geom_point(aes(color = Species))+
geom_smooth(method = 'lm')+
scale_color_brewer(palette = 'Set1')
BB <- iris %>%
ggplot(mapping = aes(x = Petal.Length, y = Sepal.Width))+
geom_point(aes(color = Species))+
geom_smooth(method = 'lm')+
scale_color_brewer(palette = 'Set1')
CC <- iris %>%
ggplot(mapping = aes(x = Petal.Width, y = Sepal.Width))+
geom_point(aes(color = Species))+
geom_smooth(method = 'lm')+
scale_color_brewer(palette = 'Set1')
DD <- iris %>%
ggplot(mapping = aes(x = Species, y = Sepal.Width))+
geom_boxplot(aes(color = Species))+
scale_color_brewer(palette = 'Set1')
gridExtra::grid.arrange(AA, BB, CC, DD, ncol = 2)## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
종속변수와 각 독립변수들의 회귀선을 추정해본 결과, Sepal.Length, Petal.Length, Petal.Width는 Sepal.width과 약간의 음의 관계를 갖는 것으로 나타난다. 하지만 이 직선은 Sepal.Width별로 구분되어 도출된 것이 아니기 때문에, 종속변수와 독립변수의 관계에 있어서 교란이 일어난 듯하다. (종속변수 내 Species의 boxplot은 마지막 그래프와 같다.)
Species를 고려하여 종속변수와 독립변수의 관계 도출하기
Species별 Sepal.Length와 Sepal.Width의 회기선
iris %>%
ggplot(mapping = aes(x = Sepal.Length, y = Sepal.Width))+
geom_point(aes(color = Species))+
geom_smooth(method = 'lm')+
facet_wrap(~Species)+
scale_color_brewer(palette = 'Set1')## `geom_smooth()` using formula 'y ~ x'
Species별 Petal.Length와 Sepal.Width의 회기선
iris %>%
ggplot(mapping = aes(x = Petal.Length, y = Sepal.Width))+
geom_point(aes(color = Species))+
geom_smooth(method = 'lm')+
facet_wrap(~Species)+
scale_color_brewer(palette = 'Set1')## `geom_smooth()` using formula 'y ~ x'
Species별 Petal.Width와 Sepal.Width의 회기선
iris %>%
ggplot(mapping = aes(x = Petal.Width, y = Sepal.Width))+
geom_point(aes(color = Species))+
geom_smooth(method = 'lm')+
facet_wrap(~Species)+
scale_color_brewer(palette = 'Set1')## `geom_smooth()` using formula 'y ~ x'
Species별로 구분해서 추정한 회귀선은 모두 우상향하는 형태를 보인다. 따라서 이 점을 반영한 종속변수와 독립변수 간의 관계는 모두 양의 상관관계임을 알 수 있다.
2. Sepal.Width를 종속변수로 하는 회귀 모형(Regression Model)을 제시하시오.
다중회기 모형 추정하기
lm_model <- lm(Sepal.Width ~ ., data = iris)
summary(lm_model)##
## Call:
## lm(formula = Sepal.Width ~ ., data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.00102 -0.14786 0.00441 0.18544 0.69719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.65716 0.25595 6.475 1.40e-09 ***
## Sepal.Length 0.37777 0.06557 5.761 4.87e-08 ***
## Petal.Length -0.18757 0.08349 -2.246 0.0262 *
## Petal.Width 0.62571 0.12338 5.072 1.20e-06 ***
## Speciesversicolor -1.16029 0.19329 -6.003 1.50e-08 ***
## Speciesvirginica -1.39825 0.27715 -5.045 1.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2678 on 144 degrees of freedom
## Multiple R-squared: 0.6352, Adjusted R-squared: 0.6225
## F-statistic: 50.14 on 5 and 144 DF, p-value: < 2.2e-16우선 종속변수와 모든 독립변수를 고려하여 다중회귀 모형을 추정한다.(traning data 사용)
효율적인 모형 추정하기
-1. 전진선택법
out1 <- step(lm_model, direction = 'forward')## Start: AIC=-389.37
## Sepal.Width ~ Sepal.Length + Petal.Length + Petal.Width + Speciessummary(out1)##
## Call:
## lm(formula = Sepal.Width ~ Sepal.Length + Petal.Length + Petal.Width +
## Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.00102 -0.14786 0.00441 0.18544 0.69719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.65716 0.25595 6.475 1.40e-09 ***
## Sepal.Length 0.37777 0.06557 5.761 4.87e-08 ***
## Petal.Length -0.18757 0.08349 -2.246 0.0262 *
## Petal.Width 0.62571 0.12338 5.072 1.20e-06 ***
## Speciesversicolor -1.16029 0.19329 -6.003 1.50e-08 ***
## Speciesvirginica -1.39825 0.27715 -5.045 1.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2678 on 144 degrees of freedom
## Multiple R-squared: 0.6352, Adjusted R-squared: 0.6225
## F-statistic: 50.14 on 5 and 144 DF, p-value: < 2.2e-16-2. 후진제거법
out2 <- step(lm_model, direction = 'backward')## Start: AIC=-389.37
## Sepal.Width ~ Sepal.Length + Petal.Length + Petal.Width + Species
##
## Df Sum of Sq RSS AIC
## <none> 10.328 -389.37
## - Petal.Length 1 0.36195 10.689 -386.21
## - Petal.Width 1 1.84465 12.172 -366.72
## - Sepal.Length 1 2.38066 12.708 -360.26
## - Species 2 3.14643 13.474 -353.48summary(out2)##
## Call:
## lm(formula = Sepal.Width ~ Sepal.Length + Petal.Length + Petal.Width +
## Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.00102 -0.14786 0.00441 0.18544 0.69719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.65716 0.25595 6.475 1.40e-09 ***
## Sepal.Length 0.37777 0.06557 5.761 4.87e-08 ***
## Petal.Length -0.18757 0.08349 -2.246 0.0262 *
## Petal.Width 0.62571 0.12338 5.072 1.20e-06 ***
## Speciesversicolor -1.16029 0.19329 -6.003 1.50e-08 ***
## Speciesvirginica -1.39825 0.27715 -5.045 1.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2678 on 144 degrees of freedom
## Multiple R-squared: 0.6352, Adjusted R-squared: 0.6225
## F-statistic: 50.14 on 5 and 144 DF, p-value: < 2.2e-16-3. 단계적 검정법
out3 <- step(lm_model, direction = 'both')## Start: AIC=-389.37
## Sepal.Width ~ Sepal.Length + Petal.Length + Petal.Width + Species
##
## Df Sum of Sq RSS AIC
## <none> 10.328 -389.37
## - Petal.Length 1 0.36195 10.689 -386.21
## - Petal.Width 1 1.84465 12.172 -366.72
## - Sepal.Length 1 2.38066 12.708 -360.26
## - Species 2 3.14643 13.474 -353.48summary(out3)##
## Call:
## lm(formula = Sepal.Width ~ Sepal.Length + Petal.Length + Petal.Width +
## Species, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.00102 -0.14786 0.00441 0.18544 0.69719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.65716 0.25595 6.475 1.40e-09 ***
## Sepal.Length 0.37777 0.06557 5.761 4.87e-08 ***
## Petal.Length -0.18757 0.08349 -2.246 0.0262 *
## Petal.Width 0.62571 0.12338 5.072 1.20e-06 ***
## Speciesversicolor -1.16029 0.19329 -6.003 1.50e-08 ***
## Speciesvirginica -1.39825 0.27715 -5.045 1.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2678 on 144 degrees of freedom
## Multiple R-squared: 0.6352, Adjusted R-squared: 0.6225
## F-statistic: 50.14 on 5 and 144 DF, p-value: < 2.2e-16회귀모형 추정결과 비교
stargazer::stargazer(out1, out2, out3,
type = 'text',
keep.stat = c('n',
'rsq',
'adj.rsq'))##
## ===============================================
## Dependent variable:
## -----------------------------
## Sepal.Width
## (1) (2) (3)
## -----------------------------------------------
## Sepal.Length 0.378*** 0.378*** 0.378***
## (0.066) (0.066) (0.066)
##
## Petal.Length -0.188** -0.188** -0.188**
## (0.083) (0.083) (0.083)
##
## Petal.Width 0.626*** 0.626*** 0.626***
## (0.123) (0.123) (0.123)
##
## Speciesversicolor -1.160*** -1.160*** -1.160***
## (0.193) (0.193) (0.193)
##
## Speciesvirginica -1.398*** -1.398*** -1.398***
## (0.277) (0.277) (0.277)
##
## Constant 1.657*** 1.657*** 1.657***
## (0.256) (0.256) (0.256)
##
## -----------------------------------------------
## Observations 150 150 150
## R2 0.635 0.635 0.635
## Adjusted R2 0.622 0.622 0.622
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01위 세 가지 과정을 통해 추정된 모형들의 통계값이 모두 일치하는 것으로 보인다. 따라서 위 세가지 방법 중 어떤 모형을 채택해도 큰 이상이 없다.(변수가 적고 다 유의하게 판단되어 모두 같게 나온다.)
VIF : 분산팽창요인 살펴보기
car::vif(out1)## GVIF Df GVIF^(1/(2*Df))
## Sepal.Length 6.124653 1 2.474804
## Petal.Length 45.132550 1 6.718076
## Petal.Width 18.373804 1 4.286468
## Species 32.701564 2 2.391344다중공선성의 문제를 방지하기 위해 분산팽창요인을 살펴본다. 보통 10이 넘으면 모형에 문제가 있다고 가정하기 때문에 문제가 없는 것 같지만, 사실 5가 넘어가는 수치 역시 작은 편이 아니기 때문에 유의할 필요가 있다. 본 과정에서는 Petal.Length의 VIF가 6.718076으로 5이상이지만, 위 모형 추정 결과 이 변수의 개별 P-value가 상대적으로 유의하다 판단하여 문제삼지 않았다.
결론 = 아래 회귀식과 같은 모형이 가장 적합한 것으로 보인다.
\[sepal.width = a_i + b_1*sepal.length + b_2*Petal.Length + b_3*Petal.Width + b_4*Species + e_i\]
training data와 test data 분리하기 ( 8 : 2 )
random <- createDataPartition(y = iris$Species,
p = 0.8,
list = F)
train <- iris[random,]
test <- iris[-random,]train : training data
test : test data
분리 데이터를 바탕으로 모형추정하기
train_model <- lm(formula = Sepal.Width ~ Sepal.Length + Petal.Length + Petal.Width +
Species, data = train)
stargazer::stargazer(train_model,
type = 'text',
keep.stat = c('n', 'rsq', 'adj.rsq'))##
## =============================================
## Dependent variable:
## ---------------------------
## Sepal.Width
## ---------------------------------------------
## Sepal.Length 0.389***
## (0.069)
##
## Petal.Length -0.166*
## (0.091)
##
## Petal.Width 0.596***
## (0.132)
##
## Speciesversicolor -1.144***
## (0.216)
##
## Speciesvirginica -1.384***
## (0.310)
##
## Constant 1.555***
## (0.265)
##
## ---------------------------------------------
## Observations 120
## R2 0.627
## Adjusted R2 0.610
## =============================================
## Note: *p<0.1; **p<0.05; ***p<0.01추정 결과 가장 적합한 것으로 보이는 모형을 trainging data를 바탕으로 학습시킨다.
예측과 평가하기
pre <- predict((train_model), newdata = test)
dat <- test %>%
select(Sepal.Width)
evaluation <- data.frame(pre, dat)
names(evaluation) <- c('predict','actual')
head(evaluation)## predict actual
## 12 3.278111 3.4
## 18 3.487734 3.5
## 20 3.471138 3.8
## 21 3.495206 3.4
## 25 3.228326 3.4
## 28 3.450500 3.5위 모형을 바탕으로 분리해둔 test data를 이용하여 예측값(predict)을 추정한다. 또한 데이터 내 실제값(actual)을 바탕으로 모형의 오차를 평가한다.
evaluation %>%
mutate(rs = (predict - actual)) %>%
summarise(sum = sum(rs),
mean = mean(rs),
var = var(rs),
sd = sd(rs))## sum mean var sd
## 1 3.291103 0.1097034 0.08559854 0.2925723이 모형의 오차 평가 결과는 다음과 같다.
sum : 총합
mean : 평균
var : 분산
sd : 표준편차