R/mmlt.R

In tram: Transformation Models

### catch constraint violations here
.log <- function(x) {
  ret <- log(pmax(.Machine$double.eps, x))
  dim(ret) <- dim(x)
  ret
}

.rbind <- function(x) {
    if (!is.list(x)) return(matrix(x, nrow = 1))
    return(do.call("rbind", x))
}

.ll <- function(dim, standardize = TRUE, args = list()) {

    if (length(dim) == 1L)
        dim <- c(dim, 0L)

    cJ <- dim[1L]
    dJ <- dim[2L]

    if (!dJ) {

        cll <- function(obs, Lambda) {

            if (dim(Lambda)[2L] > 1)
                stopifnot(!attr(Lambda, "diag"))

            if (!standardize)
                return(ldmvnorm(obs = obs, invchol = Lambda, logLik = FALSE))

            sLambda <- mvtnorm::standardize(invchol = Lambda)
            return(ldmvnorm(obs = obs, invchol = sLambda, logLik = FALSE))
        }

        csc <- function(obs, Lambda) {

            if (dim(Lambda)[2L] > 1)
                stopifnot(!attr(Lambda, "diag"))

            if (!standardize) {
                ret <- sldmvnorm(obs = obs, invchol = Lambda)
                return(list(Lambda = ret$invchol, obs = ret$obs))
            }

           chol <- solve(Lambda)
           ### START readable:
           # schol <- mvtnorm::standardize(chol = chol)
           # ret <- sldmvnorm(obs = obs, chol = schol)
           ### avoid calling solve() multiple times
           D <- sqrt(Tcrossprod(chol, diag_only = TRUE))
           sLambda <- invcholD(Lambda, D = D)
           ret <- sldmvnorm(obs = obs, invchol = sLambda)
           ret$chol <- -vectrick(sLambda, ret$invchol)
           ### END
           dobs <- ret$obs
           ret <- destandardize(chol = chol, invchol = Lambda, 
                                score_schol = ret$chol)
           return(list(Lambda = ret, obs = dobs))
       }

       return(list(logLik = cll, score = csc))
    }

    ll <- function(obs = NULL, lower, upper, Lambda) {

        if (dim(Lambda)[2L] > 1)
            stopifnot(!attr(Lambda, "diag"))

        a <- args
        a$obs <- obs
        a$mean <- 0
        a$lower <- lower
        a$upper <- upper
        a$logLik <- FALSE
        if (!standardize) {
            a$invchol <- Lambda
        } else {
            a$chol <- mvtnorm::standardize(chol = solve(Lambda))
        }
        return(do.call("ldpmvnorm", a))
    }

    sc <- function(obs = NULL, lower, upper, Lambda) {

        a <- args
        a$obs <- obs
        a$mean <- 0
        a$lower <- lower
        a$upper <- upper
        a$logLik <- TRUE
        if (!standardize) {
            a$invchol <- Lambda
            ret <- do.call("sldpmvnorm", a)
            return(list(Lambda = ret$invchol,
                        obs = ret$obs,
                        mean = ret$mean, 
                        lower = ret$lower, 
                        upper = ret$upper))
        }

        chol <- solve(Lambda)
        ### START readable:
        # a$chol <- mvtnorm::standardize(chol = chol)
        # ret <- do.call("sldpmvnorm", a)
        ### avoid calling solve() multiple times
        D <- sqrt(Tcrossprod(chol, diag_only = TRUE))
        a$invchol <- sLambda <- invcholD(Lambda, D = D)
        ret <- do.call("sldpmvnorm", a)
        ret$chol <- -vectrick(sLambda, ret$invchol)
        ### END
        smean <- ret$mean
        sobs <- ret$obs
        slower <- ret$lower
        supper <- ret$upper
        ret <- destandardize(chol = chol, invchol = Lambda, 
                             score_schol = ret$chol)
        ret <- list(Lambda = ret, mean = smean, obs = sobs, 
                    lower = slower, upper = supper)
        return(ret)
    }

    return(list(logLik = ll, score = sc))
}

