The original Vuong test

Vuong (1989) proposed a test for non-nested model. He considers two competing models, \(F_\beta = \left\{f(y|z; \beta); \beta \in \beta\right\}\) and \(G_\gamma = \left\{g(y|z; \gamma); \gamma \in \Gamma\right\}\). Denoting \(h(y | z)\) the true conditional density. The distance of \(F_\beta\) from the true model is measured by the minimum KLIC:

\[ D_f = \mbox{E}^0\left[\ln h(y\mid z)\right] - \mbox{E}^0\left[\ln f(y\mid z; \beta_*)\right] \]

where \(\mbox{E}^0\) is the expected value using the true joint distribution of \((y, X)\) and \(\beta_*\) is the pseudo-true value of \(\beta\)¹. As the true model is unobserved, denoting \(\theta^\top = (\beta ^ \top, \gamma ^ \top)\), we consider the difference of the KLIC distance to the true model of model \(G_\gamma\) and model \(F_\beta\):

\[ \Lambda(\theta) = D_g - D_f = \mbox{E}^0\left[\ln f(y\mid z; \beta_*)\right]- \mbox{E}^0\left[\ln g(y\mid z; \gamma_*)\right] = \mbox{E}^0\left[\ln \frac{f(y\mid z; \beta_*)}{g(y\mid z; \gamma_*)}\right] \]

The null hypothesis is that the distance of the two models to the true models are equal or, equivalently, that: \(\Lambda=0\). The alternative hypothesis is either \(\Lambda>0\), which means that \(F_\beta\) is better than \(G_\gamma\) or \(\Lambda<0\), which means that \(G_\gamma\) is better than \(F_\beta\). Denoting, for a given random sample of size \(N\), \(\hat{\beta}\) and \(\hat{\gamma}\) the maximum likelihood estimators of the two models and \(\ln L_f(\hat{\beta})\) and \(\ln L_g(\hat{\gamma})\) the maximum value of the log-likelihood functions of respectively \(F_\beta\) and \(G\gamma\), \(\Lambda\) can be consistently estimated by:

\[ \hat{\Lambda}_N = \frac{1}{N} \sum_{n = 1} ^ N \left(\ln f(y_n \mid x_n, \hat{\beta}) - \ln g(y_n \mid x_n, \hat{\gamma})\right) = \frac{1}{N} \left(\ln L_f(\hat{\beta}) - \ln L_g(\hat{\gamma})\right) \]

which is the likelihood ratio divided by the sample size. Note that the statistic of the standard likelihood ratio test, suitable for nested models is \(2 \left(\ln L^f(\hat{\beta}) - \ln L^g(\hat{\gamma})\right)\), which is \(2 N \hat{\Lambda}_N\).

The variance of \(\Lambda\) is:

\[ \omega^2_* = \mbox{V}^o \left[\ln \frac{f(y \mid x; \beta_*)}{g(y \mid x; \gamma_*)}\right] \]

which can be consistently estimated by:

\[ \hat{\omega}_N^2 = \frac{1}{N} \sum_{n = 1} ^ N \left(\ln f(y_n \mid x_n, \hat{\beta}) - \ln g(y_n \mid x,_n \hat{\gamma})\right) ^ 2 - \hat{\Lambda}_N ^ 2 \]

Three different cases should be considered:

• when the two models are nested, \(\omega^2_*\) is necessarily 0,
• when the two models are overlapping (which means than the models coincides for some values of the parameters), \(\omega^2_*\) may be equal to 0 or not,
• when the two models are stricly non-nested, \(\omega^2_*\) is necessarely strictly positive.

The distribution of the statistic depends on whether \(\omega^2_*\) is zero or positive. If \(\omega^2_*\) is positive, the statistic is \(\hat{T}_N = \sqrt{N}\frac{\hat{\Lambda}_N}{\hat{\omega}_N}\) and, under the null hypothesis that the two models are equivalent, follows a standard normal distribution. This is the case for two strictly non-nested models.

On the contrary, if \(\omega^2_* = 0\), the distribution is much more complicated. We need to define two matrices: \(A\) contains the expected values of the second derivates of \(\Lambda\):

