[Data Science] 데이터 종류에 따른 분석 기법

목표 : 데이터 종류에 따라 분석하는 기법들을 익혀보자.

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 모형(로지스틱) 분석은 데이터 종류에 따라 달라지는 분석 기법이다. 어떤 데이터에 어떤 분석 기법을 선택할지 파악하여 분석하는 것은 매우 중요하다.

* 참고문헌 : 실리콘밸리 데이터 과학자가 알려주는 따라하며 배우는 데이터 과학(저자 : 권재명)