.models <- function(...) {

    m <- lapply(list(...), function(x) as.mlt(x))
    nm <- abbreviate(sapply(m, function(x) x$model$response), 4)
    J <- length(m)
    Jp <- J * (J - 1) / 2
    normal <- sapply(m, function(x) x$todistr$name == "normal")
  
    w <- lapply(m, weights)
    out <- lapply(w, function(x) stopifnot(isTRUE(all.equal(x, w[[1]]))))
    w <- w[[1L]]
    if (isTRUE(all.equal(unique(w), 1))) w <- FALSE
  
    mm <- lapply(m, function(mod) {
      eY <- get("eY", environment(mod$parm))
      iY <- get("iY", environment(mod$parm))
      list(eY = eY, iY = iY)
    })

    cmod <- sapply(mm, function(x) !is.null(x$eY))  
    dmod <- sapply(mm, function(x) !is.null(x$iY))  
    stopifnot(all(xor(cmod, dmod)))
    ### continuous models first
    stopifnot(all(diff(cmod) <= 0))
    stopifnot(all(diff(dmod) >= 0))

    nobs <- unique(sapply(m, nobs))
    stopifnot(length(nobs) == 1L)
    nobs <- nobs[[1L]]

    P <- sapply(m, function(x) length(coef(x)))
    fpar <- factor(rep(1:J, P))

    parm <- function(par) {
        mpar <- par[1:sum(P)]
        split(mpar, fpar)
    }

    constr <- lapply(mm, function(m) {
        if (is.null(m$eY)) return(attr(m$iY$Yleft, "constraint"))
        return(attr(m$eY$Y, "constraint"))
    })

    ui <- do.call("bdiag", lapply(constr, function(x) x$ui))
    ci <- do.call("c", lapply(constr, function(x) x$ci))
    ui <- as(ui[is.finite(ci),,drop = FALSE], "matrix")
    ci <- ci[is.finite(ci)]

    mf <- lapply(1:J, function(j) {
        mf <- m[[j]]$data ###model.frame(m[[j]])
        if (cmod[j]) return(mf)
        yl <- m[[j]]$response$cleft
        yr <- m[[j]]$response$cright
        rp <- m[[j]]$model$response
        ml <- mr <- mf
        ml[[rp]] <- yl
        mr[[rp]] <- yr
        return(list(left = ml, right = mr))
    })

    ### needs newdata for predict
    nn <- sapply(1:J, function(j) {
        !is.null(m[[j]]$fixed) ||
        !isTRUE(all.equal(unique(m[[j]]$offset), 0)) ||
        m[[j]]$model$scale_shift
    })

    type <- lapply(1:J, function(j)
        mlt:::.type_of_response(m[[j]]$response))

    return(list(models = m, mf = mf, cont = cmod, type = type, normal = normal, 
                nobs = nobs, weights = w, nparm = P, parm = parm, 
                ui = ui, ci = ci, mm = mm, names = nm, nn = nn))
}

.model_matrix <- function(models, j = 1, newdata = NULL, prime = FALSE) {

    if (is.null(newdata)) {
        if (models$cont[j]) {
            if (prime) return(models$mm[[j]]$eY$Yprime)
            return(models$mm[[j]]$eY$Y)
        }
        stopifnot(!prime)
        Yleft <- models$mm[[j]]$iY$Yleft
        Yright <- models$mm[[j]]$iY$Yright
        return(list(Yleft = Yleft, Yright = Yright))
    }

    resp <- models$models[[j]]$model$response
    y <- R(newdata[[resp]])
    if (models$cont[j] || mlt:::.type_of_response(y) == "double") {
        if (prime) {
            drv <- 1L
            names(drv) <- models$models[[j]]$model$response
            return(model.matrix(models$models[[j]]$model, data = newdata, deriv = drv))
        } else {
            return(model.matrix(models$models[[j]]$model, data = newdata))
        }
    }
    iY <- mlt:::.mm_interval(model = models$models[[j]]$model, data = newdata, resp, y)
    Yleft <- iY$Yleft
    Yright <- iY$Yright
    return(list(Yleft = Yleft, Yright = Yright))
}

.mget <- function(models, j = 1, parm, newdata = NULL,
                  what = c("trafo", "dtrafo", "z", "zleft", 
                           "dzleft", "zright", "dzright", "zprime", 
                           "mm", "mmprime", "estfun", "scale"), ...) {

    what <- match.arg(what)

    if (length(j) > 1) {
        ret <- lapply(j, .mget, models = models, parm = parm, 
                      newdata = newdata, what = what)
        return(ret)
    }

    if (what == "scale") {
        if (models$cont[j]) {
            Y <- models$mm[[j]]$eY$Y
        } else {
            Y <- models$mm[[j]]$iY$Yleft
        }
        Ytmp <- Y
        Ytmp[!is.finite(Ytmp)] <- NA
        sc <- apply(abs(Ytmp), 2, max, na.rm = TRUE)
        lt1 <- sc < 1.1
        gt1 <- sc >= 1.1
        sc[gt1] <- 1 / sc[gt1]
        sc[lt1] <- 1
        return(sc)
    }

    prm <- models$parm(parm)[[j]]
    tmp <- models$models[[j]]
    cf <- coef(tmp)
    cf[] <- prm
    coef(tmp) <- cf

    ### check for fixed, offset, and shift_scale; go through predict if this is the case
    ### (slow but works)
    if (is.null(newdata)) {
        if (models$nn[j]) newdata <- tmp$data
    }
    if (models$cont[j]) {
        if (is.null(newdata)) {
            tr <- c(models$mm[[j]]$eY$Y %*% prm)
            trp <- c(models$mm[[j]]$eY$Yprime %*% prm)
            if (!models$normal[j])
                trd <- tmp$todistr$d(tr) * trp
        } else {
            tr <- predict(tmp, newdata = newdata, type = "trafo", ...)
            drv <- 1L
            names(drv) <- tmp$model$response
            trp <- predict(tmp, newdata = newdata, type = "trafo", 
                           deriv = drv, ...)
            if (!models$normal[j])
                trd <- predict(tmp, newdata = newdata, type = "density", ...)
        }
    } else {
        if (is.null(newdata)) {
            trl <- c(models$mm[[j]]$iY$Yleft %*% prm)
            trl[!is.finite(trl)] <- -Inf
            trr <- c(models$mm[[j]]$iY$Yright %*% prm)
            trr[!is.finite(trr)] <- Inf
        } else {
            mmj <- .model_matrix(models, j = j, newdata = newdata)
            ### continuous response for model with censoring
            if (is.matrix(mmj)) return(c(mmj %*% prm))
            trl <- c(mmj$Yleft %*% prm)
            trl[!is.finite(trl)] <- -Inf
            trr <- c(mmj$Yright %*% prm)
            trr[!is.finite(trr)] <- Inf
        }
    }

    if (what == "trafo") {
#        stopifnot(models$cont[j])
        return(tr)
    }
    if (what == "dtrafo") {
#        stopifnot(models$cont[j])
        return(tmp$todistr$d(tr))
    }
    if (what == "z") {
#        stopifnot(models$cont[j])
        if (models$normal[j]) 
            return(tr)
        return(qnorm(tmp$todistr$p(tr, log = TRUE), log.p = TRUE))
    }
    if (what == "zleft") {
#        stopifnot(!models$cont[j])
        if (models$normal[[j]])
            return(trl)
        return(qnorm(tmp$todistr$p(trl, log = TRUE), log.p = TRUE))
    }
    if (what == "dzleft") {
#        stopifnot(!models$cont[j])
        if (models$normal[[j]])
            return(rep(1, length(trl)))
        qn <- qnorm(tmp$todistr$p(trl, log = TRUE), log.p = TRUE)
        dn <- dnorm(qn)
        dn[!is.finite(dn)] <- 1
        return(tmp$todistr$d(trl) / dn)
    }
   if (what == "zright") {
#        stopifnot(!models$cont[j])
        if (models$normal[[j]])
            return(trr)
        return(qnorm(tmp$todistr$p(trr, log = TRUE), log.p = TRUE))
    }
    if (what == "dzright") {
#        stopifnot(!models$cont[j])
        if (models$normal[[j]])
            return(rep(1, length(trr)))
        qn <- qnorm(tmp$todistr$p(trr, log = TRUE), log.p = TRUE)
        dn <- dnorm(qn)
        dn[!is.finite(dn)] <- 1
        return(tmp$todistr$d(trr) / dn)
    }
    if (what == "zprime") {
#        stopifnot(models$cont[j])
        if (models$normal[[j]])
            return(trp)
        qn <- qnorm(tmp$todistr$p(tr, log = TRUE), log.p = TRUE)
        return(trd / dnorm(qn))
    }
    if (what == "estfun") {
        if (is.null(newdata))
            return(estfun(tmp))
        return(estfun(tmp, newdata = newdata))
    }
}

