logit regression의 chi-square fit과 pseudo R-square

glm으로 logit regression을 만들고, model fit을 확인하려면?

스탠포드 대학의 정치과학계산연구소가 개발한 {pscl}의 pR2() 함수를 쓰면 McFadden의 R^2를 구할 수 있어요.

다음 예제 코드를 실행시켜 보세요.
library(MASS)
data(menarche)
model <- glm(cbind(Menarche, Total-Menarche) ~ Age,family=binomial(logit), data=menarche)
if(!require(pscl)) {
  install.packages("pscl")
  library(pscl)
}
# pR2()는 {pscl}의 함수
pR2(model)["McFadden"]

결과는
 McFadden
0.9706837

97%의 pseudo R^2네요.

McFadden's의 알고리즘은 categorical data analysis 책들이나 구글 검색해서 공부해요. ㅋㅋㅋ

모델의 적합도를 보는 다른 방법은 model이 데이터를 설명하는 예측력을 따져보는 것이에요. 이럴 때는 Pearson chi-square를 계산하지요.

1 - pchisq(deviance(model), df.residual(model))

위의 계산을 보면 chi-square의 p-value를 구하는 것임이 보이지요. 이때 p-value가 유의미하지 않게 나와야 model fit이 있다고 봅니다. 왜냐구요? null hypothesis가 model이 관측치 예측을 잘 한다... 이것이니 reject하면 안되니까요.

이 계산법에는 문제가 있을 수 있어요. Grouped data를 기반으로 한 검정법이라서 covariance가 어느 정도 그룹화된 패턴이 나와야 해요.

if(!require("epiR") {
  install.packages("epiR")
  library(epiR)
}
mf <- model.frame(model)
covariation_pattern <- epi.cp(mf[-1])
nrow(covariation_pattern$cov.pattern)

이려면 결과가 25개라고 나오는데... 우리 관측치가 25개라서. ㅎㅎㅎ 그룹이 안되지요.
le Cessie – van Houwelingen – Copas – Hosmer unweighted sum of squares test for global goodness of fit 이라는 방법도 사용할 수 있지만 제가 가방끈이 짧아서 ㅜ.ㅜ


댓글

이 블로그의 인기 게시물

Bradley-Terry Model: paired comparison models

R에서 csv 파일 읽는 법

xlwings tutorial - 데이터 계산하여 붙여 넣기