이 문제는 일반선형화모형(Generalized Linear Model, GLM)에서 모수 추정을 어떤 방식으로 하느냐 하는 문제이네요.

우리가 회귀 모형에서 모수( 회귀 계수 , coefficient)를 추정할 때 , 회귀계수가 0인지 아닌지를 검정하게 됩니다

\(H_0: \beta= 0\) (여기서 베타는 회귀계수)


이 가설을 검정하는 방법은 크게 2가지 인데, 1)Wald’s test 2) Likelihood ratio test 가 있습니다

1) Wald’s test

\(Z=\frac{\beta}{S.E}\) 는 근사적으로 N(0,1)을 따른다. 혹은 \(Z^2=(\frac{\beta}{S.E})^2\) 는 근사적으로 자유도1인 \(\chi^2\) 분포를 따른다.

sov<-read.csv('sov.csv')
model1<-glm(Recalcitrant~Age+AT_PN,data = sov,family = binomial)
model1$coefficients
(Intercept)         Age       AT_PN 
-3.12391801  0.05467492 -1.33523183

여기서 왈드 검정을 실시한다

1-pchisq((model1$coef[3])^2/(coef(summary(model1))[, "Std. Error"][3])^2,1)
     AT_PN 
0.07309232

결과 p-value는 0.05보다 크며(SOL 님께서 올려주신 값과 일치하죠), 따라서 AT_PN의 계수가 0이라는 귀무가설을 기각하지 못합니다.

2)Likeihood ratio test (가능도비 검정 )

가능도비 검정은 여기서 설명하려면 길고 개념도 약간 어렵습니다. 여기서는 “복잡한 모형과 간단한 모형중에 어떤 것이 설명력이 더 좋은지 알아보는 검정” 이라고만 간단히 소개하고 넘어가겠습니다.

model2<-glm(Recalcitrant~Age,data = sov,family = binomial) # model1과 비교하기위해 AT_PN변수를 뺀  model2를 적합하였음 
library(lmtest)
lrtest(model2,model1)
Likelihood ratio test

Model 1: Recalcitrant ~ Age
Model 2: Recalcitrant ~ Age + AT_PN
  #Df  LogLik Df  Chisq Pr(>Chisq)  
1   2 -88.475                       
2   3 -86.500  1 3.9493    0.04689 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

대략 설명하면 이렇습니다. 모형을 적합할때, ‘최소한의 독립 변수로, 종속변수를 최대한 설명하고 싶어 합니다’. 즉 age+AT_PN 보다, age 하나만 가지고 reclacitrant를 설명하고 싶긴 한데, 2개의 모델을 비교해 보니 AT_PN 이라는 변수를 넣는 것이 , 변수 개수가 증가하더라도 설명력은 훨씬 좋아집니다. (P-VALUE<0.05)

즉 가능도비 검정에서는 AT_PN이 유의미한 독립변수로 나오는 것이지요.

그러면 여기서는 왜 p-value와 95% 신뢰구간 간에 불일치를 보인 것일까요?

R의 confint 함수는 가능도비 검정을 바탕으로 신뢰구간을 구하기 때문입니다. 즉 p-value 계산할때는 왈드 검정을 하고, 신뢰구간은 가능도비 검정으로 구한 것입니다. 여기서 차이가 발생하였습니다.

일반 선형화 모형(GLM)에서 모수 추정을 할때, 표본의 크기가 충분히 크기 못하면 왈드 검정과 가능도비 검정의 결과가 일치하지 않는 경우가 발생합니다. 충분히 크다는게 딱 정해져있지 않지만, 이런 경우에는가능도비 검정 결과를 따르는 것이 좋다 고 알려져 있습니다. (이건 제 대학원 교재에 적혀 있는 설명입니다. 왜 그런지 이유야 있겠지만.. 저는 통계학자가 아니니까 굳이 이유까지 알 필요는 없을 것 같네요)

아무튼 결론 적으로는, 가능도비로 구한 신뢰구간을 따르는 것이 맞으며, AT_PN은 유의미한 독립변수가 되겠네요.

만약에 95% 신뢰구간을 왈드 검정으로 구하고 싶다면 ’confint.defualt’함수를 쓰면 됩니다

confint.default(model1)
                  2.5 %      97.5 %
(Intercept) -3.79340512 -2.45443091
Age          0.03272269  0.07662715
AT_PN       -2.79540803  0.12494437
exp(confint.default(model1))
                 2.5 %     97.5 %
(Intercept) 0.02251879 0.08591207
Age         1.03326396 1.07963946
AT_PN       0.06108994 1.13308541

이것은 수기로 계산한 신뢰구간과 일치합니다.

이 질문에 대한 아이디어는 https://stats.stackexchange.com/questions/5304/why-is-there-a-difference-between-manually-calculating-a-logistic-regression-95 에서 얻었습니다.

덕분에 confint 함수가 어떻게 신뢰구간을 계산하는지 알게 되었습니다. 저도 좋은 공부가 되었습니다.