.MCw <- function(J, M, seed) {

    ### from stats:::simulate.lm
    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
        runif(1)
    if (is.null(seed)) 
        RNGstate <- get(".Random.seed", envir = .GlobalEnv)
    else {
        R.seed <- get(".Random.seed", envir = .GlobalEnv)
        set.seed(seed)
        RNGstate <- structure(seed, kind = as.list(RNGkind()))
        on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
    }

    return(matrix(runif((J - 1) * M), ncol = M))
}

.start <- function(m, xnames, fixed = NULL) {

    J <- length(m$models)
    Jp <- J * (J - 1) / 2
    Jnames <- m$names
    margin_par <- do.call("c", lapply(m$models, function(mod) coef(as.mlt(mod))))
    names(margin_par) <- paste(rep(Jnames, time = m$nparm), names(margin_par), sep = ".")

    rn <- rownames(unclass(ltMatrices(1:Jp, names = Jnames, byrow = TRUE)))
    lnames <- paste(rep(rn, each = length(xnames)),
                    rep(xnames, length(rn)), sep = ".")
    lambda_par <- rep(0, length(lnames))
    names(lambda_par) <- lnames

    start <- c(margin_par, lambda_par)

    if (!is.null(fixed))
        stopifnot(all(fixed %in% names(start)))

    return(start)
}

mmltoptim <- function(auglag = list(maxtry = 5), ...)
    mltoptim(auglag = auglag, ...)

