로지스틱 회귀 분석에서 적당한 모델을 선정했다고 해도 이것이 실제로 옳게 예측을 하는지는 별도의 검정이 필요합니다. 바로 모형의 적합도 (goodness of fit of the model) 검증이 그것입니다. 로지스틱 회귀 모형이 제대로 결과를 예측하는지 확인하는 방법 역시 여러 가지 이지만, 모형 적합도 검증을 위해 가장 흔하게 사용되는 방법은 바로 hosmer-lemeshow goodness-of-fit 혹은 hosmer-lemeshow test입니다. 아래 링크는 참고할 만한 포스팅입니다.
R에서 여러 패키지가 hosmer-lemeshow test를 지원하지만, 간단하게 코드를 이용해서 구현하는 경우가 많습니다. 아래 코드가 구글 검색시 가장 흔하게 만날 수 있는 것입니다.
hosmerlem = function(y, yhat, g=10) {
cutyhat = cut(yhat,
breaks = quantile(yhat, probs=seq(0,
1, 1/g)), include.lowest=TRUE)
obs = xtabs(cbind(1 - y, y) ~ cutyhat)
expect = xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
chisq = sum((obs - expect)^2/expect)
P = 1 - pchisq(chisq, g - 2)
return(list(chisq=chisq,p.value=P))
}
앞서 피마 인디언 예제에 적용해 보겠습니다.
#model - AIC
library(moonBook)
library(mlbench)
data(PimaIndiansDiabetes)
pima <- pimaindiansdiabetes="" span="">
pima<-subset pima="" pressure="">0)
pima<-subset mass="" pima="">0)
pima<-subset glucose="" pima="">0)
pima$obesity[pima$mass>=30]=2
pima$obesity[pima$mass<30 span="">
pima$obesity[pima$mass<25 span="">
table(pima$obesity)
out=glm(diabetes~factor(obesity)+pressure+glucose+age+pedigree,family=binomial,data=pima)
summary(out)
extractOR(out)
적용은 hosmerlem(y=, yhat=fitted())로 y에 해당하는 값은 모델에서 y (여기서는 diabetes) 이며 fitted () 안에 들어갈 것은 모델인 out입니다.
> hosmerlem(y=factor(pima$diabetes), yhat=fitted(out))
$`chisq`
[1] 724
$p.value
[1] 0
Warning message:
In Ops.factor(1, y) : ‘-’ not meaningful for factors
그런데 이렇게 하면 에러가 나옵니다. 왜일까요? 코드를 보면 그 이유를 알 수 있습니다. 이 코드는 로지스틱 회귀 모델에서 y 값이 0 혹은 1로 코딩된 경우를 염두에 두고 있습니다. 따라서 값을 그렇게 변경해주면 됩니다.
pima$DM=ifelse(pima$diabetes=="pos",1,0)
out=glm(DM~factor(obesity)+pressure+glucose+age+pedigree,family=binomial,data=pima)
summary(out)
extractOR(out)
hosmerlem(y=pima$DM, yhat=fitted(out))
> hosmerlem(y=pima$DM, yhat=fitted(out))
$`chisq`
[1] 4.727695
$p.value
[1] 0.7862444
결과값에서 chisq (카이 제곱값)이 작을 수록 모형이 적합한 것이고 P 값이 1에 가까울수록 적합한 것입니다. 판정은 P 값이 0.05 이하이면 모델에 문제가 있다고 보고 모델을 바꾸는 것입니다. 그런데 hosmer-lemeshow test는 기본적으로 표본수가 커야 결과값을 신뢰할 수 있습니다. 10개의 구간으로 나눠서 적합도를 검사하는 방법이기 때문에 만약 구간 양성 (보통 1로 코딩하는)인 표본이 5개 미만이면 에러가 있다고 평가하게 됩니다. 이 코드에는 그것이 구현되지 않았기 때문에 아래 코드를 사용해 보겠습니다.
참고로 모형 적합도를 평가하는 다른 방법인 Osius-Rojek test와 Stukel Test도 같이 코드를 넣었는데, 샘플수가 많을 때도 믿을 수 있는 결과를 주기 때문에 더 신뢰할 수 있는 방법이라는 평가도 있습니다. 다만 기본적으로 사용하는 검사법은 hosmer-lemeshow test 입니다.
#diagnostics
#H-L test
HLTest = function(obj, g) {
# first, check to see if we fed in the right kind of object
stopifnot(family(obj)$family == "binomial" && family(obj)$link == "logit")
y = obj$model[[1]]
trials = rep(1, times = nrow(obj$model))
if(any(colnames(obj$model) == "(weights)"))
trials <- model="" ncol="" obj="" span="">
# the double bracket (above) gets the index of items within an object
if (is.factor(y))
y = as.numeric(y) == 2 # Converts 1-2 factor levels to logical 0/1 values
yhat = obj$fitted.values
interval = cut(yhat, quantile(yhat, 0:g/g), include.lowest = TRUE) # Creates factor with levels 1,2,...,g
Y1 <- span="" trials="" y="">
Y0 <- -="" span="" trials="" y1="">
Y1hat <- span="" trials="" yhat="">
Y0hat <- -="" span="" trials="" y1hat="">
obs = xtabs(formula = cbind(Y0, Y1) ~ interval)
expect = xtabs(formula = cbind(Y0hat, Y1hat) ~ interval)
if (any(expect < 5))
warning("Some expected counts are less than 5. Use smaller number of groups")
pear <- -="" expect="" obs="" span="" sqrt="">
chisq = sum(pear^2)
P = 1 - pchisq(chisq, g - 2)
# by returning an object of class "htest", the function will perform like the
# built-in hypothesis tests
return(structure(list(
method = c(paste("Hosmer and Lemeshow goodness-of-fit test with", g, "bins", sep = " ")),
data.name = deparse(substitute(obj)),
statistic = c(X2 = chisq),
parameter = c(df = g-2),
p.value = P,
pear.resid = pear,
expect = expect,
observed = obs
), class = 'htest'))
}
# Osius-Rojek test
# Based on description in Hosmer and Lemeshow (2000) p. 153.
# Assumes data are aggregated into Explanatory Variable Pattern form.
o.r.test = function(obj) {
# first, check to see if we fed in the right kind of object
stopifnot(family(obj)$family == "binomial" && family(obj)$link == "logit")
mf <- model="" obj="" span="">
trials = rep(1, times = nrow(mf))
if(any(colnames(mf) == "(weights)"))
trials <- mf="" ncol="" span="">
prop = mf[[1]]
# the double bracket (above) gets the index of items within an object
if (is.factor(prop))
prop = as.numeric(prop) == 2 # Converts 1-2 factor levels to logical 0/1 values
pi.hat = obj$fitted.values
y <- prop="" span="" trials="">
yhat <- pi.hat="" span="" trials="">
nu <- pi.hat="" span="" yhat="">
pearson <- -="" nu="" span="" sum="" y="" yhat="">
c = (1 - 2*pi.hat)/nu
exclude <- c="" colnames="" mf="" span="" weights="" which="">
vars <- c="" data.frame="" exclude="" mf="" nbsp="" span="">
wlr <- .="" data="vars)</span" formula="c" lm="" weights="nu,">
rss <- nu="" residuals="" span="" sum="" wlr="">
J <- mf="" nrow="" span="">
A <- -="" 2="" span="" sum="" trials="">
z <- -="" 1="" ncol="" pearson="" rss="" span="" sqrt="" vars="">
p.value <- -="" 2="" abs="" pnorm="" span="" z="">
cat("z = ", z, "with p-value = ", p.value, "\n")
}
# Stukel Test
# Based on description in Hosmer and Lemeshow (2000) p. 155.
# Assumes data are aggregated into Explanatory Variable Pattern form.
stukel.test = function(obj) {
# first, check to see if we fed in the right kind of object
stopifnot(family(obj)$family == "binomial" && family(obj)$link == "logit")
high.prob <- fitted.values="" obj="">= 0.5)
logit2 <- linear.predictors="" obj="" span="">
z1 = 0.5*logit2*high.prob
z2 = 0.5*logit2*(1-high.prob)
mf <- model="" obj="" span="">
trials = rep(1, times = nrow(mf))
if(any(colnames(mf) == "(weights)"))
trials <- mf="" ncol="" span="">
prop = mf[[1]]
# the double bracket (above) gets the index of items within an object
if (is.factor(prop))
prop = (as.numeric(prop) == 2) # Converts 1-2 factor levels to logical 0/1 values
pi.hat = obj$fitted.values
y <- prop="" span="" trials="">
exclude <- colnames="" mf="" span="" weights="" which="">
vars <- c="" data.frame="" exclude="" mf="" span="" y="" z1="" z2="">
full <- .="" data="vars)</span" family="binomial(link" formula="y/trials" glm="" logit="" weights="trials,">
null <- .="" data="vars[,-c(1,2)])</span" family="binomial(link" formula="y/trials" glm="" logit="" weights="trials,">
LRT <- anova="" full="" null="" span="">
p.value <- -="" 1="" eviance="" f="" lrt="" pchisq="" span="">
cat("Stukel Test Stat = ", LRT$Deviance[[2]], "with p-value = ", p.value, "\n")
}
HLTest(out,g=10)
o.r.test(out)
stukel.test(out)
> HLTest(out,g=10)
Hosmer and Lemeshow goodness-of-fit test with 10 bins
data: out
X2 = 4.7277, df = 8, p-value = 0.7862
Warning message:
In HLTest(out, g = 10) :
Some expected counts are less than 5. Use smaller number of groups
> o.r.test(out)
z = -0.04210769 with p-value = 0.9664129
> stukel.test(out)
Stukel Test Stat = 8.392487 with p-value = 0.01505201
g=10 으로 했을 때 역시나 경고 메세지가 나오게 됩니다. 이 때는 5개 정도로 그룹을 작게 하면 문제를 피할 수 있습니다.
> HLTest(out,g=5)
Hosmer and Lemeshow goodness-of-fit test with 5 bins
data: out
X2 = 2.7727, df = 3, p-value = 0.428
P값이 크기 때문에 모형은 적합한 것으로 나타났습니다. 참고로 모형을 검증할때는 세 가지 방법 중 하나만 선택해야 합니다. 세 가지 방법을 사용하면 결국 드문 일도 일어날 가능성이 그 만큼 커지기 때문에 0.05를 기준으로 검증하는 의미가 희석됩니다. 보통 H-L test가 가장 무난한 방법입니다.
댓글
댓글 쓰기