¿Cómo probar la igualdad simultánea de los coeficientes elegidos en el modelo logit o probit?

¿Cómo probar la igualdad simultánea de los coeficientes elegidos en el modelo logit o probit? ¿Cuál es el enfoque estándar y cuál es el enfoque del estado del arte?

Prueba de Wald

Un enfoque estándar es la prueba de Wald . Esto es lo que hace el comando Stata test después de una regresión logit o probit. Veamos cómo funciona esto en R mirando un ejemplo:

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

Digamos que quieres probar la hipótesis vs. β g r e ≠ β g p a . Esto es equivalente a probar β g r e - β g p a = $\beta_{gre}=\beta_{gpa}$$\beta_{gre}\neq \beta_{gpa}$$\beta_{gre} - \beta_{gpa} = 0$ . La estadística de prueba de Wald es:

$$W=\frac{(\hat{\beta}-\beta_{0})}{\widehat{\operatorname{se}}(\hat{\beta})}\sim \mathcal{N}(0,1)$$

o

$$W^2 = \frac{(\hat{\theta}-\theta_{0})^2}{\operatorname{Var}(\hat{\theta})}\sim \chi_{1}^2$$

Nuestra θ aquí es β g r e - β g p un y θ 0 = 0 . Entonces, todo lo que necesitamos es el error estándar de β g r e - β g p a . Podemos calcular el error estándar con el método Delta :$\widehat{\theta}$$\beta_{gre} - \beta_{gpa}$$\theta_{0}=0$$\beta_{gre} - \beta_{gpa}$

$$\hat{se}(\beta_{gre} - \beta_{gpa})\approx \sqrt{\operatorname{Var}(\beta_{gre}) + \operatorname{Var}(\beta_{gpa}) - 2\cdot \operatorname{Cov}(\beta_{gre},\beta_{gpa})}$$

Entonces también necesitamos la covarianza de y β g p a . La matriz de varianza-covarianza se puede extraer con el$\beta_{gre}$$\beta_{gpa}$vcov comando después de ejecutar la regresión logística:

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

Finalmente, podemos calcular el error estándar:

se <- sqrt(1.196831e-06 + 0.1101040465 -2*-0.0001241775)
se
[1] 0.3321951

Entonces su valor Wald es$z$

wald.z <- (gre-gpa)/se
wald.z
[1] -2.413564

Para obtener un valor , simplemente use la distribución normal estándar:$p$

2*pnorm(-2.413564)
[1] 0.01579735

En este caso tenemos evidencia de que los coeficientes son diferentes entre sí. Este enfoque puede extenderse a más de dos coeficientes.

Utilizando multcomp

Estos cálculos bastante tediosos se pueden hacer convenientemente al Rusar el multcomppaquete. Aquí está el mismo ejemplo que el anterior pero hecho con 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)

También se puede calcular un intervalo de confianza para la diferencia de los coeficientes:

Quantile = 1.96
95% family-wise confidence level


Linear Hypotheses:
               Estimate lwr     upr    
gre - gpa == 0 -0.8018  -1.4529 -0.1507

Para ejemplos adicionales de multcomp, vea aquí o aquí .

---

Prueba de razón de probabilidad (LRT)

Los coeficientes de una regresión logística se encuentran por máxima verosimilitud. Pero debido a que la función de probabilidad involucra muchos productos, la probabilidad de registro se maximiza, lo que convierte los productos en sumas. El modelo que se ajusta mejor tiene una mayor probabilidad de registro. El modelo que involucra más variables tiene al menos la misma probabilidad que el modelo nulo. Denote la log-verosimilitud del modelo alternativo (modelo que contiene más variables) con y la log-verosimilitud del modelo nulo con L L 0 , el estadístico de prueba de razón de verosimilitud es:$LL_{a}$$LL_{0}$

$$D=2\cdot (LL_{a} - LL_{0})\sim \chi_{df1-df2}^{2}$$

$\chi^{2}$ distribución de con los grados de libertad como la diferencia en el número de variables. En nuestro caso, esto es 2.

$\beta_{gre}=\beta_{gpa}$

$$\log\left(\frac{p_{i}}{1-p_{i}}\right)=\beta_{0}+\beta_{1}\cdot \mathrm{gre} + \beta_{2}\cdot \mathrm{gpa}+\beta_{3}\cdot \mathrm{rank_{2}} + \beta_{4}\cdot \mathrm{rank_{3}}+\beta_{5}\cdot \mathrm{rank_{4}}$$ $$\log\left(\frac{p_{i}}{1-p_{i}}\right)=\beta_{0}+\beta_{1}\cdot (\mathrm{gre} + \mathrm{gpa})+\beta_{2}\cdot \mathrm{rank_{2}} + \beta_{3}\cdot \mathrm{rank_{3}}+\beta_{4}\cdot \mathrm{rank_{4}}$$

