Last updated on 2020-02-01 08:49:10 CET.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 0.1-8 | 4.40 | 190.33 | 194.73 | ERROR | |
| r-devel-linux-x86_64-debian-gcc | 0.1-8 | 3.35 | 235.52 | 238.87 | OK | |
| r-devel-linux-x86_64-fedora-clang | 0.1-8 | 234.94 | ERROR | |||
| r-devel-linux-x86_64-fedora-gcc | 0.1-8 | 232.12 | ERROR | |||
| r-devel-windows-ix86+x86_64 | 0.1-8 | 8.00 | 415.00 | 423.00 | OK | |
| r-devel-windows-ix86+x86_64-gcc8 | 0.1-8 | 13.00 | 504.00 | 517.00 | OK | |
| r-patched-linux-x86_64 | 0.1-8 | 3.90 | 301.12 | 305.02 | OK | |
| r-patched-solaris-x86 | 0.1-8 | 429.50 | OK | |||
| r-release-linux-x86_64 | 0.1-8 | 3.53 | 296.16 | 299.69 | OK | |
| r-release-windows-ix86+x86_64 | 0.1-8 | 10.00 | 493.00 | 503.00 | OK | |
| r-release-osx-x86_64 | 0.1-8 | OK | ||||
| r-oldrel-windows-ix86+x86_64 | 0.1-8 | 5.00 | 389.00 | 394.00 | OK | |
| r-oldrel-osx-x86_64 | 0.1-8 | OK |
Version: 0.1-8
Check: examples
Result: ERROR
Running examples in 'mvProbit-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: mvProbit
> ### Title: Estimation of Multivariate Probit Models
> ### Aliases: mvProbit print.mvProbit
> ### Keywords: models regression
>
> ### ** Examples
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 50
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -0.8 ),
+ c( -0.6, 1.0, -1.6 ),
+ c( 0.5, -0.6, 1.2 ) )
>
> # covariance matrix of error terms
> library( miscTools )
> sigma <- symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
>
> # estimation (BHHH optimizer and GHK algorithm)
> estResult <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ data = as.data.frame( cbind( xMat, yMat ) ), iterlim = 1, nGHK = 50 )
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
...)
where 3: fnOrig(theta, ...)
where 4: fn(start1, fixed = fixed, sumObs = TRUE, returnHessian = returnHessian,
...)
where 5: maxNRCompute(fn = function (theta, fnOrig, gradOrig = NULL, hessOrig = NULL,
fixed, sumObs = FALSE, returnHessian = TRUE, ...)
{
nParam <- length(theta)
f <- fnOrig(theta, ...)
if (any(is.na(f))) {
attr(f, "gradient") <- NA
attr(f, "hessian") <- NA
return(f)
}
gr <- attr(f, "gradient")
if (is.null(gr)) {
if (!is.null(gradOrig)) {
gr <- gradOrig(theta, ...)
}
else {
gr <- numericGradient(f = fnOrig, t0 = theta, fixed = fixed,
...)
}
}
if (is.matrix(gr)) {
if (ncol(gr) != length(theta)) {
stop(paste0("if gradient is a matrix, it must have length(parameter) colums (currently ",
length(theta), "), not ", ncol(gr)))
}
activeGr <- gr[, !fixed]
}
else {
activeGr <- gr[!fixed]
}
if (any(is.na(activeGr))) {
attr(f, "gradient") <- gr
attr(f, "hessian") <- NA
return(f)
}
if (observationGradient(gr, length(theta))) {
gr <- as.matrix(gr)
}
if (is.null(dim(gr))) {
gr[fixed] <- NA
}
else {
gr[, fixed] <- NA
}
if (isTRUE(returnHessian)) {
h <- attr(f, "hessian")
if (is.null(h)) {
if (!is.null(hessOrig)) {
h <- as.matrix(hessOrig(theta, ...))
}
else {
llFunc <- function(theta, ...) {
return(sum(fnOrig(theta, ...)))
}
if (!is.null(attr(f, "gradient"))) {
gradFunc <- function(theta, ...) {
return(sumGradients(attr(fnOrig(theta, ...),
"gradient"), nParam))
}
}
else if (!is.null(gradOrig)) {
gradFunc <- function(theta, ...) {
return(sumGradients(gradOrig(theta, ...),
nParam))
}
}
else {
gradFunc <- NULL
}
h <- numericHessian(f = llFunc, grad = gradFunc,
t0 = theta, fixed = fixed, ...)
}
}
if ((dim(h)[1] != nParam) | (dim(h)[2] != nParam)) {
stop("Wrong hessian dimension. Needed ", nParam,
"x", nParam, " but supplied ", dim(h)[1], "x",
dim(h)[2])
}
else {
h[fixed, ] <- NA
h[, fixed] <- NA
}
}
else if (tolower(returnHessian) == "bhhh") {
h <- NULL
if (is.null(dim(gr)) & any(is.na(gr[!fixed]))) {
h <- NA
}
else if (is.matrix(gr)) {
if (any(is.na(gr[, !fixed]))) {
h <- NA
}
}
if (is.null(h)) {
checkBhhhGrad(g = gr, theta = theta, analytic = (!is.null(attr(f,
"gradient")) || !is.null(gradOrig)), fixed = fixed)
h <- -crossprod(gr)
}
attr(h, "type") = "BHHH"
}
else {
h <- NULL
}
if (sumObs) {
f <- sumKeepAttr(f)
}
if (sumObs) {
gr <- sumGradients(gr, nParam)
}
if (!is.null(gradOrig) && !is.null(attr(f, "gradient"))) {
attr(f, "gradBoth") <- TRUE
}
if (!is.null(hessOrig) && !is.null(attr(f, "hessian"))) {
attr(f, "hessBoth") <- TRUE
}
attr(f, "gradient") <- gr
attr(f, "hessian") <- h
return(f)
}, fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, start = c(b_1_0 = 0.598378200054935, b_1_1 = 6.336646267785,
b_1_2 = -1.01231147815159, b_2_0 = -1.03564793824529, b_2_1 = 1.26302210717712,
b_2_2 = -1.70638777681512, b_3_0 = 0.410893287564513, b_3_1 = -0.373151392833251,
b_3_2 = 2.92074017047352, R_1_2 = 0.172617843070694, R_1_3 = -0.0160393443285098,
R_2_3 = -0.0704265674853818), finalHessian = "BHHH", bhhhHessian = TRUE,
fixed = c(b_1_0 = FALSE, b_1_1 = FALSE, b_1_2 = FALSE, b_2_0 = FALSE,
b_2_1 = FALSE, b_2_2 = FALSE, b_3_0 = FALSE, b_3_1 = FALSE,
b_3_2 = FALSE, R_1_2 = FALSE, R_1_3 = FALSE, R_2_3 = FALSE
), control = new("MaxControl", tol = 1e-08, reltol = 1.49011611938477e-08,
gradtol = 1e-06, steptol = 1e-10, lambdatol = 1e-06,
qrtol = 1e-10, qac = "stephalving", marquardt_lambda0 = 0.01,
marquardt_lambdaStep = 2, marquardt_maxLambda = 1e+12,
nm_alpha = 1, nm_beta = 0.5, nm_gamma = 2, sann_cand = NULL,
sann_temp = 10, sann_tmax = 10L, sann_randomSeed = 123L,
iterlim = 1L, printLevel = 0L), yMat = c(0, 1, 1, 1,
1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0,
1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1,
1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0), xMat = c(1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
0.584613749636069, 0.123854243844614, 0.215941568743973,
0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258,
0.448209778629426, 0.0530042267305041, 0.922267467879738,
2.05008468562714, -0.491031166056535, -2.30916887564081,
1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633,
-0.370660031792409, 0.644376548518833, -0.220486561818751,
0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
-0.627906076039371, 1.36065244853001, -0.600259587147127,
2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
), llAlgorithm = "GHK", llNGHK = 50, llIntGrad = TRUE, llOneSidedGrad = FALSE,
llEps = 1e-06, llRandom.seed = 123)
where 6: eval(as.call(cl))
where 7: eval(as.call(cl))
where 8: maxNR(fn = fn, grad = grad, hess = hess, start = start, finalHessian = finalHessian,
bhhhHessian = TRUE, ...)
where 9: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
constraints = constraints, ...)
where 10: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
llRandom.seed = random.seed, ...)
where 11: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, data = as.data.frame(cbind(xMat,
yMat)), iterlim = 1, nGHK = 50)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x2590cb0>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbit ... fn -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 0.1-8
Check: tests
Result: ERROR
Running 'mvProbitEst.R' [29s/33s]
Running 'mvProbitMargEff2.R' [2s/2s]
Running 'mvProbitTest.R' [3s/3s]
Running 'mvProbitTest5.R' [119s/130s]
Comparing 'mvProbitTest5.Rout' to 'mvProbitTest5.Rout.save' ...1604c1604
< 5 0.00 -0.01 0 0.00 0.00 0.00 0.00
---
> 5 0.00 0.00 0 0.00 0.00 0.00 0.00
1618,1620c1618,1620
< 1 0 0.00 0 0 0 0 0
< 5 0 -0.01 0 0 0 0 0
< 10 0 0.00 0 0 0 0 0
---
> 1 0 0 0 0 0 0 0
> 5 0 0 0 0 0 0 0
> 10 0 0 0 0 0 0 0
1728c1728
< 5 -0.01
---
> 5 0.00
Running 'pmvnormWrapTest.R' [2s/3s]
Running the tests in 'tests/mvProbitEst.R' failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> options( digits = 4 )
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 50
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -0.8 ),
+ c( -0.6, 1.0, -1.6 ),
+ c( 0.5, -0.6, 1.2 ) )
>
> # covariance matrix of error terms
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
>
> # create data frame
> dat <- as.data.frame( cbind( xMat, yMat ) )
>
> # estimation with the BHHH algorithm, two-sided gradients
> estResultBHHH <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz() )
> print( estResultBHHH )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
b_3_2 R_1_2 R_1_3 R_2_3
1.23283 0.86322 0.36022 -0.08406
> summary( estResultBHHH )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 1.02508 0.44548 2.301 0.0214 *
b_1_1 1.30536 0.73511 1.776 0.0758 .
b_1_2 -1.02809 0.68851 -1.493 0.1354
b_2_0 -0.78248 0.38798 -2.017 0.0437 *
b_2_1 1.14318 0.53026 2.156 0.0311 *
b_2_2 -1.66938 0.81308 -2.053 0.0401 *
b_3_0 0.58229 0.64107 0.908 0.3637
b_3_1 -0.69667 0.89818 -0.776 0.4380
b_3_2 1.23283 0.62452 1.974 0.0484 *
R_1_2 0.86322 18.60071 0.046 0.9630
R_1_3 0.36022 0.64712 0.557 0.5778
R_2_3 -0.08406 0.34909 -0.241 0.8097
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BHHH maximisation, 1 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -56.91 on 12 Df
> logLik( estResultBHHH )
'log Lik.' -56.91 (df=12)
> estResultBHHHA <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta, sigma[ lower.tri( sigma ) ] ),
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz() )
> all.equal( estResultBHHH, estResultBHHHA )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BHHH algorithm, one-sided gradients
> estResultBHHH1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBHHH1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
tol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
b_3_2 R_1_2 R_1_3 R_2_3
1.23283 0.86322 0.36022 -0.08406
> summary( estResultBHHH1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
tol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 1.02508 0.44548 2.301 0.0214 *
b_1_1 1.30536 0.73511 1.776 0.0758 .
b_1_2 -1.02809 0.68851 -1.493 0.1354
b_2_0 -0.78248 0.38797 -2.017 0.0437 *
b_2_1 1.14318 0.53026 2.156 0.0311 *
b_2_2 -1.66938 0.81307 -2.053 0.0401 *
b_3_0 0.58229 0.64106 0.908 0.3637
b_3_1 -0.69667 0.89817 -0.776 0.4380
b_3_2 1.23283 0.62452 1.974 0.0484 *
R_1_2 0.86322 18.60041 0.046 0.9630
R_1_3 0.36022 0.64712 0.557 0.5778
R_2_3 -0.08406 0.34909 -0.241 0.8097
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BHHH maximisation, 1 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -56.91 on 12 Df
> logLik( estResultBHHH1 )
'log Lik.' -56.91 (df=12)
> all.equal( estResultBHHH, estResultBHHH1, tol = 1e-5 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, two-sided gradients
> estResultBFGS <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma, method = "BFGS",
+ data = dat,
+ reltol = 0.5, algorithm = GenzBretz() )
> print( estResultBFGS )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53607 1.16101 -1.14212 -0.44952 1.05496 -1.57624 0.40751 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25958 0.29483 0.32822 -0.02382
> summary( estResultBFGS )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53607 0.82709 0.648 0.5169
b_1_1 1.16101 0.91515 1.269 0.2046
b_1_2 -1.14212 0.93847 -1.217 0.2236
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88398 1.193 0.2327
b_2_2 -1.57624 0.73187 -2.154 0.0313 *
b_3_0 0.40751 0.48897 0.833 0.4046
b_3_1 -0.85400 0.84871 -1.006 0.3143
b_3_2 1.25958 0.68040 1.851 0.0641 .
R_1_2 0.29483 1.85817 0.159 0.8739
R_1_3 0.32822 0.70923 0.463 0.6435
R_2_3 -0.02382 0.34466 -0.069 0.9449
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGS )
'log Lik.' -52.94 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients
> estResultBFGS1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma, method = "BFGS",
+ data = dat,
+ reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53607 1.16101 -1.14213 -0.44952 1.05496 -1.57624 0.40751 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25958 0.29483 0.32822 -0.02382
> summary( estResultBFGS1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53607 0.82709 0.648 0.5169
b_1_1 1.16101 0.91515 1.269 0.2046
b_1_2 -1.14213 0.93847 -1.217 0.2236
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88398 1.193 0.2327
b_2_2 -1.57624 0.73187 -2.154 0.0313 *
b_3_0 0.40751 0.48897 0.833 0.4046
b_3_1 -0.85400 0.84871 -1.006 0.3143
b_3_2 1.25958 0.68040 1.851 0.0641 .
R_1_2 0.29483 1.85817 0.159 0.8739
R_1_3 0.32822 0.70923 0.463 0.6435
R_2_3 -0.02382 0.34466 -0.069 0.9449
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGS1 )
'log Lik.' -52.94 (df=12)
> all.equal( estResultBFGS, estResultBFGS1, tol = 1e-5 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values
> estResultBFGS1a <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ data = dat, method = "BFGS",
+ reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1a )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.79353 1.88367 -1.74922 -0.40036 0.91577 -1.52106 0.97790 -1.69433
b_3_2 R_1_2 R_1_3 R_2_3
1.57625 0.75270 0.23721 -0.03453
> summary( estResultBFGS1a )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.79353 0.90412 0.878 0.3801
b_1_1 1.88367 0.90668 2.078 0.0378 *
b_1_2 -1.74922 0.97400 -1.796 0.0725 .
b_2_0 -0.40036 0.69056 -0.580 0.5621
b_2_1 0.91577 0.97998 0.934 0.3501
b_2_2 -1.52106 0.72006 -2.112 0.0347 *
b_3_0 0.97790 0.64497 1.516 0.1295
b_3_1 -1.69433 0.94788 -1.787 0.0739 .
b_3_2 1.57625 0.73399 2.148 0.0318 *
R_1_2 0.75270 4.32729 0.174 0.8619
R_1_3 0.23721 0.90535 0.262 0.7933
R_2_3 -0.03453 0.35889 -0.096 0.9234
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 3 iterations
Return code 0: successful convergence
Log-likelihood: -49.36 on 12 Df
> logLik( estResultBFGS1a )
'log Lik.' -49.36 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values for beta
> estResultBFGS1b <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ startSigma = sigma, data = dat,
+ method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1b )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1 b_3_2 R_1_2
0.8110 1.9060 -1.7627 -0.3770 0.9441 -1.5332 0.9372 -1.7498 1.6024 0.6817
R_1_3 R_2_3
0.4044 0.1663
> summary( estResultBFGS1b )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.8110 0.9297 0.872 0.3830
b_1_1 1.9060 0.9160 2.081 0.0375 *
b_1_2 -1.7627 1.0131 -1.740 0.0819 .
b_2_0 -0.3770 0.7049 -0.535 0.5928
b_2_1 0.9441 1.0345 0.913 0.3615
b_2_2 -1.5332 0.7517 -2.040 0.0414 *
b_3_0 0.9372 0.6291 1.490 0.1363
b_3_1 -1.7498 0.9996 -1.751 0.0800 .
b_3_2 1.6024 0.7848 2.042 0.0412 *
R_1_2 0.6817 3.4094 0.200 0.8415
R_1_3 0.4044 0.8024 0.504 0.6143
R_2_3 0.1663 0.3389 0.491 0.6235
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 3 iterations
Return code 0: successful convergence
Log-likelihood: -49.57 on 12 Df
> logLik( estResultBFGS1b )
'log Lik.' -49.57 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values for sigma
> estResultBFGS1s <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), data = dat,
+ method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1s )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.55814 1.14958 -1.12747 -0.44426 1.07325 -1.57226 0.37255 -0.85470
b_3_2 R_1_2 R_1_3 R_2_3
1.23011 0.34981 -0.01177 0.09435
> summary( estResultBFGS1s )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.55814 0.86615 0.644 0.5193
b_1_1 1.14958 0.90607 1.269 0.2045
b_1_2 -1.12747 0.94191 -1.197 0.2313
b_2_0 -0.44426 0.85526 -0.519 0.6035
b_2_1 1.07325 0.99314 1.081 0.2798
b_2_2 -1.57226 0.78557 -2.001 0.0453 *
b_3_0 0.37255 0.56172 0.663 0.5072
b_3_1 -0.85470 0.95214 -0.898 0.3694
b_3_2 1.23011 0.71506 1.720 0.0854 .
R_1_2 0.34981 2.02289 0.173 0.8627
R_1_3 -0.01177 0.70776 -0.017 0.9867
R_2_3 0.09435 0.33920 0.278 0.7809
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.78 on 12 Df
> logLik( estResultBFGS1s )
'log Lik.' -52.78 (df=12)
>
> # estimation with the BFGS algorithm, Miwa algorithm for obtaining integrals
> estResultBFGSm <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, method = "BFGS",
+ reltol = 0.5, algorithm = Miwa( steps = 64 ) )
There were 13 warnings (use warnings() to see them)
> print( estResultBFGSm )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53609 1.16102 -1.14213 -0.44952 1.05496 -1.57624 0.40753 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25960 0.29491 0.32819 -0.02386
> summary( estResultBFGSm )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53609 0.82702 0.648 0.5168
b_1_1 1.16102 0.91511 1.269 0.2045
b_1_2 -1.14213 0.93837 -1.217 0.2235
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88386 1.194 0.2326
b_2_2 -1.57624 0.73173 -2.154 0.0312 *
b_3_0 0.40753 0.48894 0.834 0.4046
b_3_1 -0.85400 0.84860 -1.006 0.3142
b_3_2 1.25960 0.68024 1.852 0.0641 .
R_1_2 0.29491 1.85789 0.159 0.8739
R_1_3 0.32819 0.70903 0.463 0.6435
R_2_3 -0.02386 0.34468 -0.069 0.9448
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGSm )
'log Lik.' -52.94 (df=12)
> all.equal( estResultBFGS, estResultBFGSm, tol = 1e-3 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, GHK algorithm for obtaining integrals
> estResultBFGSg <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, method = "BFGS",
+ reltol = 0.5 )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
...)
where 3: fnOrig(theta, ...)
where 4: logLikFunc(theta, fnOrig = function (param, yMat, xMat, llAlgorithm,
llNGHK, llIntGrad, llOneSidedGrad, llEps, llRandom.seed,
...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, gradOrig = ..2, hessOrig = ..3, yMat = c(1, 1, 1, 1, 1, 1,
1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0,
1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0,
0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0,
0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1,
0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,
0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1,
0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1,
1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0.253318513994755,
-0.028546755348703, -0.0428704572913161, 1.36860228401446, -0.225770985659268,
1.51647060442954, -1.54875280423022, 0.584613749636069, 0.123854243844614,
0.215941568743973, 0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258, 0.448209778629426,
0.0530042267305041, 0.922267467879738, 2.05008468562714, -0.491031166056535,
-2.30916887564081, 1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633, -0.370660031792409,
0.644376548518833, -0.220486561818751, 0.331781963915697, 1.09683901314935,
0.435181490833803, -0.325931585531227, 1.14880761845109, 0.993503855962119,
0.54839695950807, 0.238731735111441, -0.627906076039371, 1.36065244853001,
-0.600259587147127, 2.18733299301658, 1.53261062618519, -0.235700359100477,
-1.02642090030678), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
where 5: eval(f, sys.frame(sys.parent()))
where 6: eval(f, sys.frame(sys.parent()))
where 7: callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
...)
where 8: (function (theta, fName, ...)
{
return(callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
...))
})(theta = c(b_1_0 = 0.8, b_1_1 = 1.2, b_1_2 = -0.8, b_2_0 = -0.6,
b_2_1 = 1, b_2_2 = -1.6, b_3_0 = 0.5, b_3_1 = -0.6, b_3_2 = 1.2,
R_1_2 = 0.2, R_1_3 = 0.4, R_2_3 = -0.1), fName = "logLikFunc",
fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, gradOrig = NULL, hessOrig = NULL, yMat = c(1, 1, 1, 1,
1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0,
1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0,
0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0,
0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1,
1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
0.584613749636069, 0.123854243844614, 0.215941568743973,
0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258,
0.448209778629426, 0.0530042267305041, 0.922267467879738,
2.05008468562714, -0.491031166056535, -2.30916887564081,
1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633,
-0.370660031792409, 0.644376548518833, -0.220486561818751,
0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
-0.627906076039371, 1.36065244853001, -0.600259587147127,
2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
where 9: do.call(callWithoutSumt, argList)
where 10: maxOptim(fn = fn, grad = grad, hess = hess, start = start, method = "BFGS",
fixed = fixed, constraints = constraints, finalHessian = finalHessian,
parscale = parscale, control = mControl, ...)
where 11: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
constraints = constraints, ...)
where 12: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
llRandom.seed = random.seed, ...)
where 13: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, start = c(beta), startSigma = sigma,
data = dat, method = "BFGS", reltol = 0.5)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x2a96c80>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbit ... logLikFunc -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
Execution halted
Running the tests in 'tests/mvProbitMargEff2.R' failed.
Complete output:
> # I thank Mohit Batham for providing this R script that demonstrated
> # a bug in mvProbitMargEff() when called with 2 dependent variables
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> nObs <- 100
> set.seed( 123 )
> xData <- data.frame(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = as.numeric( rnorm( nObs ) > 0 ),
+ x3 = rnorm( nObs ),
+ x4 = rnorm( nObs ))
> beta <- c( 0.8, 1.2, -1.0, 1.4, -0.8,
+ -0.6, 1.0, 0.6, -1.2, -1.6)
> sigma <- symMatrix( c( 1, 0.2, 1))
> margEffCond <- try( mvProbitMargEff( ~ x1 + x2 + x3 + x4 , coef = beta,
+ sigma = sigma, data = xData, cond = TRUE ) )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = betaLower,
sigma = coef$sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
random.seed = random.seed, ...)
where 3: mvProbitMargEffInternal(yMat = yMat, xMat = xMat, coef = coef,
sigma = sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
eps = eps, dummyVars = dummyVars, random.seed = random.seed,
...)
where 4: mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
data = xData, cond = TRUE)
where 5: doTryCatch(return(expr), name, parentenv, handler)
where 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 7: tryCatchList(expr, classes, parentenv, handlers)
where 8: tryCatch(expr, error = function(e) {
call <- conditionCall(e)
if (!is.null(call)) {
if (identical(call[[1L]], quote(doTryCatch)))
call <- sys.call(-4L)
dcall <- deparse(call)[1L]
prefix <- paste("Error in", dcall, ": ")
LONG <- 75L
sm <- strsplit(conditionMessage(e), "\n")[[1L]]
w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w")
if (is.na(w))
w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L],
type = "b")
if (w > LONG)
prefix <- paste0(prefix, "\n ")
}
else prefix <- "Error : "
msg <- paste0(prefix, conditionMessage(e), "\n")
.Internal(seterrmessage(msg[1L]))
if (!silent && isTRUE(getOption("show.error.messages"))) {
cat(msg, file = outFile)
.Internal(printDeferredWarnings())
}
invisible(structure(msg, class = "try-error", condition = e))
})
where 9: try(mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
data = xData, cond = TRUE))
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x282fba8>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
> round( margEffCond, 3 )
Error in round(margEffCond, 3) :
non-numeric argument to mathematical function
Execution halted
Running the tests in 'tests/mvProbitTest.R' failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> options( digits = 4 )
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 10
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = as.numeric( rnorm( nObs ) > 0 ),
+ x3 = rnorm( nObs ),
+ x4 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -1.0, 1.4, -0.8 ),
+ c( -0.6, 1.0, 0.6, -1.2, -1.6 ),
+ c( 0.5, -0.6, -0.7, 1.1, 1.2 ) )
>
> # covariance matrix of error terms
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # all parameters in a vector
> allCoef <- c( c( beta ), sigma[ lower.tri( sigma ) ] )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
> # (yMatLin > 0 )== yMat
>
> # unconditional expectations of dependent variables
> yExp <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ) )
> round( yExp, 3 )
V1 V2 V3
1 0.021 0.725 0.194
2 0.394 0.768 0.214
3 0.125 0.788 0.196
4 0.235 0.681 0.292
5 0.680 0.435 0.579
6 0.028 0.973 0.034
7 0.958 0.186 0.784
8 0.856 0.247 0.724
9 0.061 0.968 0.034
10 0.998 0.067 0.923
> yExpA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
+ data = as.data.frame( xMat ) )
> all.equal( yExp, yExpA )
[1] TRUE
> yExp2 <- pnorm( yMatLin )
> all.equal( yExp, as.data.frame( yExp2 ) )
[1] TRUE
>
>
> # conditional expectations of dependent variables
> # (assuming that all other dependent variables are one)
> yExpCond <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz() )
> round( yExpCond, 3 )
V1 V2 V3
1 0.072 0.834 0.524
2 0.660 0.779 0.314
3 0.297 0.838 0.399
4 0.448 0.727 0.461
5 0.847 0.442 0.618
6 0.139 0.985 0.162
7 0.991 0.179 0.749
8 0.950 0.244 0.710
9 0.244 0.978 0.133
10 1.000 0.065 0.892
> yExpCondA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
+ data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz() )
> all.equal( yExpCond, yExpCondA )
[1] TRUE
> yExpCond2 <- matrix( NA, nrow = nObs, ncol = ncol( yMat ) )
> for( i in 1:nObs ) {
+ for( k in 1:ncol( yMat ) ) {
+ set.seed( 123 )
+ numerator <- pmvnorm( upper = yMatLin[ i, ], sigma = sigma )
+ set.seed( 123 )
+ denominator <- pmvnorm( upper = yMatLin[ i, -k ], sigma = sigma[ -k, -k ] )
+ yExpCond2[ i, k ] <- numerator / denominator
+ }
+ }
> all.equal( yExpCond, as.data.frame( yExpCond2 ) )
[1] TRUE
> # now with explicitly specifying the algorithm
> yExpCond3 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz )
> all.equal( yExpCond, yExpCond3 )
[1] TRUE
> identical( yExpCond, yExpCond3 )
[1] TRUE
> # now with integrals obtained by the Miwa algorithm
> yExpCond4 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = Miwa )
> all.equal( yExpCond, yExpCond4, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the Miwa algorithm, less precise
> yExpCond5 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = Miwa( steps = 32 ) )
> all.equal( yExpCond4, yExpCond5, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the TVPACK algorithm
> yExpCond6 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = TVPACK )
> all.equal( yExpCond, yExpCond6, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the TVPACK algorithm, less precise
> yExpCond7 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = TVPACK( abseps = 0.5 ) )
> all.equal( yExpCond6, yExpCond7, tol = 1e-3 )
[1] "Component \"V1\": Mean relative difference: 0.03018"
[2] "Component \"V2\": Mean relative difference: 0.04273"
[3] "Component \"V3\": Mean relative difference: 0.04242"
> # now with integrals obtained by the GHK algorithm
> yExpCond8 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = coef, sigma = sigma,
cond = cond, algorithm = algorithm, nGHK = nGHK, random.seed = random.seed,
...)
where 3: mvProbitExp(~x1 + x2 + x3 + x4, coef = c(beta), sigma = sigma,
data = as.data.frame(xMat), cond = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0xa78670>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbitExp -> mvProbitExpInternal -> pmvnormWrap
Execution halted
Running the tests in 'tests/pmvnormWrapTest.R' failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> library( "miscTools" )
>
> # covariance matrix
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> ######## only upper ##########
> upper <- c( -0.3, 0.7, -0.5 )
> # Genz + Bretz (default)
> pug <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = GenzBretz(), random.seed = 123 )
> print( pug )
[1] 0.1360783
attr(,"error")
[1] 0.000118605
attr(,"msg")
[1] "Normal Completion"
>
> # Miwa (as function)
> pum <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = Miwa, random.seed = 123 )
> print( pum )
[1] 0.136069
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> all.equal( pug, pum, check.attributes = FALSE, tol = 1e-4 )
[1] TRUE
>
> # Miwa (as object returned from function Miwa())
> pum1 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = Miwa(), random.seed = 123 )
> all.equal( pum, pum1 )
[1] TRUE
>
> # Miwa (as character string)
> pum2 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = "Miwa", random.seed = 123 )
> all.equal( pum, pum2 )
[1] TRUE
>
> # TVPACK
> put <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = TVPACK, random.seed = 123 )
> print( put )
[1] 0.136069
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> all.equal( pug, put, check.attributes = FALSE, tol = 1e-4 )
[1] TRUE
> all.equal( pum, put, check.attributes = FALSE, tol = 1e-6 )
[1] TRUE
>
> # GHK
> pughk <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = "ghk", random.seed = 123, nGHK = 1000 )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
random.seed = 123, nGHK = 1000)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
random.seed = 123, nGHK = 1000)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x46c91e0>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: <Anonymous>
Execution halted
Flavor: r-devel-linux-x86_64-debian-clang
Version: 0.1-8
Check: examples
Result: ERROR
Running examples in ‘mvProbit-Ex.R’ failed
The error most likely occurred in:
> ### Name: mvProbit
> ### Title: Estimation of Multivariate Probit Models
> ### Aliases: mvProbit print.mvProbit
> ### Keywords: models regression
>
> ### ** Examples
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 50
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -0.8 ),
+ c( -0.6, 1.0, -1.6 ),
+ c( 0.5, -0.6, 1.2 ) )
>
> # covariance matrix of error terms
> library( miscTools )
> sigma <- symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
>
> # estimation (BHHH optimizer and GHK algorithm)
> estResult <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ data = as.data.frame( cbind( xMat, yMat ) ), iterlim = 1, nGHK = 50 )
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
...)
where 3: fnOrig(theta, ...)
where 4: fn(start1, fixed = fixed, sumObs = TRUE, returnHessian = returnHessian,
...)
where 5: maxNRCompute(fn = function (theta, fnOrig, gradOrig = NULL, hessOrig = NULL,
fixed, sumObs = FALSE, returnHessian = TRUE, ...)
{
nParam <- length(theta)
f <- fnOrig(theta, ...)
if (any(is.na(f))) {
attr(f, "gradient") <- NA
attr(f, "hessian") <- NA
return(f)
}
gr <- attr(f, "gradient")
if (is.null(gr)) {
if (!is.null(gradOrig)) {
gr <- gradOrig(theta, ...)
}
else {
gr <- numericGradient(f = fnOrig, t0 = theta, fixed = fixed,
...)
}
}
if (is.matrix(gr)) {
if (ncol(gr) != length(theta)) {
stop(paste0("if gradient is a matrix, it must have length(parameter) colums (currently ",
length(theta), "), not ", ncol(gr)))
}
activeGr <- gr[, !fixed]
}
else {
activeGr <- gr[!fixed]
}
if (any(is.na(activeGr))) {
attr(f, "gradient") <- gr
attr(f, "hessian") <- NA
return(f)
}
if (observationGradient(gr, length(theta))) {
gr <- as.matrix(gr)
}
if (is.null(dim(gr))) {
gr[fixed] <- NA
}
else {
gr[, fixed] <- NA
}
if (isTRUE(returnHessian)) {
h <- attr(f, "hessian")
if (is.null(h)) {
if (!is.null(hessOrig)) {
h <- as.matrix(hessOrig(theta, ...))
}
else {
llFunc <- function(theta, ...) {
return(sum(fnOrig(theta, ...)))
}
if (!is.null(attr(f, "gradient"))) {
gradFunc <- function(theta, ...) {
return(sumGradients(attr(fnOrig(theta, ...),
"gradient"), nParam))
}
}
else if (!is.null(gradOrig)) {
gradFunc <- function(theta, ...) {
return(sumGradients(gradOrig(theta, ...),
nParam))
}
}
else {
gradFunc <- NULL
}
h <- numericHessian(f = llFunc, grad = gradFunc,
t0 = theta, fixed = fixed, ...)
}
}
if ((dim(h)[1] != nParam) | (dim(h)[2] != nParam)) {
stop("Wrong hessian dimension. Needed ", nParam,
"x", nParam, " but supplied ", dim(h)[1], "x",
dim(h)[2])
}
else {
h[fixed, ] <- NA
h[, fixed] <- NA
}
}
else if (tolower(returnHessian) == "bhhh") {
h <- NULL
if (is.null(dim(gr)) & any(is.na(gr[!fixed]))) {
h <- NA
}
else if (is.matrix(gr)) {
if (any(is.na(gr[, !fixed]))) {
h <- NA
}
}
if (is.null(h)) {
checkBhhhGrad(g = gr, theta = theta, analytic = (!is.null(attr(f,
"gradient")) || !is.null(gradOrig)), fixed = fixed)
h <- -crossprod(gr)
}
attr(h, "type") = "BHHH"
}
else {
h <- NULL
}
if (sumObs) {
f <- sumKeepAttr(f)
}
if (sumObs) {
gr <- sumGradients(gr, nParam)
}
if (!is.null(gradOrig) && !is.null(attr(f, "gradient"))) {
attr(f, "gradBoth") <- TRUE
}
if (!is.null(hessOrig) && !is.null(attr(f, "hessian"))) {
attr(f, "hessBoth") <- TRUE
}
attr(f, "gradient") <- gr
attr(f, "hessian") <- h
return(f)
}, fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, start = c(b_1_0 = 0.598378200054935, b_1_1 = 6.336646267785,
b_1_2 = -1.01231147815159, b_2_0 = -1.03564793824529, b_2_1 = 1.26302210717712,
b_2_2 = -1.70638777681512, b_3_0 = 0.410893287564513, b_3_1 = -0.373151392833251,
b_3_2 = 2.92074017047352, R_1_2 = 0.172617843070694, R_1_3 = -0.0160393443285098,
R_2_3 = -0.0704265674853818), finalHessian = "BHHH", bhhhHessian = TRUE,
fixed = c(b_1_0 = FALSE, b_1_1 = FALSE, b_1_2 = FALSE, b_2_0 = FALSE,
b_2_1 = FALSE, b_2_2 = FALSE, b_3_0 = FALSE, b_3_1 = FALSE,
b_3_2 = FALSE, R_1_2 = FALSE, R_1_3 = FALSE, R_2_3 = FALSE
), control = new("MaxControl", tol = 1e-08, reltol = 1.49011611938477e-08,
gradtol = 1e-06, steptol = 1e-10, lambdatol = 1e-06,
qrtol = 1e-10, qac = "stephalving", marquardt_lambda0 = 0.01,
marquardt_lambdaStep = 2, marquardt_maxLambda = 1e+12,
nm_alpha = 1, nm_beta = 0.5, nm_gamma = 2, sann_cand = NULL,
sann_temp = 10, sann_tmax = 10L, sann_randomSeed = 123L,
iterlim = 1L, printLevel = 0L), yMat = c(0, 1, 1, 1,
1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0,
1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1,
1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0), xMat = c(1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
0.584613749636069, 0.123854243844614, 0.215941568743973,
0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258,
0.448209778629426, 0.0530042267305041, 0.922267467879738,
2.05008468562714, -0.491031166056535, -2.30916887564081,
1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633,
-0.370660031792409, 0.644376548518833, -0.220486561818751,
0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
-0.627906076039371, 1.36065244853001, -0.600259587147127,
2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
), llAlgorithm = "GHK", llNGHK = 50, llIntGrad = TRUE, llOneSidedGrad = FALSE,
llEps = 1e-06, llRandom.seed = 123)
where 6: eval(as.call(cl))
where 7: eval(as.call(cl))
where 8: maxNR(fn = fn, grad = grad, hess = hess, start = start, finalHessian = finalHessian,
bhhhHessian = TRUE, ...)
where 9: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
constraints = constraints, ...)
where 10: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
llRandom.seed = random.seed, ...)
where 11: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, data = as.data.frame(cbind(xMat,
yMat)), iterlim = 1, nGHK = 50)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x3655878>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbit ... fn -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 0.1-8
Check: tests
Result: ERROR
Running ‘mvProbitEst.R’ [34s/36s]
Running ‘mvProbitMargEff2.R’
Running ‘mvProbitTest.R’
Running ‘mvProbitTest5.R’ [145s/148s]
Comparing ‘mvProbitTest5.Rout’ to ‘mvProbitTest5.Rout.save’ ...1604c1604
< 5 0.00 -0.01 0 0.00 0.00 0.00 0.00
---
> 5 0.00 0.00 0 0.00 0.00 0.00 0.00
1618,1620c1618,1620
< 1 0 0.00 0 0 0 0 0
< 5 0 -0.01 0 0 0 0 0
< 10 0 0.00 0 0 0 0 0
---
> 1 0 0 0 0 0 0 0
> 5 0 0 0 0 0 0 0
> 10 0 0 0 0 0 0 0
1728c1728
< 5 -0.01
---
> 5 0.00
Running ‘pmvnormWrapTest.R’
Running the tests in ‘tests/mvProbitEst.R’ failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> options( digits = 4 )
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 50
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -0.8 ),
+ c( -0.6, 1.0, -1.6 ),
+ c( 0.5, -0.6, 1.2 ) )
>
> # covariance matrix of error terms
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
>
> # create data frame
> dat <- as.data.frame( cbind( xMat, yMat ) )
>
> # estimation with the BHHH algorithm, two-sided gradients
> estResultBHHH <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz() )
> print( estResultBHHH )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
b_3_2 R_1_2 R_1_3 R_2_3
1.23283 0.86322 0.36022 -0.08406
> summary( estResultBHHH )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 1.02508 0.44548 2.301 0.0214 *
b_1_1 1.30536 0.73511 1.776 0.0758 .
b_1_2 -1.02809 0.68851 -1.493 0.1354
b_2_0 -0.78248 0.38798 -2.017 0.0437 *
b_2_1 1.14318 0.53026 2.156 0.0311 *
b_2_2 -1.66938 0.81308 -2.053 0.0401 *
b_3_0 0.58229 0.64107 0.908 0.3637
b_3_1 -0.69667 0.89818 -0.776 0.4380
b_3_2 1.23283 0.62452 1.974 0.0484 *
R_1_2 0.86322 18.60071 0.046 0.9630
R_1_3 0.36022 0.64712 0.557 0.5778
R_2_3 -0.08406 0.34909 -0.241 0.8097
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BHHH maximisation, 1 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -56.91 on 12 Df
> logLik( estResultBHHH )
'log Lik.' -56.91 (df=12)
> estResultBHHHA <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta, sigma[ lower.tri( sigma ) ] ),
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz() )
> all.equal( estResultBHHH, estResultBHHHA )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BHHH algorithm, one-sided gradients
> estResultBHHH1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBHHH1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
tol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
b_3_2 R_1_2 R_1_3 R_2_3
1.23283 0.86322 0.36022 -0.08406
> summary( estResultBHHH1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
tol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 1.02508 0.44548 2.301 0.0214 *
b_1_1 1.30536 0.73511 1.776 0.0758 .
b_1_2 -1.02809 0.68851 -1.493 0.1354
b_2_0 -0.78248 0.38797 -2.017 0.0437 *
b_2_1 1.14318 0.53026 2.156 0.0311 *
b_2_2 -1.66938 0.81307 -2.053 0.0401 *
b_3_0 0.58229 0.64106 0.908 0.3637
b_3_1 -0.69667 0.89817 -0.776 0.4380
b_3_2 1.23283 0.62452 1.974 0.0484 *
R_1_2 0.86322 18.60041 0.046 0.9630
R_1_3 0.36022 0.64712 0.557 0.5778
R_2_3 -0.08406 0.34909 -0.241 0.8097
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BHHH maximisation, 1 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -56.91 on 12 Df
> logLik( estResultBHHH1 )
'log Lik.' -56.91 (df=12)
> all.equal( estResultBHHH, estResultBHHH1, tol = 1e-5 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, two-sided gradients
> estResultBFGS <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma, method = "BFGS",
+ data = dat,
+ reltol = 0.5, algorithm = GenzBretz() )
> print( estResultBFGS )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53607 1.16101 -1.14212 -0.44952 1.05496 -1.57624 0.40751 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25958 0.29483 0.32822 -0.02382
> summary( estResultBFGS )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53607 0.82709 0.648 0.5169
b_1_1 1.16101 0.91515 1.269 0.2046
b_1_2 -1.14212 0.93847 -1.217 0.2236
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88398 1.193 0.2327
b_2_2 -1.57624 0.73187 -2.154 0.0313 *
b_3_0 0.40751 0.48897 0.833 0.4046
b_3_1 -0.85400 0.84871 -1.006 0.3143
b_3_2 1.25958 0.68040 1.851 0.0641 .
R_1_2 0.29483 1.85817 0.159 0.8739
R_1_3 0.32822 0.70923 0.463 0.6435
R_2_3 -0.02382 0.34466 -0.069 0.9449
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGS )
'log Lik.' -52.94 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients
> estResultBFGS1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma, method = "BFGS",
+ data = dat,
+ reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53607 1.16101 -1.14213 -0.44952 1.05496 -1.57624 0.40751 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25958 0.29483 0.32822 -0.02382
> summary( estResultBFGS1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53607 0.82709 0.648 0.5169
b_1_1 1.16101 0.91515 1.269 0.2046
b_1_2 -1.14213 0.93847 -1.217 0.2236
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88398 1.193 0.2327
b_2_2 -1.57624 0.73187 -2.154 0.0313 *
b_3_0 0.40751 0.48897 0.833 0.4046
b_3_1 -0.85400 0.84871 -1.006 0.3143
b_3_2 1.25958 0.68040 1.851 0.0641 .
R_1_2 0.29483 1.85817 0.159 0.8739
R_1_3 0.32822 0.70923 0.463 0.6435
R_2_3 -0.02382 0.34466 -0.069 0.9449
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGS1 )
'log Lik.' -52.94 (df=12)
> all.equal( estResultBFGS, estResultBFGS1, tol = 1e-5 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values
> estResultBFGS1a <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ data = dat, method = "BFGS",
+ reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1a )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.79353 1.88367 -1.74922 -0.40036 0.91577 -1.52106 0.97790 -1.69433
b_3_2 R_1_2 R_1_3 R_2_3
1.57625 0.75270 0.23721 -0.03453
> summary( estResultBFGS1a )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.79353 0.90412 0.878 0.3801
b_1_1 1.88367 0.90668 2.078 0.0378 *
b_1_2 -1.74922 0.97400 -1.796 0.0725 .
b_2_0 -0.40036 0.69056 -0.580 0.5621
b_2_1 0.91577 0.97998 0.934 0.3501
b_2_2 -1.52106 0.72006 -2.112 0.0347 *
b_3_0 0.97790 0.64497 1.516 0.1295
b_3_1 -1.69433 0.94788 -1.787 0.0739 .
b_3_2 1.57625 0.73399 2.148 0.0318 *
R_1_2 0.75270 4.32729 0.174 0.8619
R_1_3 0.23721 0.90535 0.262 0.7933
R_2_3 -0.03453 0.35889 -0.096 0.9234
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 3 iterations
Return code 0: successful convergence
Log-likelihood: -49.36 on 12 Df
> logLik( estResultBFGS1a )
'log Lik.' -49.36 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values for beta
> estResultBFGS1b <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ startSigma = sigma, data = dat,
+ method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1b )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1 b_3_2 R_1_2
0.8110 1.9060 -1.7627 -0.3770 0.9441 -1.5332 0.9372 -1.7498 1.6024 0.6817
R_1_3 R_2_3
0.4044 0.1663
> summary( estResultBFGS1b )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.8110 0.9297 0.872 0.3830
b_1_1 1.9060 0.9160 2.081 0.0375 *
b_1_2 -1.7627 1.0131 -1.740 0.0819 .
b_2_0 -0.3770 0.7049 -0.535 0.5928
b_2_1 0.9441 1.0345 0.913 0.3615
b_2_2 -1.5332 0.7517 -2.040 0.0414 *
b_3_0 0.9372 0.6291 1.490 0.1363
b_3_1 -1.7498 0.9996 -1.751 0.0800 .
b_3_2 1.6024 0.7848 2.042 0.0412 *
R_1_2 0.6817 3.4094 0.200 0.8415
R_1_3 0.4044 0.8024 0.504 0.6143
R_2_3 0.1663 0.3389 0.491 0.6235
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 3 iterations
Return code 0: successful convergence
Log-likelihood: -49.57 on 12 Df
> logLik( estResultBFGS1b )
'log Lik.' -49.57 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values for sigma
> estResultBFGS1s <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), data = dat,
+ method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1s )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.55814 1.14958 -1.12747 -0.44426 1.07325 -1.57226 0.37255 -0.85470
b_3_2 R_1_2 R_1_3 R_2_3
1.23011 0.34981 -0.01177 0.09435
> summary( estResultBFGS1s )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.55814 0.86615 0.644 0.5193
b_1_1 1.14958 0.90607 1.269 0.2045
b_1_2 -1.12747 0.94191 -1.197 0.2313
b_2_0 -0.44426 0.85526 -0.519 0.6035
b_2_1 1.07325 0.99314 1.081 0.2798
b_2_2 -1.57226 0.78557 -2.001 0.0453 *
b_3_0 0.37255 0.56172 0.663 0.5072
b_3_1 -0.85470 0.95214 -0.898 0.3694
b_3_2 1.23011 0.71506 1.720 0.0854 .
R_1_2 0.34981 2.02289 0.173 0.8627
R_1_3 -0.01177 0.70776 -0.017 0.9867
R_2_3 0.09435 0.33920 0.278 0.7809
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.78 on 12 Df
> logLik( estResultBFGS1s )
'log Lik.' -52.78 (df=12)
>
> # estimation with the BFGS algorithm, Miwa algorithm for obtaining integrals
> estResultBFGSm <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, method = "BFGS",
+ reltol = 0.5, algorithm = Miwa( steps = 64 ) )
There were 13 warnings (use warnings() to see them)
> print( estResultBFGSm )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53609 1.16102 -1.14213 -0.44952 1.05496 -1.57624 0.40753 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25960 0.29491 0.32819 -0.02386
> summary( estResultBFGSm )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53609 0.82702 0.648 0.5168
b_1_1 1.16102 0.91511 1.269 0.2045
b_1_2 -1.14213 0.93837 -1.217 0.2235
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88386 1.194 0.2326
b_2_2 -1.57624 0.73173 -2.154 0.0312 *
b_3_0 0.40753 0.48894 0.834 0.4046
b_3_1 -0.85400 0.84860 -1.006 0.3142
b_3_2 1.25960 0.68024 1.852 0.0641 .
R_1_2 0.29491 1.85789 0.159 0.8739
R_1_3 0.32819 0.70903 0.463 0.6435
R_2_3 -0.02386 0.34468 -0.069 0.9448
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGSm )
'log Lik.' -52.94 (df=12)
> all.equal( estResultBFGS, estResultBFGSm, tol = 1e-3 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, GHK algorithm for obtaining integrals
> estResultBFGSg <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, method = "BFGS",
+ reltol = 0.5 )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
...)
where 3: fnOrig(theta, ...)
where 4: logLikFunc(theta, fnOrig = function (param, yMat, xMat, llAlgorithm,
llNGHK, llIntGrad, llOneSidedGrad, llEps, llRandom.seed,
...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, gradOrig = ..2, hessOrig = ..3, yMat = c(1, 1, 1, 1, 1, 1,
1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0,
1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0,
0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0,
0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1,
0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,
0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1,
0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1,
1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0.253318513994755,
-0.028546755348703, -0.0428704572913161, 1.36860228401446, -0.225770985659268,
1.51647060442954, -1.54875280423022, 0.584613749636069, 0.123854243844614,
0.215941568743973, 0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258, 0.448209778629426,
0.0530042267305041, 0.922267467879738, 2.05008468562714, -0.491031166056535,
-2.30916887564081, 1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633, -0.370660031792409,
0.644376548518833, -0.220486561818751, 0.331781963915697, 1.09683901314935,
0.435181490833803, -0.325931585531227, 1.14880761845109, 0.993503855962119,
0.54839695950807, 0.238731735111441, -0.627906076039371, 1.36065244853001,
-0.600259587147127, 2.18733299301658, 1.53261062618519, -0.235700359100477,
-1.02642090030678), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
where 5: eval(f, sys.frame(sys.parent()))
where 6: eval(f, sys.frame(sys.parent()))
where 7: callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
...)
where 8: (function (theta, fName, ...)
{
return(callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
...))
})(theta = c(b_1_0 = 0.8, b_1_1 = 1.2, b_1_2 = -0.8, b_2_0 = -0.6,
b_2_1 = 1, b_2_2 = -1.6, b_3_0 = 0.5, b_3_1 = -0.6, b_3_2 = 1.2,
R_1_2 = 0.2, R_1_3 = 0.4, R_2_3 = -0.1), fName = "logLikFunc",
fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, gradOrig = NULL, hessOrig = NULL, yMat = c(1, 1, 1, 1,
1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0,
1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0,
0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0,
0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1,
1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
0.584613749636069, 0.123854243844614, 0.215941568743973,
0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258,
0.448209778629426, 0.0530042267305041, 0.922267467879738,
2.05008468562714, -0.491031166056535, -2.30916887564081,
1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633,
-0.370660031792409, 0.644376548518833, -0.220486561818751,
0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
-0.627906076039371, 1.36065244853001, -0.600259587147127,
2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
where 9: do.call(callWithoutSumt, argList)
where 10: maxOptim(fn = fn, grad = grad, hess = hess, start = start, method = "BFGS",
fixed = fixed, constraints = constraints, finalHessian = finalHessian,
parscale = parscale, control = mControl, ...)
where 11: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
constraints = constraints, ...)
where 12: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
llRandom.seed = random.seed, ...)
where 13: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, start = c(beta), startSigma = sigma,
data = dat, method = "BFGS", reltol = 0.5)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x1074890>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbit ... logLikFunc -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
Execution halted
Running the tests in ‘tests/mvProbitMargEff2.R’ failed.
Complete output:
> # I thank Mohit Batham for providing this R script that demonstrated
> # a bug in mvProbitMargEff() when called with 2 dependent variables
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> nObs <- 100
> set.seed( 123 )
> xData <- data.frame(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = as.numeric( rnorm( nObs ) > 0 ),
+ x3 = rnorm( nObs ),
+ x4 = rnorm( nObs ))
> beta <- c( 0.8, 1.2, -1.0, 1.4, -0.8,
+ -0.6, 1.0, 0.6, -1.2, -1.6)
> sigma <- symMatrix( c( 1, 0.2, 1))
> margEffCond <- try( mvProbitMargEff( ~ x1 + x2 + x3 + x4 , coef = beta,
+ sigma = sigma, data = xData, cond = TRUE ) )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = betaLower,
sigma = coef$sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
random.seed = random.seed, ...)
where 3: mvProbitMargEffInternal(yMat = yMat, xMat = xMat, coef = coef,
sigma = sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
eps = eps, dummyVars = dummyVars, random.seed = random.seed,
...)
where 4: mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
data = xData, cond = TRUE)
where 5: doTryCatch(return(expr), name, parentenv, handler)
where 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 7: tryCatchList(expr, classes, parentenv, handlers)
where 8: tryCatch(expr, error = function(e) {
call <- conditionCall(e)
if (!is.null(call)) {
if (identical(call[[1L]], quote(doTryCatch)))
call <- sys.call(-4L)
dcall <- deparse(call)[1L]
prefix <- paste("Error in", dcall, ": ")
LONG <- 75L
sm <- strsplit(conditionMessage(e), "\n")[[1L]]
w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w")
if (is.na(w))
w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L],
type = "b")
if (w > LONG)
prefix <- paste0(prefix, "\n ")
}
else prefix <- "Error : "
msg <- paste0(prefix, conditionMessage(e), "\n")
.Internal(seterrmessage(msg[1L]))
if (!silent && isTRUE(getOption("show.error.messages"))) {
cat(msg, file = outFile)
.Internal(printDeferredWarnings())
}
invisible(structure(msg, class = "try-error", condition = e))
})
where 9: try(mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
data = xData, cond = TRUE))
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x29ccd30>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
> round( margEffCond, 3 )
Error in round(margEffCond, 3) :
non-numeric argument to mathematical function
Execution halted
Running the tests in ‘tests/mvProbitTest.R’ failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> options( digits = 4 )
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 10
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = as.numeric( rnorm( nObs ) > 0 ),
+ x3 = rnorm( nObs ),
+ x4 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -1.0, 1.4, -0.8 ),
+ c( -0.6, 1.0, 0.6, -1.2, -1.6 ),
+ c( 0.5, -0.6, -0.7, 1.1, 1.2 ) )
>
> # covariance matrix of error terms
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # all parameters in a vector
> allCoef <- c( c( beta ), sigma[ lower.tri( sigma ) ] )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
> # (yMatLin > 0 )== yMat
>
> # unconditional expectations of dependent variables
> yExp <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ) )
> round( yExp, 3 )
V1 V2 V3
1 0.021 0.725 0.194
2 0.394 0.768 0.214
3 0.125 0.788 0.196
4 0.235 0.681 0.292
5 0.680 0.435 0.579
6 0.028 0.973 0.034
7 0.958 0.186 0.784
8 0.856 0.247 0.724
9 0.061 0.968 0.034
10 0.998 0.067 0.923
> yExpA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
+ data = as.data.frame( xMat ) )
> all.equal( yExp, yExpA )
[1] TRUE
> yExp2 <- pnorm( yMatLin )
> all.equal( yExp, as.data.frame( yExp2 ) )
[1] TRUE
>
>
> # conditional expectations of dependent variables
> # (assuming that all other dependent variables are one)
> yExpCond <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz() )
> round( yExpCond, 3 )
V1 V2 V3
1 0.072 0.834 0.524
2 0.660 0.779 0.314
3 0.297 0.838 0.399
4 0.448 0.727 0.461
5 0.847 0.442 0.618
6 0.139 0.985 0.162
7 0.991 0.179 0.749
8 0.950 0.244 0.710
9 0.244 0.978 0.133
10 1.000 0.065 0.892
> yExpCondA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
+ data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz() )
> all.equal( yExpCond, yExpCondA )
[1] TRUE
> yExpCond2 <- matrix( NA, nrow = nObs, ncol = ncol( yMat ) )
> for( i in 1:nObs ) {
+ for( k in 1:ncol( yMat ) ) {
+ set.seed( 123 )
+ numerator <- pmvnorm( upper = yMatLin[ i, ], sigma = sigma )
+ set.seed( 123 )
+ denominator <- pmvnorm( upper = yMatLin[ i, -k ], sigma = sigma[ -k, -k ] )
+ yExpCond2[ i, k ] <- numerator / denominator
+ }
+ }
> all.equal( yExpCond, as.data.frame( yExpCond2 ) )
[1] TRUE
> # now with explicitly specifying the algorithm
> yExpCond3 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz )
> all.equal( yExpCond, yExpCond3 )
[1] TRUE
> identical( yExpCond, yExpCond3 )
[1] TRUE
> # now with integrals obtained by the Miwa algorithm
> yExpCond4 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = Miwa )
> all.equal( yExpCond, yExpCond4, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the Miwa algorithm, less precise
> yExpCond5 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = Miwa( steps = 32 ) )
> all.equal( yExpCond4, yExpCond5, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the TVPACK algorithm
> yExpCond6 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = TVPACK )
> all.equal( yExpCond, yExpCond6, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the TVPACK algorithm, less precise
> yExpCond7 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = TVPACK( abseps = 0.5 ) )
> all.equal( yExpCond6, yExpCond7, tol = 1e-3 )
[1] "Component \"V1\": Mean relative difference: 0.03018"
[2] "Component \"V2\": Mean relative difference: 0.04273"
[3] "Component \"V3\": Mean relative difference: 0.04242"
> # now with integrals obtained by the GHK algorithm
> yExpCond8 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = coef, sigma = sigma,
cond = cond, algorithm = algorithm, nGHK = nGHK, random.seed = random.seed,
...)
where 3: mvProbitExp(~x1 + x2 + x3 + x4, coef = c(beta), sigma = sigma,
data = as.data.frame(xMat), cond = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x1b424a0>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbitExp -> mvProbitExpInternal -> pmvnormWrap
Execution halted
Running the tests in ‘tests/pmvnormWrapTest.R’ failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> library( "miscTools" )
>
> # covariance matrix
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> ######## only upper ##########
> upper <- c( -0.3, 0.7, -0.5 )
> # Genz + Bretz (default)
> pug <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = GenzBretz(), random.seed = 123 )
> print( pug )
[1] 0.1360783
attr(,"error")
[1] 0.000118605
attr(,"msg")
[1] "Normal Completion"
>
> # Miwa (as function)
> pum <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = Miwa, random.seed = 123 )
> print( pum )
[1] 0.136069
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> all.equal( pug, pum, check.attributes = FALSE, tol = 1e-4 )
[1] TRUE
>
> # Miwa (as object returned from function Miwa())
> pum1 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = Miwa(), random.seed = 123 )
> all.equal( pum, pum1 )
[1] TRUE
>
> # Miwa (as character string)
> pum2 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = "Miwa", random.seed = 123 )
> all.equal( pum, pum2 )
[1] TRUE
>
> # TVPACK
> put <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = TVPACK, random.seed = 123 )
> print( put )
[1] 0.136069
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> all.equal( pug, put, check.attributes = FALSE, tol = 1e-4 )
[1] TRUE
> all.equal( pum, put, check.attributes = FALSE, tol = 1e-6 )
[1] TRUE
>
> # GHK
> pughk <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = "ghk", random.seed = 123, nGHK = 1000 )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
random.seed = 123, nGHK = 1000)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
random.seed = 123, nGHK = 1000)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x3dcc488>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: <Anonymous>
Execution halted
Flavor: r-devel-linux-x86_64-fedora-clang
Version: 0.1-8
Check: examples
Result: ERROR
Running examples in ‘mvProbit-Ex.R’ failed
The error most likely occurred in:
> ### Name: mvProbit
> ### Title: Estimation of Multivariate Probit Models
> ### Aliases: mvProbit print.mvProbit
> ### Keywords: models regression
>
> ### ** Examples
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 50
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -0.8 ),
+ c( -0.6, 1.0, -1.6 ),
+ c( 0.5, -0.6, 1.2 ) )
>
> # covariance matrix of error terms
> library( miscTools )
> sigma <- symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
>
> # estimation (BHHH optimizer and GHK algorithm)
> estResult <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ data = as.data.frame( cbind( xMat, yMat ) ), iterlim = 1, nGHK = 50 )
Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
...)
where 3: fnOrig(theta, ...)
where 4: fn(start1, fixed = fixed, sumObs = TRUE, returnHessian = returnHessian,
...)
where 5: maxNRCompute(fn = function (theta, fnOrig, gradOrig = NULL, hessOrig = NULL,
fixed, sumObs = FALSE, returnHessian = TRUE, ...)
{
nParam <- length(theta)
f <- fnOrig(theta, ...)
if (any(is.na(f))) {
attr(f, "gradient") <- NA
attr(f, "hessian") <- NA
return(f)
}
gr <- attr(f, "gradient")
if (is.null(gr)) {
if (!is.null(gradOrig)) {
gr <- gradOrig(theta, ...)
}
else {
gr <- numericGradient(f = fnOrig, t0 = theta, fixed = fixed,
...)
}
}
if (is.matrix(gr)) {
if (ncol(gr) != length(theta)) {
stop(paste0("if gradient is a matrix, it must have length(parameter) colums (currently ",
length(theta), "), not ", ncol(gr)))
}
activeGr <- gr[, !fixed]
}
else {
activeGr <- gr[!fixed]
}
if (any(is.na(activeGr))) {
attr(f, "gradient") <- gr
attr(f, "hessian") <- NA
return(f)
}
if (observationGradient(gr, length(theta))) {
gr <- as.matrix(gr)
}
if (is.null(dim(gr))) {
gr[fixed] <- NA
}
else {
gr[, fixed] <- NA
}
if (isTRUE(returnHessian)) {
h <- attr(f, "hessian")
if (is.null(h)) {
if (!is.null(hessOrig)) {
h <- as.matrix(hessOrig(theta, ...))
}
else {
llFunc <- function(theta, ...) {
return(sum(fnOrig(theta, ...)))
}
if (!is.null(attr(f, "gradient"))) {
gradFunc <- function(theta, ...) {
return(sumGradients(attr(fnOrig(theta, ...),
"gradient"), nParam))
}
}
else if (!is.null(gradOrig)) {
gradFunc <- function(theta, ...) {
return(sumGradients(gradOrig(theta, ...),
nParam))
}
}
else {
gradFunc <- NULL
}
h <- numericHessian(f = llFunc, grad = gradFunc,
t0 = theta, fixed = fixed, ...)
}
}
if ((dim(h)[1] != nParam) | (dim(h)[2] != nParam)) {
stop("Wrong hessian dimension. Needed ", nParam,
"x", nParam, " but supplied ", dim(h)[1], "x",
dim(h)[2])
}
else {
h[fixed, ] <- NA
h[, fixed] <- NA
}
}
else if (tolower(returnHessian) == "bhhh") {
h <- NULL
if (is.null(dim(gr)) & any(is.na(gr[!fixed]))) {
h <- NA
}
else if (is.matrix(gr)) {
if (any(is.na(gr[, !fixed]))) {
h <- NA
}
}
if (is.null(h)) {
checkBhhhGrad(g = gr, theta = theta, analytic = (!is.null(attr(f,
"gradient")) || !is.null(gradOrig)), fixed = fixed)
h <- -crossprod(gr)
}
attr(h, "type") = "BHHH"
}
else {
h <- NULL
}
if (sumObs) {
f <- sumKeepAttr(f)
}
if (sumObs) {
gr <- sumGradients(gr, nParam)
}
if (!is.null(gradOrig) && !is.null(attr(f, "gradient"))) {
attr(f, "gradBoth") <- TRUE
}
if (!is.null(hessOrig) && !is.null(attr(f, "hessian"))) {
attr(f, "hessBoth") <- TRUE
}
attr(f, "gradient") <- gr
attr(f, "hessian") <- h
return(f)
}, fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, start = c(b_1_0 = 0.598378200054935, b_1_1 = 6.336646267785,
b_1_2 = -1.01231147815159, b_2_0 = -1.03564793824529, b_2_1 = 1.26302210717712,
b_2_2 = -1.70638777681512, b_3_0 = 0.410893287564513, b_3_1 = -0.373151392833251,
b_3_2 = 2.92074017047352, R_1_2 = 0.172617843070694, R_1_3 = -0.0160393443285098,
R_2_3 = -0.0704265674853818), finalHessian = "BHHH", bhhhHessian = TRUE,
fixed = c(b_1_0 = FALSE, b_1_1 = FALSE, b_1_2 = FALSE, b_2_0 = FALSE,
b_2_1 = FALSE, b_2_2 = FALSE, b_3_0 = FALSE, b_3_1 = FALSE,
b_3_2 = FALSE, R_1_2 = FALSE, R_1_3 = FALSE, R_2_3 = FALSE
), control = new("MaxControl", tol = 1e-08, reltol = 1.49011611938477e-08,
gradtol = 1e-06, steptol = 1e-10, lambdatol = 1e-06,
qrtol = 1e-10, qac = "stephalving", marquardt_lambda0 = 0.01,
marquardt_lambdaStep = 2, marquardt_maxLambda = 1e+12,
nm_alpha = 1, nm_beta = 0.5, nm_gamma = 2, sann_cand = NULL,
sann_temp = 10, sann_tmax = 10L, sann_randomSeed = 123L,
iterlim = 1L, printLevel = 0L), yMat = c(0, 1, 1, 1,
1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1,
1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1,
1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1,
0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0,
1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1,
1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0), xMat = c(1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
0.584613749636069, 0.123854243844614, 0.215941568743973,
0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258,
0.448209778629426, 0.0530042267305041, 0.922267467879738,
2.05008468562714, -0.491031166056535, -2.30916887564081,
1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633,
-0.370660031792409, 0.644376548518833, -0.220486561818751,
0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
-0.627906076039371, 1.36065244853001, -0.600259587147127,
2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
), llAlgorithm = "GHK", llNGHK = 50, llIntGrad = TRUE, llOneSidedGrad = FALSE,
llEps = 1e-06, llRandom.seed = 123)
where 6: eval(as.call(cl))
where 7: eval(as.call(cl))
where 8: maxNR(fn = fn, grad = grad, hess = hess, start = start, finalHessian = finalHessian,
bhhhHessian = TRUE, ...)
where 9: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
constraints = constraints, ...)
where 10: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
llRandom.seed = random.seed, ...)
where 11: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, data = as.data.frame(cbind(xMat,
yMat)), iterlim = 1, nGHK = 50)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x36e8a08>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbit ... fn -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc
Version: 0.1-8
Check: tests
Result: ERROR
Running ‘mvProbitEst.R’ [33s/50s]
Running ‘mvProbitMargEff2.R’
Running ‘mvProbitTest.R’
Running ‘mvProbitTest5.R’ [144s/204s]
Comparing ‘mvProbitTest5.Rout’ to ‘mvProbitTest5.Rout.save’ ...1604c1604
< 5 0.00 -0.01 0 0.00 0.00 0.00 0.00
---
> 5 0.00 0.00 0 0.00 0.00 0.00 0.00
1618,1620c1618,1620
< 1 0 0.00 0 0 0 0 0
< 5 0 -0.01 0 0 0 0 0
< 10 0 0.00 0 0 0 0 0
---
> 1 0 0 0 0 0 0 0
> 5 0 0 0 0 0 0 0
> 10 0 0 0 0 0 0 0
1728c1728
< 5 -0.01
---
> 5 0.00
Running ‘pmvnormWrapTest.R’
Running the tests in ‘tests/mvProbitEst.R’ failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> options( digits = 4 )
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 50
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -0.8 ),
+ c( -0.6, 1.0, -1.6 ),
+ c( 0.5, -0.6, 1.2 ) )
>
> # covariance matrix of error terms
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
>
> # create data frame
> dat <- as.data.frame( cbind( xMat, yMat ) )
>
> # estimation with the BHHH algorithm, two-sided gradients
> estResultBHHH <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz() )
> print( estResultBHHH )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
b_3_2 R_1_2 R_1_3 R_2_3
1.23283 0.86322 0.36022 -0.08406
> summary( estResultBHHH )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), tol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 1.02508 0.44548 2.301 0.0214 *
b_1_1 1.30536 0.73511 1.776 0.0758 .
b_1_2 -1.02809 0.68851 -1.493 0.1354
b_2_0 -0.78248 0.38798 -2.017 0.0437 *
b_2_1 1.14318 0.53026 2.156 0.0311 *
b_2_2 -1.66938 0.81308 -2.053 0.0401 *
b_3_0 0.58229 0.64107 0.908 0.3637
b_3_1 -0.69667 0.89818 -0.776 0.4380
b_3_2 1.23283 0.62452 1.974 0.0484 *
R_1_2 0.86322 18.60071 0.046 0.9630
R_1_3 0.36022 0.64712 0.557 0.5778
R_2_3 -0.08406 0.34909 -0.241 0.8097
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BHHH maximisation, 1 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -56.91 on 12 Df
> logLik( estResultBHHH )
'log Lik.' -56.91 (df=12)
> estResultBHHHA <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta, sigma[ lower.tri( sigma ) ] ),
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz() )
> all.equal( estResultBHHH, estResultBHHHA )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BHHH algorithm, one-sided gradients
> estResultBHHH1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, tol = 0.5,
+ algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBHHH1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
tol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
1.02508 1.30536 -1.02809 -0.78248 1.14318 -1.66938 0.58229 -0.69667
b_3_2 R_1_2 R_1_3 R_2_3
1.23283 0.86322 0.36022 -0.08406
> summary( estResultBHHH1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, algorithm = GenzBretz(), oneSidedGrad = TRUE,
tol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 1.02508 0.44548 2.301 0.0214 *
b_1_1 1.30536 0.73511 1.776 0.0758 .
b_1_2 -1.02809 0.68851 -1.493 0.1354
b_2_0 -0.78248 0.38797 -2.017 0.0437 *
b_2_1 1.14318 0.53026 2.156 0.0311 *
b_2_2 -1.66938 0.81307 -2.053 0.0401 *
b_3_0 0.58229 0.64106 0.908 0.3637
b_3_1 -0.69667 0.89817 -0.776 0.4380
b_3_2 1.23283 0.62452 1.974 0.0484 *
R_1_2 0.86322 18.60041 0.046 0.9630
R_1_3 0.36022 0.64712 0.557 0.5778
R_2_3 -0.08406 0.34909 -0.241 0.8097
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BHHH maximisation, 1 iterations
Return code 2: successive function values within tolerance limit
Log-likelihood: -56.91 on 12 Df
> logLik( estResultBHHH1 )
'log Lik.' -56.91 (df=12)
> all.equal( estResultBHHH, estResultBHHH1, tol = 1e-5 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, two-sided gradients
> estResultBFGS <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma, method = "BFGS",
+ data = dat,
+ reltol = 0.5, algorithm = GenzBretz() )
> print( estResultBFGS )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53607 1.16101 -1.14212 -0.44952 1.05496 -1.57624 0.40751 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25958 0.29483 0.32822 -0.02382
> summary( estResultBFGS )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53607 0.82709 0.648 0.5169
b_1_1 1.16101 0.91515 1.269 0.2046
b_1_2 -1.14212 0.93847 -1.217 0.2236
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88398 1.193 0.2327
b_2_2 -1.57624 0.73187 -2.154 0.0313 *
b_3_0 0.40751 0.48897 0.833 0.4046
b_3_1 -0.85400 0.84871 -1.006 0.3143
b_3_2 1.25958 0.68040 1.851 0.0641 .
R_1_2 0.29483 1.85817 0.159 0.8739
R_1_3 0.32822 0.70923 0.463 0.6435
R_2_3 -0.02382 0.34466 -0.069 0.9449
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGS )
'log Lik.' -52.94 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients
> estResultBFGS1 <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma, method = "BFGS",
+ data = dat,
+ reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53607 1.16101 -1.14213 -0.44952 1.05496 -1.57624 0.40751 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25958 0.29483 0.32822 -0.02382
> summary( estResultBFGS1 )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = GenzBretz(),
oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53607 0.82709 0.648 0.5169
b_1_1 1.16101 0.91515 1.269 0.2046
b_1_2 -1.14213 0.93847 -1.217 0.2236
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88398 1.193 0.2327
b_2_2 -1.57624 0.73187 -2.154 0.0313 *
b_3_0 0.40751 0.48897 0.833 0.4046
b_3_1 -0.85400 0.84871 -1.006 0.3143
b_3_2 1.25958 0.68040 1.851 0.0641 .
R_1_2 0.29483 1.85817 0.159 0.8739
R_1_3 0.32822 0.70923 0.463 0.6435
R_2_3 -0.02382 0.34466 -0.069 0.9449
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGS1 )
'log Lik.' -52.94 (df=12)
> all.equal( estResultBFGS, estResultBFGS1, tol = 1e-5 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values
> estResultBFGS1a <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ data = dat, method = "BFGS",
+ reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1a )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.79353 1.88367 -1.74922 -0.40036 0.91577 -1.52106 0.97790 -1.69433
b_3_2 R_1_2 R_1_3 R_2_3
1.57625 0.75270 0.23721 -0.03453
> summary( estResultBFGS1a )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, method = "BFGS",
algorithm = GenzBretz(), oneSidedGrad = TRUE, reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.79353 0.90412 0.878 0.3801
b_1_1 1.88367 0.90668 2.078 0.0378 *
b_1_2 -1.74922 0.97400 -1.796 0.0725 .
b_2_0 -0.40036 0.69056 -0.580 0.5621
b_2_1 0.91577 0.97998 0.934 0.3501
b_2_2 -1.52106 0.72006 -2.112 0.0347 *
b_3_0 0.97790 0.64497 1.516 0.1295
b_3_1 -1.69433 0.94788 -1.787 0.0739 .
b_3_2 1.57625 0.73399 2.148 0.0318 *
R_1_2 0.75270 4.32729 0.174 0.8619
R_1_3 0.23721 0.90535 0.262 0.7933
R_2_3 -0.03453 0.35889 -0.096 0.9234
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 3 iterations
Return code 0: successful convergence
Log-likelihood: -49.36 on 12 Df
> logLik( estResultBFGS1a )
'log Lik.' -49.36 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values for beta
> estResultBFGS1b <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ startSigma = sigma, data = dat,
+ method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1b )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1 b_3_2 R_1_2
0.8110 1.9060 -1.7627 -0.3770 0.9441 -1.5332 0.9372 -1.7498 1.6024 0.6817
R_1_3 R_2_3
0.4044 0.1663
> summary( estResultBFGS1b )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, startSigma = sigma,
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.8110 0.9297 0.872 0.3830
b_1_1 1.9060 0.9160 2.081 0.0375 *
b_1_2 -1.7627 1.0131 -1.740 0.0819 .
b_2_0 -0.3770 0.7049 -0.535 0.5928
b_2_1 0.9441 1.0345 0.913 0.3615
b_2_2 -1.5332 0.7517 -2.040 0.0414 *
b_3_0 0.9372 0.6291 1.490 0.1363
b_3_1 -1.7498 0.9996 -1.751 0.0800 .
b_3_2 1.6024 0.7848 2.042 0.0412 *
R_1_2 0.6817 3.4094 0.200 0.8415
R_1_3 0.4044 0.8024 0.504 0.6143
R_2_3 0.1663 0.3389 0.491 0.6235
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 3 iterations
Return code 0: successful convergence
Log-likelihood: -49.57 on 12 Df
> logLik( estResultBFGS1b )
'log Lik.' -49.57 (df=12)
>
> # estimation with the BFGS algorithm, one-sided gradients, no starting values for sigma
> estResultBFGS1s <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), data = dat,
+ method = "BFGS", reltol = 0.5, algorithm = GenzBretz(), oneSidedGrad = TRUE )
> print( estResultBFGS1s )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.55814 1.14958 -1.12747 -0.44426 1.07325 -1.57226 0.37255 -0.85470
b_3_2 R_1_2 R_1_3 R_2_3
1.23011 0.34981 -0.01177 0.09435
> summary( estResultBFGS1s )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
method = "BFGS", algorithm = GenzBretz(), oneSidedGrad = TRUE,
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.55814 0.86615 0.644 0.5193
b_1_1 1.14958 0.90607 1.269 0.2045
b_1_2 -1.12747 0.94191 -1.197 0.2313
b_2_0 -0.44426 0.85526 -0.519 0.6035
b_2_1 1.07325 0.99314 1.081 0.2798
b_2_2 -1.57226 0.78557 -2.001 0.0453 *
b_3_0 0.37255 0.56172 0.663 0.5072
b_3_1 -0.85470 0.95214 -0.898 0.3694
b_3_2 1.23011 0.71506 1.720 0.0854 .
R_1_2 0.34981 2.02289 0.173 0.8627
R_1_3 -0.01177 0.70776 -0.017 0.9867
R_2_3 0.09435 0.33920 0.278 0.7809
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.78 on 12 Df
> logLik( estResultBFGS1s )
'log Lik.' -52.78 (df=12)
>
> # estimation with the BFGS algorithm, Miwa algorithm for obtaining integrals
> estResultBFGSm <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, method = "BFGS",
+ reltol = 0.5, algorithm = Miwa( steps = 64 ) )
There were 13 warnings (use warnings() to see them)
> print( estResultBFGSm )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
reltol = 0.5)
Coefficients:
b_1_0 b_1_1 b_1_2 b_2_0 b_2_1 b_2_2 b_3_0 b_3_1
0.53609 1.16102 -1.14213 -0.44952 1.05496 -1.57624 0.40753 -0.85400
b_3_2 R_1_2 R_1_3 R_2_3
1.25960 0.29491 0.32819 -0.02386
> summary( estResultBFGSm )
Call:
mvProbit(formula = cbind(y1, y2, y3) ~ x1 + x2, data = dat, start = c(beta),
startSigma = sigma, method = "BFGS", algorithm = Miwa(steps = 64),
reltol = 0.5)
Coefficients:
Estimate Std. error t value Pr(> t)
b_1_0 0.53609 0.82702 0.648 0.5168
b_1_1 1.16102 0.91511 1.269 0.2045
b_1_2 -1.14213 0.93837 -1.217 0.2235
b_2_0 -0.44952 0.79215 -0.567 0.5704
b_2_1 1.05496 0.88386 1.194 0.2326
b_2_2 -1.57624 0.73173 -2.154 0.0312 *
b_3_0 0.40753 0.48894 0.834 0.4046
b_3_1 -0.85400 0.84860 -1.006 0.3142
b_3_2 1.25960 0.68024 1.852 0.0641 .
R_1_2 0.29491 1.85789 0.159 0.8739
R_1_3 0.32819 0.70903 0.463 0.6435
R_2_3 -0.02386 0.34468 -0.069 0.9448
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
BFGS maximization, 4 iterations
Return code 0: successful convergence
Log-likelihood: -52.94 on 12 Df
> logLik( estResultBFGSm )
'log Lik.' -52.94 (df=12)
> all.equal( estResultBFGS, estResultBFGSm, tol = 1e-3 )
[1] "Component \"call\": target, current do not match when deparsed"
>
> # estimation with the BFGS algorithm, GHK algorithm for obtaining integrals
> estResultBFGSg <- mvProbit( cbind( y1, y2, y3 ) ~ x1 + x2,
+ start = c( beta ), startSigma = sigma,
+ data = dat, method = "BFGS",
+ reltol = 0.5 )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBetaTmp, sigma = sigmaTmp, algorithm = algorithm,
nGHK = nGHK, random.seed = randomSeed, ...)
where 2: mvProbitLogLikInternal(yMat = yMat, xMat = xMat, coef = param,
sigma = NULL, algorithm = llAlgorithm, nGHK = llNGHK, returnGrad = llIntGrad,
oneSidedGrad = llOneSidedGrad, eps = llEps, randomSeed = llRandom.seed,
...)
where 3: fnOrig(theta, ...)
where 4: logLikFunc(theta, fnOrig = function (param, yMat, xMat, llAlgorithm,
llNGHK, llIntGrad, llOneSidedGrad, llEps, llRandom.seed,
...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, gradOrig = ..2, hessOrig = ..3, yMat = c(1, 1, 1, 1, 1, 1,
1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1,
1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0,
1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0,
0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0,
0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1,
0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1,
0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1,
0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1,
1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0.253318513994755,
-0.028546755348703, -0.0428704572913161, 1.36860228401446, -0.225770985659268,
1.51647060442954, -1.54875280423022, 0.584613749636069, 0.123854243844614,
0.215941568743973, 0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258, 0.448209778629426,
0.0530042267305041, 0.922267467879738, 2.05008468562714, -0.491031166056535,
-2.30916887564081, 1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633, -0.370660031792409,
0.644376548518833, -0.220486561818751, 0.331781963915697, 1.09683901314935,
0.435181490833803, -0.325931585531227, 1.14880761845109, 0.993503855962119,
0.54839695950807, 0.238731735111441, -0.627906076039371, 1.36065244853001,
-0.600259587147127, 2.18733299301658, 1.53261062618519, -0.235700359100477,
-1.02642090030678), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
where 5: eval(f, sys.frame(sys.parent()))
where 6: eval(f, sys.frame(sys.parent()))
where 7: callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
...)
where 8: (function (theta, fName, ...)
{
return(callWithoutArgs(theta, fName = fName, args = names(formals(sumt)),
...))
})(theta = c(b_1_0 = 0.8, b_1_1 = 1.2, b_1_2 = -0.8, b_2_0 = -0.6,
b_2_1 = 1, b_2_2 = -1.6, b_3_0 = 0.5, b_3_1 = -0.6, b_3_2 = 1.2,
R_1_2 = 0.2, R_1_3 = 0.4, R_2_3 = -0.1), fName = "logLikFunc",
fnOrig = function (param, yMat, xMat, llAlgorithm, llNGHK,
llIntGrad, llOneSidedGrad, llEps, llRandom.seed, ...)
{
logLikVal <- mvProbitLogLikInternal(yMat = yMat, xMat = xMat,
coef = param, sigma = NULL, algorithm = llAlgorithm,
nGHK = llNGHK, returnGrad = llIntGrad, oneSidedGrad = llOneSidedGrad,
eps = llEps, randomSeed = llRandom.seed, ...)
return(logLikVal)
}, gradOrig = NULL, hessOrig = NULL, yMat = c(1, 1, 1, 1,
1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 0,
1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0,
0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1,
1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1,
0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0,
0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1,
1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0), xMat = c(1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0,
1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0,
1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0,
1, 0, 0.253318513994755, -0.028546755348703, -0.0428704572913161,
1.36860228401446, -0.225770985659268, 1.51647060442954, -1.54875280423022,
0.584613749636069, 0.123854243844614, 0.215941568743973,
0.379639482759882, -0.502323453109302, -0.33320738366942,
-1.01857538310709, -1.07179122647558, 0.303528641404258,
0.448209778629426, 0.0530042267305041, 0.922267467879738,
2.05008468562714, -0.491031166056535, -2.30916887564081,
1.00573852446226, -0.709200762582393, -0.688008616467358,
1.0255713696967, -0.284773007051009, -1.22071771225454, 0.18130347974915,
-0.138891362439045, 0.00576418589988693, 0.38528040112633,
-0.370660031792409, 0.644376548518833, -0.220486561818751,
0.331781963915697, 1.09683901314935, 0.435181490833803, -0.325931585531227,
1.14880761845109, 0.993503855962119, 0.54839695950807, 0.238731735111441,
-0.627906076039371, 1.36065244853001, -0.600259587147127,
2.18733299301658, 1.53261062618519, -0.235700359100477, -1.02642090030678
), llAlgorithm = "GHK", llNGHK = 1000, llIntGrad = TRUE,
llOneSidedGrad = FALSE, llEps = 1e-06, llRandom.seed = 123)
where 9: do.call(callWithoutSumt, argList)
where 10: maxOptim(fn = fn, grad = grad, hess = hess, start = start, method = "BFGS",
fixed = fixed, constraints = constraints, finalHessian = finalHessian,
parscale = parscale, control = mControl, ...)
where 11: maxRoutine(fn = logLik, grad = grad, hess = hess, start = start,
constraints = constraints, ...)
where 12: maxLik(logLik = logLik, start = start, method = method, finalHessian = finalHessian,
yMat = yMat, xMat = xMat, llAlgorithm = algorithm, llNGHK = nGHK,
llIntGrad = intGrad, llOneSidedGrad = oneSidedGrad, llEps = eps,
llRandom.seed = random.seed, ...)
where 13: mvProbit(cbind(y1, y2, y3) ~ x1 + x2, start = c(beta), startSigma = sigma,
data = dat, method = "BFGS", reltol = 0.5)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x2797898>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbit ... logLikFunc -> fnOrig -> mvProbitLogLikInternal -> pmvnormWrap
Execution halted
Running the tests in ‘tests/mvProbitMargEff2.R’ failed.
Complete output:
> # I thank Mohit Batham for providing this R script that demonstrated
> # a bug in mvProbitMargEff() when called with 2 dependent variables
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> nObs <- 100
> set.seed( 123 )
> xData <- data.frame(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = as.numeric( rnorm( nObs ) > 0 ),
+ x3 = rnorm( nObs ),
+ x4 = rnorm( nObs ))
> beta <- c( 0.8, 1.2, -1.0, 1.4, -0.8,
+ -0.6, 1.0, 0.6, -1.2, -1.6)
> sigma <- symMatrix( c( 1, 0.2, 1))
> margEffCond <- try( mvProbitMargEff( ~ x1 + x2 + x3 + x4 , coef = beta,
+ sigma = sigma, data = xData, cond = TRUE ) )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = betaLower,
sigma = coef$sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
random.seed = random.seed, ...)
where 3: mvProbitMargEffInternal(yMat = yMat, xMat = xMat, coef = coef,
sigma = sigma, cond = cond, algorithm = algorithm, nGHK = nGHK,
eps = eps, dummyVars = dummyVars, random.seed = random.seed,
...)
where 4: mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
data = xData, cond = TRUE)
where 5: doTryCatch(return(expr), name, parentenv, handler)
where 6: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 7: tryCatchList(expr, classes, parentenv, handlers)
where 8: tryCatch(expr, error = function(e) {
call <- conditionCall(e)
if (!is.null(call)) {
if (identical(call[[1L]], quote(doTryCatch)))
call <- sys.call(-4L)
dcall <- deparse(call)[1L]
prefix <- paste("Error in", dcall, ": ")
LONG <- 75L
sm <- strsplit(conditionMessage(e), "\n")[[1L]]
w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w")
if (is.na(w))
w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L],
type = "b")
if (w > LONG)
prefix <- paste0(prefix, "\n ")
}
else prefix <- "Error : "
msg <- paste0(prefix, conditionMessage(e), "\n")
.Internal(seterrmessage(msg[1L]))
if (!silent && isTRUE(getOption("show.error.messages"))) {
cat(msg, file = outFile)
.Internal(printDeferredWarnings())
}
invisible(structure(msg, class = "try-error", condition = e))
})
where 9: try(mvProbitMargEff(~x1 + x2 + x3 + x4, coef = beta, sigma = sigma,
data = xData, cond = TRUE))
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x37a32c0>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
> round( margEffCond, 3 )
Error in round(margEffCond, 3) :
non-numeric argument to mathematical function
Execution halted
Running the tests in ‘tests/mvProbitTest.R’ failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> options( digits = 4 )
>
> ## generate a simulated data set
> set.seed( 123 )
> # number of observations
> nObs <- 10
>
> # generate explanatory variables
> xMat <- cbind(
+ const = rep( 1, nObs ),
+ x1 = as.numeric( rnorm( nObs ) > 0 ),
+ x2 = as.numeric( rnorm( nObs ) > 0 ),
+ x3 = rnorm( nObs ),
+ x4 = rnorm( nObs ) )
>
> # model coefficients
> beta <- cbind( c( 0.8, 1.2, -1.0, 1.4, -0.8 ),
+ c( -0.6, 1.0, 0.6, -1.2, -1.6 ),
+ c( 0.5, -0.6, -0.7, 1.1, 1.2 ) )
>
> # covariance matrix of error terms
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> # all parameters in a vector
> allCoef <- c( c( beta ), sigma[ lower.tri( sigma ) ] )
>
> # generate dependent variables
> yMatLin <- xMat %*% beta
> yMat <- ( yMatLin + rmvnorm( nObs, sigma = sigma, pre0.9_9994 = TRUE ) ) > 0
> colnames( yMat ) <- paste( "y", 1:3, sep = "" )
> # (yMatLin > 0 )== yMat
>
> # unconditional expectations of dependent variables
> yExp <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ) )
> round( yExp, 3 )
V1 V2 V3
1 0.021 0.725 0.194
2 0.394 0.768 0.214
3 0.125 0.788 0.196
4 0.235 0.681 0.292
5 0.680 0.435 0.579
6 0.028 0.973 0.034
7 0.958 0.186 0.784
8 0.856 0.247 0.724
9 0.061 0.968 0.034
10 0.998 0.067 0.923
> yExpA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
+ data = as.data.frame( xMat ) )
> all.equal( yExp, yExpA )
[1] TRUE
> yExp2 <- pnorm( yMatLin )
> all.equal( yExp, as.data.frame( yExp2 ) )
[1] TRUE
>
>
> # conditional expectations of dependent variables
> # (assuming that all other dependent variables are one)
> yExpCond <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz() )
> round( yExpCond, 3 )
V1 V2 V3
1 0.072 0.834 0.524
2 0.660 0.779 0.314
3 0.297 0.838 0.399
4 0.448 0.727 0.461
5 0.847 0.442 0.618
6 0.139 0.985 0.162
7 0.991 0.179 0.749
8 0.950 0.244 0.710
9 0.244 0.978 0.133
10 1.000 0.065 0.892
> yExpCondA <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = allCoef,
+ data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz() )
> all.equal( yExpCond, yExpCondA )
[1] TRUE
> yExpCond2 <- matrix( NA, nrow = nObs, ncol = ncol( yMat ) )
> for( i in 1:nObs ) {
+ for( k in 1:ncol( yMat ) ) {
+ set.seed( 123 )
+ numerator <- pmvnorm( upper = yMatLin[ i, ], sigma = sigma )
+ set.seed( 123 )
+ denominator <- pmvnorm( upper = yMatLin[ i, -k ], sigma = sigma[ -k, -k ] )
+ yExpCond2[ i, k ] <- numerator / denominator
+ }
+ }
> all.equal( yExpCond, as.data.frame( yExpCond2 ) )
[1] TRUE
> # now with explicitly specifying the algorithm
> yExpCond3 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = GenzBretz )
> all.equal( yExpCond, yExpCond3 )
[1] TRUE
> identical( yExpCond, yExpCond3 )
[1] TRUE
> # now with integrals obtained by the Miwa algorithm
> yExpCond4 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = Miwa )
> all.equal( yExpCond, yExpCond4, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the Miwa algorithm, less precise
> yExpCond5 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = Miwa( steps = 32 ) )
> all.equal( yExpCond4, yExpCond5, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the TVPACK algorithm
> yExpCond6 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = TVPACK )
> all.equal( yExpCond, yExpCond6, tol = 1e-3 )
[1] TRUE
> # now with integrals obtained by the TVPACK algorithm, less precise
> yExpCond7 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE,
+ algorithm = TVPACK( abseps = 0.5 ) )
> all.equal( yExpCond6, yExpCond7, tol = 1e-3 )
[1] "Component \"V1\": Mean relative difference: 0.03018"
[2] "Component \"V2\": Mean relative difference: 0.04273"
[3] "Component \"V3\": Mean relative difference: 0.04242"
> # now with integrals obtained by the GHK algorithm
> yExpCond8 <- mvProbitExp( ~ x1 + x2 + x3 + x4, coef = c( beta ),
+ sigma = sigma, data = as.data.frame( xMat ), cond = TRUE )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: pmvnormWrap(upper = xBeta[i, ], sigma = coef$sigma, algorithm = algorithm,
nGHK = nGHK, random.seed = random.seed, ...)
where 2: mvProbitExpInternal(yMat = yMat, xMat = xMat, coef = coef, sigma = sigma,
cond = cond, algorithm = algorithm, nGHK = nGHK, random.seed = random.seed,
...)
where 3: mvProbitExp(~x1 + x2 + x3 + x4, coef = c(beta), sigma = sigma,
data = as.data.frame(xMat), cond = TRUE)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x20ade30>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: mvProbitExp -> mvProbitExpInternal -> pmvnormWrap
Execution halted
Running the tests in ‘tests/pmvnormWrapTest.R’ failed.
Complete output:
> library( "mvProbit" )
Loading required package: mvtnorm
Loading required package: maxLik
Loading required package: miscTools
Please cite the 'maxLik' package as:
Henningsen, Arne and Toomet, Ott (2011). maxLik: A package for maximum likelihood estimation in R. Computational Statistics 26(3), 443-458. DOI 10.1007/s00180-010-0217-1.
If you have questions, suggestions, or comments regarding the 'maxLik' package, please use a forum or 'tracker' at maxLik's R-Forge site:
https://r-forge.r-project.org/projects/maxlik/
Loading required package: abind
> library( "miscTools" )
>
> # covariance matrix
> sigma <- miscTools::symMatrix( c( 1, 0.2, 0.4, 1, -0.1, 1 ) )
>
> ######## only upper ##########
> upper <- c( -0.3, 0.7, -0.5 )
> # Genz + Bretz (default)
> pug <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = GenzBretz(), random.seed = 123 )
> print( pug )
[1] 0.1360783
attr(,"error")
[1] 0.000118605
attr(,"msg")
[1] "Normal Completion"
>
> # Miwa (as function)
> pum <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = Miwa, random.seed = 123 )
> print( pum )
[1] 0.136069
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"
> all.equal( pug, pum, check.attributes = FALSE, tol = 1e-4 )
[1] TRUE
>
> # Miwa (as object returned from function Miwa())
> pum1 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = Miwa(), random.seed = 123 )
> all.equal( pum, pum1 )
[1] TRUE
>
> # Miwa (as character string)
> pum2 <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = "Miwa", random.seed = 123 )
> all.equal( pum, pum2 )
[1] TRUE
>
> # TVPACK
> put <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = TVPACK, random.seed = 123 )
> print( put )
[1] 0.136069
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"
> all.equal( pug, put, check.attributes = FALSE, tol = 1e-4 )
[1] TRUE
> all.equal( pum, put, check.attributes = FALSE, tol = 1e-6 )
[1] TRUE
>
> # GHK
> pughk <- mvProbit:::pmvnormWrap( upper = upper, sigma = sigma,
+ algorithm = "ghk", random.seed = 123, nGHK = 1000 )
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
mvProbit
--- call from context ---
mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
random.seed = 123, nGHK = 1000)
--- call from argument ---
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
--- R stacktrace ---
where 1: mvProbit:::pmvnormWrap(upper = upper, sigma = sigma, algorithm = "ghk",
random.seed = 123, nGHK = 1000)
--- value of length: 2 type: logical ---
[1] FALSE FALSE
--- function from context ---
function (lower = -Inf, upper = Inf, sigma, algorithm, random.seed,
nGHK = NULL, ...)
{
if (!is.matrix(sigma)) {
stop("argument 'sigma' must be a matrix")
}
else if (nrow(sigma) != ncol(sigma)) {
stop("argument 'sigma' must be a square matrix")
}
else if (!all(is.numeric(sigma))) {
stop("argument 'sigma must be a numeric matrix")
}
if (!all(is.numeric(lower))) {
stop("argument 'lower' must be numeric")
}
else if (length(lower) == 1) {
lower <- rep(lower, nrow(sigma))
}
else if (length(lower) != nrow(sigma)) {
stop("argument 'lower' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
if (!all(is.numeric(upper))) {
stop("argument 'upper' must be numeric")
}
else if (length(upper) == 1) {
upper <- rep(upper, nrow(sigma))
}
else if (length(upper) != nrow(sigma)) {
stop("argument 'upper' must either be a single numeric value",
" or a vector with length equal to the number of rows/columns",
" of argument 'sigma'")
}
algOkay <- TRUE
ghk <- FALSE
if (!is.list(algorithm) && length(algorithm) != 1) {
stop("argument 'algorithm' must be a single function",
" or a single character string")
}
else if (is.function(algorithm)) {
algResult <- do.call(algorithm, list())
if (!class(algResult)[1] %in% c("GenzBretz", "Miwa",
"TVPACK")) {
algOkay <- FALSE
}
}
else if (is.character(algorithm)) {
if (tolower(algorithm) == "ghk") {
ghk <- TRUE
}
else if (!algorithm %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
}
else if (!class(algorithm) %in% c("GenzBretz", "Miwa", "TVPACK")) {
algOkay <- FALSE
}
if (!algOkay) {
stop("argument 'algorithm' must be either one of the functions",
" 'GenzBretz()', 'Miwa()', or 'TVPACK()'", " or one of the character strings",
" \"GenzBretz\", \"Miwa\", or \"TVPACK\"")
}
if (length(random.seed) != 1) {
stop("argument 'random.seed' must be a single numerical value")
}
else if (!is.numeric(random.seed)) {
stop("argument 'random.seed' must be numerical")
}
if (exists(".Random.seed")) {
savedSeed <- .Random.seed
}
set.seed(random.seed)
if (exists("savedSeed")) {
on.exit(assign(".Random.seed", savedSeed, envir = sys.frame()))
}
else {
on.exit(rm(.Random.seed, envir = sys.frame()))
}
if (ghk) {
if (is.null(nGHK)) {
stop("if the GHK algorithm is used,", " argument 'nGHK' must be specified")
}
else if (length(nGHK) != 1) {
stop("argument 'nGHK' must be a single integer value")
}
else if (!is.numeric(nGHK)) {
stop("argument 'nGHK' must be numeric")
}
else if (nGHK <= 0) {
stop("argument 'nGHK' must be positive")
}
L <- try(t(chol(sigma)), silent = TRUE)
if (class(L) == "try-error") {
warning("the correlation matrix is not positive definite")
return(NA)
}
trunpt <- rep(NA, length(lower))
above <- rep(NA, length(lower))
for (i in 1:length(lower)) {
if (lower[i] == -Inf) {
trunpt[i] <- upper[i]
above[i] <- 1
}
else if (upper[i] == Inf) {
trunpt[i] <- lower[i]
above[i] <- 0
}
else {
stop("if algorithm 'GHK' is used,", " either the lower truncation point must be '-Inf'",
" or the upper truncation point must be 'Inf'")
}
}
sink(tempfile())
on.exit(sink(), add = TRUE)
result <- ghkvec(L = L, trunpt = trunpt, above = above,
r = nGHK)
result <- drop(result)
}
else {
result <- pmvnorm(lower = lower, upper = upper, sigma = sigma,
algorithm = algorithm, ...)
}
return(result)
}
<bytecode: 0x39b0a58>
<environment: namespace:mvProbit>
--- function search by body ---
Function pmvnormWrap in namespace mvProbit has this body.
----------- END OF FAILURE REPORT --------------
Error in if (class(L) == "try-error") { : the condition has length > 1
Calls: <Anonymous>
Execution halted
Flavor: r-devel-linux-x86_64-fedora-gcc