mmlt <- function(..., formula = ~ 1, data, conditional = FALSE, 
                 theta = NULL, fixed = NULL, scale = FALSE,
                 optim = mmltoptim(), args = list(seed = 1, M = 1000), 
                 dofit = TRUE, domargins = TRUE)
{
  
    call <- match.call()

    m <- .models(...)

    if (conditional && !all(m$normal))
        stop("Conditional models only available for marginal probit-type models.")

    if (conditional && !domargins)
        stop("Conditional models must fit marginal and joint parameters.")

    ### compute starting values for lambda
    if (is.null(theta) && dofit && domargins) {
        cl <- match.call()
        cl$conditional <- FALSE
        cl$domargins <- FALSE
        sm <- eval(cl, parent.frame())
        if (!is.null(sm)) {
            theta <- coef(sm, type = "all")
            if (conditional) {
                ### theta are conditional parameters, scale with sigma
                class(sm)[1] <- "cmmlt" ### do NOT standardize Lambda
                d <- rowMeans(diagonals(coef(sm, newdata = data, type = "Sigma")))
                theta[1:sum(m$nparm)] <- theta[1:sum(m$nparm)] * rep(sqrt(d), times = m$nparm)
            }
        } else {
            theta <- do.call("c", lapply(m$models, function(mod) coef(mod)))
        }
    } 

    cJ <- sum(m$cont)
    dJ <- sum(!m$cont)
    J <- cJ + dJ
    Jp <- J * (J - 1) / 2
    llsc <- .ll(c(cJ, dJ), standardize = !conditional, args)

    if (dJ && is.null(args$w))
        args$w <- .MCw(J = dJ, M = args$M, seed = args$seed)

    if (isTRUE(all.equal(formula, ~ 1))) {
        lX <- matrix(1)
        colnames(lX) <- "(Intercept)"
        bx <- NULL
    } else {
        bx <- formula
        if (inherits(formula, "formula")) {
            bx <- as.basis(formula, data)
        } 
        lX <- model.matrix(bx, data = data)
        if (conditional)
            warning("Conditional models with covariate-dependent correlations are order-dependent")
    }
    .Xparm <- function(parm) {
        parm <- parm[-(1:sum(m$nparm))]
        return(matrix(parm, nrow = ncol(lX)))
    }

    start <- .start(m, colnames(lX), names(fixed))
    parnames <- eparnames <- names(start)
    lparnames <- names(start)[-(1:sum(m$nparm))]
    if (!is.null(fixed)) eparnames <- eparnames[!eparnames %in% names(fixed)]
    if (!is.null(fixed)) lparnames <- lparnames[!lparnames %in% names(fixed)]

    if (!is.null(theta)) {
        if (!is.null(fixed)) theta <- theta[!names(theta) %in% names(fixed)]
        stopifnot(length(theta) == length(eparnames))
        names(theta) <- eparnames
    }


    if (cJ) {
        mm <- lapply(1:cJ, function(j) .model_matrix(m, j = j))
        mmp <- lapply(1:cJ, function(j) .model_matrix(m, j = j, prime = TRUE))
    }
    if (dJ) {
        dmm <- lapply(cJ + 1:dJ, function(j) .model_matrix(m, j = j))
        dmm <- lapply(dmm, function(x) {
            x$Yleft[!is.finite(x$Yleft[,1]),] <- 0
            x$Yright[!is.finite(x$Yright[,1]),] <- 0
            x
        })
    }


    ### note: estfun() already has weights already multiplied to scores
    weights <- m$weights
    if (weights) {
        if (cJ) {
            mm <- lapply(mm, function(x) x * weights)
            mmp <- lapply(mmp, function(x) x * weights)
        }
        if (dJ) {
            dmm <- lapply(dmm, function(x) list(Yleft = x$Yleft * weights,
                                                Yright = x$Yright * weights))
        }
    }

    LAMBDA <- ltMatrices(matrix(0, nrow = Jp, ncol = nrow(lX)),
                         byrow = TRUE, diag = FALSE, names = names(m$models))

    ll <- function(parm, newdata = NULL) {

        if (!is.null(newdata) && !isTRUE(all.equal(formula, ~ 1))) 
            lX <- model.matrix(bx, data = newdata)

        # Lambda <- ltMatrices(t(lX %*% .Xparm(parm)), byrow = TRUE, diag = FALSE, 
        #                      names = names(m$models))
        # saves time in ltMatrices
        Lambda <- LAMBDA
        Lambda[] <- t(lX %*% .Xparm(parm))
        ret <- 0
        if (cJ) {
            z <- .rbind(.mget(m, j = which(m$cont), parm = parm, what = "z", newdata = newdata))
            zp <- .rbind(.mget(m, j = which(m$cont), parm = parm, what = "zprime", newdata = newdata))
            ret <- colSums(.log(zp))
            if (!dJ) return(ret + llsc$logLik(obs = z, Lambda = Lambda))
        }
        if (dJ) {
            lower <- .rbind(.mget(m, j = which(!m$cont), parm = parm, what = "zleft", newdata = newdata))
            upper <- .rbind(.mget(m, j = which(!m$cont), parm = parm, what = "zright", newdata = newdata))
            if (!cJ)
                return(llsc$logLik(lower = lower, upper = upper, Lambda = Lambda))
        }
        return(ret + llsc$logLik(obs = z, lower = lower, upper = upper, Lambda = Lambda))
    }

    sc <- function(parm, newdata = NULL, scores = FALSE) {

        if (scores) {
            RS <- CS <- function(x) x
        } else {
            RS <- function(x) rowSums(x, na.rm = TRUE)
            CS <- function(x) colSums(x, na.rm = TRUE)
        }

        if (!is.null(newdata) && !isTRUE(all.equal(formula, ~ 1))) 
            lX <- model.matrix(bx, data = newdata)

        # Lambda <- ltMatrices(t(lX %*% .Xparm(parm)), byrow = TRUE, diag = FALSE, 
        #                      names = names(m$models))
        # saves time in ltMatrices
        Lambda <- LAMBDA
        Lambda[] <- t(lX %*% .Xparm(parm))

        if (cJ) {
            z <- .rbind(.mget(m, j = which(m$cont), parm = parm, what = "z", newdata = newdata))
            if (!dJ)
                sc <- llsc$score(obs = z, Lambda = Lambda)
        }
        if (dJ) {
            lower <- .rbind(.mget(m, j = which(!m$cont), parm = parm, what = "zleft", newdata = newdata))
            upper <- .rbind(.mget(m, j = which(!m$cont), parm = parm, what = "zright", newdata = newdata))
            if (!cJ)
                sc <- llsc$score(lower = lower, upper = upper, Lambda = Lambda)
        }
        if (cJ && dJ)
            sc <- llsc$score(obs = z, lower = lower, upper = upper, Lambda = Lambda)

        Lmat <- Lower_tri(sc$Lambda)[rep(1:Jp, each = ncol(lX)), , drop = FALSE]
        if (identical(c(lX), 1)) {
            scL <- RS(Lmat) ### NaN might appear in scores
        } else {
            scL <- RS(Lmat * t(lX[,rep(1:ncol(lX), Jp), drop = FALSE]))
        }
      
        scp <- vector(mode = "list", length = cJ + dJ)

        if (cJ) {
            if (all(m$normal)) {
                zp <- .rbind(.mget(m, j = which(m$cont), parm = parm, what = "zprime", newdata = newdata))
                scp[1:cJ] <- lapply(1:cJ, function(j) {
                    CS(mm[[j]] * c(sc$obs[j,])) + CS(mmp[[j]] / c(zp[j,]))
                })
            } else {
                dz <- .rbind(.mget(m, j = which(m$cont), parm = parm, what = "dtrafo", newdata = newdata))
                ef <- lapply(which(m$cont), function(j) .mget(m, j = j, parm = parm, what = "estfun", newdata = newdata))
                scp[1:cJ] <- lapply(1:cJ, function(j) {
                    CS(mm[[j]] * c(sc$obs[j,] + z[j,]) / c(dnorm(z[j,])) * c(dz[j,])) - CS(ef[[j]])
                })
            }
        }

        if (dJ) {
            if (all(m$normal)) {
                scp[cJ + 1:dJ] <- lapply(1:dJ, function(j) {
                    CS(dmm[[j]]$Yleft * c(sc$lower[j,])) +
                    CS(dmm[[j]]$Yright * c(sc$upper[j,]))
                })
            } else {
                dzl <- .rbind(.mget(m, j = which(!m$cont), parm = parm, what = "dzleft", newdata = newdata))
                dzl[!is.finite(dzl)] <- 0
                dzr <- .rbind(.mget(m, j = which(!m$cont), parm = parm, what = "dzright", newdata = newdata))
                dzr[!is.finite(dzr)] <- 0
                scp[cJ + 1:dJ] <- lapply(1:dJ, function(j) {
                    return(CS(dmm[[j]]$Yleft * c(dzl[j,]) * c(sc$lower[j,])) +
                           CS(dmm[[j]]$Yright * c(dzr[j,]) * c(sc$upper[j,])))
                })
            }
        }
        
        if (!scores) {
            ret <- c(do.call("c", scp), c(scL))
            names(ret) <- parnames
            return(ret)
        }
        ret <- cbind(do.call("cbind", scp), t(scL))
        colnames(ret) <- parnames
        return(ret)
    }

    ### scale
    if (scale) {
        scl <- rep(apply(abs(lX), 2, max, na.rm = TRUE), times = Jp)
        lt1 <- scl < 1.1
        gt1 <- scl >= 1.1
        scl[gt1] <- 1 / scl[gt1]
        scl[lt1] <- 1
        scl <- c(do.call("c", .mget(m, j = 1:J, parm = NULL, what = "scale")), scl)
        names(scl) <- parnames
    } else {
        scl <- numeric(length(parnames))
        scl[] <- 1
        names(scl) <- parnames
    }

    if (!weights) 
        weights <- 1

    f <- function(par, scl, ...) {
        if (!is.null(fixed)) {
            p <- par
            names(p) <- eparnames
            p[names(fixed)] <- fixed
            par <- p[parnames]
        }
        return(-sum(weights * ll(par * scl, ...)))
    }

    ### note: x * weights was already computed
    g <- function(par, scl, ...) {
        if (!is.null(fixed)) {
            p <- par
            names(p) <- eparnames
            p[names(fixed)] <- fixed
            par <- p[parnames]
        }
        ret <- -sc(par * scl, ...) * scl
        if (is.null(fixed)) return(ret)
        if (is.matrix(ret))
            return(ret[, !parnames %in% names(fixed)])
        return(ret[!parnames %in% names(fixed)])
    }

    if (!domargins) {

        stopifnot(!conditional)

        start <- start[1:sum(m$nparm)]

        cll <- function(cpar) f(c(start / scl[names(start)], cpar), scl = scl)
        csc <- function(cpar) {
            ret <- g(c(start / scl[names(start)], cpar), scl = scl)
            return(ret[names(lambdastart)])
        }

        if (is.null(theta)) {
            lambdastart <- rep(0, length(lparnames))
            names(lambdastart) <- lparnames
        } else {
            lambdastart <- theta[lparnames]
        }

        if (length(lambdastart)) {
            for (i in 1:length(optim)) {
                op <- optim[[i]](theta = lambdastart, f = cll, g = csc)
                if (op$convergence == 0) break()
            }
            names(op$par) <- names(lambdastart)
            ret <- c(start, op$par * scl[names(lambdastart)])
            ### note: We throw away optim_hessian. vcov.mmlt
            ### uses numDeriv::hessian INCLUDING marginal parameters
            ### (which is probably the right thing to do)
            ret <- list(par = ret, value = -op$value)
        } else {
            ### no parameters to optimise over
            return(NULL)
        }
    } else {

        ui <- m$ui
        ui <- cbind(ui, matrix(0, nrow = nrow(ui), 
                               ncol = length(parnames) - ncol(ui)))
        if (!is.null(fixed)) 
            ui <- ui[, !parnames %in% names(fixed), drop = FALSE]
        ci <- m$ci

        if (is.null(theta) && !dofit) 
            return(list(ll = function(...) f(..., scl = 1), 
                        score = function(...) g(..., scl = 1), 
                        ui = ui, ci = ci))

        start <- theta / scl[eparnames]
        ui <- t(t(ui) * scl[eparnames])
  
        if (dofit) {
            for (i in 1:length(optim)) {
                ret <- optim[[i]](theta = start, f = function(par) f(par, scl = scl), 
                                                 g = function(par) g(par, scl = scl), 
                                  ui = ui, ci = ci)
                if (ret$convergence == 0) break()
            }
            if (ret$convergence != 0)
                warning("Optimisation did not converge")
        } else {
            ret <- list(par = start, value = f(theta, scl = 1), convergence = NA,
                        optim_hessian = NA)
        }
        names(ret$par) <- eparnames
        ret$par[eparnames] <- ret$par[eparnames] * scl[eparnames]
    }
  
    ret$ll <- function(...) f(..., scl = 1)
    ret$score <- function(...) g(..., scl = 1)
    ret$args <- args
    ret$logLik <- -ret$value
    ret$models <- m
    ret$formula <- formula
    ret$bx <- bx
    ret$parm <- function(par) {
        if (!is.null(fixed)) 
            par <- c(par, fixed)[parnames]
        return(c(m$parm(par), list(.Xparm(par))))
    }
    if (!missing(data))
        ret$data <- data
    ret$names <- m$names
    ret$call <- match.call()
    class(ret) <- c(ifelse(conditional, "cmmlt", "mmmlt"), "mmlt")
    ret$mmlt <- "Multivariate Conditional Transformation Model"
    ret
}


