36.일반선형모델 - 포아송회귀분석.Rmd
0.02MB
36.일반선형모델---포아송회귀분석.html
1.95MB

 

https://www.youtube.com/watch?v=-Za5MrWUohg&list=PLY0OaF78qqGAxKX91WuRigHpwBU0C2SB_&index=36

  • 목표: 포아송회귀분석 이해

1. 포아송회귀분석 개요

  • 포아송회기분석은 결과변수가 특정기간동안의 사건발생 횟수인 경우에 적용할 수 있음
  • 1시간동안 걸려오는 상담전화수, 하루동안 발생하는 범죄횟수, 한 달동안 발생하는 교통사고횟수, 한해동안 일어나는 지진건수 등은 모두 포아송회귀분석으로 사건발생 횟수를 예측 가능
  • 사건발생 횟수는 일반적으로 제한된 값을 가지며 음수일 수 없음
  • 사건발생 횟수는 정규분포를 따르지 않기 때문에 표준적인 선형회귀분석으로 예측할 수는 없음

2. 포아송회귀모델

  • 포아송회귀모델은 결과변수 y가 포아송분포를 따른다고 가정하며 아래 선형모델로 정의 됨

\[ln(\lambda) = \beta_0 + \beta_0x_2 + \dots + + \beta_mx_m\]

  • 여기서 \(\lambda\)는 결과변수 y의 평균을 나타냄
  • 결과변수 \(\lambda\)에 로그를 취한 \(ln(\lambda)\)는 예측변수들의 선형결합으로 나타낼 수 있음
  • 포아송회귀모델을 일반선형모델 틀로 보면, 링크함수는 \(ln(\lambda\),) 결과변수 y의 확률분포는 포아송분포임

3. 포아송회기분석 수행 절차

  • 분석대상 데이터(breslow.dat) 설명
    • 여기서는 robust패키지 breslow.dat를 사용함
    • breslow.dat데이터셋은 뇌전증환자의 치료제 투약 전 8주간의 발작일수와 투여 후 8주간의 발작횟수가 기록돼 있음
    • 뇌전증 치료제가 투약 후 8주동안 발생하는 발작횟수에 미치는 영향을 포아송회귀분석으로 분석할 것임
library(robust)
data("breslow.dat")
str(breslow.dat)
## 'data.frame':    59 obs. of  12 variables:
##  $ ID   : int  104 106 107 114 116 118 123 126 130 135 ...
##  $ Y1   : int  5 3 2 4 7 5 6 40 5 14 ...
##  $ Y2   : int  3 5 4 4 18 2 4 20 6 13 ...
##  $ Y3   : int  3 3 0 1 9 8 0 23 6 6 ...
##  $ Y4   : int  3 3 5 4 21 7 2 12 5 0 ...
##  $ Base : int  11 11 6 8 66 27 12 52 23 10 ...
##  $ Age  : int  31 30 25 36 22 29 31 42 37 28 ...
##  $ Trt  : Factor w/ 2 levels "placebo","progabide": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Ysum : int  14 14 11 13 55 22 12 95 22 33 ...
##  $ sumY : int  14 14 11 13 55 22 12 95 22 33 ...
##  $ Age10: num  3.1 3 2.5 3.6 2.2 2.9 3.1 4.2 3.7 2.8 ...
##  $ Base4: num  2.75 2.75 1.5 2 16.5 6.75 3 13 5.75 2.5 ...
  • breslow.dat데이터셋 변수 중 분석에 필요한 4개 변수만 추출한 별도 데이터셋으로 분석수행
seizure <- breslow.dat[c("Base", "Age", "Trt", "sumY")]
  • 앞의 세 개 변수, Base·Age·Trt변수는 예측변수, 네번째 변수 sumY는 결과변수로서 사용함
  • 변수설명
    • 결과변수로서 사용할 sumY은 치료 시작 후 8주 동안의 발작횟수를 나타냄
    • Base변수는 치료 시작 전 8주간의 발작횟수를 나타냄
    • Age변수는 환자의 나이를 나타내고
    • Trt변수는 치료제로서 2개의 범주를 가짐. placebo은 위약, progabide은 진약을 나타냄
    • Base·Age·Trt 이 세 개의 예측변수가 결과변수인 sumY, 즉 치료 시작 후 8주 동안의 발작횟수에 미치는 영향을 포아송회귀분석으로 분석할 것임
  • 먼저 결과변수와 예측변수에 요약통계량을 살펴볼 것임