\[ A(\theta_*) = \mbox{E}^0\left[\frac{\partial^2 \Lambda}{\partial \theta \partial \theta ^ \top}\right] = \mbox{E}^0\left[\begin{array}{cc} \frac{\partial^2 \ln f}{\partial \beta \partial \beta ^ \top} & 0 \\ 0 & -\frac{\partial^2 \ln g}{\partial \beta \partial \beta ^ \top} \end{array}\right] = \left[ \begin{array}{cc} A_f(\beta_*) & 0 \\ 0 & - A_g(\gamma_*) \end{array} \right] \]

and \(B\) the variance of its first derivatives:

\[ B(\theta_*) = \mbox{E}^0\left[\frac{\partial \Lambda}{\partial \theta}\frac{\partial \Lambda}{\partial \theta ^ \top}\right]= \mbox{E}^0\left[ \left(\frac{\partial \ln f}{\partial \beta}, - \frac{\partial \ln g}{\partial \gamma} \right) \left(\frac{\partial \ln f}{\partial \beta ^ \top}, - \frac{\partial \ln g}{\partial \gamma ^ \top} \right) \right] = \mbox{E}^0\left[ \begin{array}{cc} \frac{\partial \ln f}{\partial \beta} \frac{\partial \ln f}{\partial \beta^\top} & - \frac{\partial \ln f}{\partial \beta} \frac{\partial \ln g}{\partial \gamma ^ \top} \\ - \frac{\partial \ln g}{\partial \gamma} \frac{\partial \ln f}{\partial \beta^\top} & \frac{\partial \ln g}{\partial \gamma} \frac{\partial \ln g}{\partial \gamma^\top} \end{array} \right] \]

or:

\[ B(\theta_*) = \left[ \begin{array}{cc} B_f(\beta_*) & - B_{fg}(\beta_*, \gamma_*) \\ - B_{gf}(\beta_*, \gamma_*) & B_g(\gamma_*) \end{array} \right] \]

Then:

\[ W(\theta_*) = B(\theta_*) \left[-A(\theta_*)\right] ^ {-1}= \left[ \begin{array}{cc} -B_f(\beta_*) A^{-1}_f(\beta_*) & - B_{fg}(\beta_*, \gamma_*) A^{-1}_g(\gamma_*) \\ B_{gf}(\gamma_*, \beta_*) A^{-1}_f(\beta_*) & B_g(\gamma_*) A^{-1}_g(\gamma_*) \end{array} \right] \]

Denote \(\lambda_*\) the eigen values of \(W\). When \(\omega_*^2 = 0\) (which is always the case for nested models), the statistic is the one used in the standard likelihood ratio test: \(2 (\ln L_f - \ln L_g) = 2 N \hat{\Lambda}_N\) which, under the null, follows a weighted \(\chi ^ 2\) distribution with weights equal to \(\lambda_*\). The Vuong test can be seen in this context as a more robust version of the standard likelihood ratio test, because it doesn’t assume, under the null, that the larger model is correctly specified.

Note that, if the larger model is correctly specified, the information matrix equality implies that \(B_f(\theta_*)-A_f(\theta_*)\). In this case, the two matrices on the diagonal of \(W\) reduce to \(-I_{K_f}\) and \(I_{K_g}\), the trace of \(W\) to \(K_g - K_f\) and the distribution of the statistic under the null reduce to a \(\chi^2\) with \(K_g - K_f\) degrees of freedom.

The \(W\) matrice can be consistently estimate by computing the first and the second derivatives of the likelihood functions of the two models for \(\hat{\theta}\). For example,

\[ \hat{A}_f(\hat{\beta}) = \frac{1}{N} \sum_{n= 1} ^ N \frac{\partial^2 \ln f}{\partial \beta \partial \beta ^ \top}(\hat{\beta}, x_n, y_n) \]

\[ \hat{B}_{fg}(\hat{\theta})= \frac{1}{N} \sum_{n=1}^N \frac{\partial \ln f}{\partial \beta}(\hat{\beta}, x_n, y_n) \frac{\partial \ln g}{\partial \gamma^\top}(\hat{\gamma}, x_n, y_n) \]

For the overlapping case, the test should be performed in two steps:

• the first step consists on testing whether \(\omega_*^*\) is 0 or not. This hypothesis is based on the statistic \(N \hat{\omega} ^ 2\) which, under the null (\(\omega_*^2=0\)) follows a weighted \(\chi ^ 2\) distributions with weights equal to \(\lambda_* ^ 2\). If the null hypothesis is not rejected, the test stops at this step and the conclusion is that the two models are equivalent,
• if the null hypothesis is reject, the second step consists on applying the test for non-nested models previously described.

