R/mhurdle.R

In mhurdle: Multiple Hurdle Tobit Models

#' Estimation of limited dependent variable models
#' 
#' mhurdle fits a large set of models relevant when the dependent variable is 0
#' for a part of the sample.
#' 
#' 
#' `mhurdle` fits models for which the dependent variable is zero for
#' a part of the sample. Null values of the dependent variable may
#' occurs because of one or several mechanisms : good rejection, lack
#' of ressources and purchase infrequency. The model is described
#' using a three-parts formula : the first part describes the
#' selection process if any, the second part the regression equation
#' and the third part the purchase infrequency process.  `y ~ 0 | x1 +
#' x2 | z1 + z2` means that there is no selection process.  `y ~ w1 +
#' w2 | x1 + x2 | 0` and `y ~ w1 + w2 | x1 + x2` describe the same
#' model with no purchase infrequency process. The second part is
#' mandatory, it explains the positive values of the dependant
#' variable. The `dist` argument indicates the distribution of the
#' error term. If `dist = "n"`, the error term is normal and (at least
#' part of) the zero observations are also explained by the second
#' part as the result of a corner solution. Several models described
#' in the litterature are obtained as special cases :
#' 
#' A model with a formula like `y~0|x1+x2` and `dist="n"` is the Tobit
#' model proposed by \insertCite{TOBIN/58}{mhurdle}.
#' 
#' `y~w1+w2|x1+x2` and `dist="l"` or `dist="t"` is the single hurdle
#' model proposed by \insertCite{CRAGG/71}{mhurdle}. With `dist="n"`,
#' the double hurdle model also proposed by
#' \insertCite{CRAGG/71}{mhurdle} is obtained. With `corr="h1"` we get
#' the correlated version of this model described by
#' \insertCite{BLUNDELL/87}{mhurdle}.
#' 
#' `y~0|x1+x2|z1+z2` is the P-Tobit model of
#' \insertCite{DEATO/IRISH/84}{mhurdle}, which can be a single hurdle
#' model if `dist="t"` or `dist="l"` or a double hurdle model if
#' `dist="n"`.
#' 
#' @name mhurdle
#' @aliases mhurdle
#' @param formula a symbolic description of the model to be fitted,
#' @param data a `data.frame`,
#' @param subset see [stats::lm()],
#' @param weights see [stats::lm()],
#' @param na.action see [stats::lm()],
#' @param start starting values,
#' @param dist the distribution of the error of the consumption
#'     equation: one of `"n"` (normal), `"ln"` (log-normal) `"bc"`
#'     (box-cox normal) and `"ihs"` (inverse hyperbolic sinus
#'     transformation),
#' @param h2 if `TRUE` the second hurdle is effective, it is not
#'     otherwise,
#' @param scaled if `TRUE`, the dependent variable is divided by its
#'     geometric mean,
#' @param corr a boolean indicating whether the errors of the
#'     different equations are correlated or not,
#' @param robust transformation of the structural parameters in order
#'     to avoid numerical problems,
#' @param check_gradient if `TRUE`, a matrix containing the analytical and
#'     the numerical gradient for the starting values are returned,
#' @param \dots further arguments.
#' @return
#' #' an object of class `c("mhurdle", "maxLik")`.
#' 
#' A `mhurdle` object has the following elements :
#' 
#' - coefficients: the vector of coefficients,
#' - vcov: the covariance matrix of the coefficients,
#' - fitted.values: a matrix of fitted.values, the first column being
#' the probability of 0 and the second one the mean values for the
#' positive observations,
#' - logLik: the log-likelihood,
#' - gradient: the gradient at convergence,
#' - model: a data.frame containing the variables used for the estimation,
#' - coef.names: a list containing the names of the coefficients in the
#' selection equation, the regression equation, the infrequency of purchase
#' equation and the other coefficients (the standard deviation of the error
#' term and the coefficient of correlation if `corr = TRUE`,
#' - formula: the model formula, an object of class `Formula`
#' - call: the call,
#' - rho: the lagrange multiplier test of no correlation.
#' 
#' @references
#'
#' \insertRef{BLUNDELL/87}{mhurdle}
#'
#' \insertRef{CRAGG/71}{mhurdle}
#'
#' \insertRef{DEATO/IRISH/84}{mhurdle}
#'
#' \insertRef{TOBIN/58}{mhurdle}
#'
#' @keywords regression
#' @examples
#' 
#' data("Interview", package = "mhurdle")
#' 
#' # independent double hurdle model
#' idhm <- mhurdle(vacations ~ car + size | linc + linc2 | 0, Interview,
#'               dist = "ln", h2 = TRUE, method = "bfgs")
#' 
#' # dependent double hurdle model
#' ddhm <- mhurdle(vacations ~ car + size | linc + linc2  | 0, Interview,
#'               dist = "ln", h2 = TRUE, method = "bfgs", corr = TRUE)
#' 
#' # a double hurdle p-tobit model
#' ptm <- mhurdle(vacations ~ 0 | linc + linc2 | car + size, Interview,
#'               dist = "ln", h2 = TRUE, method = "bfgs", corr = TRUE)
#' @importFrom Formula Formula
#' @importFrom survival survreg Surv
#' @importFrom truncreg truncreg
#' @importFrom stats binomial cor df.residual dnorm ecdf formula glm
#'     integrate lm lm.fit model.frame model.matrix model.response
#'     optimize pchisq pnorm printCoefmat qnorm quantile rnorm terms
#'     uniroot var .getXlevels deviance
#' @importFrom maxLik maxLik activePar
#' @export
mhurdle <- function(formula, data, subset, weights, na.action,
                    start = NULL, dist = c("ln", "n", "bc", "ihs"), h2 = FALSE,
                    scaled = TRUE, corr = FALSE, robust = TRUE,
                    check_gradient = FALSE, ...){
    fitted = TRUE
    dots <- list(...)
    oldoptions <- options(warn = -1)
    on.exit(options(oldoptions))
    cl <- match.call()
    posT <- as.list(cl) == "T" ; posF <- as.list(cl) == "F"
    cl[posT] <- TRUE           ; cl[posF] <- FALSE
    dist <- match.arg(dist)

    if (dist == "ln" & h2) dist <- "ln2"
    if (dist == "n" & ! h2) dist <- "tn"
    if (dist == "bc" & h2) dist <- "bc2"
    # 1. Compute the model.frame and the model.matrix

    if (! inherits(formula, "Formula")) formula <- Formula(formula)
    if (length(formula)[2] > 4) stop("at most 4 rhs should be provided in the formula")
    mf <- match.call(expand.dots = FALSE)
    mc <- match.call(expand.dots = TRUE)
    
    m <- match(c("formula", "data", "subset", "na.action", "weights"),
               names(mc), 0L)
    mf <- mc[c(1L, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1L]] <- as.name("model.frame")
    mf$formula <- formula
    mf <- eval(mf, parent.frame())
    .terms <- terms(mf)
    X1 <- model.matrix(formula, data = mf, rhs = 1)
    X2 <- model.matrix(formula, data = mf, rhs = 2)
    X3 <- model.matrix(formula, data = mf, rhs = 3)
    X4 <- model.matrix(formula, data = mf, rhs = 4)
    # Remove the intercept if any for the heteroscedastic model
    if (ncol(X4)){
        int4 <- attr(terms(formula, data = mf, rhs = 4), "intercept")
        X4 <- X4[, - 1, drop = FALSE]
    }
    y <- model.response(mf)
    if (scaled){
        geomean <- exp(mean(log(y[y > 0])))
        y <- y / geomean
        attr(y, "geomean") <- attr(mf, "geomean") <- geomean
    }
    N <- length(y)
    if (length(X1) == 0) X1 <- NULL
    if (length(X2) == 0) stop("the second hurdle (consumption equation) is mandatory")
    if (length(X3) == 0) X3 <- NULL
    if (length(X4) == 0) X4 <- NULL
    h1 <- ! is.null(X1)
    h3 <- ! is.null(X3)
    
    #  2. One equation models
    if (! h1 && ! h3 && is.null(X4)){
        result <- one_equation_model(X2, y, dist, check_gradient = check_gradient, start = start)
        if (check_gradient){
            result$call <- cl
            return(result)
        }
        result$naive <- NULL#naive
        result$call <- cl
        result$model <- mf
        result$formula <- formula
        result$terms <- .terms
        result$R2 <- c(null = NA, positive = NA)
        return(result)
    }
    
    # 3. Selection single hurdle models without correlation that can
    # be estimated simply in two parts using seperate.mhurdle()
    if (h1 && !h3  && ! corr && !h2){
        result <- two_parts_model(X1, X2, y, dist = dist, check_gradient = check_gradient, start = start)
        if (check_gradient){
            result$call <- cl
            return(result)
        }
        result$naive <- NULL#naive
        result$call <- cl
        result$model <- mf
        result$formula <- formula
        result$R2 <- c(null = NA, positive = NA)
        result$terms <- .terms
        return(result)
    }

    # 5. Compute the starting values if not provided (use the linear
    # specification as the starting value for ihs and the log-linear
    # specification for Box-Cox)

    if (is.null(start)){
        dist.start <- dist
        if (dist %in% c("bc", "bc2", "ln2")) dist.start <- "ln"
        if (dist == "ihs") dist.start <- "n"
        start <- start.mhurdle(X1, X2, X3, y, dist.start)
        # in case of heteroscedasctic model, add K4 zeros to the start
        # vector and the intercept should be ln(sigma_o) (not sigma_o)
        # because of the exp form
        sd.pos <- getindex(X1, X2, X3, X4, corr, dist, which = "sd")
        if (robust) start[sd.pos] <- log(start[sd.pos])
        if (! is.null(X4)) start <- c(start[1:sd.pos], rep(0, ncol(X4)))
        
        # add shape and/or scale parameters
        if (corr){
            if (robust) rhoinit <- tan(0.1 * pi / 2) else rhoinit <- 0.1
            if (h1 + h3 == 2) start <- c(start, rho12 = rhoinit, rho13 = rhoinit, rho23 = rhoinit)
            else start <- c(start, rho = rhoinit)
        }
        if (dist %in% c("bc", "bc2", "ihs")) start <- c(start, tr =  0.1)
        if (dist %in% c("bc2", "ln2")){
            pos = 1
            if (robust) pos <- max(log(pos), -100)
            start <- c(start, pos = pos)
        }
    }
    else{
        if (robust){
            sd.pos <- getindex(X1, X2, X3, X4, corr, dist, which = "sd")
            if (dist %in% c("ln2", "bc2")) mu.pos <- getindex(X1, X2, X3, X4, corr, dist, which = "pos")
            rho.pos <- getindex(X1, X2, X3, X4, corr, dist, which = "corr")
            start[sd.pos] <- log(start[sd.pos])
            if (dist %in% c("ln2", "bc2")) start[mu.pos] <- max(log(start[mu.pos]), - 50)
            if (corr) start[rho.pos] <- tan(start[rho.pos] * pi / 2)
        }
    }

    result <- mhurdle.fit(start, X1, X2, X3, X4, y, gradient = TRUE,
                          dist = dist, corr = corr, robust = robust,
                          fitted = fitted,
                          check_gradient = check_gradient,
                          ...)
    if (check_gradient){
        result$call <- cl
        return(result)
    }
    
    # 3. Compute the naive model

    Pnull <- mean(y == 0)
    dist.naive <- dist
    if (dist %in% c("ihs")) dist.naive <- "n"
    if (dist %in% c("bc", "bc2", "ln2")) dist.naive <- "ln"
    if (dist.naive != "ln"){
        Ec <- mean(y[y > 0])
        Vc <- var(y[y > 0])}
    else{
        Ec <- mean(log(y[y > 0]))
        Vc <- var(log(y[y > 0]))
    }

    start.naive <- c(rep(0.1, 1 + h1 + h3), 1)
    moments <- c(Pnull, Ec, Vc)
    naive <- maxLik::maxLik(lnl.naive, start = start.naive,
                    dist = dist.naive, moments = moments,
                    h1 = h1, h3 = h3)
    coef.naive <- naive$est
    parts <- attr(naive$max, "parts") * c(Pnull, 1 - Pnull, 1 - Pnull) * N
    logLik.naive <- structure(as.numeric(naive$max) * N, nobs = length(y),
                              parts = parts, 
                              df = length(coef.naive), class = "logLik")
    naive <- list(coefficients = coef.naive, logLik = logLik.naive, code = naive$code)

    result$naive <- naive
    result$call <- cl
    result$formula <- formula
    result$model <- mf
    result$terms <- .terms
    result$xlevels = .getXlevels(.terms, mf)
    
    lnL1c <- attr( naive$logLik, "parts")[1] + attr( naive$logLik, "parts")[2]
    lnL1u <- attr(result$logLik, "parts")[1] + attr(result$logLik, "parts")[2]

    lnL2c <- attr( naive$logLik, "parts")[3] - attr( naive$logLik, "parts")[2]
    lnL2u <- attr(result$logLik, "parts")[3] - attr(result$logLik, "parts")[2]

    R1 <- 1 - (lnL1c / lnL1u) ^ (2 / N * lnL1c)
    R2old <- 1 - (exp(lnL2c) / exp(lnL2u)) ^ (2 / (N * (1 - Pnull)))
    R2 <- 1 - exp(- 2 * (lnL2u - lnL2c) / (N * (1 - Pnull)))
    result$R2 <- c(null = unname(R1), positive = unname(R2), old = unname(R2old))
    result$dpar <- - sum(diag(solve(result$hessian, crossprod(result$gradient))))
    result
}