summary(seizure)
##       Base             Age               Trt          sumY       
##  Min.   :  6.00   Min.   :18.00   placebo  :28   Min.   :  0.00  
##  1st Qu.: 12.00   1st Qu.:23.00   progabide:31   1st Qu.: 11.50  
##  Median : 22.00   Median :28.00                  Median : 16.00  
##  Mean   : 31.22   Mean   :28.34                  Mean   : 33.05  
##  3rd Qu.: 41.00   3rd Qu.:32.00                  3rd Qu.: 36.00  
##  Max.   :151.00   Max.   :42.00                  Max.   :302.00
  • 발작횟수를 나타내는 Base변수와 sumY변수 요약통계량에서 중위수가 평균에 비해 모두 작음
  • 중위수가 평균에 비해 작기 때문에 발작횟수는 오른쪽으로 꼬리가 긴 분포를 갖는 것으로 예상할 수 있음
  • 히스토그램으로 그러한 분포를 갔는지 시각적으로 확인할 것임
hist(seizure$sumY, breaks = 20, col="cornflowerblue", 
     xlab="Seizure Count", 
     main="Distribution of Seizures")
  • 결과
    • 발작횟수는 실제 오른쪽으로 꼬리가 긴 편향된 분포를 가짐
    • 이러한 오른쪽으로 꼬리가 긴 편향된 분포를 갖고 있기 때문에 결과변수가 정규분포를 따를 것을 요구하는 표준적인 선형회귀분석을 적용할 수 없음
    • 이처럼 결과변수가 정규분포에서 요구하는 형태를 따르지 않을 때 일반선형모델 범주에 속하는 포아송 회귀분석은 좋은 대안이 될 수 있음
  • 포아송회귀모델은 일반선형모델의 범주에 속하기 때문에 로지스틱회귀분석을 수행하는 glm()를 포아송회귀분석에서도 동일하게 이용할 수 있음
  • glm() 함수 인수설명
    • 1st: 결과변수와 예측변수 간의 관계를 포뮬러 형식으로 지정함
    • data: 예측모델 개발에 사용할 데이터셋을 지정함
    • family: 결과변수의 확률분포함수를 지정함. 여기서는 결과변수가 포아송분포를 따르는 것으로 가정하기 때문에 poisson함수를 지정함
    • 분석결과는 summary()를 로 확인할 수 있음

3-1. 회귀모델 구축 + 회귀계수 유의성 검정

seizure.poisson <- glm(sumY ~ Base + Age + Trt, 
                       data=seizure, family=poisson())
summary(seizure.poisson)
## 
## Call:
## glm(formula = sumY ~ Base + Age + Trt, family = poisson(), data = seizure)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.0569  -2.0433  -0.9397   0.7929  11.0061  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.9488259  0.1356191  14.370  < 2e-16 ***
## Base          0.0226517  0.0005093  44.476  < 2e-16 ***
## Age           0.0227401  0.0040240   5.651 1.59e-08 ***
## Trtprogabide -0.1527009  0.0478051  -3.194   0.0014 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2122.73  on 58  degrees of freedom
## Residual deviance:  559.44  on 55  degrees of freedom
## AIC: 850.71
## 
## Number of Fisher Scoring iterations: 5
  • 결과
    • 포아송회귀분석의 결과로써 Estimate(회귀계수), Std. Error(표준오차), Pr(>|z|)(회귀계수의 유의성 검정결과가 출력됨)
    • 여기서 위의 세 예측변수는 유의수준 0.05에서 모두 통계적으로 유의함