.coef.mmlt <- function(object, newdata,
                      type = c("all", "Lambda", "Lambdainv", "Precision", 
                               "PartialCorr", "Sigma", "Corr", "Spearman", "Kendall"), 
                      ...)
{
  
    type <- match.arg(type)
    if (type == "all") return(object$par)

    if (type == "Spearman")
        return(6 * asin(coef(object, newdata = newdata, type = "Cor") / 2) / pi)
  
    if (type == "Kendall")
        return(2 * asin(coef(object, newdata = newdata, type = "Cor")) / pi)

    prm <- object$parm(object$par)
    prm <- prm[[length(prm)]]

    if (missing(newdata) || is.null(object$bx)) {
        if (NROW(prm) > 1L && type != "Lambda")
            stop("newdata not specified")
        ret <- ltMatrices(t(prm), byrow = TRUE, diag = FALSE, names = object$names)
    } else {
        X <- model.matrix(object$bx, data = newdata)
        ret <- ltMatrices(t(X %*% prm), byrow = TRUE, diag = FALSE, names = object$names)
    }

    if (inherits(object, "mmmlt")) ret <- invcholD(ret)

    ret <- switch(type, "Lambda" = ret,
                        "Lambdainv" = solve(ret),
                        "Precision" = invchol2pre(ret),
                        "PartialCorr" = invchol2pc(ret),
                        "Sigma" = invchol2cov(ret),
                        "Corr" = invchol2cor(ret))
    return(ret)
}

