-
[Data Science] 데이터 종류에 따른 분석 기법Data Science/Data Science in R 2018. 10. 25. 02:28
목표 : 데이터 종류에 따라 분석하는 기법들을 익혀보자.
1. 데이터형, 분석 기법, R함수에 따른 분류
데이터 형태
분석 기법과 R함수
0. 모든 데이터
데이터 내용, 구조파악(glimpse)
요약 통계량(summary)
단순 시각화
1. 수량형 변수
분포 시각(hist, boxplot, density)
요약 통계량(mean, median)
t-검정 t.test()
2. 범주형 변수(성공 - 실패)
도수 분포 table(), xtabs()
바그래프 barplot()
이항검정 binom.test()
3. 수량형 x, 수량형 y
산점도 plot()
상관계수 cor()
단순회귀 lm()
로버스트 회귀 lqs()
비모수회귀
4. 범주형 x, 수량형 y
병렬상자그림 boxplot()
분산분석(ANOVA) lm(y~x)
5. 수량형 x, 범주형 y(성공 - 실패)
산점도/병렬상자그림 plot(), boxplot()
로지스틱 회귀분석 glm(family = 'binom')
2. 모든 데이터에 행해야 할 분석
각각의 변수형을 다루기 전에, 어떤 변수든 기본적으로 행할 분석들이 있다.
1) 데이터 내용, 구조, 타입 파악
- dplyr::glimpse()
2) 요약 통계량 파악
- summary()
3) 결측치 확인
- summary()
4) 무작정 시각화
- plot(), pairs() 등
- 데이터의 관측치가 많을 경우 실행시간이 길으니 dplyr::sample_n()을 사용할 것을 권유
다음 예로써 mpg 데이터를 통해 살펴보자.
> library(dplyr)
> library(ggplot2)
> glimpse(mpg)
Observations: 234
Variables: 11
$ manufacturer <chr> "audi", "audi", "audi", "audi", "audi", "audi", "audi", "audi", "audi", "audi", "...
$ model <chr> "a4", "a4", "a4", "a4", "a4", "a4", "a4", "a4 quattro", "a4 quattro", "a4 quattro...
$ displ <dbl> 1.8, 1.8, 2.0, 2.0, 2.8, 2.8, 3.1, 1.8, 1.8, 2.0, 2.0, 2.8, 2.8, 3.1, 3.1, 2.8, 3...
$ year <int> 1999, 1999, 2008, 2008, 1999, 1999, 2008, 1999, 1999, 2008, 2008, 1999, 1999, 200...
$ cyl <int> 4, 4, 4, 4, 6, 6, 6, 4, 4, 4, 4, 6, 6, 6, 6, 6, 6, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, ...
$ trans <chr> "auto(l5)", "manual(m5)", "manual(m6)", "auto(av)", "auto(l5)", "manual(m5)", "au...
$ drv <chr> "f", "f", "f", "f", "f", "f", "f", "4", "4", "4", "4", "4", "4", "4", "4", "4", "...
$ cty <int> 18, 21, 20, 21, 16, 18, 18, 18, 16, 20, 19, 15, 17, 17, 15, 15, 17, 16, 14, 11, 1...
$ hwy <int> 29, 29, 31, 30, 26, 26, 27, 26, 25, 28, 27, 25, 25, 25, 25, 24, 25, 23, 20, 15, 2...
$ fl <chr> "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "p", "...
$ class <chr> "compact", "compact", "compact", "compact", "compact", "compact", "compact", "com...
> summary(mpg)
manufacturer model displ year cyl trans
Length:234 Length:234 Min. :1.600 Min. :1999 Min. :4.000 Length:234
Class :character Class :character 1st Qu.:2.400 1st Qu.:1999 1st Qu.:4.000 Class :character
Mode :character Mode :character Median :3.300 Median :2004 Median :6.000 Mode :character
Mean :3.472 Mean :2004 Mean :5.889
3rd Qu.:4.600 3rd Qu.:2008 3rd Qu.:8.000
Max. :7.000 Max. :2008 Max. :8.000
drv cty hwy fl class
Length:234 Min. : 9.00 Min. :12.00 Length:234 Length:234
Class :character 1st Qu.:14.00 1st Qu.:18.00 Class :character Class :character
Mode :character Median :17.00 Median :24.00 Mode :character Mode :character
Mean :16.86 Mean :23.44
3rd Qu.:19.00 3rd Qu.:27.00
Max. :35.00 Max. :44.00
3. 수량형 변수의 분석
일변량 변수에 대한 통계량 분석은 크게 다음과 같다.
1) 데이터 분포 시각화
- ggplot() + geom_{histogram, desity}() 등
2) 요약 통계량 계산
- summary(), mean(), median(), var(), sd(), mad(), quantile() 등
3) 데이터 정규성 검사 : 분포가 정규분포와 얼마나 유사한지 검사
- qqplot, qqline() 등
4) 가설검정과 신뢰구간 : t-검정(평균에 대한 검정)과 신뢰구간
- t.test()
5) 이상점 찾아보기 : 로버스트 통계량 계산
1) 데이터 분포 시각화 및 요약 통계량 계산 : 다음 예로써 mpg$hwy 데이터의 구조와 내용을 파악해보자
> summary(mpg$hwy)
Min. 1st Qu. Median Mean 3rd Qu. Max.
12.00 18.00 24.00 23.44 27.00 44.00
> quantile(mpg$hwy)
0% 25% 50% 75% 100%
12 18 24 27 44
> opar = par(mfrow = c(2,2))
> hist(mpg$hwy)
> boxplot(mpg$hwy)
> qqnorm(mpg$hwy)
> qqline(mpg$hwy)
> par(opar)
★2) 가설검정과 신뢰구간 : 다음은 mpg$hwy 데이터를 통해 일변량 연속형 변수에 흔하게 사용되는 통계 추정 절차인 t-검정을 진행할 것이다.
* t-검정이란, 가설을 통해 귀무가설(효과가 없는 가설)인지 혹은 대립가설(효과가 있는 가설)인지 P값을 통해 통계적으로 판단하는 것이다.
따라서, 다음과 같은 단측 검정을 시행하고자 한다.
H0 : mu <= 22.9 vs. H1 : mu > 22.9
위의 가설검정을 t.test()를 통해 알아보자.
> hwy = mpg$hwy
> n=length(hwy)
> mu0 = 22.9
> t.test(hwy, mu=mu0, alternative = "greater")
One Sample t-test
data: hwy
t = 1.3877, df = 233, p-value = 0.08328
alternative hypothesis: true mean is greater than 22.9
95 percent confidence interval:
22.79733 Inf
sample estimates:
mean of x
23.44017
위의 파란색 글씨가 유의깊게 봐야할 부분이다.
1. 결과는 p값은 0.083이다.
-> 만약 모 평균 고속도로 연비가 22.9라면 우리가 관측한 것만큼 큰 표본 평균값과 t통계량(t=1.3877)이 관측될 확률이 8.3%라는 것이다.
따라서, 유의수준 alpha = 10%라면 고속도로 연비가 22.9보다 크다고 결론지을 수 있지만, 유의수준이 5%라면 고속도로 연비가 22.9보다 크다고
결론지을 만한 증거가 충분하지 않다고 할 수 있다.(즉, 귀무가설에 대한 증거가 불충분하므로 위 가설은 통계적으로 유의미하다.)
2. 신뢰구간은 어떨까?
-> t.test()에 alternative = "greater" 나 "less" 없이 실행하면 95% 신뢰구간을 구해준다.
> t.test(hwy)
One Sample t-test
data: hwy
t = 60.216, df = 233, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
22.67324 24.20710
sample estimates:
mean of x
23.44017
위의 출력결과에서 알 수 있듯 신뢰구간은 [22.6, 24.2]인 것으로 보아 다시 한번 22.9는 통계적으로 유의미한 가설인 것을 판단할 수 있다.
3) 이상점과 로버스트 통계 방법 : 다음은 이상점을 찾아내는 것이다. (나타나는 흔한 이유 중 하나는 데이터입력 오류이다.)
* 이상점이란? 여러 이유로 다른 관측치와 매우 다른 값을 가진 관측치. 즉, 다시 말해 군집 상태의 데이터와는 동 떨어진 데이터를 말한다.
이상점을 찾는 방법 중 하나는 상자 그림을 살펴보는 것이다. 상자 그림의 상자는 25% 백분위 수(Q1)와 75% 백분위 수(Q2)를 나타낸다. 상자 그림에서 이상점은 [Q1 - 1.5 * IQR, Q3 + 1.5* IQR] 구간 바깥의 값들이다.
* IQR : inter-quartile range(=Q3 - Q1)을 나타낸다.
로버스트 통계 방법은 이상점의 영향을 적게 받는 절차이다. 기본적으로 mean 대신 median, SD 대신 MAD(Median Absolute Deviance)를 사용하면 된다.
> c(mean(hwy), sd(hwy))
[1] 23.440171 5.954643
> c(median(hwy), mad(hwy))
[1] 24.000 7.413
4. 성공-실패값 범주형 변수의 분석
범주형 변수의 대표 틀은 '성공' or ' 실패' 범주만을 가진 데이터이므로 해당 범주만을 가진 데이터의 분석을 다룰 것이다.( ex) 여론조사의 '찬성' or '반대' )
이러한 데이터를 접하면 다음과 같이 분석해주면 된다.
1) 요약 통계량 계산
- table() : 도수 분포 계산
- xtabs() : 도수 분포 계산 및 클래스 속성을 가지게 됨
- prop.table() : 도수를 상대도수로 바꿔줌
2) 데이터 분포의 시각화
- barplot()
3) '성공률'에 대한 검정과 신뢰구간 계산
- binom.test()
다음 예로써 일상생활에서 가장 흔하게 접할 수 있는 여론조사를 예로 들자.
다음 데이터는 모집단, 즉 전국의 성인이 천만 명이라 할 때 랜덤하게 선정한 n=100명의 사람에게 '지지' or '지지하지 않음'을 택하도록 한 것이며, 다음 처럼 수집 될 것이다.
응답자 번호 응답
1 지지
... ...
n 지지하지 않음
참 지지율이 p=50%라 하고, 지지 여부 관측치는 0 = 'no', 1='yes'레벨을 가진 factor 변수로 생성하였다.
> set.seed(1606)
> n=100
> p=0.5
> x=rbinom(n,1,p)
> x=factor(x, levels = c(0,1), labels = c("no", "yes"))
> x
[1] yes yes yes yes no no yes no no yes yes no no yes no no yes yes yes no yes yes no no yes
[26] no no yes no no yes no yes yes no no yes yes no yes no no no yes yes no yes yes yes yes
[51] no yes no yes yes no no yes yes yes yes yes yes no no no yes no no yes no no yes no yes
[76] yes yes no yes no no no no yes yes yes yes no yes yes no no no no yes yes yes yes no yes
Levels: no yes
> table(x)
x
no yes
46 54
> prop.table(table(x)) x
no yes
0.46 0.54
> barplot(table(x)) 우리는 실제 '지지율'이 50%인 것을 알고 있지만, 그것을 모른다고 가정하고 다음 가설을 검정한다고 해보자.
H0 : P = 0.5 vs. H1 : P! =/= 0.5
R에서 binom.test()로 쉽게 신뢰구간과 가설검정을 시행 할 수 있다.
> binom.test(x=length(x[x=='yes']), n=length(x), p=0.5)
Exact binomial test
data: length(x[x == "yes"]) and length(x)
number of successes = 54, number of trials = 100, p-value = 0.4841
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.4374116 0.6401566
sample estimates:
probability of success
0.54
결과에서 보듯이 95% 신뢰구간은 [0.437 0.640]이다. 그리고 P값은 0.48로 주어졌다.
즉, 귀무가설 H0 : P = 0.5 가 참 일 때, 주어진 데이터만큼 극단적 데이터를 관측할 확률은 48%로 아주 높은 편이다.
따라서, 귀무가설을 기각할 증거는 희박하다고 할 수 있다.
5. 수량형 X(설명변수 X), 수량형 Y(반응변수 Y)의 분석
수량형 X 변수와 수량형 Y변수의 관계를 연구하는 방법의 절차는 다음을 따를 것을 권장한다.
1) 산점도를 통한 관계 모양 파악
- plot(), ggplot2의 geom_point()
- 중복치가 많을 땐 jitter 사용
- 데이터가 많을 땐 alpha = 옵션을 통해 표본화
- 관계가 선형인지, 강한지, 이상치가 있는지 파악
2) 상관계수 계산
- 선형 관계의 강도만 재는 것에 염두(선형적인 연관성이 있는지만)
★3) 선형 모형을 적합
- 모형의 적합도와 모수의 의미 파악
- 잔차의 분포 파악(잔차에 이상점이 있는지, 잔차가 예측변수에 따라 다른 분산을 갖지는 않는지)
4) 이상치가 있을 땐, 로버스트 회귀분석 사용
5) 비선형 데이터에는 LOESS 등의 평활법 사용(세부 내용 생략)
1) 산점도를 통한 관계 모양 파악 : 두 연속형 변수 X와 Y가 있을 때의 시각화는 산점도가 기본이다.
> ggplot(mpg, aes(cty, hwy)) + geom_jitter() + geom_smooth(method = "lm")
2) 상관계수 : cor()함수는 상관계수를 계산해주며, 기본적으로 피어슨 상관관계를 계산해준다
* 피어슨 상관관계 : https://namu.wiki/w/%EC%83%81%EA%B4%80%20%EA%B3%84%EC%88%98
피어슨 상관계수는 두 변량의 '선형'관계의 강도를 -1(완벽한 선형 감소 관계)에서 1(완벽한 선형 증가 관계) 사이의 숫자로 나타낸다. (0은 선형 관계가 없음을 나타냄)
하지만, 산점도를 그리지 않고 상관계수만 보는 것은 위험하다. 하단 그림의 가운데 줄에서 보듯 상관계수는 경사를 얘기해주지 않는다. 또한, 가장 하단의 그림은 다양한 비선형적 관계와 데이터의 군집성 등의 패턴은 전혀 잡아내지 못한다.
※ 상관계수는 이상치의 영향을 많이 받는다!!
상관계수는 이상치의 영향을 많이 받으므로 로버스트한 방법인 켄달의 타우나 스피어맨의 로 통계량을 계산하는 것도 좋은 방법이다.
* 켄달의 타우 상관계수 : https://thebook.io/006723/ch07/06/03/
* 스피어맨의 로 상관계수 : https://thebook.io/006723/ch07/06/02/
켄달의 타우와 스피어맨의 로를 계산하려면 다음처럼 cor() 함수의 method = "kendall" or "spearman" 옵션을 지정하면 된다.
> cor(mpg$cty, mpg$hwy)
[1] 0.9559159
> with(mpg, cor(cty,hwy))
[1] 0.9559159
> with(mpg, cor(cty, hwy, method = "kendall"))
[1] 0.8628045
> with(mpg, cor(cty, hwy, method = "spearman"))
[1] 0.9542104
★3) 선형회귀 모형 적합 : 하나의 설명변수와 하나의 반응변수로 이루어진 모형을 단순 회귀분석 모형 이라고 한다.
lm()과 summary.lm()함수는 선형 모형을 최소제곱법으로 추정한다. 즉, 잔차의 제곱합을 최소화하는 문제를 풀어서 추정치를 구한다.
lm.summary()은 각 추정치와 더불어 각 모수값이 0인지 아닌지에 대한 가설검정 결과를 보여준다.
> hwy_lm = lm(hwy ~ cty, data = mpg)
> summary(hwy_lm)
Call:
lm(formula = hwy ~ cty, data = mpg)
Residuals:
Min 1Q Median 3Q Max
-5.3408 -1.2790 0.0214 1.0338 4.0461
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.89204 0.46895 1.902 0.0584 .
cty 1.33746 0.02697 49.585 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.752 on 232 degrees of freedom
Multiple R-squared: 0.9138, Adjusted R-squared: 0.9134
F-statistic: 2459 on 1 and 232 DF, p-value: < 2.2e-16
* Estimate Std. : 추정값, Pr(>ltl) : P값, R-squared : 같은 데이터를 가지고 여러가지 유형을 다룰 때 유용한 값
즉, 적합된 모형은 다음과 같다
hwy = 0.892 + 1.337 * cty
hwy값에 대한 cty의 선형 효과는 통계적으로 유의하다.
추정값은 1.337, 표준편차는 0.027, t-값은 49.5, P값은 2e-16 즉, 실질적으로 0이다.
만약 hwy가 cty와 선형 관계가 전혀 없다면 49.5라는 큰 t-값을 관측하는 것은 실질적으로 불가능하다. 따라서 이 데이터는 귀무가설이 참이 아니라는 증거가 된다.
4) 선형회귀모형 예측 : 선형회귀모형을 예측해보자.
predict()함수는 반응변수의 예측값을 계산한다. 디폴트는 적합에 사용한 데이터로 계산하며, 다른 데이터를 제공하려면 newdata = 옵션을 사용한다.
resid()함수는 잔차를 계산한다.
* 잔차 : 잔차는 관측값과 최확값(참값을 대신하는 값)의 차이
다음 예제는 데이터의 추정량과 잔차뿐 아니라 새로운 데이터 cty=10, 20, 30에 대한 값에 대한 예측값, 그리고 예측오차를 계산한다.
> predict(hwy_lm)
결과 생략
> resid(hwy_lm)
결과 생략
> predict(hwy_lm, newdata = data.frame(cty=c(10,20,30)))
1 2 3
14.26660 27.64115 41.01571
> predict(hwy_lm, newdata = data.frame(cty=c(10,20,30)),se.fit = TRUE) #se.fit : standard error
$`fit`
1 2 3
14.26660 27.64115 41.01571
$se.fit
1 2 3
0.2176003 0.1424778 0.3725052
$df
[1] 232
$residual.scale
[1] 1.752289
5) 선형회귀 모형의 가정 진단 : 선형회귀 결과가 의미 있으려면 다음 조건이 충족되어야 한다.
① x와 y의 관계가 선형이다
② 오차의 분포가 독립이다.
③ 잔차의 분포가 동일하다.
④ 잔차의 분포가 N(0,s^2)이다.
R에서 lm클래스 객체에 plot()을 실행하면 내부적으로 plot.lm()이 실행된다.
이것은 선형 모형 진단 플롯들을 그려주며, 예제는 다음과 같다.
> class(hwy_lm)
> opar = par(mfrow = c(2,2), oma = c(0,0,1.1,0))
> plot(hwy_lm, las = 1) #Residuals, Fitted...
> par(opar)
1번 그림 : 잔차의 분포가 예측값과 독립적인지 확인하는데 사용
2번 그림 : 표준화된 잔차의 분포가 정규분포에 가까운지 확인
3번 그림 : 잔차의 절댓값의 제곱근과 예측값 사이의 관계
4번 그림 : 표준화된 잔차와 Leverage 간의 관계
* Leverage : x 값이 전체적으로 형성된 그룹에서 멀리 떨어진 것을 말함
6) 로버스트 선형회귀분석 : 목표는 같지만 이상치가 있다면 "피해가자" 라는 분석 기법
선형회귀분석의 모수의 추정값은 잔차의 분포에서 이상치가 있을 때 지나치게 민감하게 반응한다.(이유는 모수를 추정할 때 최소제곱 방법을 사용하기 때문)
이상치에 민감하지 않은 추정 방법이 필요할 때는 로버스트 회귀분석 방법론을 사용한다.
다양한 로버스트 선형회귀분석 방법 중 한 가지 예로 MASS 패키지의 lqs() 함수가 있다.
MASS::lqs()는 데이터 중 "좋은" 관측지만 적합에 사용한다.(n관측치와 p 독립변수가 있을 때 "lqs"는 제곱오차의 floor((n+p+11)/2) quantule을 최소화한다.)
실행 예는 ?lqs의 예에서 따왔다.
사용한 데이터는 암모니아를 질산으로 산화하는 화학공정에 21가지 다른 냉각기류, 냉각수 온도, 산농도 값에 따른 에너지 손실로 이루어진 "stackloss" 데이터이다.
> ibrary(MASS)
> set.seed(123) #make reproducible
> lqs(stack.loss ~ ., data = stackloss) #로버스트
Call:
lqs.formula(formula = stack.loss ~ ., data = stackloss)
Coefficients:
(Intercept) Air.Flow Water.Temp Acid.Conc.
-3.631e+01 7.292e-01 4.167e-01 -8.131e-17
Scale estimates 0.9149 1.0148
로버스트 방법인 lqs()로부터 얻은 위의 결과를 최소제곱법 방법인 lm()으로 계산된 아래의 결과와 비교해보자
> lm(stack.loss ~ ., data = stackloss) #보통 선형모형
Call:
lm(formula = stack.loss ~ ., data = stackloss)
Coefficients:
(Intercept) Air.Flow Water.Temp Acid.Conc.
-39.9197 0.7156 1.2953 -0.1521
두 방법이 상당이 큰 차이를 보이는 것을 알 수 있다. 특히, lm()의 결과에 비해 lqs()방법에 의한 결과에서는 산농도 변수의 효과는 거의 없다.
6. 범주형 x, 수량형 y
mpg 데이터에서 5개의 차종 간에 연비의 차이가 있는지
3개의 혈압약 간에 혈압의 감소량이 차이가 있는지 등의 문제를 생각해보면
범주형 설명변수는 차종과 3개의 혈압약이 되고
수량형 반응변수는 연비, 혈압감소 등이 될 것이다.
이러한 데이터를 분석 할 때는 다음 방법들을 사용하면 된다.
1) 병렬상자그림을 이용하여 데이터 시각화
2) lm()함수로 ANOVA 선형 모형을 적합
3) plot.lm()으로 잔차의 분포를 확인
★ 1) 분산 분석 : 설명변수가 범주형이고, 반응변수가 수량형일 경우 분산분석(ANOVA)를 사용한다.
mpg 데이터를 통해 예를 살펴보자.
> mpg %>% ggplot(aes(class, hwy)) + geom_boxplot()
분산분석은 lm()함수로 간단히 실행 가능하다.
> hwy_lm2 = lm(hwy ~ class, data = mpg)
> summary(hwy_lm2)
Call:
lm(formula = hwy ~ class, data = mpg)
Residuals:
Min 1Q Median 3Q Max
-8.1429 -1.8788 -0.2927 1.1803 15.8571
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.800 1.507 16.454 < 2e-16 ***
classcompact 3.498 1.585 2.206 0.0284 *
classmidsize 2.493 1.596 1.561 0.1198
classminivan -2.436 1.818 -1.340 0.1815
classpickup -7.921 1.617 -4.898 1.84e-06 ***
classsubcompact 3.343 1.611 2.075 0.0391 *
classsuv -6.671 1.567 -4.258 3.03e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.37 on 227 degrees of freedom
Multiple R-squared: 0.6879, Adjusted R-squared: 0.6797
F-statistic: 83.39 on 6 and 227 DF, p-value: < 2.2e-16
위의 분산분석의 결과는 앞서 회귀분석에서 lm()함수를 실행했을 때와 결과가 비슷하다 우선 모수의 추정치로부터, 적합된 모형은 다음과 같다.
★ yhat = 24.8 + 3.50 * 1(class=compact) + 2.49 * 1(class=midsize)+...-6.67*1(class=suv) ★
P-값이 0이면 ***, 0.001보다 작으면 **, 0.1보다 작으면 .으로 나타내는 식이다. 쉽게, classpickup과 classsuv가 다른 집단보다는 유의하게 다른 평균 연비를 가지고 있음을 알 수 있다.
위 결과에서 classpickup 집단의 효과를 설명해보자
1. 모수 추정값 -7.92
- class=pickup이면 연비가 평균 7.92만큼 감소
2. 표준오차 1.62
- 95% 신뢰구간은 -7.92 + -1.96 * 1.62, 즉[-11.1, -4.74] 정도
3. H0 : βpickup = 0 vs H1 not H0 의 가설에 대한 P값은 거의 0
4. Multiple R-squared(모형의 적합도) : 0.6879
- 연비의 총 변동량 중 자동차 클래스로 설명되는 변동량의 비율이 69%
5. Adjusted R-squared : 0.6797
- 총 6개의 모수를 사용했음을 감안하여 그 비율을 68%로 줄여서 보고한다는 뜻
6. F-statistic: 83.39 & p-value: < 2.2e-16
- H0 : 연비에 클래스 집단의 효과가 전혀 없다 vs. H1: not H0 에 대한 검정결과
새로운 x값에 대한 예측/추정값은 predict() 함수로 얻을 수 있다.
> predict(hwy_lm2, newdata = data.frame(class="pickup"))
1
16.87879
2) 분산분석의 진단 : 분산분석 결과가 의미 있기 위해서는 다음 가정이 충족되어야 한다.
① 잔차의 분포가 독립이다.
② 잔차의 분산이 동일이다.
③ 잔차의 분포가 N(0,s^2)이다.
회귀분석과 마찬가지로 가장 중요한 것은 분포의 독립성과 이상치 유무이다.
> opar = par(mfrow = c(2,2), oma = c(0,0,1.1,0))
> plot(hwy_lm2, las=1)
> par(opar)
7. 수량형 x, 범주형 y(성공 - 실패)
이러한 데이터를 분석할 때는 다음의 방법들을 사용하면 된다.
1) X와(jitter된) Y 변수의 산점도를 그리고 Y변수의 그룹별로 X 변수의 병렬 상자그림을 그려본다.
- Y변수 값에 따라 X 변수의 분포가 차이가 있는지 확인
- 두 변수 간에 어떤 관계가 있는지 확인
- 이상치가 존재하는지 확인
- 표본 로그오즈와 x의 산점도에 선형 패턴이 있는 확인
2) glm() 함수로 일반화 선형 모형을 적합
- summary.glm()
3) plot.glm()으로 잔차 분포 확인
- 이상점은 없는지 확인
- 모형의 가정은 만족하는지 확인
성공-실패 범주형 y 변수와 수량형 설명변수를 가진 데이터는 전통적인 선형 모형으로 다룰 수 없다.
그 이유는, 전통적인 선형 모형은 반응변수 y의 범위가 무한대이기 때문이다.
대신, 일반화 선형 모형(Generalized Linear Model, GLM)을 사용하여야 한다.
1) 챌린저 데이터 분석
이 데이터 분석 기법에 대해서는 챌린저 O링 데이터를 사용할 것이다. 1986년 쏘아올린 챌린저 호 우주선 폭발 사건의 데이터인데, O링 부품이 망가진 것이 폭발의 원인이라 성공-실패 여부를 데이터로 만들어 놓은 것이다.
> chall <- read.csv('https://raw.githubusercontent.com/stedy/Machine-Learning-with-R-datasets/master/challenger.csv')
> chall <- tbl_df(chall)
> glimpse(chall)
Observations: 23
Variables: 5
$ o_ring_ct <int> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6
$ distress_ct <int> 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 1
$ temperature <int> 66, 70, 69, 68, 67, 72, 73, 70, 57, 63, 70, 78, 67, 53, 67, 75, 70, 81, 76, 79, 75...
$ pressure <int> 50, 50, 50, 50, 50, 50, 100, 100, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200...
$ launch_id <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23
여러 변수들 중 온도와 O링 실패변수만 고려하자. 이 데이터를 시각화하면 온도와 실패한 O링 개수 간의 의미 있는 관계를 관찰할 수 있다.
> library(gridExtra)
> p1= chall %>% ggplot(aes(temperature, distress_ct)) + geom_point()
> p2= chall %>% ggplot(aes(factor(distress_ct), temperature)) + geom_boxplot()
> grid.arrange(p1,p2,nrow=1,ncol=2)
이제 glm() 함수로 모형을 추정해보자. glm()함수로 로지스틱 모형을 적합하는 코드는 다음과 같다.
> chall_glm = glm(cbind(distress_ct, o_ring_ct - distress_ct)~ + temperature, data = chall, family = 'binomial')
> chall_glm
Call: glm(formula = cbind(distress_ct, o_ring_ct - distress_ct) ~ +temperature,
family = "binomial", data = chall)
Coefficients:
(Intercept) temperature
8.8169 -0.1795
Degrees of Freedom: 22 Total (i.e. Null); 21 Residual
Null Deviance: 20.71
Residual Deviance: 9.527 AIC: 24.87
> summary(chall_glm)
Call:
glm(formula = cbind(distress_ct, o_ring_ct - distress_ct) ~ +temperature,
family = "binomial", data = chall)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.7526 -0.5533 -0.3388 -0.1901 1.5388
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 8.81692 3.60697 2.444 0.01451 *
temperature -0.17949 0.05822 -3.083 0.00205 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 20.706 on 22 degrees of freedom
Residual deviance: 9.527 on 21 degrees of freedom
AIC: 24.865
Number of Fisher Scoring iterations: 6
temperature 변수의 효과는 p값이 0.002로 유의미한 것으로 밝혀졌다.
temperature의 효과는 -0.179는 온도가 1도 상승할 때 속칭 로그 오즈비가 0.179만큼 감소한다는 것이다. (무슨 말인지 이해가 안가면 다음에 나올 시각화 파트에서 집중하기 바란다.)
2) 로지스틱 모형 예측, 링크와 반응 변수
만약 temperature = 30(화씨) 이라면 성공확률은 어떨까? predict.glm()으로 예측해보자
> predict(chall_glm, data.frame(temperature=30))
1
3.432159
결과는 predict.glm()은 type=c("link", "response", "terms") 옵션이 있고, 디폴트인 "link"는 선형예측값 Xβhat을 출력하기 때문이다.
따라서, 위의 결과는 β0hat + β1hatx = 8.82 - 0.179 * 30 값이다.
선형예측값이 아닌 확률값을 얻으려면 이 값을 로지스틱 변환 or predict.glm()에서 type="response" 옵션을 사용하면 된다.
> predict(chall_glm, data.frame(temperature = 30), type = 'response')
1
0.9686946
3) 로지스틱 모형 적합결과의 시각화 : 시각화를 위해서는 베이스 그래픽을 사용 할 수 있다.
> logistic <- function(x){exp(x)/(exp(x)+1)}
> plot(c(20,85), c(0,1), type = "n", xlab = "temperature", ylab = "prob")
> tp <- seq(20, 85, 1)
> chall_glm_pred <- predict(chall_glm, data.frame(temperature = tp), se.fit = TRUE)
> lines(tp, logistic(chall_glm_pred$fit))
> lines(tp, logistic(chall_glm_pred$fit - 1.96 * chall_glm_pred$se.fit), lty=2)
> lines(tp, logistic(chall_glm_pred$fit + 1.96 * chall_glm_pred$se.fit), lty=2)
> abline(v=30, lty=2, col='blue')
※ 선형 분석, 분산 분석, GLM 모형(로지스틱) 분석은 데이터 종류에 따라 달라지는 분석 기법이다. 어떤 데이터에 어떤 분석 기법을 선택할지 파악하여 분석하는 것은 매우 중요하다.
* 참고문헌 : 실리콘밸리 데이터 과학자가 알려주는 따라하며 배우는 데이터 과학(저자 : 권재명)
'Data Science > Data Science in R' 카테고리의 다른 글
[Data Science] Adult 데이터로 알아보는 분류분석 모형 개념 (1) 2018.12.15 [Data Science] 분류 분석의 기본 개념 (0) 2018.12.15 [Data Science 기초] IMDB 영화 정보 데이터 시각화 및 분석 (0) 2018.10.24 [Data Science 기초] 데이터 시각화 응용 - 변수의 종류에 따른 시각화 (0) 2018.10.23 [R 기초] 데이터 시각화 기초 - Base R package vs. ggplot2 (0) 2018.10.23