3-2. 회귀모델 유의성 검정

  • 회귀모델 자체 유의성 검정은 출력 결과의 아랫부분, deviance를 이용해 할 수 있음
  • deviance(이탈도) 설명
    • deviance는 이탈도라고 번역해서 얘기하기도 함
    • 이탈도는 예측모델의 적합도 정도를 나타낸 지표로, 엄밀히는 비적합도를 나타낸 지표임
    • 왜냐하면 이탈도가 클수록 예측모델은 적합도가 좋지 않다는 것을 나타내기 때문임
    • Null deviance는 예측모델의 상수항만을 포함하는 모델인 Null 모델의 이탈도를 나타냄
    • Residual deviance는 현재 구축한 3개의 예측변수를 포함하고 있는 예측모델의 이탈도를 나타냄
    • 회귀모델의 유의성 검정은 Null 모델 이탈도와 구축한 현재 모델 이탈도 간의 차이가 통계적으로 유의한 지를 검정함으로써 수행할 수 있음
    • 이탈도는 \(\chi^2\)분포를 따르기 때문에 \(\chi^2\)분포 하에서 이탈도의 차이가 통계적으로 유의한 지 검정함
    • 이는 pchisq()로 검정할 수 있음
    • pchisq() 인수설명
      • q: 2개의 이탈도의 차이값 지정
      • df: 대응되는 자유도의 차이값 지정
      • lower.tail: \(\chi^2\)분포 상 오른 쪽 꼬리 부분 면적이 유의확률이기 때문에 FALSE로 지정
pchisq(q=2122.73-559.44, df= 55-53, 
       lower.tail=FALSE)
## [1] 0
  • 결과
-   유의확률은 0 으로 산출 됨

-   상수항만을 포함하고 있는 Null 모델의 이탈도와 추가로 3개의 예측변수를 투입하여 현재 구축한 예측모델의 이 이탈도 간 차이는 통계적으로 유의함

-   따라서 추가로 투입한 3개의 예측변수는 모델의 적합도 개선에 통계적으로 의미있는 기여를 한다고 볼 수 있고, 이는 3개의 예측변수를 포함하고 있는 이 모델이 통계적으로 유의하다는 것을 의미

3-3. 포아송회귀모델에서 결과변수 해석 - 로그값을 지수함수로 변환해 해석

  • 포아송회귀모델에서 결과변수는 로그값으로 표현됨
    • 나이(Age) 회귀계수 0.023은 나이 1살 증가할 때 발작횟수의 로그값은 0.023만큼 커진다는 의미임
    • 그런데 로그값은 원래 관측값과 척도가 달라 해석상 어려움이 있음
    • 따라서 로그값보다는 일반적으로 결과변수의 원래 척도로 예측변수의 회귀계수를 해석함
      • 회귀계수에 지수함수를 취하면 원래 척도를 바탕으로 회귀계수를 해석할 수 있음
      • 지수함수로 변환한 나이의 회귀계수를 계산하면 아래와 같음
exp(coef(seizure.poisson))
##  (Intercept)         Base          Age Trtprogabide 
##    7.0204403    1.0229102    1.0230007    0.8583864
  • 결과
    • 지수함수로 변환한 나이의 회귀계수는 1.0 이상임
    • 나이 1살 증가는 발작횟수를 1.023배 증가시킨다는 의미임
    • 나이 1살 증가는 발작횟수를 2.3% 증가시킴
    • 이는 나이가 많아질수록 발작횟수도 증가한다는 것을 의미함
    • 치료제 유형 1단위의 증가는 발작횟수를 0.858배 증가시킴
    • 치료제 유형 1단위의 증가는
      • Trt변수값이 기준범주인 1에서 2로 변화한다는 것을 나타냄
      • 이는 위약에서 진약으로의 이동을 의미함
      • 즉 뇌전증 치료제를 처방받은 환자집단은 위약복용 환자집단에 비해 발작횟수가 14.2% 감소함

4. 포아송회귀분석에서 과산포 문제

