출처: 시작하세요! 데이터 분석 with R : R로 배우는 기초 통계와 데이터 분석 기법
Github Code: Jun4871 Github
Overview
대부분의 물건에는 수명이 존재하고, 더 오랫동안 문제없이 사용하기 위한 사후관리가 필요하다. 이러한 배경에서 A/S 라는 것이 탄생하였을 것이다. 물론 A/S가 발생하는 원인 중에는 초기생산 불량이나 사용부주의에 의한 조기 파손 등 다양한 이유로 발생하겠지만, 이러한 케이스를 제외하고 본다면 A/S는 언제, 얼마나 발생하게 될까? 이는 기업의 관점에서 중요한 문제일 것이고 곧 인력배치나 재쟁편성과도 직결된다고 볼 수 있다. 제공된 데이터를 활용하여 유사상황에서의 솔루션을 만들어보자.
1. 데이터 로드 및 라이브러리 활성화
2. 데이터 탐색
데이터 프레임 구조이며 총 3개의 컬럼과 31개의 행으로 이루어져 있고, 각각 전월 정수기 총 대여 수, 전월 10년 이상 노후 정수기 총 대여 수, 당월 AS에 소요된 시간을 의미한다.
## 'data.frame': 31 obs. of 3 variables:
## $ purifier : num 168750 171450 172800 174000 174810 ...
## $ old_purifier: num 33750 42863 31104 40020 26222 ...
## $ as_time : num 20453 21850 20214 21660 19659 ...## purifier old_purifier as_time
## Min. :168750 Min. :26145 Min. :19659
## 1st Qu.:179708 1st Qu.:31805 1st Qu.:21200
## Median :194145 Median :36888 Median :22994
## Mean :192377 Mean :37215 Mean :22896
## 3rd Qu.:204000 3rd Qu.:43058 3rd Qu.:24316
## Max. :216375 Max. :51552 Max. :270703. 데이터 시각화
“총 정수기 대여 대수”와 “노후 정수기 대여 대수”에 대해 “AS 소요시간”과의 관계를 산점도로 표현해보도록 하자.
1) 총 정수기 대여 수와 AS시간
그래프의 결과와 같이 “양의 상관관계”를 가지고 있음을 확인할 수 있다. 조금 더 정확하기 확인하기 위해 AS시간과의 상관계수를 산출해보도록 하자. cor()함수를 통해 0.9102497의 수치로 실제 상관관계가 매우 높음이 확인되었다.
#plot(purifier_df$purifier, purifier_df$as_time, xlab = "총정수기대여수", ylab = "AS시간")
ggplot(data = purifier_df) +
geom_point(aes(x=purifier_df$purifier, y=purifier_df$as_time,
color = purifier_df$purifier), size = 3) +
theme_classic() +
theme(axis.text.x = element_text(size=10)) +
theme(axis.text.y = element_text(size=10)) +
xlab("총 정수기 대여 대수(전월)") +
ylab("A/S 시간(당월)")
## [1] 0.91024972) 노후 정수기 대여 수와 AS시간
0.8795908의 수치로 노후 정수기 역시 AS 소요시간과 양의 상관관계를 가지고 있다는 것을 확인할 수 있다.
#plot(purifier_df$old_purifier, purifier_df$as_time, xlab = "노후 정수기대여수", ylab = "AS시간")
ggplot(data = purifier_df) +
geom_point(aes(x=purifier_df$old_purifier, y=purifier_df$as_time,
color = purifier_df$purifier), size = 3) +
theme_classic() +
theme(axis.text.x = element_text(size=10)) +
theme(axis.text.y = element_text(size=10)) +
xlab("노후 정수기 대여 대수(전월)") +
ylab("A/S 시간(당월)")
## [1] 0.87959083) 의문점
“총 정수기 대여 대수”, “노후 정수기 대여 대수” 모두 “AS 소요시간”과 “양의 상관관계”를 가지고 있다는 사실만으로 아래 질문들에 대한 답이 가능한가 ?
정수기 대여 수가 1대 증가할 때마다 AS 소요시간은 얼마나 늘어날까?
노후 정수기가 1대 증가할 때마다 AS 소요시간은 또 얼마나 늘어날까?
정수기 대여 수와 노후 정수기 대여 수 중 AS 소요시간에 더 영향을 주는 항목은 무엇일까?
이번 달 정수기 대여 수와 노후 정수기 대여 수를 알고 있다면 익월의 AS 시간을 예측할 수 있을까?
4. 회귀분석 실험
R 내에 내장되어 있는 “cars” 데이터로 위 의문점들을 해결해보자.
1) 데이터 로드
데이터 프레임 형식에 2개의 컬럼을 가지고 있고 50개의 관측값을 가지고 있다. lm() 함수를 통해 회귀분석을 수행할 수 있는데 이 결과를 lm_result에 할당해주자. 여기서 ~의 의미는 종속변수와 독립변수를 구분짓는 기호라고 보면 된다.
## 'data.frame': 50 obs. of 2 variables:
## $ speed: num 4 4 7 7 8 9 10 10 10 11 ...
## $ dist : num 2 10 4 22 16 10 18 26 34 17 ...2) 결괏값 확인하기
summary의 결괏값을 확인해보자.
Adjusted R-squared:0.6438 모델의 정확도를 의미하고 도출된 회귀모델을 가지고 현 데이터를 64.38% 설명할 수 있다는 의미이다.
p-value: 1.49e-12 모델의 통계적 유의성을 의미하고, 1.49*10-12의 값으로 0.05보다 매우 작은 값이며, 보통 유의 수준에서 95%에서 0.05보다 작다면 도출된 회귀모델이 통계적으로 유의하다고 인정한다.
이 결과를 기준으로 회귀모델 방정식은 다음과 같이 도출된다. dist(제동거리 feet)=3.9324 * speed(차속도 mi/h) - 17.5791 각 독립변수의 p-value 또한 모두 0.05 미만이므로 통계적으로 유의하다.
##
## Call:
## lm(formula = dist ~ speed, data = cars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -29.069 -9.525 -2.272 9.215 43.201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -17.5791 6.7584 -2.601 0.0123 *
## speed 3.9324 0.4155 9.464 1.49e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438
## F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12## (Intercept) speed
## -17.579095 3.932409## 2.5 % 97.5 %
## (Intercept) -31.167850 -3.990340
## speed 3.096964 4.767853## [1] 11353.52## 1 2 3 4 5 6 7 8
## -1.849460 -1.849460 9.947766 9.947766 13.880175 17.812584 21.744993 21.744993
## 9 10 11 12 13 14 15 16
## 21.744993 25.677401 25.677401 29.609810 29.609810 29.609810 29.609810 33.542219
## 17 18 19 20 21 22 23 24
## 33.542219 33.542219 33.542219 37.474628 37.474628 37.474628 37.474628 41.407036
## 25 26 27 28 29 30 31 32
## 41.407036 41.407036 45.339445 45.339445 49.271854 49.271854 49.271854 53.204263
## 33 34 35 36 37 38 39 40
## 53.204263 53.204263 53.204263 57.136672 57.136672 57.136672 61.069080 61.069080
## 41 42 43 44 45 46 47 48
## 61.069080 61.069080 61.069080 68.933898 72.866307 76.798715 76.798715 76.798715
## 49 50
## 76.798715 80.731124## 1 2 3 4 5 6 7
## 3.849460 11.849460 -5.947766 12.052234 2.119825 -7.812584 -3.744993
## 8 9 10 11 12 13 14
## 4.255007 12.255007 -8.677401 2.322599 -15.609810 -9.609810 -5.609810
## 15 16 17 18 19 20 21
## -1.609810 -7.542219 0.457781 0.457781 12.457781 -11.474628 -1.474628
## 22 23 24 25 26 27 28
## 22.525372 42.525372 -21.407036 -15.407036 12.592964 -13.339445 -5.339445
## 29 30 31 32 33 34 35
## -17.271854 -9.271854 0.728146 -11.204263 2.795737 22.795737 30.795737
## 36 37 38 39 40 41 42
## -21.136672 -11.136672 10.863328 -29.069080 -13.069080 -9.069080 -5.069080
## 43 44 45 46 47 48 49
## 2.930920 -2.933898 -18.866307 -6.798715 15.201285 16.201285 43.201285
## 50
## 4.2688763) 결괏값 시각화
이번에는 실제로 관측된 독립변수와 종속변수 간의 산점도 위에 회귀모델 방정식의 그래프를 추가해 회귀분석이 실제 값을 얼마나 설명하는지 확인해 보도록 하자. 종속변수인 제동거리를 y축, 독립변수인 속도를 x축으로 하는 산점도를 그려보자.
ggplot(data = cars) +
geom_point(aes(cars$speed, cars$dist, color = cars$speed), size = 2) +
theme_light() +
theme(axis.text.x = element_text(size=10)) +
theme(axis.text.y = element_text(size=10)) +
xlab("속도") +
ylab("제동거리") 
이제 종속변수와 독립변수가 표현된 산점도 위에 회귀모델 방정식의 선을 추가해 보자. 기존 그래프 위에 선을 추가하려면 abline() 함수를 사용하면 된다.
Residuals vs. Fitted : 실제 값과 예측 값을 나타내는 그래프로서 빨간 선의 기울기가 0에 가까울수록(가로 점선) 좋은 모델이라고 볼 수 있다.
Normal Q-Q : 실제 값과 예측 값의 차이(잔차)가 정규분포를 따르는지를 확인하는 그래프로스, 직선에 점들이 밀집돼 있으면 좋은 모델이라고 불 수 있다.
Scale-Location : 실제 값과 예측 값의 차이가 등분산성을 따르는지 확인하는 그래프로서 Residuals vs. Fitted 그래프처럼 빨간 선의 기울기가 0에 가까울수록 좋은 모델이라고 볼 수 있다.
Residuals vs. Leverage : 이상치를 표현하는 그래프로서 빨간 점선 안에 점들이 들어있지 않을수록 좋은 모델이라고 볼 수 있다.
# "차속도"에 따른 "제동거리" 회귀분석
lm_result_1 <- lm(formula = dist ~ speed, data = cars)
# 회귀모델 방정식 선 그리기
ggplot(data = lm_result_1) +
geom_point(aes(x=lm_result_1$model$speed, y=lm_result_1$model$dist, color = lm_result_1$model$speed)) +
geom_abline(intercept=-17.579095, slope=3.932409) +
theme_light() +
theme(axis.text.x = element_text(size=10)) +
theme(axis.text.y = element_text(size=10)) +
xlab("속도") +
ylab("제동거리")
# Resi_Fitt <- ggplot(data = lm_result_1) +
# geom_point(aes( x = lm_result_1$fitted.values, y = lm_result_1$residuals, color = lm_result_1$residuals)) +
# theme_light() +
# theme(axis.text.x = element_text(size=10)) +
# theme(axis.text.y = element_text(size=10)) +
# xlab("속도") +
# ylab("제동거리")
# Theore_Stdliz <- ggplot(data = lm_result_1) +
# geom_point(aes( x = lm_result_1$, y = lm_result_1$residuals, color = lm_result_1$residuals)) +
# theme_light() +
# theme(axis.text.x = element_text(size=10)) +
# theme(axis.text.y = element_text(size=10)) +
# xlab("속도") +
# ylab("제동거리")
# 가로 2, 세로 2개의 그래프를 한 번에 그리도록 설정
plot(lm_result_1)
4) 예측
새로운 독립변수를 회귀모델 방정식에 대입해 종속변수를 예측할 수 있다. 예를 들어 속도 70mi/h를 도출한 회귀모델 방정식에 대입해보면 예상되는 제동거리는 257.68 feet가 된다. (3.9324 * 70 - 17.5791 = 257.6889)
점 추정과 구간 추정
이처럼 종속변수 값을 특정 값 하나로 예측해 주는 것을 “점 추정(point esimation)” 방식이라 한다. 이러한 방식은 “가장 좋은 단일 예측 값”을 제시함으로써 간단명료하게 표현할 수 있다는 장점이 있지만 이 예측 값의 불확실성을 표현하지 못한다는 단점도 있다. 실제로 회귀모델 방정식만 보더라도 다음 두 가지의 불 확실성을 고려해 볼 수 있다.
- 회귀모델 방정식의 계수에 대한 불확실성
- 회귀모델 방정식을 통해 나온 결괏값에 대한 불확실성(결과오차범위)
이러한 불확실성을 감안해 종속변수를 하나의 값이 아닌 범위 값으로 제시하는 방식을 “구간 추정”이라고 한다.예를 들어, 70 mi/h 의 제동거리를 “구간 추정” 방식으로 표현한다면 다음과 같이 표현할 수 있다.
“제동거리는 202 feet ~ 312 feet (신뢰구간) 사이일 확률이 95%(신뢰수준)다” 이때 신뢰수준이 높고, 신뢰구간이 좁을수록 좋은 모델이다. 이제 predict() 함수를 통해 회귀모델의 예측 값을 산출해보자. 먼저 50, 60, 70, 80, 90 100 mi/h 속도일 때 제동거리를 예측해보자.
먼저 점 추정 방식으로 제동거리를 예측해 봤는데, predict함수의 결과가 독립변수 없이 예측 값만 출력하다 보니 가독성이 떨어진다. 결과값에 독립변수를 붙여 함께 확인해보자.
# "차속도"에 따른 "제동거리" 회귀분석
lm_result_1 <- lm(formula = dist ~ speed, data =cars)
# 예측할 독립변수 데이터프레임 생성
# 데이터프레임을 생성할 때는 회구분석 시 사용한 독립변수명과 동일하게 칼럼명을 생성
speed <- c(50, 60, 70, 80, 90, 100)
df_input <- data.frame(speed)
# 입력값 확인
df_input## 1 2 3 4 5 6
## 179.0413 218.3654 257.6895 297.0136 336.3377 375.661860 mi/h 속도라면 제동거리가 218.3654 feet이 된다는 것을 예측할 수 있다.
# predict 결괏값 SET
predict_dist <- predict(lm_result_1, df_input)
# predict 결괏값 구조 확인 -> 숫자로 된 1차 배열
str(predict_dist)## Named num [1:6] 179 218 258 297 336 ...
## - attr(*, "names")= chr [1:6] "1" "2" "3" "4" ...이번에는 구간 추정으로 예측을 해보자. fit은 점 추정 예측 값이며, lwr은 구간 최솟값, upr은 구간 최댓값을 의미한다. 아례의 예측 결과를 토대로 60 mi/h 속도인 경우 제동거리를 예측하면 다음과 같다.
- 점 추정방식: 218.3654 feet
- 구간 추정방식 : 180.8489 feet ~ 2565.8820 feet 사이일 확률 95%
# 모델계수에 대한 불확실성을 감안한 구간 추정(confidence) / 신뢰구간 95%
predict_dist <- predict(lm_result_1, df_input, interval = "confidence", level = 0.95)
# 예측 결과
predict_dist## fit lwr upr
## 1 179.0413 149.8060 208.2766
## 2 218.3654 180.8489 255.8820
## 3 257.6895 211.8651 303.5139
## 4 297.0136 242.8670 351.1602
## 5 336.3377 273.8603 398.8151
## 6 375.6618 304.8480 446.4755이번에는 interval에 “prediction”을 설정해 모델계수의 불확실성과 결과에 대한 오차까지 감안한 구간을 추정해보도록 하자. interval을 “prediction”으로 설정하면 도출된 회귀분석 모델 결과의 오차까지도 감안하기 때문에 “confidence”보다 더 넓은 구간으로 제동거리를 예측하게 된다.
- confidence 구간추정 : ‘180.8489 feet’ ~ ‘255,8820 feet’ 사이일 확률이 95%
- prediction 구간추정 : ‘169.7474 feet’ ~ ‘266.9834 feet’ 사이일 확률이 95%
# 모델계수에 대한 불확실성 + 결괏값에 대한 오차를 감안한 구간 추정 / 신뢰구간 95%
predict_dist <- predict(lm_result_1, df_input, interval = "prediction", level = 0.95)
# 입력값과 함께 보기
cbind(df_input, predict_dist)5. 정수기 A/S 기사는 몇명이 적당할까?
이제 회귀분석을 활용해 대여된 정수기 대수를 기준으로 추후 소요될 AS시간을 예측하고 AS기사를 몇 명 고용할지 산출해보자.
purifier : 전월 정수기 총 대여 수
old_purifier : 전월 10년 이상 노후 정수기 총 대여 수
as_time: 당월 AS에 소요된 시간
## purifier old_purifier as_time
## Min. :168750 Min. :26145 Min. :19659
## 1st Qu.:179708 1st Qu.:31805 1st Qu.:21200
## Median :194145 Median :36888 Median :22994
## Mean :192377 Mean :37215 Mean :22896
## 3rd Qu.:204000 3rd Qu.:43058 3rd Qu.:24316
## Max. :216375 Max. :51552 Max. :27070위의 데이터프레임을 통해 회귀분석에 필요한 독립변수와 종속변수를 간단히 나눠 보면 아래와 같다.
독립변수 : 전월 기준 정수기 총 대여 수, 전월 기준 10년 이상 노후 정수기 총 대여 수
종속변수 : 당월 기준 AS에 소요된 시간
독립변수를 선정할 때는 독립변수 간의 상관성, 포함성이 있는지 확인해야 하며 존재한다면 제거해 회귀분석 모델의 왜곡을 최소화해야 한다. “전월 기준 정수기 총 대여 수”에는 “전월 기준 10년 이상 노후 정수기 대여 수”가 포함돼 있고, 이로 인해 “전월 10년 이상 노후 정수기 대여 수”가 중복으로 회귀분석에 반영됨으로써 모델상 왜곡이 발생하게 된다. 이 경우에는 총 대여 수를 “10년 미만 정수기”와 “10년 이상 정수기”로 나눠 독립변수에 적용하는 것이 바람직 할 것이다. 또 아래 상관성 분석을 통해서도 확인 가능하다.
## [1] 0.6042548상관계수가 0.6 이상으로 “전월 정수기 총 대여 수”와 “10년 이상 정수기 대여 수” 간에는 상관성이 매우 크다. 그렇다면 “10년 미만 정수기”와 “10년 이상 정수기”로 나눠 변수 간의 상관성을 확인해 보도록 하자.
# 10년 미만 정수기 = 총 대여 수 - 10년 이상 정수기
# 10년 미만 정수기 vs 10년 이상 정수기 상관성 분석
cor((purifier_df$purifier - purifier_df$old_purifier), purifier_df$old_purifier)## [1] 0.1151678상관계수가 0.1 이하로 상관성이 많이 제거된 것을 확인할 수 있다. 위의 결과를 근거로 독립변수를 다시 정의해 보자.
- 독립변수 : 전월 기준 10년 미만 정수기 대여 수, 전월 기준 10년 이상 노후 정수기 대여 수
- 종속변수 : 당월 기준 AS에 소요된 시간
이제 새로운 독립변수인 “10년 미만 정수기 대여 수”를 데이터프레임에 추가해 보도록 하자.
## 'data.frame': 31 obs. of 3 variables:
## $ purifier : num 168750 171450 172800 174000 174810 ...
## $ old_purifier: num 33750 42863 31104 40020 26222 ...
## $ as_time : num 20453 21850 20214 21660 19659 ...# 10년 미안 정수기 대여 수 (new_purifier) 추가 (총 대여 수 - 10년 이상 정수기)
purifier_df$new_purifier <- purifier_df$purifier - purifier_df$old_purifier
#정수기 데이터 프레임 확인
str(purifier_df)## 'data.frame': 31 obs. of 4 variables:
## $ purifier : num 168750 171450 172800 174000 174810 ...
## $ old_purifier: num 33750 42863 31104 40020 26222 ...
## $ as_time : num 20453 21850 20214 21660 19659 ...
## $ new_purifier: num 135000 128587 141696 133980 148588 ...데이터프레임에 ‘10년 미만 정수기 대여 수’ 컬럼이 추가 되었고, 이를 포함하여 다시 회귀분석을 한다.
# 회귀분석 수행(lm)
# 종속변수 : AS 시간 (as_time)
# 독립변수 : 10년 미만 정수기 (new_purifier), 10년 이상 정수기(old_purifier)
lm_result <- lm(as_time ~ new_purifier + old_purifier, data = purifier_df)아래는 회귀분석 결과이다. ’도출한 모델은 얼마나 적절한가’를 보면 도출된 회귀분석 모델은 p-value가 0.000000000000000022로서 0.05 보다 매우 작고, R-squared 값 또한 0.9996으로 정확성이 높다고 볼 수 있다. 그 다음 독립변수를 보면 10년 미만 정수기(new_purifier), 10년 이상 정수기(old_purifier)의 p-value가 2.2 * 10-16미만으로서 통계적으로 유의한 독립변수임을 알 수 있다. 위 결과를 기반으로 회귀모델 방정식을 구하면 아래와 같다.
AS시간 = 0.0881 * (10년 미만 정수기 대여 수) + 0.2398 * (10년 이상 정수기 대여 수) + 193.7
AS 시간은 10년 미만 정수기 1대당 0.0881 시간이 소요되고, 10년 이상 정수기는 1대당 0.2398 시간이 소요되는 것으로 보인다. 즉 “10년 미만 정수기 대여 수”와 “10년 이상 정수기 대여 수” 모두 “AS시간”에 양의 상관성을 가지고 있으며 그 중 “10년 이상 정수기 대여 수”가 “AS 시간”에 더 영향을 주는 것으로 분석된다.
##
## Call:
## lm(formula = as_time ~ new_purifier + old_purifier, data = purifier_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -75.122 -14.427 -2.473 10.416 178.170
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.937e+02 1.077e+02 1.799 0.0828 .
## new_purifier 8.881e-02 6.742e-04 131.713 <2e-16 ***
## old_purifier 2.398e-01 1.057e-03 226.933 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 40.97 on 28 degrees of freedom
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9996
## F-statistic: 3.837e+04 on 2 and 28 DF, p-value: < 2.2e-16월말 최종 대여 수를 보니 총 대여 수가 300,000대, 그 중 10년 이상 노후 정수기 대수가 70,000대로 집계되었다. 회귀모델 방정식을 통해 예측되는 익월 AS 시간은 43.619시간이 되며 이로 인해 필요한 AS기사는 약 273명으로 예측된다. 물론 구간 추정 방식으로도 확인이 가능하다.
# 예측할 독립변수 값 설정(데이터프레임)
input_predict <- data.frame(new_purifier = 300000, old_purifier = 70000)
# 예측값 저장
predict_as_time <- predict(lm_result, input_predict)
# 예측값 출력
predict_as_time## 1
## 43619.54## 1
## 272.6221# 구간 추정
predict_as_time <- predict(lm_result, input_predict, interval = "confidence", level = 0.95)
# 익월 AS시간이 '43,414시간' ~ '43,824시간' 이 될 가능성을 95%로 예측
predict_as_time## fit lwr upr
## 1 43619.54 43414.58 43824.5