coef.cmmlt <- function(object, newdata,
                       type = c("all", "conditional", "Lambda", "Lambdainv", 
                                "Precision", "PartialCorr", "Sigma", "Corr", 
                                "Spearman", "Kendall"), 
                       ...)
{

    type <- match.arg(type)
    if (type == "conditional") {
        prm <- object$parm(object$par)
        return(prm[-length(prm)])
    }
    return(.coef.mmlt(object = object, newdata = newdata, type = type, ...))
}

coef.mmmlt <- function(object, newdata,
                       type = c("all", "marginal", "Lambda", "Lambdainv", 
                                "Precision", "PartialCorr", "Sigma", "Corr", 
                                "Spearman", "Kendall"), 
                       ...)
{

    type <- match.arg(type)
    if (type == "marginal") {
        prm <- object$parm(object$par)
        return(prm[-length(prm)])
    }
    return(.coef.mmlt(object = object, newdata = newdata, type = type, ...))
}


vcov.mmlt <- function(object, ...) {
    step <- 0
    lam <- 1e-6
    H <- object$optim_hessian
    if (is.null(H)) {
        if (requireNamespace("numDeriv")) {
            H <- numDeriv::hessian(object$ll, object$par)
        } else {
            stop("Hessian not available")
        }
    }
    while((step <- step + 1) <= 3) {
          ret <- try(solve(H + (step - 1) * lam * diag(nrow(H))))
          if (!inherits(ret, "try-error")) break
    }
    if (inherits(ret, "try-error"))
        stop("Hessian is not invertible")
    if (step > 1)
        warning("Hessian is not invertible, an approximation is used")
    rownames(ret) <- colnames(ret) <- names(coef(object))
    ret
}

logLik.mmlt <- function (object, parm = coef(object), newdata = NULL, ...) 
{
    args <- list(...)
    if (length(args) > 0) 
        warning("Arguments ", names(args), " are ignored")
    ret <- -object$ll(parm, newdata = newdata)
    attr(ret, "df") <- length(object$par)
    class(ret) <- "logLik"
    ret
}

estfun.mmlt <- function(x, parm = coef(x, type = "all"), 
                        newdata = NULL, ...) {
    args <- list(...)
    if (length(args) > 0)
        warning("Arguments ", names(args), " are ignored")
    return(x$score(parm, newdata = newdata, scores = TRUE))
}