4-1. 포아송회귀분석에서 과산포 문제 설명

  • 포아송분포에서는 평균과 분산이 같음
  • 포아송회귀분석에서 과산포 문제 발생 조건
    • 결과변수의 실제 관측된 분산이 포아송분포의 의해 기대되는 분산보다 더 클 때 발생함
    • 즉 분산 대 평균의 비율이 1보다 클 때 발생함
    • 사건발생 횟수데이터에서는 과산포가 빈번히 발생하고 과산포는 잘못된 분석결과를 유도할 수 있기 때문에 주의 깊게 다룰 필요가 있음
  • 과산포 발생으로 인한 표준오차, t값, 유의확률의 왜곡
    • 데이터셋에 과산포가 존재하면 표준오차가 매우 작아지게 됨
    • 회귀계수의 유의성 검정을 위해 사용하는 검정통계량 t값은 회귀계수를 표준오차로 나눠 계산하기 때문에 표준오차가 매우 작아지게 되면 t값은 매우 커지게 되고
    • 이는 유의확률을 지나치게 작게 만들어 회귀계수가 0이라는 귀무가설을 쉽게 기각하게 만듦
  • 과산포는 종종 상태의존성이라는 현상으로 인해 발생함
    • 포아송분포는 각 사건이 독립적이고 사건발생률이 일정하다고 가정함
      • 현재 사용하고 있는 뇌전증 발작데이터에서 이는 이미 환자에 대해서 발생하는 발작은
        • 서로 독립적이고(즉, 발작의 발생은 다음에 발생하는 다른 발작들과 관련성이 없고),
        • 발작이 발생하는 간격은 일정하다는 것을 의미함
    • 하지만 이러한 가정은 종종 현실적이지 못함
      • 이는 어떤 한 환자에 대해 1번째 발작 발생확률은 이미 20번째 발작을 겪은 상태에서 21번째 발작이 일어날 확률과 같을 가능성은 없기 때문임
      • 따라서 관측값의 변동성은 기대한 것보다 더욱 커지게 됨
  • 로지스틱회귀분석에서처럼 잔차이탈도 대 잔차자유도의 비율이 1을 크게 상회하면 과산포를 의심함
    • 아래는 뇌전증 발작데이터에서의 과산포 비율을 계산한 것임
deviance(seizure.poisson) / df.residual(seizure.poisson)
## [1] 10.1717
  • 결과
    • 과산포 비율이 10.17로 1보다 많이 큰 것으로 나타남
    • 따라서 과산포 가능성은 매우 높은 것으로 생각됨

4-2. 포아송회귀분석에서 과산포가능성의 통계적 검정

  • qcc패키지 qcc.overdispersion.test()로 과산포 가능성을 통계적으로 검정할 수 있음
  • qcc.overdispersion.test() 인수설명
    • 1st: 과산포를 검정하고자 하는 데이터를 지정함. 여기서는 결과변수인 발작횟수를 지정함
    • type: 확률분포를 지정함. 여기서는 포아송을 지정함
    • qcc.overdispersion.test()는 ‘데이터가 과산포가 아니다’ 라는 귀무가설을 검정함
library(qcc)
qcc.overdispersion.test(seizure$sumY, type="poisson")
##                    
## Overdispersion test Obs.Var/Theor.Var Statistic p-value
##        poisson data          62.87013  3646.468       0
  • 결과
    • 유의확률 0으로서 유의수준 0.05보다 작기 때문에 귀무가설을 기각하고
    • 이 데이터는 과산포문제를 갖고 있는 것으로 판단해 볼 수 있음

4-3. 과산포 문제를 지닌 데이터에 대한 포아송회귀분석 수행

  • 과산포 문제가 존재하더라도 여전히 glm()로 포아송회귀분석을 수행할 수 있음
  • family인수에 poisson()을 지정하는 대신에 quasipoisson()을 지정함
seizure.qpoisson <- glm(sumY ~ Base + Age + Trt, 
                       data=seizure, family=quasipoisson())