The non-degenerate Vuong test

Shi (2015) proposed a non-degenerate version of the Vuong (1989) test. She shows that the Vuong test has size distortion, leading to subsequent overrejection. The cause of this problem is that the distribution of \(\hat{\Lambda}\) is discontinuous in the \(\omega^2\) parameter (namely a normal distribution if \(\omega^2 > 0\) and a distribution related to a weight \(\chi^2\) distribution if \(\omega^2 = 0\)). Especially in small samples, it may be difficult to distinguish a positive versus a zero value of \(\omega ^ 2\) because of sampling error. To solve this problem, using local asymptotic theory, Shi (2015) showed that, rewriting the Vuong statistic as:

\[ \hat{T} = \frac{N \hat{\Lambda}_N}{\sqrt{N \hat{\omega} ^ 2_N}} \]

the asymptotic distribution of the numerator and of the square of the denominator of the Vuong statistic is the same as:

\[ \left( \begin{array}{cc} N \hat{\Lambda}_N \\ N \hat{\omega} ^ 2 _ N \end{array} \right) \rightarrow^d \left( \begin{array}{cc} J_\Lambda \\ J_\omega \end{array} \right) = \left( \begin{array}{cc} \sigma z_\omega - z_\theta ^ \top V z_\theta / 2 \\ \sigma ^ 2 - 2 \sigma \rho_* ^ \top V z_\theta + z_\theta ^ \top V ^ 2 z_\theta \end{array} \right) \]

where:

\[ \left(\begin{array}{c}z_\omega \\ z_\theta \end{array}\right) \sim N \left(0, \left(\begin{array}{cc} 1 & \rho_* ^ \top \\ \rho_* & I \end{array}\right) \right), \]

\(\rho_*\) is a vector of length \(K_f + K_g\), \(\sigma\) a positive scalar and V is the diagonal matrix containing the eigen values of \(B ^ {\frac{1}{2}} A ^ {-1} B ^ {\frac{1}{2}}\).

Based on this result, Shi (2015) showed:

• that the expected value of the numerator is \(-\mbox{trace}(V) / 2\), the classical Vuong statistic is therefore biased and this bias can be severe in small samples and when the degree of parametrization of the two models are very different²,
• that the denominator, being random, can take values close to zero with a significant probability, which can generate fat tails in the distribution of the statistic.

Shi (2015) therefore proposed to modify the numerator of the Vuong statistic:

\[\hat{\Lambda}^{\mbox{mod}}_N = \hat{\Lambda}_N + \frac{\mbox{tr}(V)}{2 N}\]

and to add a constant to the denominator, so that:

\[ \left(\hat{\omega}^{\mbox{mod}}(c)\right) ^ 2 = \hat{\omega} ^ 2 + c \; \mbox{tr}(V) ^ 2 / N \]

The non-degenarate Vuong test is then:

\[ T_N^{\mbox{mod}} = \frac{\hat{\Lambda}^{\mbox{mod}}_N}{\hat{\omega}^{\mbox{mod}}}= \sqrt{N}\frac{\hat{\Lambda}_N + \mbox{tr}(V) / 2N}{\sqrt{\hat{\omega} ^ 2 + c \;\mbox{tr}(V) ^ 2 / N}} \]

The distribution of the modified Vuong statistic can be estimated by simulations: drawing in the distribution of \((z_\omega, z_\theta^\top)\), we compute for every draw \(J_\Lambda\), \(J_\omega\) and \(J_\Lambda / \sqrt{J_\omega}\). As \(\sigma\) and \(\rho_*\) can’t be estimated consistently, the supremum other these parameters are taken, and Shi (2015) indicates that \(\rho_*\) should be in this case a vector where all the elements are zero except for the one that coincides with the highest absolute value of \(V\) which is set to 1.

The Shi test is then computed as follow:

1. start with a given size for the test, say \(\alpha = 0.05\),
2. for a given value of \(c\), choose \(\sigma\) which maximize the simulated critical value for \(c\) and \(\alpha\),
3. adjust \(c\) so that this critical value equals the normal critical value, up to a small disperency (say 0.1); for example, if the size is 5%, the target is \(v_{1 - \alpha / 2} = 1.96 + 0.1 = 2.06\),
4. compute \(\hat{T}_N^{\mbox{mod}}\) for the given values of \(c\) and \(\sigma\) ; if \(\hat{T}_N^{\mbox{mod}} > v_{1 - \alpha / 2}\), reject the null hypothesis at the \(\alpha\) level,
5. to get a p-value, if \(\hat{T}_N^{\mbox{mod}} > v_{1 - \alpha / 2}\) increase \(\alpha\) and repeat the previous steps until a new value of \(\alpha\) is obtained so that \(\hat{T}_N^{\mbox{mod}} = v_{1 - \alpha^* / 2}\), \(\alpha^*\) being the p-value of the test.

Simulations

Shi (2015) provides an example of simulations of non-nested linear models that shows that the distribution of the Vuong statistic can be very different from a standard normal. The data generating process used for the simulations is:

\[ y = 1 + \sum_{k = 1} ^ {K_f} z^f_k + \sum_{k = 1} ^ {K_g} z^g_k + \epsilon \]

where \(z^f\) is the set of \(K_f\) covariates that are used in the first model and \(z^g\) the set of \(K_g\) covariates used in the second model and \(\epsilon \sim N(0, 1 - a ^ 2)\). \(z^f_k \sim N(0, a / \sqrt{K_f})\) and \(z^g_k \sim N(0, a / \sqrt{K_g})\), so that the explained variance explained by the two competing models is the same (equal to \(a ^ 2\)) and the null hypothesis of the Vuong test is true. The vuong_sim unables to simulate values of the Vuong test. As in Shi (2015), we use a very different degree of parametrization for the two models, with \(K_f = 15\) and \(K_G = 1\).

library(micsr)
Vuong <- vuong_sim(N = 100, R = 1000, Kf = 15, Kg = 1, a = 0.5)
head(Vuong)

## [1]  1.5588591  2.5625929 -1.0904569  0.7207765  2.7263322  1.1310495

mean(Vuong)

## [1] 1.119354

mean(abs(Vuong) > 1.96)

## [1] 0.207

We can see that the the mean of the statistic for the 1000 replications is far away from 0, which means that the numerator of the Vuong statistic is seriously biased. 20.7% of the values of the statistic are greater than the critical value so that the Vuong test will lead in such context a noticeable overrejection. The empirical pdf is shown in the following figure, along with the normal pdf.

library("ggplot2")
ggplot(data = data.frame(Vuong = Vuong)) + geom_density(aes(x = Vuong)) +
    geom_function(fun = dnorm, linetype = "dotted")

Implementation of the non-degenarate Vuong test

The micsr package provides a ndvuong function that implements the classical Vuong test. It has a nest argument (that is FALSE by default but can be set to TRUE to get the nested version of the Vuong test). This package also provide a llcont generic which returns a vector of length \(N\) containing the contribution of every observation to the log-likelihood.