summary.mmlt <- function(object, ...) {
    ret <- list(call = object$call,
                #                tram = object$tram,
                test = cftest(object, parm = names(coef(object, with_baseline = FALSE))),
                ll = logLik(object))
    class(ret) <- "summary.mmlt"
    ret
}

print.summary.mmlt <- function(x, digits = max(3L, getOption("digits") - 3L), ...) {
    cat("\n", x$mmlt, "\n")
    cat("\nCall:\n")
    print(x$call)
    cat("\nCoefficients:\n")
    pq <- x$test$test
    mtests <- cbind(pq$coefficients, pq$sigma, pq$tstat, pq$pvalues)
    colnames(mtests) <- c("Estimate", "Std. Error", "z value", "Pr(>|z|)")
    sig <- .Machine$double.eps
    printCoefmat(mtests, digits = digits, has.Pvalue = TRUE, 
                 P.values = TRUE, eps.Pvalue = sig)
    cat("\nLog-Likelihood:\n ", x$ll, " (df = ", attr(x$ll, "df"), 
          ")", sep = "")
    cat("\n\n")
    invisible(x)
}

print.mmlt <- function(x, ...) {
    cat("\n", "Multivariate conditional transformation model", "\n")
    cat("\nCall:\n")
    print(x$call)
    cat("\nCoefficients:\n")
    print(coef(x))
    invisible(x)
}


predict.mmlt <- function (object, newdata, margins = 1:J, 
    type = c("trafo", "distribution", "survivor", "density", "hazard"), log = FALSE, 
    args = object$args, ...) 
{
    J <- length(object$models$models)
    margins <- sort(margins)
    stopifnot(all(margins %in% 1:J))

    if (length(margins) == 1L) {
        ### ... may carry q = something
        tmp <- object$models$models[[margins]]
        cf <- coef(tmp)
        cf[] <- object$models$parm(coef(object))[[margins]]
        coef(tmp) <- cf
        ### marginal models
        if (!inherits(object, "cmmlt")) {
            ret <- predict(tmp, newdata = newdata, type = type, log = log, ...)
            return(ret)
        }
        ### conditional models
        mcov <- coef(object, newdata = newdata, type = "Sigma")
        msd <- sqrt(diagonals(mcov)[margins,])
        if (length(unique(msd)) == 1L && 
            !"bscaling" %in% names(tmp$model$model)) { ### no stram model
            cf <- cf / unique(msd)
            coef(tmp) <- cf
            ret <- predict(tmp, newdata = newdata, type = type, log = log, ...)
            return(ret)
        }
        type <- match.arg(type)
        tr <- predict(tmp, newdata = newdata, type = "trafo", ...) 
        msd <- matrix(msd, nrow = nrow(tr), ncol = ncol(tr), byrow = TRUE)
        tr <- tr / msd
        switch(type, "trafo" = return(tr),
                     "distribution" = return(pnorm(tr, log.p = log)),
                     "survivor" = return(pnorm(tr, log.p = log, lower.tail = FALSE)),
                     "density" = {
                         dx <- 1
                         names(dx) <- variable.names(tmp)[1L]
                         dtr <- predict(tmp, newdata = newdata, type = "trafo", deriv = dx, ...)
                         ret <- dnorm(tr, log = TRUE) - .log(msd) + .log(dtr)
                         if (log) return(ret)
                         return(exp(ret))
                     },
                     stop("not yet implemented"))
    }

    type <- match.arg(type)
    ### don't feed ...
    z <- .mget(object$models, margins, parm = coef(object, type = "all"),
               newdata = newdata, what = "z")
    z <- .rbind(z)

    if (type == "trafo") {
        stopifnot(!log)
        L <- coef(object, newdata = newdata, type = "Lambda")
        if (length(margins) != J) 
            L <- marg_mvnorm(invchol = L, which = margins)$invchol
        return(Mult(L, z))
    }
    if (type == "distribution") {
        lower <- matrix(-Inf, ncol = ncol(z), nrow = nrow(z))
        upper <- z
        Linv <- coef(object, newdata = newdata, type = "Lambdainv")
        if (length(margins) != J) 
            Linv <- marg_mvnorm(chol = Linv, which = margins)$chol
        a <- args
        a$lower <- lower
        a$upper <- upper
        a$logLik <- FALSE
        a$chol <- Linv
        ret <- do.call("lpmvnorm", a)
        if (log) return(ret)
        return(exp(ret))
    }
    if (type == "survivor") {
        lower <- z 
        upper <- matrix(Inf, ncol = ncol(z), nrow = nrow(z))
        Linv <- coef(object, newdata = newdata, type = "Lambdainv")
        if (length(margins) != J) 
            Linv <- marg_mvnorm(chol = Linv, which = margins)$chol
        a <- args
        a$lower <- lower
        a$upper <- upper
        a$logLik <- FALSE
        a$chol <- Linv
        ret <- do.call("lpmvnorm", a)
        if (log) return(ret)
        return(exp(ret))
    }
    stopifnot(type == "density")
    stopifnot(all(object$models$cont))
    zprime <- .mget(object$models, margins, parm = coef(object, type = "all"),
                    newdata = newdata, what = "zprime")
    if (length(margins) > 1L) {
        zprime <- .rbind(zprime)
    } else {
        zprime <- matrix(zprime, nrow = 1)
    }
    L <- coef(object, newdata = newdata, type = "Lambda")
        if (length(margins) != J) 
            L <- marg_mvnorm(invchol = L, which = margins)$invchol
    ret <- ldmvnorm(obs = z, invchol = L, logLik = FALSE)
    ret <- ret + colSums(.log(zprime))
    if (log) return(ret)
    return(exp(ret))
}