summary(seizure.poisson)
## 
## Call:
## glm(formula = sumY ~ Base + Age + Trt, family = poisson(), data = seizure)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.0569  -2.0433  -0.9397   0.7929  11.0061  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   1.9488259  0.1356191  14.370  < 2e-16 ***
## Base          0.0226517  0.0005093  44.476  < 2e-16 ***
## Age           0.0227401  0.0040240   5.651 1.59e-08 ***
## Trtprogabide -0.1527009  0.0478051  -3.194   0.0014 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2122.73  on 58  degrees of freedom
## Residual deviance:  559.44  on 55  degrees of freedom
## AIC: 850.71
## 
## Number of Fisher Scoring iterations: 5
summary(seizure.qpoisson)
## 
## Call:
## glm(formula = sumY ~ Base + Age + Trt, family = quasipoisson(), 
##     data = seizure)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -6.0569  -2.0433  -0.9397   0.7929  11.0061  
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.948826   0.465091   4.190 0.000102 ***
## Base          0.022652   0.001747  12.969  < 2e-16 ***
## Age           0.022740   0.013800   1.648 0.105085    
## Trtprogabide -0.152701   0.163943  -0.931 0.355702    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 11.76075)
## 
##     Null deviance: 2122.73  on 58  degrees of freedom
## Residual deviance:  559.44  on 55  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
  • 결과:
    • 과산포를 고려한 포아송회귀 분석결과를 살펴보고 회귀계수/표준오차를 종전모델과 비교할 것임
    • 회귀계수: 종전의 모델과 동일함
    • 표준오차: 종전보다 매우 커졌음
      • 따라서 이 경우는 더 커진 표준오차로 인해
      • 나이와 치료제에 대한 유의확률이 유의수준 0.05보다 커져서
      • 회귀계수가 0이라는 귀무가설을 기각하지 못함
    • 검정통계량 t값은 회귀계수를 표준오차로 나눠 구할 수 있음. 그래서 지금처럼 표준오차가 커지면 t값이 작아지게 되고 이는 유의확률 p값을 증가시킴
    • 과산포를 고려할 때,
      • 뇌전증 치료제 복용 환자집단이 위약을 처방받은 환자집단보다 발작횟수가 더 적다고 할 만한 증거가 충분하지 못함
      • 나이 또한 발작횟수의 유의확률에 영향을 미치지 못함
    • 따라서 여기에서는 과거의 발작횟수만이 향후 발작횟수를 예측하는데 유의한 영향을 미치는 것으로 나타남

5. 관측값 시간 간격이 다른 경우의 포아송회귀분석

  • 지금까지 다른 포아송회귀분석에서는 결과변수로 사용한 사건발생횟수를 하나의 고정된 기간동안 측정하는 것으로 제한하였음
    • 뇌전증 발작데이터에서 각 관측값은 모두 8주 동안의 동일한 기간에서의 발작횟수로 측정됨
    • 발작횟수를 측정한 시간의 길이는 모든 환자에 대해서 동일함
  • 하지만 많은 경우, 각 관측값에 대한 시간의 길이가 서로 다를 수 있음
    • 어떤 환자는 8주간의 기간동안 발작횟수를 측정할 수 있고
    • 어떤 환자는 12주간의 기간동안 발작횟수를 기록할 수 있음
  • 기존 포아송회귀분석을 확장해 이처럼 각 관측값 시간간격이 다른 경우의 통계분석을 수행할 수 있음
    • 이 경우 결과변수는 사건발생횟수를 단위시간으로 나눈 사건발생률이 됨
    • 즉, 단위시간당 발생횟수가 분석 대상이 됨
  • 사건발생률을 분석하기 위해 각 관측값에 대해 얼마의 기간동안 사건발생횟수를 측정했는지를 나타내는 시간변수가 포아송회귀모델에 포함되어야 함
    • 시간이 포함된 포아송회귀모델은 이와 같은 수식으로 표현될 수 있음\[ln(\lambda) = \beta_0 + \beta_0x_2 + \dots + + \beta_mx_m\]
    • \[ln(\frac{\lambda}{time}) = \beta_0 + \beta_0x_2 + \dots + + \beta_mx_m\]
    • 기존의 포아송회귀모델에, 시간변수에 로그를 취한 값이 새로운 항으로 추가됨

5-1. 관측값 시간 간격이 다른 경우의 포아송회귀분석 수행

  • 사건발생횟수의 측정기간이 서로 달라 가변적인 측정기간을 갖는 포아송회귀분석을 수행할 것임
  • MASS패키지로 ships데이터셋 설명
    • 파도로 인한 선박의 손상횟수와 월단위로 측정된 사용기간이 기록돼 있음
    • 선박 손상횟수는 incidents변수에 저장되어 있음
    • 선박의 사용기간은 service변수에 기록되어 있음
    • 선박 종류인 type변수, 건조연도인 year변수, 운행기간인 period변수가 포함돼 있음
