로짓 또는 프로 빗 모델에서 선택한 계수의 동시 동등성을 테스트하는 방법은 무엇입니까? 표준 접근법은 무엇이며 최첨단 접근법은 무엇입니까?
답변
왈드 테스트
하나의 표준 접근법은 Wald 테스트 입니다. 이것은 무엇 STATA 명령이 test
로짓 또는 프로 빗 회귀 후 않습니다. 예제를 통해 R에서 어떻게 작동하는지 봅시다 :
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv") # Load dataset from the web
mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial") # calculate the logistic regression
summary(mylogit)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.989979 1.139951 -3.500 0.000465 ***
gre 0.002264 0.001094 2.070 0.038465 *
gpa 0.804038 0.331819 2.423 0.015388 *
rank2 -0.675443 0.316490 -2.134 0.032829 *
rank3 -1.340204 0.345306 -3.881 0.000104 ***
rank4 -1.551464 0.417832 -3.713 0.000205 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
예를 들어 대 가설을 테스트하려고합니다 . 이는 테스트와 동일 합니다. Wald 검정 통계량은 다음과 같습니다.
또는
우리의 여기 와 . 따라서 필요한 것은 표준 오류 입니다. 델타 방법으로 표준 오차를 계산할 수 있습니다 . βgRE–βgPθ0=0βgRE–βgP
따라서 및 의 공분산도 필요합니다 . 로지스틱 회귀 분석을 실행 한 후 명령을 사용하여 분산 공분산 행렬을 추출 할 수 있습니다 . β g p a
vcov
var.mat <- vcov(mylogit)[c("gre", "gpa"),c("gre", "gpa")]
colnames(var.mat) <- rownames(var.mat) <- c("gre", "gpa")
gre gpa
gre 1.196831e-06 -0.0001241775
gpa -1.241775e-04 0.1101040465
마지막으로 표준 오차를 계산할 수 있습니다.
se <- sqrt(1.196831e-06 + 0.1101040465 -2*-0.0001241775)
se
[1] 0.3321951
따라서 Wald 값은
wald.z <- (gre-gpa)/se
wald.z
[1] -2.413564
값 을 얻으려면 표준 정규 분포를 사용하십시오.
2*pnorm(-2.413564)
[1] 0.01579735
이 경우 계수가 서로 다르다는 증거가 있습니다. 이 접근법은 두 개 이상의 계수로 확장 될 수 있습니다.
사용 multcomp
이 다소 지루한 계산은 패키지 를 R
사용하여 편리하게 수행 할 수 있습니다 multcomp
. 위와 같은 예이지만 다음과 같이 수행됩니다 multcomp
.
library(multcomp)
glht.mod <- glht(mylogit, linfct = c("gre - gpa = 0"))
summary(glht.mod)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
gre - gpa == 0 -0.8018 0.3322 -2.414 0.0158 *
confint(glht.mod)
계수 차이에 대한 신뢰 구간도 계산할 수 있습니다.
Quantile = 1.96
95% family-wise confidence level
Linear Hypotheses:
Estimate lwr upr
gre - gpa == 0 -0.8018 -1.4529 -0.1507
에 대한 추가 예는 여기 또는 여기를multcomp
참조 하십시오 .
우도 비 검정 (LRT)
로지스틱 회귀의 계수는 최대 가능성으로 구합니다. 그러나 우도 기능에는 많은 제품이 포함되므로 로그 우도는 최대화되어 제품을 합산합니다. 더 잘 맞는 모델은 로그 가능성이 더 높습니다. 더 많은 변수를 포함하는 모델은 null 모델과 적어도 같은 가능성을 갖습니다. 의 대체 모델 (더 많은 변수를 포함하는 모델) 의 로그 우도와 의 null 모델의 로그 우도를 나타내는 우도 비율 검정 통계량은 다음과 같습니다. L L 0
우도 비율 검정 통계량 은 자유도가 변수 개수의 차이 인 분포를 따릅니다 . 우리의 경우 이것은 2입니다.
우도 비 검정을 수행하려면 두 우도를 비교할 수 있도록 제약 조건을 사용하여 모형을 적합시켜야합니다 . 전체 모델은 . 제약 조건 모델은 . 로그 ( P의 난
mylogit2 <- glm(admit ~ I(gre + gpa) + rank, data = mydata, family = "binomial")
이 경우 logLik
로지스틱 회귀 분석 후 두 모델의 로그 우도를 추출하는 데 사용할 수 있습니다 .
L1 <- logLik(mylogit)
L1
'log Lik.' -229.2587 (df=6)
L2 <- logLik(mylogit2)
L2
'log Lik.' -232.2416 (df=5)
에 제약을 포함하는 모델 gre
과는 gpa
전체 모델 (-229.26)에 비해 약간 더 높은 로그 우도 (-232.24)가 있습니다. 우리의 우도 비 검정 통계량은 다음과 같습니다.
D <- 2*(L1 - L2)
D
[1] 16.44923
이제 의 CDF를 사용 하여 값 을 계산할 수 있습니다 .
1-pchisq(D, df=1)
[1] 0.01458625
-value의 계수가 다른 것을 나타내는 매우 작다.
R에는 우도 비 검정이 내장되어 있습니다. 이 anova
함수를 사용하여 우도 비 검정을 계산할 수 있습니다 .
anova(mylogit2, mylogit, test="LRT")
Analysis of Deviance Table
Model 1: admit ~ I(gre + gpa) + rank
Model 2: admit ~ gre + gpa + rank
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 395 464.48
2 394 458.52 1 5.9658 0.01459 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
다시 말하지만, 계수 gre
와 계수 gpa
가 크게 다르다는 강력한 증거가 있습니다 .
점수 시험 (일명 Rao의 점수 시험, 일명 라그랑주 승수 시험)
점수 함수 로그 – 우도 함수의 유도체 ( ) 여기서 파라미터이고 데이터 (단 변량 경우를 돕기 위해 여기에 도시 목적) :
이것은 기본적으로 로그 우도 함수의 기울기입니다. 또한하자 될 피셔 정보 행렬 에 대한 상기 로그 – 우도 함수의 이차 미분의 음 기대이다 . 점수 시험 통계는 다음과 같습니다.
점수 테스트는 다음을 사용하여 계산할 수도 있습니다 anova
(점수 테스트 통계는 “라오”라고 함).
anova(mylogit2, mylogit, test="Rao")
Analysis of Deviance Table
Model 1: admit ~ I(gre + gpa) + rank
Model 2: admit ~ gre + gpa + rank
Resid. Df Resid. Dev Df Deviance Rao Pr(>Chi)
1 395 464.48
2 394 458.52 1 5.9658 5.9144 0.01502 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
결론은 이전과 동일합니다.
노트
모형이 선형 일 때 다른 검정 통계 간의 흥미로운 관계는 (Johnston and DiNardo (1997) : Econometric Methods ) : Wald LR Score입니다.
답변
변수가 이진이거나 다른 변수 인 경우 변수를 지정하지 않았습니다. 이진 변수에 대해 이야기한다고 생각합니다. 다항식 버전의 프로 빗 및 로짓 모델도 있습니다.
일반적으로 테스트 접근법의 완전한 삼위 일체를 사용할 수 있습니다.
우도 비 테스트
LM 테스트
왈드 테스트
각 테스트는 서로 다른 테스트 통계를 사용합니다. 표준 접근법은 세 가지 테스트 중 하나를 수행하는 것입니다. 세 가지를 모두 공동 테스트에 사용할 수 있습니다.
LR 테스트는 제한된 모델과 무제한 모델의 로그 가능성의 차이를 사용합니다. 따라서 제한된 모델은 지정된 계수가 0으로 설정된 모델입니다. 무제한은 “정상”모델입니다. Wald 검정은 무조건 모델 만 추정한다는 장점이 있습니다. 기본적으로 제한되지 않은 MLE에서 평가되는 경우 제한이 거의 충족되는지 묻습니다. Lagrange-Multiplier 테스트의 경우 제한된 모델 만 추정하면됩니다. 제한된 ML 추정기는 무제한 모델의 점수를 계산하는 데 사용됩니다. 이 점수는 일반적으로 0이 아니므로이 불일치는 LR 테스트의 기초입니다. 컨텍스트에서 LM-Test를 사용하여 이분산성을 테스트 할 수 있습니다.
답변
표준 접근법은 Wald 검정, 우도 비율 검정 및 점수 검정입니다. 무의식적으로 그들은 동일해야합니다. 내 경험상 우도 비율 검정은 유한 표본에 대한 시뮬레이션에서 약간 더 나은 성능을 보이는 경향이 있지만,이 문제가 매우 극단적 인 (작은 표본) 시나리오에있을 경우 이러한 모든 검정을 대략적인 근사치로만 사용합니다. 그러나 모형 (공변량의 수, 교호 작용 효과의 존재) 및 데이터 (다중 선형성, 종속 변수의 한계 분포)에 따라, “엄청난 증상이없는 무증상 왕국”은 놀라 울 정도로 적은 수의 관측으로 근사치가 될 수 있습니다.
아래는 단지 150 개의 관측치 샘플에서 Wald, 우도 비 및 점수 테스트를 사용한 Stata에서의 시뮬레이션의 예입니다. 이러한 작은 표본에서도 3 개의 검정은 상당히 유사한 p- 값을 생성하고 귀무 가설이 참일 때 p- 값의 샘플링 분포는 균일 한 분포를 따르는 것처럼 보입니다 (또는 적어도 균일 한 분포와의 편차) 몬테카를로 실험에서 무작위 상속으로 인해 예상되는 것보다 크지 않습니다.
clear all
set more off
// data preparation
sysuse nlsw88, clear
gen byte edcat = cond(grade < 12, 1, ///
cond(grade == 12, 2, 3)) ///
if grade < .
label define edcat 1 "less than high school" ///
2 "high school" ///
3 "more than high school"
label value edcat edcat
label variable edcat "education in categories"
// create cascading dummies, i.e.
// edcat2 compares high school with less than high school
// edcat3 compares more than high school with high school
gen byte edcat2 = (edcat >= 2) if edcat < .
gen byte edcat3 = (edcat >= 3) if edcat < .
keep union edcat2 edcat3 race south
bsample 150 if !missing(union, edcat2, edcat3, race, south)
// constraining edcat2 = edcat3 is equivalent to adding
// a linear effect (in the log odds) of edcat
constraint define 1 edcat2 = edcat3
// estimate the constrained model
logit union edcat2 edcat3 i.race i.south, constraint(1)
// predict the probabilities
predict pr
gen byte ysim = .
gen w = .
program define sim, rclass
// create a dependent variable such that the null hypothesis is true
replace ysim = runiform() < pr
// estimate the constrained model
logit ysim edcat2 edcat3 i.race i.south, constraint(1)
est store constr
// score test
tempname b0
matrix `b0' = e(b)
logit ysim edcat2 edcat3 i.race i.south, from(`b0') iter(0)
matrix chi = e(gradient)*e(V)*e(gradient)'
return scalar p_score = chi2tail(1,chi[1,1])
// estimate unconstrained model
logit ysim edcat2 edcat3 i.race i.south
est store full
// Wald test
test edcat2 = edcat3
return scalar p_Wald = r(p)
// likelihood ratio test
lrtest full constr
return scalar p_lr = r(p)
end
simulate p_score=r(p_score) p_Wald=r(p_Wald) p_lr=r(p_lr), reps(2000) : sim
simpplot p*, overall reps(20000) scheme(s2color) ylab(,angle(horizontal))