simulate.mmlt <- function(object, nsim = 1L, seed = NULL, newdata, K = 50, ...) {

    ### from stats:::simulate.lm
    if (!exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) 
        runif(1)
    if (is.null(seed)) 
        RNGstate <- get(".Random.seed", envir = .GlobalEnv)
    else {
        R.seed <- get(".Random.seed", envir = .GlobalEnv)
        set.seed(seed)
        RNGstate <- structure(seed, kind = as.list(RNGkind()))
        on.exit(assign(".Random.seed", R.seed, envir = .GlobalEnv))
    }

    if (!is.data.frame(newdata))
        stop("not yet implemented")

    args <- list(...)
    if (length(args) > 0L)
        stop("argument(s)", paste(names(args), collapse = ", "), "ignored")

    if (nsim > 1L) 
        return(replicate(nsim, simulate(object, newdata = newdata, K = K, ...), 
                         simplify = FALSE))

    J <- length(object$models$models)
    L <- coef(object, newdata = newdata, type = "Lambda")
    N <- nrow(newdata)

    Z <- matrix(rnorm(J * N), ncol = N)
    Ztilde <- solve(L, Z)

    ret <- matrix(0.0, nrow = N, ncol = J)

    if (inherits(object, "cmmlt")) {
        for (j in 1:J) {
            q <- mkgrid(object$models$models[[j]], n = K)[[1L]]
            pr <- predict(object$models$models[[j]], newdata = newdata, type = "trafo", q = q)
            if (!is.matrix(pr)) pr <- matrix(pr, nrow = length(pr), ncol = NROW(newdata))
            ret[,j] <- as.double(mlt:::.invf(object$models$models[[j]], f = t(pr), 
                                             q = q, z = t(Ztilde[j,,drop = FALSE])))
        }
    } else {
        Ztilde <- pnorm(Ztilde, log.p = TRUE)
        for (j in 1:J) {
            q <- mkgrid(object$models$models[[j]], n = K)[[1L]]
            pr <- predict(object$models$models[[j]], newdata = newdata, type = "logdistribution", q = q)
            if (!is.matrix(pr)) pr <- matrix(pr, nrow = length(pr), ncol = NROW(newdata))
            ret[,j] <- as.double(mlt:::.invf(object$models$models[[j]], f = t(pr), 
                                             q = q, z = t(Ztilde[j,,drop = FALSE])))
        }
    }
    colnames(ret) <- variable.names(object, response_only = TRUE)
    return(ret)
}

variable.names.mmlt <- function(object, response_only = FALSE, ...) {

    if (response_only)
        return(sapply(object$models$models, function(x) variable.names(x)[1L]))
    vn <- unique(c(sapply(object$models$models, function(x) variable.names(x)), 
                 all.vars(object$formula)))
    return(vn)
}
    
confregion <- function(object, level = .95, ...)
    UseMethod("confregion")

confregion.mmlt <- function(object, level = .95, newdata, K = 250, ...) {

    if (!missing(newdata)) stopifnot(nrow(newdata) == 1)

    Linv <- coef(object, newdata = newdata, type = "Lambdainv")
    Linv <- as.array(Linv)[,,1]
    J <- nrow(Linv)

    q <- qchisq(level, df = J)

    if (J == 2) {
        angle <- seq(0, 2 * pi, length = K)
        x <- cbind(cos(angle), sin(angle))
    } else {
        x <- matrix(rnorm(K * J), nrow = K, ncol = J)
        x <- x / sqrt(rowSums(x^2))
    }
    x <- sqrt(q) * x
    a <- x %*% t(Linv)

    nd <- if (missing(newdata)) data.frame(1) else newdata

    ret <- lapply(1:J, function(j) {
        prb <- object$models$models[[j]]$todistr$p(a[,j])
        predict(object$models$models[[j]], newdata = nd, type = "quantile", prob = prb)
    })
    
    ret <- do.call("cbind", ret)
    return(ret)
}

HDR <- function(object, level = .95, ...)
    UseMethod("HDR")

HDR.mmlt <- function(object, level = .95, newdata, nsim = 1000L, K = 25, ...) {

    if (!missing(newdata)) {
        stopifnot(nrow(newdata) == 1)
    } else {
        newdata <- data.frame(1)
    }

    ### https://doi.org/10.2307/2684423 Section 3.2
    y <- simulate(object, newdata = newdata[rep(1, nsim),,drop = FALSE])
    y <- cbind(y, newdata)
    d <- predict(object, newdata = y, type = "density")

    ret <- do.call("expand.grid", lapply(object$models$models, function(x) mkgrid(x, n = K)[[1L]]))
    colnames(ret) <- variable.names(object, response_only = TRUE)
    ret <- cbind(ret, newdata)
    ret$density <- predict(object, newdata = ret, type = "density")
    attr(ret, "cuts") <- quantile(d, prob = 1 - level)
    ret
}

mkgrid.mmlt <- function(object, ...) {

    lx <- mkgrid(as.basis(object$formula, data = object$data), ...)
    grd <- do.call("c", lapply(object$models$models, mkgrid, ...))
    grd <- c(grd, lx)
    do.call("expand.grid", grd[unique(names(grd))])
}