Running examples in 'SemiParSampleSel-Ex.R' failed
The error most likely occurred in:

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: SemiParSampleSel
> ### Title: Semiparametric Sample Selection Modelling with Continuous or
> ###   Discrete Response
> ### Aliases: SemiParSampleSel
> ### Keywords: copula sample selection semiparametric sample selection
> ###   modelling sample selection model smooth regression spline shrinkage
> ###   smoother variable selection
> 
> ### ** Examples
> 
> 
> library(SemiParSampleSel)
> 
> ######################################################################
> ## Generate data
> ## Correlation between the two equations and covariate correlation 0.5 
> ## Sample size 2000 
> ######################################################################
> 
> set.seed(0)
> 
> n <- 2000
> 
> rhC <- rhU <- 0.5      
> 
> SigmaU <- matrix(c(1, rhU, rhU, 1), 2, 2)
> U      <- rmvnorm(n, rep(0,2), SigmaU)
> 
> SigmaC <- matrix(rhC, 3, 3); diag(SigmaC) <- 1
> 
> cov    <- rmvnorm(n, rep(0,3), SigmaC, method = "svd")
> cov    <- pnorm(cov)
> 
> bi <- round(cov[,1]); x1 <- cov[,2]; x2 <- cov[,3]
>   
> f11 <- function(x) -0.7*(4*x + 2.5*x^2 + 0.7*sin(5*x) + cos(7.5*x))
> f12 <- function(x) -0.4*( -0.3 - 1.6*x + sin(5*x))  
> f21 <- function(x) 0.6*(exp(x) + sin(2.9*x)) 
> 
> ys <-  0.58 + 2.5*bi + f11(x1) + f12(x2) + U[, 1] > 0
> y  <- -0.68 - 1.5*bi + f21(x1) +         + U[, 2]
> yo <- y*(ys > 0)
>   
> 
> dataSim <- data.frame(ys, yo, bi, x1, x2)
> 
> ## CLASSIC SAMPLE SELECTION MODEL
> ## the first equation MUST be the selection equation
> 
> out <- SemiParSampleSel(list(ys ~ bi + x1 + x2, 
+                              yo ~ bi + x1), 
+                         data = dataSim)
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
SemiParSampleSel
 --- call from context --- 
SemiParSampleSel(list(ys ~ bi + x1 + x2, yo ~ bi + x1), data = dataSim)
 --- call from argument --- 
if (class(X2s) == "try-error") stop("Check that the range of the covariate values of the selected sample is the same as that of the complete dataset.")
 --- R stacktrace ---
where 1: SemiParSampleSel(list(ys ~ bi + x1 + x2, yo ~ bi + x1), data = dataSim)

 --- value of length: 2 type: logical ---
[1] FALSE FALSE
 --- function from context --- 