library(MASS)
str(ships)
## 'data.frame':    40 obs. of  5 variables:
##  $ type     : Factor w/ 5 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 2 2 ...
##  $ year     : int  60 60 65 65 70 70 75 75 60 60 ...
##  $ period   : int  60 75 60 75 60 75 60 75 60 75 ...
##  $ service  : int  127 63 1095 1095 1512 3353 0 2244 44882 17176 ...
##  $ incidents: int  0 0 3 4 6 18 0 11 39 29 ...
  • 선박의 손상횟수를 측정한 기간이 각 관측값마다 다르기 때문에 여기서는 사용기간 당 손상횟수를 결과변수로 하는 포아송회귀모델을 적용하는 것이 적절해 보임
  • 이를 위해 먼저 간단한 전처리 작업을 수행할 것임
    • 먼저 사용기간이 0개월인 선박은 분석해서 제외함
    • 그리고 year변수와 period변수는 범주형 변수로 변환할 것임
shipsinc <- subset(ships, service > 0)
shipsinc$year <- factor(shipsinc$year)
shipsinc$period <- factor(shipsinc$period)
levels(shipsinc$year)
## [1] "60" "65" "70" "75"
levels(shipsinc$period)
## [1] "60" "75"
  • 결과: 4개의 건조시점으로 구분된 year변수와 2개의 운행시점으로 구분된 period변수는 범주형 변수인 factor로 변환됨
  • 포아송회귀분석 수행 - glm()
    • 종전과 동일하게 주에 glm()을 이용할 수 있음
  • glm() 인수설명
    • 1st: 결과변수와 예측변수 간의 관계를 포뮬러 형식으로 지정함
    • data: 모델 생성을 위한 데이터셋을 지정함
    • family: poisson() 함수를 지정함
    • offset: 시간변수를 지정함. 여기서는 선박의 사용기간을 지정함
shipsinc.poisson <- glm(incidents ~ type + year + period, 
                        data=shipsinc, 
                        family=poisson(), 
                        offset=log(service))
summary(shipsinc.poisson)
## 
## Call:
## glm(formula = incidents ~ type + year + period, family = poisson(), 
##     data = shipsinc, offset = log(service))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6768  -0.8293  -0.4370   0.5058   2.7912  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -6.40590    0.21744 -29.460  < 2e-16 ***
## typeB       -0.54334    0.17759  -3.060  0.00222 ** 
## typeC       -0.68740    0.32904  -2.089  0.03670 *  
## typeD       -0.07596    0.29058  -0.261  0.79377    
## typeE        0.32558    0.23588   1.380  0.16750    
## year65       0.69714    0.14964   4.659 3.18e-06 ***
## year70       0.81843    0.16977   4.821 1.43e-06 ***
## year75       0.45343    0.23317   1.945  0.05182 .  
## period75     0.38447    0.11827   3.251  0.00115 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 146.328  on 33  degrees of freedom
## Residual deviance:  38.695  on 25  degrees of freedom
## AIC: 154.56
## 
## Number of Fisher Scoring iterations: 5
  • 분석결과를 보기 전 먼저 과산포 여부를 점검할 것임
deviance(shipsinc.poisson) / df.residual(shipsinc.poisson)
## [1] 1.547802
qcc.overdispersion.test(shipsinc$incidents , type="poisson")
##                    
## Overdispersion test Obs.Var/Theor.Var Statistic p-value
##        poisson data          23.64624  780.3258       0
  • 결과
    • 과산포 비율이 1보다 크고 qcc.overdispersion.test()의의 통계적 검정결과 또한 과산포를 지지함
    • 따라서 family 인수에 포아송을 지정하는 대신에 quasipoisson()을 지정해 poisson회귀모델을 다시 추정해 볼 것임
  • 그런데 update()를 이용하면 기존 생성모델을 기반으로 새로운 모델을 쉽게 추정할 수 있음
shipsinc.poisson <- update(shipsinc.poisson, 
                           family=quasipoisson())
summary(shipsinc.poisson)
## 
## Call:
## glm(formula = incidents ~ type + year + period, family = quasipoisson(), 
##     data = shipsinc, offset = log(service))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6768  -0.8293  -0.4370   0.5058   2.7912  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -6.40590    0.28276 -22.655  < 2e-16 ***
## typeB       -0.54334    0.23094  -2.353  0.02681 *  
## typeC       -0.68740    0.42789  -1.607  0.12072    
## typeD       -0.07596    0.37787  -0.201  0.84230    
## typeE        0.32558    0.30674   1.061  0.29864    
## year65       0.69714    0.19459   3.583  0.00143 ** 
## year70       0.81843    0.22077   3.707  0.00105 ** 
## year75       0.45343    0.30321   1.495  0.14733    
## period75     0.38447    0.15380   2.500  0.01935 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.691028)
## 
##     Null deviance: 146.328  on 33  degrees of freedom
## Residual deviance:  38.695  on 25  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
  • 회귀계수에 지수함수를 취해서 원래의 척도를 바탕으로 회귀계수를 해석할 수 있도록 변환할 것임