mhurdle.fit <- function(start, X1, X2, X3, X4, y, gradient = FALSE, fit = FALSE,
                        dist = c("ln", "n", "tn", "bc", "ihs", "bc2", "ln2"),
                        corr = FALSE, robust = TRUE,  fitted = FALSE,
                        check_gradient = FALSE, ...){
    start.time <- proc.time()
    h1 <- ! is.null(X1)
    h3 <- ! is.null(X3)
    h4 <- ! is.null(X4)
    
    # fancy coefficients names
    sd.names <- "sd"
    
    if (corr){
        if (h1 & h3) rho.names <- c("corr12", "corr13", "corr23")
        else rho.names <- ifelse(h1, "corr12", "corr23")
    }
    else rho.names <- NULL
    if (dist %in% c("bc", "bc2", "ihs")) tr.names <- "tr" else tr.names <- NULL
    if (dist %in% c("ln2", "bc2")) mu.names <- "pos" else mu.names <- NULL

    coef.names <- list(h1   = colnames(X1),
                       h2   = colnames(X2),
                       h3   = colnames(X3),
                       sd   = sd.names,
                       h4   = colnames(X4),
                       corr = rho.names,
                       tr   = tr.names,
                       pos   = mu.names)

    start.names <- coef.names
    if (h1) start.names$h1 <- paste("h1", start.names$h1, sep = ".")
    start.names$h2 <- paste("h2", start.names$h2, sep = ".")
    if (h3) start.names$h3 <- paste("h3", start.names$h3, sep = ".")
    if (h4) start.names$h4 <- paste("h4", start.names$h4, sep = ".")
    names(start) <- Reduce("c", start.names)
    f <- function(param) mhurdle.lnl(param, X1 = X1, X2 = X2, X3 = X3, X4 = X4, y = y,
                                     gradient = TRUE, fitted = FALSE,
                                     dist = dist, corr = corr,
                                     robust = robust)

    if (check_gradient) return(compare_gradient(f, start))

    maxl <- maxLik::maxLik(f, start = start, control = list(lambdatol = 1E-20),  ...)
    nb.iter <- maxl$iterations
    convergence.OK <- maxl$code <= 2
    coefficients <- maxl$estimate

    if (fitted) fitted.values <- attr(mhurdle.lnl(coefficients, X1 = X1, X2 = X2, X3 = X3, X4 = X4, y = y,
                                                  gradient = FALSE, fitted = TRUE, robust = robust,
                                                  dist = dist, corr = corr), "fitted")
    else fitted.values <- NULL

    # contribution of every single observation to the likelihood and
    # its gradient (as an attribute)
    logLik <- f(coefficients)
    gradi <- attr(logLik, "gradi")
    logLik <- structure(as.numeric(logLik), df = length(coefficients),
                        parts = attr(logLik, "parts"),
                        nobs = length(y), class = "logLik")
    hessian <- maxl$hessian
    ev <- eigen(- hessian)$values
    if (any(ev < 0)) cat("The hessian is not negative definite\n")
    else (if (any(abs(ev) < 1E-07)) cat("the hessian is singular\n"))
    opg <- crossprod(gradi)
    elaps.time <- proc.time() - start.time
    
#    eps <- with(maxl, gradient %*% solve(- hessian) %*% gradient)
    eps <- with(maxl, gradient %*% solve(opg) %*% gradient)
    est.stat <- list(elaps.time = elaps.time,
                     nb.iter = nb.iter,
                     eps = eps,
                     method = maxl$type,
                     message = maxl$message
                     )
    class(est.stat) <- "est.stat"
    gtheta <- rep(1, length(coefficients))
    if (robust){
        if (corr){
            poscor <- sub.mhurdle(coef.names, "corr")
            gtheta[poscor] <- 2 / pi / (1 + coefficients[poscor] ^ 2)
            coefficients[poscor] <- atan(coefficients[poscor]) * 2 / pi
        }
        if (dist %in% c("bc2", "ln2")){
            posmu <- sub.mhurdle(coef.names, "pos")
            gtheta[posmu] <- exp(coefficients[posmu])
            coefficients[posmu] <- exp(coefficients[posmu])
        }
        possd <- sub.mhurdle(coef.names, "sd")
        gtheta[possd] <- exp(coefficients[possd])
        coefficients[possd] <- exp(coefficients[possd])
    }
    result <- list(coefficients  = coefficients,
                   fitted.values = fitted.values,
                   logLik        = logLik,
                   gradient      = gradi,
                   hessian       = hessian,
                   formula       = NULL,
                   terms         = NULL,
                   xlevels       = NULL,
                   model         = NULL,
                   coef.names    = coef.names,
                   call          = NULL,
                   est.stat      = est.stat,
                   naive         = NULL
                   )
    class(result) <- c("mhurdle", class(result))
    result
    
}


sanitize <- function(x, myeps = 1E-07, mymax = 1E02, string = c("", ""), replace = TRUE, verbal = TRUE){
    string <- paste("of", string[1], "in", string[2])
    if (replace){
        if (any(is.na(x))){
            if (verbal) cat(paste(sum(is.na(x)), "NA values", string, "replaced by 0\n"))
            x[is.na(x)] <- 0
        }
        if ( any(x > 0 & x < myeps)){
            if (verbal) cat(paste(sum(x > 0 & x < myeps), "values", string, "lower than",  myeps,"replaced by", myeps, "\n"))
            x[x > 0 & x < myeps] <- myeps
        }
        if (any(x < - mymax)){
            if (verbal) cat(paste(sum(x < - mymax), "values", string, "lower than", - mymax, "replaced by", - mymax, "\n"))
            x[x < - mymax] <- - mymax
        }
        if (any(x > mymax)){
            if (verbal) cat(paste(sum(x >  mymax), "values", string, "greater than",  mymax, "replaced by", mymax, "\n"))
            x[x > mymax] <- mymax
        }
    }
    x
}