The ndvuong package provides the ndvuong function. As for the vuongtest function, the two main arguments are two fitted models (say model1 and model2). The \(\hat{\Lambda}_n\) vector is obtained using llcont(model1) - llcont(model2). The relevant matrices \(A_i\) and \(B_i\) are computed from the fitted models using the estfun and the meat functions from the sandwich package. More precisely, \(A^{-1}\) is bdiag(-bread(model1), bread(model2)³ and \(B\) is crossprod(estfun(model1), - estfun(model2)) / N, where N is the sample size. Therefore, the ndvuong function can be used with any models for which a llcont, a estfun and a bread method is available.

Applications

library("micsr")
test <- ndvuong(turnout$group, turnout$intens, size = 0.05, pval = FALSE)
test

## 
##  Non-degenerate Vuong test for non-nested models
## 
## data:  turnout$group-turnout$intens
## z = 1.7759, size = 0.050000, vuong_stat = 2.084528, constant =
## 0.381107, crit-value = 2.059963, sum e.v. = -10.997224, vuong_p.value =
## 0.018556
## alternative hypothesis: different models

ndvuong(turnout$group, turnout$intens, nd = FALSE)

## 
##  Vuong test for non-nested models
## 
## data:  turnout$group-turnout$intens
## z = 2.0845, p-value = 0.01856
## alternative hypothesis: different models

test <- ndvuong(turnout$group, turnout$intens, pval = TRUE)
test

## 
##  Non-degenerate Vuong test for non-nested models
## 
## data:  turnout$group-turnout$intens
## z = 1.8125, vuong_stat = 2.084528, constant = 0.000000, sum e.v. =
## -10.997224, vuong_p.value = 0.018556, p-value = 0.0864
## alternative hypothesis: different models

library("mlogit")

## Le chargement a nécessité le package : dfidx

## 
## Attachement du package : 'dfidx'

## L'objet suivant est masqué depuis 'package:stats':
## 
##     filter

## 
## Attachement du package : 'mlogit'

## Les objets suivants sont masqués depuis 'package:micsr':
## 
##     rg, scoretest, stdev

data("ModeCanada")
MC <- mlogit.data(ModeCanada, subset = noalt == 4, chid.var = "case",
                  alt.var = "alt", drop.index = TRUE)

ml <- mlogit(choice ~ freq + cost + ivt + ovt | urban + income, MC, 
             reflevel = 'car', alt.subset = c("car", "train", "air"))

hl <- mlogit(choice ~ freq + cost + ivt + ovt | urban + income, MC, 
             reflevel = 'car', alt.subset = c("car", "train", "air"),
             heterosc = TRUE)
coef(summary(hl))

##                       Estimate  Std. Error   z-value     Pr(>|z|)
## (Intercept):train  0.678393435 0.332762598  2.038671 4.148288e-02
## (Intercept):air    0.656754399 0.468163091  1.402832 1.606668e-01
## freq               0.063924677 0.004916769 13.001360 0.000000e+00
## cost              -0.026961457 0.004283139 -6.294789 3.078178e-10
## ivt               -0.009680773 0.001053874 -9.185892 0.000000e+00
## ovt               -0.032165526 0.003593007 -8.952258 0.000000e+00
## urban:train        0.797131578 0.120739176  6.602096 4.053868e-11
## urban:air          0.445472634 0.082160945  5.421951 5.895197e-08
## income:train      -0.012597857 0.003994180 -3.154053 1.610196e-03
## income:air         0.018859983 0.003215926  5.864558 4.503294e-09
## sp.train           1.237182865 0.110460959 11.200182 0.000000e+00
## sp.air             0.540323852 0.111835294  4.831425 1.355592e-06

lmtest::waldtest(hl, heterosc = FALSE)

## 
##  Wald test
## 
## data:  homoscedasticity
## chisq = 25.196, df = 2, p-value = 3.38e-06

scoretest(ml, heterosc = TRUE)

## 
##  score test
## 
## data:  heterosc = TRUE
## chisq = 9.4883, df = 2, p-value = 0.008703
## alternative hypothesis: heteroscedastic model

lmtest::lrtest(hl, ml)

## Likelihood ratio test
## 
## Model 1: choice ~ freq + cost + ivt + ovt | urban + income
## Model 2: choice ~ freq + cost + ivt + ovt | urban + income
##   #Df  LogLik Df  Chisq Pr(>Chisq)  
## 1  12 -1838.1                       
## 2  10 -1841.6 -2 6.8882    0.03193 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ndvuong(hl, ml, nested = TRUE)

## 
##  Non-degenerate Vuong test for nested models
## 
## data:  hl-ml
## z = 0.4554, vuong_stat = 6.888241, vuong_p.value = 0.047927, p-value =
## 0.211
## alternative hypothesis: different models

library("mlogit")
data("RiskyTransport")
RT <- mlogit.data(RiskyTransport, shape = "long", choice = "choice",
                  chid.var = "chid", alt.var = "mode", id.var = "id")

ml <- mlogit(choice ~ cost, data = RT)
hl <- mlogit(choice ~ cost , data = RT, heterosc = TRUE)
nl <- mlogit(formula = choice ~ cost, data = RT,
             nests = list(fast = c("Helicopter", "Hovercraft"),
                          slow = c("WaterTaxi", "Ferry")),
             un.nest.el = TRUE)
modelsummary::msummary(list(multinomial = ml, heteroscedastic = hl, nested =  nl))

ndvuong(nl, hl, vartest = TRUE)

## 
##  variance test
## 
## data:  nl-hl
## w2 = 0.047327, p-value < 2.2e-16
## alternative hypothesis: positive variance

ndvuong(hl, nl)

## 
##  Non-degenerate Vuong test for non-nested models
## 
## data:  hl-nl
## z = 1.7298, vuong_stat = 1.829208, constant = 0.000000, sum e.v. =
## -1.832021, vuong_p.value = 0.033684, p-value = 0.0975
## alternative hypothesis: different models