function (formula, data = list(), weights = NULL, subset = NULL, 
    start.v = NULL, start.theta = NULL, BivD = "N", margins = c("probit", 
        "N"), fp = FALSE, infl.fac = 1, rinit = 1, rmax = 100, 
    iterlimsp = 50, pr.tolsp = 1e-06, bd = NULL, parscale) 
{
    if (margins[2] == "GA") 
        margins[2] <- "G"
    if (margins[1] == "probit") 
        margins[1] <- "N"
    marg2 <- c("N", "G", "P", "NB", "D", "PIG", "S", "BB", "BI", 
        "GEOM", "LG", "NBII", "WARING", "YULE", "ZIBB", "ZABB", 
        "ZABI", "ZIBI", "ZALG", "ZANBI", "ZINBI", "ZAP", "ZIP", 
        "ZIP2", "ZIPIG")
    marg2D <- c("P", "NB", "D", "PIG", "S", "BB", "BI", "GEOM", 
        "LG", "NBII", "WARING", "YULE", "ZIBB", "ZABB", "ZABI", 
        "ZIBI", "ZALG", "ZANBI", "ZINBI", "ZAP", "ZIP", "ZIP2", 
        "ZIPIG")
    copul <- c("N", "FGM", "F", "AMH", "C0", "C90", "C180", "C270", 
        "J0", "J90", "J180", "J270", "G0", "G90", "G180", "G270")
    if (!(margins[1] %in% c("N", "probit"))) 
        stop("Error in selection equation's margin name. It can currently only be probit or equivalently N.")
    if (!(margins[2] %in% marg2)) 
        stop("Error in outcome equation's margin name. It should be one of: N, GA, P, NB, D, PIG, S, BB, BI, GEOM, LG, NBII, WARING, YULE, ZIBB, ZABB, ZABI, ZIBI, ZALG, ZANBI, ZINBI, ZAP, ZIP, ZIP2, ZIPIG.")
    if (!(BivD %in% copul)) 
        stop("Error in parameter BivD value. It should be one of: N, FGM, F, AMH, C0, C90, C180, C270, J0, J90, J180, J270, G0, G90, G180, G270.")
    if (length(formula) > 2 && margins[2] %in% c("N", "G", "NB", 
        "PIG", "BB", "NBII", "WARING", "ZABI", "ZIBI", "ZALG", 
        "ZAP", "ZIP", "ZIP2")) {
        if (length(formula) != 4) 
            stop("You need to specify four equations.")
    }
    if (length(formula) > 2 && margins[2] %in% c("P", "BI", "GEOM", 
        "LG", "YULE")) {
        if (length(formula) != 3) 
            stop("You need to specify three equations.")
    }
    if (length(formula) > 2 && margins[2] %in% c("D", "S", "ZIBB", 
        "ZABB", "ZANBI", "ZINBI", "ZIPIG")) {
        if (length(formula) != 5) 
            stop("You need to specify five equations.")
    }
    qu.mag <- sp <- nu <- gam2.1 <- KendTau <- gamlss.fit <- edf <- edf1 <- NULL
    sp3 <- gp3 <- gam3 <- X3 <- sp.theta <- NULL
    sp4 <- gp4 <- gam4 <- X4 <- NULL
    sp5 <- gp5 <- gam5 <- X5 <- NULL
    X3.d2 <- X4.d2 <- X5.d2 <- l.sp3 <- l.sp4 <- l.sp5 <- 0
    if (length(formula) > 2) {
        f3t <- try(formula[[3]][[3]], silent = TRUE)
        if (class(f3t) != "try-error") 
            stop("The third equation does not require a response.")
        if (length(formula) > 3) {
            f4t <- try(formula[[4]][[3]], silent = TRUE)
            if (class(f4t) != "try-error") 
                stop("The fourth equation does not require a response.")
            if (length(formula) > 4) {
                f5t <- try(formula[[5]][[3]], silent = TRUE)
                if (class(f5t) != "try-error") 
                  stop("The fifth equation does not require a response.")
            }
        }
    }
    ig <- interpret.gam(formula)
    mf <- match.call(expand.dots = FALSE)
    if (length(formula) == 2) 
        pred.n <- union(ig[[1]]$pred.names, c(ig[[2]]$pred.names, 
            ig[[2]]$response))
    if (length(formula) == 3) 
        pred.n <- union(ig[[1]]$pred.names, c(ig[[2]]$pred.names, 
            ig[[3]]$pred.names, ig[[2]]$response))
    if (length(formula) == 4) 
        pred.n <- union(ig[[1]]$pred.names, c(ig[[2]]$pred.names, 
            ig[[3]]$pred.names, ig[[4]]$pred.names, ig[[2]]$response))
    if (length(formula) == 5) 
        pred.n <- union(ig[[1]]$pred.names, c(ig[[2]]$pred.names, 
            ig[[3]]$pred.names, ig[[4]]$pred.names, ig[[5]]$pred.names, 
            ig[[2]]$response))
    if (length(formula) == 6) 
        pred.n <- union(ig[[1]]$pred.names, c(ig[[2]]$pred.names, 
            ig[[3]]$pred.names, ig[[4]]$pred.names, ig[[5]]$pred.names, 
            ig[[6]]$pred.names, ig[[2]]$response))
    fake.formula <- paste(ig[[1]]$response, "~", paste(pred.n, 
        collapse = " + "))
    environment(fake.formula) <- environment(ig$fake.formula)
    mf$formula <- fake.formula
    mf$start.v <- mf$method <- mf$start.theta <- mf$BivD <- mf$margins <- mf$infl.fac <- mf$fp <- mf$rinit <- mf$rmax <- mf$iterlimsp <- mf$pr.tolsp <- mf$parscale <- mf$bd <- NULL
    mf$drop.unused.levels <- TRUE
    mf$na.action <- na.pass
    mf[[1]] <- as.name("model.frame")
    data <- eval(mf, parent.frame())
    indS <- as.logical(data[, ig[[1]]$response]) == FALSE
    indS <- ifelse(is.na(indS), FALSE, indS)
    data[indS, ig[[2]]$response] <- ifelse(is.na(data[indS, ig[[2]]$response]), 
        0, data[indS, ig[[2]]$response])
    data <- na.omit(data)
    if (is.null(weights)) {
        weights <- rep(1, dim(data)[1])
        data$weights <- weights
        names(data)[length(names(data))] <- "(weights)"
    }
    else weights <- data[, "(weights)"]
    formula.eq1 <- formula[[1]]
    formula.eq2 <- formula[[2]]
    gam1 <- eval(substitute(gam(formula.eq1, binomial(link = "probit"), 
        gamma = infl.fac, weights = weights, data = data), list(weights = weights)))
    y1 <- gam1$y
    inde <- y1 > 0
    if (margins[2] == "N") 
        gam2 <- eval(substitute(gam(formula.eq2, gamma = infl.fac, 
            data = data, weights = weights, subset = inde), list(weights = weights, 
            inde = inde)))
    if (margins[2] == "G") 
        gam2 <- eval(substitute(gam(formula.eq2, gamma = infl.fac, 
            data = data, weights = weights, subset = inde, family = Gamma(link = "log")), 
            list(weights = weights, inde = inde)))
    if (margins[2] %in% c("P", "NB", "D", "PIG", "S", "GEOM", 
        "NBII", "WARING", "YULE", "ZANBI", "ZINBI", "ZAP", "ZIP", 
        "ZIP2", "ZIPIG", "LG", "ZALG", "BB", "BI", "ZIBB", "ZABB", 
        "ZABI", "ZIBI")) {
        gam2 <- eval(substitute(gam(formula.eq2, gamma = infl.fac, 
            data = data, weights = weights, subset = inde, family = poisson(link = "log")), 
            list(weights = weights, inde = inde)))
    }
    n <- length(y1)
    X1 <- model.matrix(gam1)
    X1.d2 <- dim(X1)[2]
    X2.d2 <- length(coef(gam2))
    X2 <- matrix(0, n, X2.d2, dimnames = list(c(1:n), c(names(coef(gam2)))))
    X2[inde, ] <- model.matrix(gam2)
    if (margins[2] == "LG") {
        y2 <- rep(1, n)
    }
    else {
        y2 <- rep(0, n)
    }
    y2[inde] <- gam2$y
    l.sp1 <- length(gam1$smooth)
    if (l.sp1 != 0) 
        sp1 <- gam1$sp
    else sp1 <- NULL
    l.sp2 <- length(gam2$smooth)
    if (l.sp2 != 0) 
        sp2 <- gam2$sp
    else sp2 <- NULL
    gp1 <- gam1$nsdf
    gp2 <- gam2$nsdf
    if (length(formula) == 3) 
        theta <- rnorm(n)
    nad1 <- "theta"
    if (length(formula) == 4) {
        sigma <- rnorm(n)
        nad1 <- "sigma"
        theta <- rnorm(n)
        nad2 <- "theta"
    }
    if (length(formula) == 5) {
        sigma <- rnorm(n)
        nad1 <- "sigma"
        nu <- rnorm(n)
        nad2 <- "nu"
        theta <- rnorm(n)
        nad3 <- "theta"
    }
    if (length(formula) > 2) {
        formula.eq3 <- formula[[3]]
        formula.eq3 <- as.formula(paste(nad1, "~", formula.eq3[2], 
            sep = ""))
        gam3 <- eval(substitute(gam(formula.eq3, weights = weights, 
            data = data)))
        X3 <- model.matrix(gam3)
        X3.d2 <- llM <- dim(X3)[2]
        l.sp3 <- length(gam3$sp)
        if (l.sp3 != 0) 
            sp3 <- sp.theta <- gam3$sp
        else sp3 <- NULL
        environment(gam3$formula) <- environment(gam2$formula)
        gp3 <- gam3$nsdf
        a.theta <- coef(gam3)
    }
    if (length(formula) > 3) {
        formula.eq4 <- formula[[4]]
        formula.eq4 <- as.formula(paste(nad2, "~", formula.eq4[2], 
            sep = ""))
        gam4 <- eval(substitute(gam(formula.eq4, weights = weights, 
            data = data)))
        X4 <- model.matrix(gam4)
        X4.d2 <- llM <- dim(X4)[2]
        l.sp4 <- length(gam4$sp)
        if (l.sp4 != 0) 
            sp4 <- sp.theta <- gam4$sp
        else sp4 <- NULL
        environment(gam4$formula) <- environment(gam2$formula)
        gp4 <- gam4$nsdf
        a.theta <- coef(gam4)
    }
    if (length(formula) > 4) {
        formula.eq5 <- formula[[5]]
        formula.eq5 <- as.formula(paste(nad3, "~", formula.eq5[2], 
            sep = ""))
        gam5 <- eval(substitute(gam(formula.eq5, weights = weights, 
            data = data)))
        X5 <- model.matrix(gam5)
        X5.d2 <- llM <- dim(X5)[2]
        l.sp5 <- length(gam5$sp)
        if (l.sp5 != 0) 
            sp5 <- sp.theta <- gam5$sp
        else sp5 <- NULL
        environment(gam5$formula) <- environment(gam2$formula)
        gp5 <- gam5$nsdf
        a.theta <- coef(gam5)
    }
    if ((l.sp1 != 0 || l.sp2 != 0 || l.sp3 != 0 || l.sp4 != 0 || 
        l.sp5 != 0) && fp == FALSE) {
        qu.mag <- S.m(gam1, gam2, gam3, gam4, gam5, l.sp1, l.sp2, 
            l.sp3, l.sp4, l.sp5, margins, eq1 = "yes")
        sp <- c(sp1, sp2, sp3, sp4, sp5)
        qu.mag1 <- S.m(gam1, gam2, gam3, gam4, gam5, l.sp1, l.sp2, 
            l.sp3, l.sp4, l.sp5, margins, eq1 = "no")
        sp1 <- c(sp2, sp3, sp4, sp5)
    }
    if (length(formula) < 3) {
        p.g1 <- predict(gam1)
        imr <- data$imr <- dnorm(p.g1)/pnorm(p.g1)
        data$bd <- bd
        formula.eq2.1 <- update.formula(formula.eq2, ~. + imr)
        if (margins[2] %in% c("P", "NB", "D", "PIG", "S", "NBII", 
            "ZANBI", "ZINBI", "ZAP", "ZIP", "ZIP2", "ZIPIG")) {
            formula.eq2.1 <- update.formula(formula.eq2, round((. + 
                mean(.[. >= 0]))/2) ~ . + imr)
        }
        else if (margins[2] %in% c("BB", "BI")) {
            formula.eq2.1 <- update.formula(formula.eq2, cbind(round(bd * 
                (. + 0.5)/(bd + 1)), bd - round(bd * (. + 0.5)/(bd + 
                1))) ~ . + imr)
        }
        else if (margins[2] %in% c("GEOM", "YULE", "WARING", 
            "LG", "ZALG")) {
            formula.eq2.1 <- update.formula(formula.eq2, ~. + 
                imr)
        }
        else if (margins[2] %in% c("ZABI", "ZIBI", "ZIBB", "ZABB")) {
            formula.eq2.1 <- update.formula(formula.eq2, cbind(., 
                bd - .) ~ . + imr)
        }
        if (margins[2] == "N") 
            gam2.1 <- eval(substitute(gam(formula.eq2.1, gamma = infl.fac, 
                data = data, weights = weights, subset = inde), 
                list(weights = weights, inde = inde)))
        if (margins[2] == "G") 
            gam2.1 <- eval(substitute(gam(formula.eq2.1, gamma = infl.fac, 
                data = data, weights = weights, subset = inde, 
                family = Gamma(link = "log")), list(weights = weights, 
                inde = inde)))
        if (margins[2] %in% c("P", "NB", "D", "PIG", "S", "GEOM", 
            "NBII", "WARING", "YULE", "ZANBI", "ZINBI", "ZAP", 
            "ZIP", "ZIP2", "ZIPIG", "LG", "ZALG")) {
            gam2.1 <- eval(substitute(gam(formula.eq2.1, gamma = infl.fac, 
                data = data, weights = weights, subset = inde, 
                family = poisson(link = "log")), list(weights = weights, 
                inde = inde)))
        }
        if (margins[2] %in% c("BB", "BI", "ZIBB", "ZABB", "ZABI", 
            "ZIBI")) {
            gam2.1 <- eval(substitute(gam(formula.eq2.1, gamma = infl.fac, 
                data = data, weights = weights, subset = inde, 
                family = binomial(link = "logit")), list(weights = weights, 
                inde = inde)))
        }
        sigma <- sqrt(mean(residuals(gam2.1, type = "deviance")^2) + 
            mean(imr[inde] * (imr[inde] + p.g1[inde])) * gam2.1$coef["imr"]^2)[[1]]
        names(sigma) <- "sigma"
        co <- (gam2.1$coef["imr"]/sigma)[[1]]
        if (margins[2] == "G") {
            sigma <- k <- (summary(gam2)$dispersion)^(-1)
            names(k) <- "sigma"
        }
        a.theta <- st.theta.star(start.theta, co, BivD)
        names(a.theta) <- "theta.star"
        if (margins[2] %in% c("N", "NB", "PIG", "WARING", "BB", 
            "NBII")) {
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], log(sigma), a.theta)
            names(start.v)[length(start.v) - 1] <- "log.sigma"
        }
        if (margins[2] %in% c("ZIP", "ZAP", "ZIP2", "ZALG", "ZIBI", 
            "ZABI")) {
            sigma.star.start1 <- sum(gam2.1$y == 0)/length(gam2.1$y)
            sigma.star.start1 <- ifelse(sigma.star.start1 > 0.999, 
                0.999, sigma.star.start1)
            sigma.star.start1 <- ifelse(sigma.star.start1 < 0.001, 
                0.001, sigma.star.start1)
            sigma.star.start1 <- qlogis(sigma.star.start1)
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], sigma.star.start1, a.theta)
            names(start.v)[length(start.v) - 1] <- "logi.sigma"
        }
        if (margins[2] == "G") {
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], log(k), a.theta)
            names(start.v)[length(start.v) - 1] <- "log.shape"
        }
        if (margins[2] %in% c("P", "YULE", "GEOM", "BI", "LG")) {
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], a.theta)
        }
        if (margins[2] %in% c("D")) {
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], log(sigma), qlogis(0.51), a.theta)
            names(start.v)[length(start.v) - 2] <- "log.sigma"
            names(start.v)[length(start.v) - 1] <- "logi.nu"
        }
        if (margins[2] %in% c("ZINBI")) {
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], log(sigma), qlogis(((sum(gam2$y == 
                  0)/length(gam2$y)) + 0.01)/2), a.theta)
            names(start.v)[length(start.v) - 2] <- "log.sigma"
            names(start.v)[length(start.v) - 1] <- "logi.nu"
        }
        if (margins[2] %in% c("ZANBI")) {
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], log(sigma), qlogis(max(sum(gam2$y == 
                  0)/length(gam2$y), 0.01)), a.theta)
            names(start.v)[length(start.v) - 2] <- "log.sigma"
            names(start.v)[length(start.v) - 1] <- "logi.nu"
        }
        if (margins[2] %in% c("ZIPIG")) {
            nu.star.start1 <- sum(gam2.1$y == 0)/length(gam2.1$y)
            nu.star.start1 <- ifelse(nu.star.start1 > 0.999, 
                0.999, nu.star.start1)
            nu.star.start1 <- ifelse(nu.star.start1 < 0.001, 
                0.001, nu.star.start1)
            nu.star.start1 <- qlogis(nu.star.start1)
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], log(sigma), nu.star.start1, a.theta)
            names(start.v)[length(start.v) - 2] <- "log.sigma"
            names(start.v)[length(start.v) - 1] <- "logi.nu"
        }
        if (margins[2] %in% c("ZIBB", "ZABB")) {
            nu.star.start1 <- sum(gam2.1$y == 0)/length(gam2.1$y)
            nu.star.start1 <- ifelse(nu.star.start1 > 0.999, 
                0.999, nu.star.start1)
            nu.star.start1 <- ifelse(nu.star.start1 < 0.001, 
                0.001, nu.star.start1)
            nu.star.start1 <- qlogis(nu.star.start1)
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], log(sigma), nu.star.start1, a.theta)
            names(start.v)[length(start.v) - 2] <- "log.sigma"
            names(start.v)[length(start.v) - 1] <- "logi.nu"
        }
        if (margins[2] == "S") {
            if (is.null(start.v)) 
                start.v <- c(coef(gam1), coef(gam2.1)[names(coef(gam2.1)) != 
                  "imr"], log(sigma), -0.5, a.theta)
            names(start.v)[length(start.v) - 2] <- "log.sigma"
            names(start.v)[length(start.v) - 1] <- "nu"
        }
        start.v1 <- start.v[-c(1:X1.d2, length(start.v))]
    }
    else {
        start.v <- start.v <- c(coef(gam1), coef(gam2), coef(gam3), 
            coef(gam4), coef(gam5))
        start.v1 <- start.v[-c(1:X1.d2, (length(start.v) - (llM - 
            1)):length(start.v))]
    }
    if (missing(parscale)) 
        parscale <- 1
    if (margins[2] %in% c("N", "G")) 
        funcop <- ghss
    else funcop <- ghssD
    respvec <- cbind(y1, y2)
    VC <- list(X1 = X1, X2 = X2, X3 = X3, X4 = X4, X5 = X5, X1.d2 = X1.d2, 
        X2.d2 = X2.d2, X3.d2 = X3.d2, X4.d2 = X4.d2, X5.d2 = X5.d2, 
        gp1 = gp1, gp2 = gp2, gp3 = gp3, gp4 = gp4, gp5 = gp5, 
        n = n, bd = bd, l.sp1 = l.sp1, l.sp2 = l.sp2, l.sp3 = l.sp3, 
        l.sp4 = l.sp4, l.sp5 = l.sp5, infl.fac = infl.fac, weights = weights, 
        fp = fp, BivD = BivD, margins = margins, parscale = parscale)
    if (margins[2] %in% marg2D) {
        gamlss.fit <- fit.SemiParSampleSel(funcop = ghssDuniv, 
            start.v = start.v1, dat = respvec, rinit = rinit, 
            rmax = rmax, VC = VC, qu.mag = qu.mag1, sp = sp1, 
            iterlimsp = iterlimsp, pr.tolsp = pr.tolsp, univariate = TRUE)
        start.v <- c(coef(gam1), gamlss.fit$fit$argument, a.theta)
        sp <- c(sp1, gamlss.fit$sp, sp.theta)
        gam2$coefficients <- gamlss.fit$fit$argument[1:X2.d2]
        gam2$Vp <- gamlss.fit$magpp$Vb[1:X2.d2, 1:X2.d2]
        gam2$sig2 <- 1
        gam2$edf <- gamlss.fit$magpp$edf[1:X2.d2]
        if (length(formula) > 3) {
            if (margins[2] %in% c("NB", "PIG", "D", "S", "BB", 
                "NBII", "WARING", "ZABI", "ZIBI", "ZALG", "ZAP", 
                "ZIP", "ZIP2", "ZIBB", "ZABB", "ZANBI", "ZINBI", 
                "ZIPIG")) {
                gam3$coefficients <- gamlss.fit$fit$argument[X2.d2 + 
                  (1:X3.d2)]
                gam3$Vp <- gamlss.fit$magpp$Vb[X2.d2 + (1:X3.d2), 
                  X2.d2 + (1:X3.d2)]
                gam3$sig2 <- 1
                gam3$edf <- gamlss.fit$magpp$edf[X2.d2 + (1:X3.d2)]
            }
            if (margins[2] %in% c("D", "S", "ZIBB", "ZABB", "ZANBI", 
                "ZINBI", "ZIPIG")) {
                gam4$coefficients <- gamlss.fit$fit$argument[X2.d2 + 
                  X3.d2 + (1:X4.d2)]
                gam4$Vp <- gamlss.fit$magpp$Vb[X2.d2 + X3.d2 + 
                  (1:X4.d2), X2.d2 + X3.d2 + (1:X4.d2)]
                gam4$sig2 <- 1
                gam4$edf <- gamlss.fit$magpp$edf[X2.d2 + X3.d2 + 
                  (1:X4.d2)]
            }
        }
    }
    fit <- fit.SemiParSampleSel(funcop = funcop, start.v = start.v, 
        dat = respvec, rinit = rinit, rmax = rmax, VC = VC, qu.mag = qu.mag, 
        sp = sp, iterlimsp = iterlimsp, pr.tolsp = pr.tolsp, 
        univariate = FALSE)
    epsilon <- sqrt(.Machine$double.eps)
    He <- fit$fit$hessian
    He.eig <- eigen(He, symmetric = TRUE)
    if (min(He.eig$values) < epsilon) 
        He.eig$values[which(He.eig$values < epsilon)] <- 1e-07
    Vb <- He.eig$vectors %*% tcrossprod(diag(1/He.eig$values), 
        He.eig$vectors)
    if ((l.sp1 != 0 || l.sp2 != 0 || l.sp3 != 0 || l.sp4 != 0 || 
        l.sp5 != 0) && fp == FALSE) {
        HeSh <- He - fit$fit$S.h
        F <- Vb %*% HeSh
        F1 <- fit$magpp$edf1
        R <- fit$bs.mgfit$R
        Ve <- F %*% Vb
    }
    else {
        HeSh <- He
        Ve <- Vb
        F <- F1 <- diag(rep(1, dim(Vb)[1]))
        R <- fit$bs.mgfit$R
    }
    t.edf <- sum(diag(F))
    param <- fit$fit$argument[-c(1:X1.d2)]
    param <- param[1:X2.d2]
    X2s <- try(predict.gam(gam2, newdata = data, type = "lpmatrix"), 
        silent = TRUE)
    if (class(X2s) == "try-error") 
        stop("Check that the range of the covariate values of the selected sample is the same as that of the complete dataset.")
    eta2 <- X2s %*% param
    dimnames(fit$fit$hessian)[[1]] <- dimnames(fit$fit$hessian)[[2]] <- dimnames(Vb)[[1]] <- dimnames(Vb)[[2]] <- dimnames(HeSh)[[1]] <- dimnames(HeSh)[[2]] <- dimnames(F)[[1]] <- dimnames(F)[[2]] <- dimnames(He)[[1]] <- dimnames(He)[[2]] <- names(fit$fit$argument)
    if ((l.sp1 != 0 || l.sp2 != 0 || l.sp3 != 0 || l.sp4 != 0 || 
        l.sp5 != 0)) {
        edf <- edf1 <- list(0, 0, 0, 0, 0)
        for (i in 1:5) {
            if (i == 1) {
                mm <- l.sp1
                if (mm == 0) 
                  next
            }
            if (i == 2) {
                mm <- l.sp2
                if (mm == 0) 
                  next
            }
            if (i == 3) {
                mm <- l.sp3
                if (mm == 0) 
                  next
            }
            if (i == 4) {
                mm <- l.sp4
                if (mm == 0) 
                  next
            }
            if (i == 5) {
                mm <- l.sp5
                if (mm == 0) 
                  break
            }
            for (k in 1:mm) {
                if (i == 1) {
                  gam <- gam1
                  ind <- (gam$smooth[[k]]$first.para):(gam$smooth[[k]]$last.para)
                }
                if (i == 2) {
                  gam <- gam2
                  ind <- (gam$smooth[[k]]$first.para:gam$smooth[[k]]$last.para) + 
                    X1.d2
                }
                if (i == 3) {
                  gam <- gam3
                  ind <- (gam$smooth[[k]]$first.para:gam$smooth[[k]]$last.para) + 
                    X1.d2 + X2.d2
                }
                if (i == 4) {
                  gam <- gam4
                  ind <- (gam$smooth[[k]]$first.para:gam$smooth[[k]]$last.para) + 
                    X1.d2 + X2.d2 + X3.d2
                }
                if (i == 5) {
                  gam <- gam5
                  ind <- (gam$smooth[[k]]$first.para:gam$smooth[[k]]$last.para) + 
                    X1.d2 + X2.d2 + X3.d2 + X4.d2
                }
                edf[[i]][k] <- sum(diag(F)[ind])
                edf1[[i]][k] <- sum(F1[ind])
            }
        }
        if (l.sp1 != 0) 
            names(edf[[1]]) <- names(edf1[[1]]) <- names(gam1$sp)
        if (l.sp2 != 0) 
            names(edf[[2]]) <- names(edf1[[2]]) <- names(gam2$sp)
        if (l.sp3 != 0) 
            names(edf[[3]]) <- names(edf1[[3]]) <- names(gam3$sp)
        if (l.sp4 != 0) 
            names(edf[[4]]) <- names(edf1[[4]]) <- names(gam4$sp)
        if (l.sp5 != 0) 
            names(edf[[5]]) <- names(edf1[[5]]) <- names(gam5$sp)
    }
    KendTau <- theta <- sigma <- k <- phi <- nu <- NULL
    if (length(formula) < 3) {
        theta <- theta.tau(BivD = BivD, theta.star = fit$fit$argument["theta.star"])
        KendTau <- theta[, 2]
        theta <- theta[, 1]
        names(theta) <- names(KendTau) <- NULL
        if (margins[2] %in% c("N", "NB", "PIG", "D", "S", "BB", 
            "NBII", "WARING", "ZIBB", "ZABB", "ZANBI", "ZINBI", 
            "ZIPIG")) {
            sigma <- exp(fit$fit$argument["log.sigma"])
            names(sigma) <- k <- NULL
            phi <- sigma^2
        }
        if (margins[2] %in% c("ZABI", "ZIBI", "ZALG", "ZAP", 
            "ZIP", "ZIP2")) {
            sigma <- plogis(fit$fit$argument["logi.sigma"])
            names(sigma) <- k <- NULL
            phi <- sigma^2
        }
        if (margins[2] == "G") {
            sigma <- k <- exp(fit$fit$argument["log.shape"])
            names(k) <- names(sigma) <- NULL
            phi <- k^{
                -1
            }
        }
        if (margins[2] == "S") {
            nu <- fit$fit$argument["nu"]
            names(nu) <- NULL
        }
        if (margins[2] %in% c("D", "ZIBB", "ZABB", "ZANBI", "ZINBI", 
            "ZIPIG")) {
            nu <- plogis(fit$fit$argument["logi.nu"])
            names(nu) <- NULL
        }
        theta.a <- theta
        KendTau.a <- KendTau
        sigma.a <- sigma
        phi.a <- phi
        nu.a <- nu
    }
    if (length(formula) > 2) {
        sigma <- k <- phi <- nu <- NULL
        etatheta <- fit$fit$etatheta
        etasqv <- fit$fit$etasqv
        etak <- fit$fit$etak
        etak <- fit$fit$etak
        theta <- theta.tau(BivD = BivD, theta.star = fit$fit$etatheta)
        KendTau <- theta[, 2]
        theta <- as.matrix(theta[, 1])
        if (margins[2] %in% c("N", "NB", "PIG", "D", "S", "BB", 
            "NBII", "WARING", "ZIBB", "ZABB", "ZANBI", "ZINBI", 
            "ZIPIG")) {
            sigma <- exp(fit$fit$etasqv)
            phi <- sigma^2
        }
        if (margins[2] %in% c("ZABI", "ZIBI", "ZALG", "ZAP", 
            "ZIP", "ZIP2")) {
            sigma <- plogis(fit$fit$etasqv)
            phi <- sigma^2
        }
        if (margins[2] == "G") {
            sigma <- k <- exp(fit$fit$etak)
            phi <- k^{
                -1
            }
        }
        if (margins[2] == "S") 
            nu <- fit$fit$etanu
        if (margins[2] %in% c("D", "ZIBB", "ZABB", "ZANBI", "ZINBI", 
            "ZIPIG")) 
            nu <- plogis(fit$fit$etanu)
        theta.a <- mean(as.numeric(theta))
        KendTau.a <- mean(as.numeric(KendTau))
        sigma.a <- mean(as.numeric(sigma))
        phi.a <- mean(as.numeric(phi))
        nu.a <- mean(as.numeric(nu))
    }
    L <- list(fit = fit$fit, gam1 = gam1, gam2 = gam2, gam3 = gam3, 
        gam4 = gam4, gam5 = gam5, gam2.1 = gam2.1, coefficients = fit$fit$argument, 
        weights = weights, sp = fit$sp, iter.if = fit$iter.if, 
        iter.sp = fit$iter.sp, iter.inner = fit$iter.inner, start.v = start.v, 
        formula = formula, phi = phi, sigma = sigma, nu = nu, 
        theta = theta, tau = KendTau, n = n, n.sel = length(gam2$y), 
        bd = bd, X1 = X1, X2 = X2, X1.d2 = X1.d2, X2.d2 = X2.d2, 
        X3 = X3, X3.d2 = X3.d2, X4 = X4, X4.d2 = X4.d2, X5 = X5, 
        X5.d2 = X5.d2, l.sp1 = l.sp1, l.sp2 = l.sp2, l.sp3 = l.sp3, 
        l.sp4 = l.sp4, l.sp5 = l.sp5, He = He, HeSh = HeSh, Vb = Vb, 
        F = F, BivD = BivD, margins = margins, t.edf = t.edf, 
        bs.mgfit = fit$bs.mgfit, conv.sp = fit$conv.sp, wor.c = fit$wor.c, 
        eta1 = fit$fit$eta1, eta2 = eta2, y1 = y1, y2 = y2, logLik = -fit$fit$l, 
        fp = fp, gp1 = gp1, gp2 = gp2, gp3 = gp3, gp4 = gp4, 
        gp5 = gp5, X2s = X2s, magpp = fit$magpp, sigma.a = sigma.a, 
        phi.a = phi.a, theta.a = theta.a, nu.a = nu.a, tau.a = KendTau.a, 
        gamlss.fit = gamlss.fit, edf = edf, edf11 = edf1, R = R, 
        Ve = Ve)
    class(L) <- "SemiParSampleSel"
    L
}
<bytecode: 0xb92a600>
<environment: namespace:SemiParSampleSel>
 --- function search by body ---
Function SemiParSampleSel in namespace SemiParSampleSel has this body.
 ----------- END OF FAILURE REPORT -------------- 
Fatal error: the condition has length > 1