mylogit2 <- glm(admit ~ I(gre + gpa) + rank, data = mydata, family = "binomial")

En nuestro caso, podemos usar logLikpara extraer la probabilidad logarítmica de los dos modelos después de una regresión logística:

L1 <- logLik(mylogit)
L1
'log Lik.' -229.2587 (df=6)

L2 <- logLik(mylogit2)
L2
'log Lik.' -232.2416 (df=5)

El modelo que contiene la restricción grey gpatiene una probabilidad de registro ligeramente mayor (-232.24) en comparación con el modelo completo (-229.26). Nuestro estadístico de prueba de razón de probabilidad es:

D <- 2*(L1 - L2)
D
[1] 16.44923

$\chi^{2}_{2}$$p$

1-pchisq(D, df=1)
[1] 0.01458625

$p$

R tiene la prueba de razón de probabilidad incorporada; podemos usar la anovafunción para calcular la prueba de razón de probabilidad:

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

Nuevamente, tenemos una fuerte evidencia de que los coeficientes grey gpason significativamente diferentes entre sí.

---

Prueba de puntuación (también conocida como prueba de puntuación de Rao, también conocida como prueba multiplicadora de Lagrange)

La función de puntuación $U(\theta)$ es la derivada de la función log-verosimilitud ($\text{log} L(\theta|x)$) dónde $\theta$ son los parámetros y $x$ los datos (el caso univariante se muestra aquí con fines ilustrativos):

$$U(\theta) = \frac{\partial \text{log} L(\theta|x)}{\partial \theta}$$

Esta es básicamente la pendiente de la función log-verosimilitud. Además, deja$I(\theta)$ser la matriz de información de Fisher, que es la expectativa negativa de la segunda derivada de la función de log-verosimilitud con respecto a$\theta$. Las estadísticas de la prueba de puntuación son:

$$S(\theta_{0})=\frac{U(\theta_{0}^{2})}{I(\theta_{0})}\sim\chi^{2}_{1}$$

La prueba de puntaje también se puede calcular usando anova(las estadísticas de la prueba de puntaje se llaman "Rao"):

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

La conclusión es la misma que antes.

---

Nota

An interesting relationship between the different test statistics when the model is linear is (Johnston and DiNardo (1997): Econometric Methods): Wald $\geq$ LR $\geq$ Score.

No especificó sus variables, si son binarias u otra cosa. Creo que hablas de variables binarias. También existen versiones multinomiales del modelo probit y logit.

En general, puede usar la trinidad completa de los enfoques de prueba, es decir

Prueba de razón de verosimilitud

Prueba LM

Prueba de Wald

Cada prueba usa diferentes estadísticas de prueba. El enfoque estándar sería tomar una de las tres pruebas. Los tres pueden usarse para hacer pruebas conjuntas.

La prueba LR utiliza la diferencia de la probabilidad logarítmica de un modelo restringido y no restringido. Por lo tanto, el modelo restringido es el modelo, en el que los coeficientes especificados se establecen en cero. El sin restricciones es el modelo "normal". La prueba de Wald tiene la ventaja de que solo se estima el modelo sin restricciones. Básicamente pregunta si la restricción está casi satisfecha si se evalúa en el MLE sin restricciones. En el caso de la prueba de multiplicador de Lagrange, solo se debe estimar el modelo restringido. El estimador de NM restringido se utiliza para calcular la puntuación del modelo sin restricciones. Este puntaje generalmente no será cero, por lo que esta discrepancia es la base de la prueba LR. La prueba LM también se puede utilizar en su contexto para evaluar la heterocedasticidad.

Los enfoques estándar son la prueba de Wald, la prueba de razón de probabilidad y la prueba de puntaje. Asintóticamente deberían ser lo mismo. En mi experiencia, las pruebas de razón de probabilidad tienden a funcionar ligeramente mejor en simulaciones en muestras finitas, pero los casos en los que esto importa sería en escenarios muy extremos (muestra pequeña) donde tomaría todas estas pruebas solo como una aproximación aproximada. Sin embargo, dependiendo de su modelo (número de covariables, presencia de efectos de interacción) y sus datos (multicolinealidad, la distribución marginal de su variable dependiente), el "maravilloso reino de Asymptotia" puede aproximarse bien por un número sorprendentemente pequeño de observaciones.

A continuación se muestra un ejemplo de tal simulación en Stata utilizando el Wald, la razón de probabilidad y la prueba de puntaje en una muestra de solo 150 observaciones. Incluso en una muestra tan pequeña, las tres pruebas producen valores p bastante similares y la distribución de muestreo de los valores p cuando la hipótesis nula es verdadera parece seguir una distribución uniforme como debería (o al menos las desviaciones de la distribución uniforme no son más grandes de lo que cabría esperar debido a la aleatoriedad heredada en un experimento de Monte Carlo).

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))