exp(coef(shipsinc.poisson))
## (Intercept)       typeB       typeC       typeD       typeE      year65 
##  0.00165178  0.58080262  0.50288104  0.92685196  1.38483286  2.00800246 
##      year70      year75    period75 
##  2.26693019  1.57369544  1.46883116

5-2. 관측값 시간 간격이 다른 경우의 포아송회귀분석 결과분석

  • 지수함수를 이용하여 변환한 회귀계수를 바탕으로 아래사항들을 확인할 수 있음

https://youtu.be/-Za5MrWUohg?list=PLY0OaF78qqGAxKX91WuRigHpwBU0C2SB_&t=1564

  • 선박 종류에 따른 분석
    • 선박 종류에서 A~E형 선박까지 5가지 종류가 있고, 기준범주는 A형 선박임
    • 따라서 여기 회귀계수는 기준범주인 A형 선박 대비 손상위험 정도를 나타냄
    • A형 선박 대비 B형은 0.58배, C형은 0.50배, D형은 0.92배, E형은 1.38배 더 손상 위험이 높음.
    • 따라서 E형의 손상위험이 가장 높고 C형의 손상위험이 가장 낮음
    • 기준범주인 A형 선박은 두번째로 손상위험이 높음
    • 가장 손상위험이 높은 E형 대비 가장 손상위험이 낮은 C형의 비는 2.75임
    • 따라서 E형 선박은 C형 선박에 비해 손상 위험이 2.75베 더 높음
  • 건조시점에 따른 분석
    • 건조시점은 4개의 범주를 가짐
      • 1960~64년인 선박, 65~69년, 70~74년, 75~79년인 선박으로 나눠짐
      • 그 중 기준범주는 1960~64년 사이 건조된 선박임
    • 회귀계수를 보면,
      • 65~69년 건조선박은 기준범주인 60~64년 대비 2배 더 손상위험이 높음
      • 70~74년 건조선박은 기준범주 대비 2.27배 더 손상위험이 높음
      • 75~79년 건조선박은 기준범주 대비 1.57배 더 손상 위험이 높음
      • 따라서 70~74년 건조선박이 가장 손상에 취약함
      • 기준범주인 60~64년 건조선박은 가장 안전함
      • 70~74년 건조선박은 60~64년 건조선박에 비해 2.27배 손상위험이 더 큼
  • 운행기간에 따른 분석
    • 운행기간은 2개의 범주를 가짐
    • 60~74년, 75~79년 사이 운행선박임
    • 그 중 첫 번째 범주인 60~74년 운행선박이 기준범주임
    • 기준범주 대비 75~79년 운행선박은 손상위험이 1.47배 더 높음

 

 

'R programming > 통계데이터분석' 카테고리의 다른 글

[R/통계] 38. 차원분석 - 요인분석(FA) 🔑 factor analysis | 공통요인 | 요인적재값(factor loading) | 요인점수(factor score) | 회전  (0) 2022.11.22
[R/통계] 37.차원분석 - 주성분분석(PCA) 🔑 principal component analysis | 차원축소 | 성분 | 주성분은 변수들의 선형결합으로 표현  (0) 2022.11.06
[R/통계] 34.일반선형모델 - 다항 로지스틱회귀분석 🔑 multinomial logistic regression analysis | 세 개 이상의 범주 | 로그오즈  (0) 2022.10.13
[R/통계] 34.일반선형모델 - 페널티 로지스틱회귀분석(ridge, lasso, elasticnet) 🔑 penalized regression analysis | 라소 | 릿지  (0) 2022.10.11
[R/통계] 32.일반선형모델 - 이항 로지스틱회귀분석 🔑 binomial logistic regression analysis | 이분형 범주 | 로그오즈(log odds)  (0) 2022.10.11

+ Recent posts

티스토리툴바