# R/ndvuongtest.R In mhurdle: Multiple Hurdle Tobit Models

#### Documented in ndvuongtest

```gbdiag <- function(...){
if (nargs() == 1)
x <- as.list(...)
else
x <- list(...)
n <- length(x)
if(n==0) return(NULL)
x <- lapply(x, function(y) if(length(y)) as.matrix(y) else
stop("Zero-length component in x"))
d <- array(unlist(lapply(x, dim)), c(2, n))
rr <- d[1,]
cc <- d[2,]
rsum <- sum(rr)
csum <- sum(cc)
out <- array(0, c(rsum, csum))
ind <- array(0, c(4, n))
rcum <- cumsum(rr)
ccum <- cumsum(cc)
ind[1,-1] <- rcum[-n]
ind[2,] <- rcum
ind[3,-1] <- ccum[-n]
ind[4,] <- ccum
imat <- array(1:(rsum * csum), c(rsum, csum))
iuse <- apply(ind, 2, function(y, imat) imat[(y[1]+1):y[2],
(y[3]+1):y[4]], imat=imat)
iuse <- as.vector(unlist(iuse))
out[iuse] <- unlist(x)
return(out)
}

sqrtm <- function(x){
eigenx <- eigen(x)
#    if (any(eigenx\$values < 0)) stop("negative eigen values")
#    eigenx\$vectors %*% diag(sqrt(abs(eigenx\$values))) %*% t(eigenx\$vectors)
eigenx\$vectors %*% diag(sqrt(eigenx\$values)) %*% t(eigenx\$vectors)
}

#' Vuoung test for non-nested models
#'
#' The Vuong test is suitable to discriminate between two non-nested models.
#'
#'
#' @aliases ndvuongtest
#' @param x a first fitted model of class \code{"mhurdle"},
#' @param y a second fitted model of class \code{"mhurdle"},
#' @param size the size of the test,
#' @param pval should the p-value be computed ?
#' @param type the kind of test to be computed,
#' @param ndraws the number of draws for the simulations,
#' @param diffnorm a creuser,
#' @param seed the seed,
#' @param print.level the level of details to be printed.
#' @return an object of class \code{"htest"}
#' @seealso \code{vuong} in package \code{pscl}.
#' @references Vuong Q.H. (1989) Likelihood ratio tests for model selection and
#' non-nested hypothesis, Econometrica, vol.57(2), pp.307-33.
#' @keywords htest
#' @examples
#'
#' data("Interview", package = "mhurdle")
#' # dependent double hurdle model
#' dhm <- mhurdle(vacations ~ car + size | linc + linc2 | 0, Interview,
#'               dist = "ln", h2 = TRUE, method = "bhhh", corr = TRUE)
#'
#' # a double hurdle p-tobit model
#' ptm <- mhurdle(vacations ~ 0 | linc + linc2 | car + size, Interview,
#'               dist = "ln", h2 = TRUE, method = "bhhh", corr = TRUE)
#' vuongtest(dhm, ptm)
#' @export
ndvuongtest <- function(x, y, size = 0.05, pval = TRUE,
type = c("non-nested", "nested", "overlapping"),
ndraws = 1E04, diffnorm = 0.1, seed = 1,
print.level = 0){
type <- match.arg(type)
eigentol <- 1E-12
data.name <- c(
paste(deparse(substitute(x))),
paste(deparse(substitute(y)))
)
data.name <- paste(data.name, collapse = "-")
set.seed(seed)

if (!is.null(x\$lnl)) x\$logLik <- x\$lnl
if (!is.null(y\$lnl)) y\$logLik <- y\$lnl

fmodel <- x                    ;  gmodel <- y
hessian.f <- fmodel\$hessian    ;  hessian.g <- gmodel\$hessian
# this version is uncorrect where ther
hessian.f <- - solve(vcov(fmodel))    ;  hessian.g <- - solve(vcov(gmodel))

lnl.f <- fmodel\$logLik         ;  lnl.g <- gmodel\$logLik
K <- K.f + K.g
A <- bdiag(hessian.f, - hessian.g) / N
solveA <- bdiag(solve(hessian.f / N), - solve(hessian.g / N))
B <- crossprod(grad.tot) / N + diag(rep(eigentol, K))
sqrtB <- sqrtm(B)
V <- sqrtB %*% solveA %*% sqrtB
V <- eigen(V, symmetric = TRUE)\$values
V <- round(V, 5)
Vf <- - sum(eigen(solve(hessian.f / N) %*%
only.values = TRUE)\$values)
Vg <- - sum(eigen(solve(hessian.g / N) %*%
only.values = TRUE)\$values)
if (type != "nested"){
rho <- as.numeric((abs(V) - max(abs(V))) == 0)
rho <- rho / sqrt(sum(rho))
}
else rho <- rep(0, length(V))
SIGMA <- rbind(c(1 , rho),
cbind(rho, diag(rep(1, K))))
Z <- matrix(rnorm(ndraws * (K + 1)), ndraws)
#    Z <- Z %*% sqrtm(SIGMA + diag(rep(eigentol, K + 1)))
Z <- Z %*% sqrtm(SIGMA)
ZL <- Z[,   1]
ZP <- Z[, - 1]
trVsq <- sum(V ^ 2)
Vnmlzd <- V / sqrt(trVsq)
A1 <- ZL
A2 <- as.numeric(apply(ZP, 1, function(x) sum(Vnmlzd * x ^ 2) / 2) -
sum(Vnmlzd) / 2)
A3 <- as.numeric(apply(ZP, 1, function(x) sum(Vnmlzd * rho * x)))
A4 <- as.numeric(apply(ZP, 1, function(x) sum(x ^ 2 * Vnmlzd ^ 2)))

LR <- mean(lnl.f - lnl.g)
w2 <- mean( (lnl.f - lnl.g) ^ 2) - LR ^ 2
Tvuong <- sqrt(N) * LR / sqrt(w2)

Tmod <- function(sigma, cst){
# R draws in the distribution of the Vuong statistic
num <- sigma * A1 - A2
denom <- sigma ^ 2 - 2 * sigma * A3 + A4 + cst
num / sqrt(denom)
}

quant <- function(sigma, cst, size)
# compute the empirical quantile of level 1 - size in the
# distribution of the Vuong statistic
as.numeric(quantile(abs(Tmod(sigma, cst)), 1 - size))

sigstar <- function(cst, size)
# compute the value of sigma which maximize the empirical
# quantile of level 1 - size of the distribution of the Vuong
# statistic
optimize(function(x) quant(x, cst, size), c(0, 5), maximum = TRUE)\$maximum

seekpval <- function(size){
# for a given size, compute the constant so that the critical
# value equals the target
c.value <- 0
cv.value <- quant(sigstar(0, size), 0, size)
cv.normal <- qnorm(1 - size / 2)
cv.target <- cv.normal + diffnorm

if (cv.value < cv.target){
cv.value <- max(cv.value, cv.normal)
}
else {
froot <- function(cst) quant(sigstar(cst, size), cst, size) - cv.target
zo <- uniroot(froot, c(0, 10))
c.value <- zo\$root
cv.value <- quant(sigstar(c.value, size), c.value, size)
}
LRmod <- LR + sum(V) / (2 * N)
w2mod <- w2 + c.value * sum(V ^ 2) / N
Tnd <- sqrt(N) * (LRmod) / sqrt(w2mod)
pvalue <- 1 - ecdf(abs(Tmod(sigstar(c.value, size), c.value)))(abs(Tnd))
list(pvalue = pvalue, stat = Tnd,
constant = c.value, cv = cv.value)#,
}

if (type == "nested"){
LRmod <- LR + sum(V) / (2 * N)
Tnd <- sqrt(N) * (LRmod) / sqrt(w2)
if (Tnd < 0) pvalue <- ecdf(Tmod(0, 0))(Tnd)
else pvalue <- 1 - ecdf(Tmod(0, 0))(Tnd)
TndVuong <- structure(list(statistic = c(z = Tnd),
method = "Non-degenerate Vuong test",
p.value = pvalue,
data.name = data.name,
alternative = "different models"),
class = "htest")
trad <- vuongtest(x, y, type = "nested")
class = "ndvuongtest")
}
else{
if (pval){
froot <- function(alpha) seekpval(alpha)\$pvalue - alpha
pvalue <- uniroot(froot, c(1E-7, 1-1E-7))
results <- seekpval(pvalue\$root)
Tndvuong <- structure(list(statistic = c(z = as.numeric(results\$stat)),
method = "Non-degenerate Vuong test",
p.value = results\$pvalue,
data.name = data.name,
alternative = "different models",
parameters = c(constant = results\$constant,
'crit-value' = results\$cv,
'sum e.v.' = sum(V),
'tracef' = Vf,
'traceg' = Vg)),
class = "htest")
}
else{
results <- seekpval(size)
Tndvuong <- structure(list(statistic = c(z = as.numeric(results\$stat)),
method = "Non-degenerate Vuong test",
p.value = results\$pvalue,
data.name = data.name,
size = size,
alternative = "different models",
parameters = c(constant = results\$constant,
'crit. value' = results\$cv,
'sum e.v.' = sum(V))),
class = "htest")
}
Tvuong <- structure(list(statistic = c(z = Tvuong),
method = "Vuong test",
p.value = 2 * pnorm(abs(Tvuong), lower.tail = FALSE),
data.name = data.name,
alternative = "different models"),
class = "htest")

structure(list(trad = Tvuong, nd = Tndvuong), class = "ndvuongtest")
}
}

print.ndvuongtest <- function(x, model = c("nd", "trad", "both"), ...){
model <- match.arg(model)
if (model != "nd") print(x\$trad, ...)
if (model != "trad") print(x\$nd, ...)
}
```

## Try the mhurdle package in your browser

Any scripts or data that you put into this service are public.

mhurdle documentation built on Dec. 11, 2021, 9:21 a.m.