¿Cuál es la diferencia entre estas dos pruebas Breusch-Pagan?

Question

Usando R en algunos datos y tratando de ver si mis datos son heteroscedásticos o no, he encontrado dos implementaciones de la prueba de Breusch-Pagan, bptest (paquete lmtest) y ncvTest (coche de paquete). Sin embargo, estos producen resultados diferentes. ¿Cuál es la diferencia entre ambos? ¿Cuándo se debe optar por uno u otro?

> model <- lm(y ~ x)
> bp <- bptest(model)
> bp
studentized Breusch-Pagan test

data:  model 
BP = 3.3596, df = 1, p-value = 0.06681

> ncvTest(model)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 3.858704    Df = 1     p = 0.04948855 

Este ejemplo muestra que, según las pruebas, mis datos son en un caso heteroscedásticos y en el otro homoscedásticos. He encontrado esta pregunta aquí por lo que bptest podría ser estudiado y ncvTest podría no serlo, sin embargo, ¿qué significa esto entonces?

Answer 1

Tu suposición es correcta, ncvTest realiza la versión original de la prueba Breusch-Pagan. En realidad, esto puede verificarse comparándolo con bptest(model, studentize = FALSE) . (Como señaló @Helix123, dos funciones también difieren en otros aspectos como los argumentos por defecto, uno debe comprobar los manuales del paquete de lmtest y car para más detalles).

El test de Breusch-Pagan estudiado fue propuesto por R. Koenker en su artículo de 1981 Nota sobre el estudio de una prueba de heteroscedasticidad . La diferencia más evidente entre ambos es que utilizan estadísticas de prueba diferentes. En concreto, dejemos que $\xi^\ast$ sean los estadísticos de la prueba de estudio y $\hat{\xi}$ sea el original, $$\newcommand{\Var}{\operatorname{Var}}\hat{\xi}=\lambda\xi^\ast,\qquad\lambda=\frac{\Var(\varepsilon^2)}{2\Var(\varepsilon)^2}.$$

Aquí hay un fragmento de código que demuestra lo que acabo de escribir (datos tomados de faraway paquete):

> mdl = lm(final ~ midterm, data = stat500)
> bptest(mdl)

    studentized Breusch-Pagan test

data:  mdl
BP = 0.86813, df = 1, p-value = 0.3515

> bptest(mdl, studentize = FALSE)

    Breusch-Pagan test

data:  mdl
BP = 0.67017, df = 1, p-value = 0.413

> ncvTest(mdl)
Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.6701721    Df = 1     p = 0.4129916 
> 
> n = nrow(stat500)
> e = residuals(mdl)
> bpmdl = lm(e^2 ~ midterm, data = stat500)
> lambda = (n - 1) / n * var(e^2) / (2 * ((n - 1) / n * var(e))^2)
> Studentized_bp = n * summary(bpmdl)$r.squared
> Original_bp = Studentized_bp * lambda
> 
> Studentized_bp
[1] 0.8681335
> Original_bp
[1] 0.6701721

En cuanto a la razón por la que se quiere estudiar la prueba original de BP, una cita directa del artículo de R. Koenker puede ser útil:

... De este análisis se desprenden dos conclusiones:

1. La potencia asintótica de la prueba de Breusch y Pagan es extremadamente sensible a la curtosis de la distribución de $\varepsilon$ y
2. el tamaño asintótico de la prueba es correcto sólo en el caso especial de la curtosis gaussiana.

La primera conclusión se amplía en Koenker y Bassett (1981) donde se sugieren pruebas alternativas y robustas para la heteroscedasticidad. La segunda conclusión implica que los niveles de significación sugeridos sugeridos por Breusch y Pagan serán correctos sólo bajo condiciones gaussianas en $\varepsilon$ . Dado que estas condiciones se asumen generalmente por fe ciega y son notoriamente difíciles de verificar, se sugiere una modificación de la prueba de Breusch y Pagan, que "estudian" correctamente el estadístico de la prueba y conduce a niveles de significación asintóticamente correctos para una clase razonablemente amplia de distribuciones para $\varepsilon$ .

En resumen, la prueba de la PA estudiada es más robusta que la original.

Answer 2

En términos prácticos, ncvTest utiliza el lado izquierdo de la ecuación y bptest utiliza el lado derecho, por defecto.

Esto significa que en un caso de Y ~ X ambas pruebas proporcionarán los mismos resultados (en cuanto a la studentize = F opción de bptest ). Pero en un análisis multivariante como Y ~ X1 + X2 los resultados serán diferentes. (Como señaló @Helix123)

En el archivo de ayuda de ncvTest : var.formula : " una fórmula unilateral para la varianza del error; si se omite, la varianza del error depende de los valores ajustados ." Lo que significa que, por defecto, se utilizarán los valores ajustados, pero también permite utilizar una combinación lineal de las variables independientes (X1 + X2).

En el archivo de ayuda de bptest : varformula : "Por defecto las mismas variables explicativas se toman como en el modelo de regresión principal".

Continuando con el mismo ejemplo de @Francis (datos stat500 , de faraway paquete):

> mdl_t = lm(final ~ midterm + hw, data = stat500)

Prueba de PA, utilizando valores ajustados:

> ncvTest(mdl_t) # Default

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 0.6509135    Df = 1     p = 0.4197863 

> bptest(mdl_t, varformula = ~ fitted.values(mdl_t), studentize = F)

Breusch-Pagan test

data:  mdl_t
BP = 0.65091, df = 1, p-value = 0.4198

Prueba de PA, utilizando una combinación lineal de predictores:

> ncvTest(mdl_t, var.formula = ~ midterm + hw)
Non-constant Variance Score Test 
Variance formula: ~ midterm + hw 
Chisquare = 0.7689743    Df = 2     p = 0.6807997 

> bptest(mdl_t, studentize = F) # Default

Breusch-Pagan test

data:  mdl_t
BP = 0.76897, df = 2, p-value = 0.6808

La "opción de combinación lineal" permite investigar la heteroscedasticidad asociada a la dependencia lineal de una variable independiente específica. Por ejemplo, sólo la hw variable:

> ncvTest(mdl_t, var.formula = ~ hw)
Non-constant Variance Score Test 
Variance formula: ~ hw 
Chisquare = 0.04445789    Df = 1     p = 0.833004 

> bptest(mdl_t, varformula = ~ stat500$hw, studentize = F)

Breusch-Pagan test

data:  mdl_t
BP = 0.044458, df = 1, p-value = 0.833

Por último, como ha resumido @Francis, "En resumen, la prueba de PA estudiada es más robusta que la original", yo suelo ir con bptest con studentize = TRUE (por defecto) y varformula = ~ fitted.values(my.lm) como opciones, para una primera aproximación a la homoscedasticidad.