R/clmm2.R

Defines functions update.clmm2 plot.profile.clmm2 confint.profile.clmm2 profile.clmm2 print.summary.clmm2 summary.clmm2 vcov.clmm2 print.clmm2 finalizeRhoM getNGHQinC getNGHQinR getNAGQinC tmpAGQ getNAGQinR getNLA2 clmm2 update.u2.v3 .negLogLikMfast .negLogLikBase clmm2.control

Documented in clmm2 clmm2.control confint.profile.clmm2 plot.profile.clmm2 profile.clmm2 profile.clmm2 update.clmm2

#############################################################################
##    Copyright (c) 2010-2022 Rune Haubo Bojesen Christensen
##
##    This file is part of the ordinal package for R (*ordinal*)
##
##    *ordinal* is free software: you can redistribute it and/or modify
##    it under the terms of the GNU General Public License as published by
##    the Free Software Foundation, either version 2 of the License, or
##    (at your option) any later version.
##
##    *ordinal* is distributed in the hope that it will be useful,
##    but WITHOUT ANY WARRANTY; without even the implied warranty of
##    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##    GNU General Public License for more details.
##
##    A copy of the GNU General Public License is available at
##    <https://www.r-project.org/Licenses/> and/or
##    <http://www.gnu.org/licenses/>.
#############################################################################
## This file contains:
## The main clmm2 function and some related auxiliary functions.

clmm2.control <-
    function(method = c("ucminf", "nlminb", "model.frame"),
             ..., trace = 0, maxIter = 50, gradTol = 1e-4,
             maxLineIter = 50,
             innerCtrl = c("warnOnly", "noWarn", "giveError"))
{
    method <- match.arg(method)
    innerCtrl <- match.arg(innerCtrl)
    ctrl <- list(trace=ifelse(trace < 0, 1, 0),
                 maxIter=maxIter,
                 gradTol=gradTol,
                 maxLineIter=maxLineIter,
                 innerCtrl=innerCtrl)
    optCtrl <- list(trace = abs(trace), ...)

    if(!is.numeric(unlist(ctrl[-5])))
        stop("maxIter, gradTol, maxLineIter and trace should all be numeric")
    if(any(ctrl[-c(1, 5)] <= 0))
       stop("maxIter, gradTol and maxLineIter have to be > 0")
    if(method == "ucminf" && !"grtol" %in% names(optCtrl))
        optCtrl$grtol <- 1e-5
    if(method == "ucminf" && !"grad" %in% names(optCtrl))
        optCtrl$grad <- "central"

    list(method = method, ctrl = ctrl, optCtrl = optCtrl)
}


.negLogLikBase <- function(rho) {
### Update stDev, sigma, eta1Fix, and eta2Fix given new par:
    with(rho, {
        if(estimLambda > 0)
            lambda <- par[nxi + p + k + 1:estimLambda]
        if(estimStDev)
            stDev <- exp(par[p+nxi+k+estimLambda+ 1:s])
        sigma <-
            if(k > 0) expSoffset * exp(drop(Z %*% par[nxi+p + 1:k]))
            else expSoffset
        eta1Fix <- drop(B1 %*% par[1:(nxi + p)])
        eta2Fix <- drop(B2 %*% par[1:(nxi + p)])
    })
    return(invisible())
}

.negLogLikMfast <- function(rho) { ## negative log-likelihood
    fit <- with(rho, {
        .C("nll",
           as.double(u),
           length(u),
           as.integer(grFac),
           as.double(stDev),
           as.double(o1),
           as.double(o2),
           length(o1),
           eta1 = as.double(eta1),
           eta2 = as.double(eta2),
           as.double(eta1Fix),
           as.double(eta2Fix),
           as.double(sigma),
           pr = as.double(pr),
           as.double(weights),
           as.double(lambda),
           as.integer(linkInt),
           nll = double(1)
           )[c("eta1", "eta2", "pr", "nll")]
    })
    rho$eta1 <- fit$eta1
    rho$eta2 <- fit$eta2
    rho$pr <- fit$pr
    fit$nll
}

update.u2.v3 <- function(rho) {
### third version: C-implementation of NR-algorithm.
    .negLogLikBase(rho) ## update: par, stDev, eta1Fix, eta2Fix eta2Fix, sigma
    fit <- with(rho,
                .C("NRalgv3",
                   as.integer(ctrl$trace),
                   as.integer(ctrl$maxIter),
                   as.double(ctrl$gradTol),
                   as.integer(ctrl$maxLineIter),
                   as.integer(grFac),
                   as.double(stDev),
                   as.double(o1),
                   as.double(o2),
                   as.double(eta1Fix),
                   as.double(eta2Fix),
                   as.double(sigma),
                   as.integer(linkInt),
                   as.double(weights),
                   u = as.double(uStart),
                   pr = as.double(pr),
                   funValue = double(1),
                   gradValues = as.double(uStart),
                   hessValues = as.double(rep(1, length(uStart))),
                   length(pr),
                   length(uStart),
                   maxGrad = double(1),
                   conv = 0L,
                   as.double(lambda),
                   Niter = as.integer(Niter)
                   )[c("u", "funValue", "gradValues",
                       "hessValues", "maxGrad", "conv", "Niter")] )
    ## Get message:
    message <- switch(as.character(fit$conv),
                      "1" = "max|gradient| < tol, so current iterate is probably solution",
                      "0" = "Non finite negative log-likelihood",
                      "-1" = "iteration limit reached when updating the random effects",
                      "-2" = "step factor reduced below minimum when updating the random effects")
    ## Check for convergence and report warning/error:
    if(rho$ctrl$trace > 0 && fit$conv == 1)
        cat("\nOptimizer converged! ", "max|grad|:",
            fit$maxGrad, message, fill = TRUE)
    if(fit$conv != 1 && rho$ctrl$innerCtrl == "warnOnly")
        warning(message, "\n  at iteration ", rho$Niter)
    else if(fit$conv != 1 && rho$ctrl$innerCtrl == "giveError")
        stop(message, "\n  at iteration ", rho$Niter)
    ## Store values and return:
    rho$Niter <- fit$Niter
    rho$u <- fit$u
    rho$D <- fit$hessValue
    rho$gradient <- fit$gradValue
    if(!is.finite(rho$negLogLik <- fit$funValue))
        return(FALSE)
    return(TRUE)
}

clmm2 <-
  function(location, scale, nominal, random, data, weights, start, subset,
           na.action, contrasts, Hess = FALSE, model = TRUE, sdFixed,
           link = c("logistic", "probit", "cloglog", "loglog",
           "cauchit", "Aranda-Ordaz", "log-gamma"), lambda,
           doFit = TRUE, control, nAGQ = 1,
           threshold = c("flexible", "symmetric", "equidistant"), ...)
    ## Handle if model = FALSE
### Marginal fitted values? (pr | u = 0) or (pr | u = u.tilde) ?
### How can we (should we?) get u.tilde and var(u.tilde) with GHQ?
### Make safeStart function if !is.finite(negLogLik)
### Write test suite for doFit argument
{
    R <- match.call(expand.dots = FALSE)
    Call <- match.call()
    if(missing(random)) {
        Call[[1]] <- as.name("clm2")
        return(eval.parent(Call))
    }
    if(missing(lambda)) lambda <- NULL
    if(missing(contrasts)) contrasts <- NULL
    if(missing(control)) control <- clmm2.control(...)
    if(!setequal(names(control), c("method", "ctrl", "optCtrl")))
       stop("specify 'control' via clmm2.control()")
    if (missing(data)) data <- environment(location)
    if (is.matrix(eval.parent(R$data)))
        R$data <- as.data.frame(data)
### Collect all variables in a single formula and evaluate to handle
### missing values correctly:
    m <- match(c("location", "scale", "nominal"), names(R), 0)
    F <- lapply(as.list(R[m]), eval.parent) ## evaluate in parent
    varNames <- unique(unlist(lapply(F, all.vars)))
    longFormula <-
        eval(parse(text = paste("~", paste(varNames, collapse = "+")))[1])
    m <- match(c("location", "data", "subset", "weights", "random",
                 "na.action"), names(R), 0)
    R <- R[c(1, m)]
    R$location <- longFormula
    R$drop.unused.levels <- TRUE
    R[[1]] <- as.name("model.frame")
    names(R)[names(R) == "location"] <- "formula"
    R <- eval.parent(R)
    nonNA <- rownames(R)
### Append nonNA index to Call$subset to get the right design matrices
### from clm2:
    Call$subset <-
        if(is.null(Call$subset)) nonNA
        else c(paste(deparse(Call$subset), "&"), nonNA)
    Call$start <-
        if(is.null(Call$start) || !is.null(Call$sdFixed)) Call$start
        else start[-length(start)]
    Call$random <- Call$control <- Call$nAGQ <- Call$sdFixed <-
        Call$innerCtrl <- NULL
    Call$method <- control$method
    Call$doFit <- Call$Hess <- FALSE
    Call[[1]] <- as.name("clm2")
    rhoM <- eval.parent(Call)
    if(control$method == "model.frame")
        return(rhoM)
    rhoM$call <- match.call()
    rhoM$randomName <- deparse(rhoM$call$random)
### Set grouping factor and stDev parameter:
    rhoM$grFac <- R[,"(random)"]
    if(!missing(sdFixed) && !is.null(sdFixed)) {
        stopifnot(length(sdFixed) == 1 && sdFixed > 0)
        rhoM$estimStDev <- FALSE
        rhoM$stDev <- sdFixed
    }
    else
        rhoM$estimStDev <- TRUE
    with(rhoM, {
        r <- nlevels(grFac) ## no. random effects
        grFac <- as.integer(unclass(grFac))
        if(r <= 2) stop("Grouping factor must have 3 or more levels")
        s <- ifelse(estimStDev, 1L, 0L) ## no. variance parameters
        Niter <- 0L
    })
### set starting values:
    if(missing(start)) {
        suppressWarnings(fitCLM(rhoM))
        if(rhoM$estimStDev) rhoM$start <- rhoM$par <- c(rhoM$par, log(1))
        else rhoM$start <- rhoM$par
    } else
        rhoM$start <- rhoM$par <- start
    rhoM$uStart <- rhoM$u <- rep(0, rhoM$r)
### Test starting values:
    if(length(rhoM$start) !=
       with(rhoM, nxi + p + k + estimLambda + estimStDev))
        stop("'start' is ", length(rhoM$start),
             " long, but should be ", with(rhoM, nxi + p + k + estimLambda + estimStDev),
             " long")
    if(rhoM$ncolXX == 0) {
        if(!all(diff(c(rhoM$tJac %*% rhoM$start[1:rhoM$nalpha])) > 0))
            stop("Threshold starting values are not of increasing size")
    }
### Change the lower limit if lambda is estimated with the
### Aranda-Ordaz link and sdFixed is not supplied:
    if(rhoM$estimLambda > 0 && rhoM$link == "Aranda-Ordaz" &&
       is.null(rhoM$call$sdFixed))
        rhoM$limitLow <- c(rep(-Inf, length(rhoM$par)-2), 1e-5, -Inf)
### This should hardly ever be the case:
    .negLogLikBase(rhoM) ## set lambda, stDev, sigma, eta1Fix and eta2Fix
    if(!is.finite(.negLogLikMfast(rhoM)))
        stop("Non-finite integrand at starting values")
    rhoM$ctrl <- control$ctrl
    rhoM$optCtrl <- control$optCtrl
    if(rhoM$method == "nlminb") {
        m <- match(names(rhoM$optCtrl), c("grad","grtol"), 0)
        rhoM$optCtrl <- rhoM$optCtrl[!m]
    }
### Match doFit:
    if(is.logical(doFit) || is.numeric(doFit)) {
        if(doFit) doFit <- "C"
        else doFit <- "no"
    }
    else if(!is.character(doFit) || !doFit %in% c("no", "R", "C"))
        stop("argument 'doFit' not recognized. 'doFit' should be\n
numeric, logical or one of c('no', 'R', 'C')")

### Set ObjFun parameters:
    ObjFun <- getNLA2 ## same for "R" and "C"
    rhoM$updateU <-
        if(doFit == "R") update.u2
        else update.u2.v3
    rhoM$nAGQ <- as.integer(nAGQ)
    if(rhoM$nAGQ >= 2) {
        ghq <- gauss.hermite(rhoM$nAGQ)
        rhoM$ghqns <- ghq$nodes
        rhoM$ghqws <- ghq$weights
        if(doFit == "R") {
            ObjFun <- getNAGQinR
            rhoM$PRnn <- array(0, dim=c(rhoM$n, rhoM$nAGQ))
            rhoM$PRrn <- array(0, dim=c(rhoM$r, rhoM$nAGQ))
            rhoM$ghqws <- ghq$weights * exp(rhoM$ghqns^2)
        }
        else
            ObjFun <- getNAGQinC
    }
    if(rhoM$nAGQ <= -1) {
        ghq <- gauss.hermite(abs(rhoM$nAGQ))
        rhoM$ghqns <- ghq$nodes
        rhoM$ghqws <- ghq$weights * exp((ghq$nodes^2)/2)
        if(doFit == "R"){
            ObjFun <- getNGHQinR
        }
        else {
            ObjFun <- getNGHQinC
            rhoM$ghqws <- log(ghq$weights) + (ghq$nodes^2)/2
        }
    }
    stopifnot(rhoM$nAGQ != 0) ## test needed?

### Fit the model:
    if(!doFit %in% c("C", "R"))
        return(rhoM)
    if(rhoM$nAGQ > -1)
        rhoM$updateU(rhoM) # Try updating the random effects
    rhoM$optRes <- switch(rhoM$method,
                       "ucminf" = ucminf(rhoM$start, function(x)
                       ObjFun(rhoM, x), control=rhoM$optCtrl),
                       "nlminb" = nlminb(rhoM$start, function(x)
                       ObjFun(rhoM, x), control=rhoM$optCtrl,
                       lower = rhoM$limitLow, upper = rhoM$limitUp))
    rhoM$par <- rhoM$optRes[[1]]
    if(Hess) {
        if(rhoM$link == "Aranda-Ordaz" &&
           rhoM$estimLambda > 0 && rhoM$lambda < 1e-3)
            message("Cannot get Hessian because lambda = ",rhoM$lambda
                    ," is too close to boundary.\n",
                    " Fit model with link == 'logistic' to get Hessian")
        else {
            rhoM$Hessian <- myhess(function(x) ObjFun(rhoM, x),
                                    rhoM$par)
            rhoM$par <- rhoM$optRes[[1]]
        }
    }
    .negLogLikMfast(rhoM) ## update pr
    ## if(rhoM$nAGQ > -1)
    rhoM$updateU(rhoM) # Makes sure ranef's are evaluated at the optimum
### Post processing:
    res <- finalizeRhoM(rhoM)
    res$call <- match.call()
    res$na.action <- attr(R, "na.action")
    res$contrasts <- contrasts
    class(res) <- c("clmm2", "clm2")
    res
}

getNLA2 <- function(rho, par) {
### negative log-likelihood by the Laplace approximation
### (with update.u2 in C or R):
    if(!missing(par)) rho$par <- par
    if(!rho$updateU(rho)) return(Inf)
    if(any(rho$D < 0)) return(Inf)
    ## logDetD <- sum(log(rho$D/(2*pi)))
    logDetD <- sum(log(rho$D)) - rho$r * log(2*pi)
    rho$negLogLik + logDetD / 2
}

getNAGQinR <- function(rho, par) {
### negative log-likelihood by adaptive Gauss-Hermite quadrature
### implemented in R:
    if(!missing(par))
        rho$par <- par
    if(!rho$updateU(rho)) return(Inf)
    if(any(rho$D < 0)) return(Inf)
    with(rho, {
        K <- sqrt(2/D)
        agqws <- K %*% t(ghqws)
        agqns <- apply(K %*% t(ghqns), 2, function(x) x + u)
        ranNew <- apply(agqns, 2, function(x) x[grFac] * stDev)

        eta1Tmp <- (eta1Fix + o1 - ranNew) / sigma
        eta2Tmp <- (eta2Fix + o2 - ranNew) / sigma
        if(nlambda)
            ## PRnn <- (pfun(eta1Tmp, lambda) - pfun(eta2Tmp, lambda))^weights
            ## This is likely a computationally more safe solution:
          PRnn <- exp(weights * log(pfun(eta1Tmp, lambda) -
                                    pfun(eta2Tmp, lambda)))
        else
            ## PRnn <- (pfun(eta1Tmp) - pfun(eta2Tmp))^weights
            PRnn <- exp(weights * log(pfun(eta1Tmp) - pfun(eta2Tmp)))
### OPTION: The fitted values could be evaluated with getFittedC for
### better precision.
        for(i in 1:r)
            ## PRrn[i,] <- apply(PRnn[grFac == i, ], 2, prod)
            PRrn[i,] <- apply(PRnn[grFac == i, ,drop = FALSE], 2, prod)
        PRrn <- PRrn * agqws * dnorm(x=agqns, mean=0, sd=1)
### OPTION: Could this be optimized by essentially computing dnorm 'by hand'?
    })
    -sum(log(rowSums(rho$PRrn)))
}

## tmpAGQ(rho)

tmpAGQ <- function(rho, par) {
    if(!missing(par))
        rho$par <- par
    with(rho, {
        ls()
        stDev <- exp(ST[[1]][1, 1])
        nlambda <- 0
        K <- sqrt(2/D)
        agqws <- K %*% t(ghqws)
        agqns <- apply(K %*% t(ghqns), 2, function(x) x + u)
        grFac <- unclass(grFac)
        ranNew <- apply(agqns, 2, function(x) x[grFac] * stDev)
        eta1Tmp <- (eta1Fix + o1 - ranNew) / sigma
        eta2Tmp <- (eta2Fix + o2 - ranNew) / sigma
        if(nlambda)
            PRnn <- exp(weights * log(pfun(eta1Tmp, lambda) -
                                      pfun(eta2Tmp, lambda)))
        else
            PRnn <- exp(wts * log(pfun(eta1Tmp) - pfun(eta2Tmp)))

        dim(eta1Tmp)


        exp(wts[IND] * log(pfun(eta1Tmp[IND, ]) - pfun(eta2Tmp[IND, ])))

        PRrn <- do.call(rbind, lapply(1:dims$q, function(i) {
            apply(PRnn[grFac == i, ,drop = FALSE], 2, prod)
        }))
        head(PRrn)

        PRrn <- do.call(rbind, lapply(1:dims$q, function(i) {
            apply(PRnn[grFac == i, ,drop = FALSE], 2, function(x) sum(log(x)))
        }))
        head(PRrn)
        ## Could we do something like
        PRnn <- wts * log(pfun(eta1Tmp) - pfun(eta2Tmp))
        PRrn <- do.call(rbind, lapply(1:dims$q, function(i) {
            apply(PRnn[grFac == i, ,drop = FALSE], 2, function(x) sum(x))
        }))
        head(PRrn, 20)
        ## to avoid first exp()ing and then log()ing?
        head(exp(PRrn), 20)
        range(PRrn)
        exp(range(PRrn))

        out <- PRrn + log(agqws) + log(dnorm(x=agqns, mean=0, sd=1))


        log(2 * 3)
        log(2) + log(3)

        PRnn[grFac == 12, , drop=FALSE]
        IND <- which(grFac == 12)
        cbind(IND, wts[IND], PRnn[IND, ])

        dim(PRrn)
        ## There seems to be underfloow allready in the computations
        ## in PRnn, which propagates to PRrn
        PRrn <- PRrn * agqws * dnorm(x=agqns, mean=0, sd=1)
    })
    -sum(log(rowSums(rho$PRrn)))
}

getNAGQinC <- function(rho, par) {
### negative log-likelihood by adaptive Gauss-Hermite quadrature
### implemented in C:
    if(!missing(par))
        rho$par <- par
    if(!rho$updateU(rho)) return(Inf)
    if(any(rho$D < 0)) return(Inf)
    with(rho, {
        .C("getNAGQ",
           nll = double(1), ## nll
           as.integer(grFac), ## grFac
           as.double(stDev), ## stDev
           as.double(eta1Fix),
           as.double(eta2Fix),
           as.double(o1),
           as.double(o2),
           as.double(sigma), ## Sigma
           as.double(weights),
           length(sigma), ## nx - no. obs
           length(uStart), ## nu - no. re
           as.double(ghqns),
           as.double(log(ghqws)), ## lghqws
           as.double(ghqns^2), ## ghqns2
           as.double(u),
           as.double(D),
           as.integer(abs(nAGQ)),
           as.integer(linkInt),
           as.double(lambda))$nll
    })
}

getNGHQinR <- function(rho, par) {
### negative log-likelihood by standard Gauss-Hermite quadrature
### implemented in R:
  if(!missing(par))
    rho$par <- par
  .negLogLikBase(rho) ## Update lambda, stDev, sigma and eta*Fix
  with(rho, {
    ns <- ghqns * stDev
    SS <- numeric(r) ## summed likelihood
    for(i in 1:r) {
      ind <- grFac == i
      eta1Fi <- eta1Fix[ind]
      eta2Fi <- eta2Fix[ind]
      o1i <- o1[ind]
      o2i <- o2[ind]
      si <- sigma[ind]
      wt <- weights[ind]
      for(h in 1:abs(nAGQ)) {
        eta1s <- (eta1Fi + o1i - ns[h]) / si
        eta2s <- (eta2Fi + o2i - ns[h]) / si
        ## SS[i] <- exp(sum(wt * log(pfun(eta1s) - pfun(eta2s)))) *
        ##     ghqws[h] * exp(ghqns[h]^2) * dnorm(x=ghqns[h]) + SS[i]
        SS[i] <- exp(sum(wt * log(pfun(eta1s) - pfun(eta2s)))) *
          ghqws[h] + SS[i]
### OPTION: The fitted values could be evaluated with getFittedC for
### better precision.
      }
    }
    -sum(log(SS)) + r * log(2*pi)/2
  })
}

getNGHQinC <- function(rho, par) {
### negative log-likelihood by standard Gauss-Hermite quadrature
### implemented in C:
    if(!missing(par))
        rho$par <- par
    .negLogLikBase(rho) ## Update lambda, stDev, sigma and eta*Fix
    with(rho, {
        .C("getNGHQ_C",
           nll = double(1),
           as.integer(grFac),
           as.double(stDev),
           as.double(eta1Fix),
           as.double(eta2Fix),
           as.double(o1),
           as.double(o2),
           as.double(sigma),
           as.double(weights),
           length(sigma),
           length(uStart),
           as.double(ghqns),
           as.double(ghqws),
           as.integer(abs(nAGQ)),
           as.integer(linkInt),
           as.double(ghqns * stDev),
           as.double(lambda))$nll
    })
}

finalizeRhoM <- function(rhoM) {
    if(rhoM$method == "ucminf") {
        if(rhoM$optRes$info[1] > rhoM$optCtrl[["grtol"]])
            warning("clmm2 may not have converged:\n  optimizer 'ucminf' terminated with max|gradient|: ",
                    rhoM$optRes$info[1], call.=FALSE)
        rhoM$convergence <-
            ifelse(rhoM$optRes$info[1] > rhoM$optCtrl[["grtol"]], FALSE, TRUE)
    }
    if(rhoM$method == "nlminb") {
        rhoM$convergence <- ifelse(rhoM$optRes$convergence == 0, TRUE, FALSE)
        if(!rhoM$convergence)
            warning("clmm2 may not have converged:\n  optimizer 'nlminb' terminated with message: ",
                    rhoM$optRes$message, call.=FALSE)
    }
    if(rhoM$ctrl$gradTol < max(abs(rhoM$gradient)))
        warning("Inner loop did not converge at termination:\n  max|gradient| = ",
                max(abs(rhoM$gradient)))
    with(rhoM, {
        if(nxi > 0) {
            xi <- par[1:nxi]
            names(xi) <- xiNames
            thetaNames <- paste(lev[-length(lev)], lev[-1], sep="|")
            Alpha <- Theta <- matrix(par[1:nxi], nrow=ncolXX, byrow=TRUE)
            Theta <- t(apply(Theta, 1, function(x) c(tJac %*% x)))
            if(ncolXX > 1){
                dimnames(Theta) <- list(dnXX[[2]], thetaNames)
                dimnames(Alpha) <- list(dnXX[[2]], alphaNames)
            }
            else {
                Theta <- c(Theta)
                Alpha <- c(Alpha)
                names(Theta) <- thetaNames
                names(Alpha) <- alphaNames
            }
            coefficients <- xi
        }
        else coefficients <- numeric(0)
        if(p > 0) {
            beta <- par[nxi + 1:p]
            names(beta) <- dnX[[2]]
            coefficients <- c(coefficients, beta)
        }
        if(k > 0) {
            zeta <- par[nxi+p + 1:k]
            names(zeta) <- dnZ[[2]]
            coefficients <- c(coefficients, zeta)
        }
        if(estimLambda > 0) {
            names(lambda) <- "lambda"
            coefficients <- c(coefficients, lambda)
        }
        if(s > 0) {
            stDev <- exp(par[nxi+p+k + estimLambda + 1:s])
            coefficients <- c(coefficients, stDev)
        }
        names(stDev) <- randomName
        if(exists("Hessian", inherits=FALSE))
            dimnames(Hessian) <- list(names(coefficients),
                                      names(coefficients))
        edf <- p + nxi + k + estimLambda + s
        nobs <- sum(weights)
        fitted.values <- pr
        df.residual = nobs - edf
        ranef <- u * stDev
        condVar <- 1/D * stDev^2
        logLik <- -optRes[[2]]
    })
    res <- as.list(rhoM)
    keepNames <-
        c("ranef", "df.residual", "fitted.values", "edf", "start",
          "stDev", "beta", "coefficients", "zeta", "Alpha", "Theta",
          "xi", "lambda", "convergence", "Hessian",
          "gradient", "optRes", "logLik", "Niter", "grFac", "call",
          "scale", "location", "nominal", "method", "y", "lev",
          "nobs", "threshold", "estimLambda", "link", "nAGQ",
          "condVar", "contrasts", "na.action")
    m <- match(keepNames, names(res), 0)
    res <- res[m]
    res
}

anova.clmm2 <- function (object, ..., test = c("Chisq", "none"))
{
    anova.clm2(object, ..., test = c("Chisq", "none"))
}

print.clmm2 <- function(x, ...)
{
  if(x$nAGQ >= 2)
    cat(paste("Cumulative Link Mixed Model fitted with the adaptive",
              "Gauss-Hermite \nquadrature approximation with",
              x$nAGQ ,"quadrature points"), "\n\n")
  else if(x$nAGQ <= -1)
    cat(paste("Cumulative Link Mixed Model fitted with the",
              "Gauss-Hermite \nquadrature approximation with",
              abs(x$nAGQ) ,"quadrature points"), "\n\n")
  else
    cat("Cumulative Link Mixed Model fitted with the Laplace approximation\n",
        fill=TRUE)
    if(!is.null(cl <- x$call)) {
        cat("Call:\n")
        dput(cl, control=NULL)
    }
    if(length(x$stDev)) {
        cat("\nRandom effects:\n")
        varMat <- matrix(c(x$stDev^2, x$stDev), nrow =
                         length(x$stDev), ncol=2)
        rownames(varMat) <- names(x$stDev)
        colnames(varMat) <- c("Var", "Std.Dev")
        print(varMat, ...)
    } else {
        cat("\nNo random effects\n")
    }
    if(length(x$beta)) {
        cat("\nLocation coefficients:\n")
        print(x$beta, ...)
    } else {
        cat("\nNo location coefficients\n")
    }
    if(length(x$zeta)) {
        cat("\nScale coefficients:\n")
        print(x$zeta, ...)
    } else {
        cat("\nNo Scale coefficients\n")
    }
    if(x$estimLambda > 0) {
        cat("\nLink coefficient:\n")
        print(x$lambda)
    }
    if(length(x$xi) > 0) {
        cat("\nThreshold coefficients:\n")
        print(x$Alpha, ...)
        if(x$threshold != "flexible") {
            cat("\nThresholds:\n")
            print(x$Theta, ...)
        }
    }
    cat("\nlog-likelihood:", format(x$logLik, nsmall=2), "\n")
    cat("AIC:", format(-2*x$logLik + 2*x$edf, nsmall=2), "\n")
    if(nzchar(mess <- naprint(x$na.action))) cat("(", mess, ")\n", sep="")
    invisible(x)
}

vcov.clmm2 <- function(object, ...)
{
    if(is.null(object$Hessian)) {
        stop("Model needs to be fitted with Hess = TRUE")
    }
    dn <- names(object$coefficients)
    structure(solve(object$Hessian), dimnames = list(dn, dn))
}

summary.clmm2 <- function(object, digits = max(3, .Options$digits - 3),
                         correlation = FALSE, ...)
{
    estimStDev <- !("sdFixed" %in% names(as.list(object$call)))
    edf <- object$edf
    coef <- with(object,
                 matrix(0, edf-estimStDev, 4,
                        dimnames =
                         list(names(coefficients[seq_len(edf-estimStDev)]),
                        c("Estimate", "Std. Error", "z value", "Pr(>|z|)"))))
    coef[, 1] <- object$coefficients[seq_len(edf-estimStDev)]
    if(is.null(object$Hessian)) {
      stop("Model needs to be fitted with Hess = TRUE")
    }
    vc <- try(vcov(object), silent = TRUE)
    if(inherits(vc, "try-error")) {
        warning("Variance-covariance matrix of the parameters is not defined")
        coef[, 2:4] <- NaN
        if(correlation) warning("Correlation matrix is unavailable")
        object$condHess <- NaN
    } else {
        sd <- sqrt(diag(vc))
        coef[, 2] <- sd[seq_len(edf - estimStDev)]
        object$condHess <-
            with(eigen(object$Hessian, only.values = TRUE),
                 abs(max(values) / min(values)))
        coef[, 3] <- coef[, 1]/coef[, 2]
        coef[, 4] <- 2*pnorm(abs(coef[, 3]), lower.tail=FALSE)
        if(correlation)
            object$correlation <- (vc/sd)/rep(sd, rep(object$edf, object$edf))
    }
    object$coefficients <- coef
    object$digits <- digits
    varMat <- matrix(c(object$stDev^2, object$stDev),
                     nrow = length(object$stDev), ncol=2)
    rownames(varMat) <- names(object$stDev)
    colnames(varMat) <- c("Var", "Std.Dev")
    object$varMat <- varMat
    class(object) <- "summary.clmm2"
    object
}

print.summary.clmm2 <- function(x, digits = x$digits, signif.stars =
                              getOption("show.signif.stars"), ...)
{
    if(x$nAGQ >=2)
        cat(paste("Cumulative Link Mixed Model fitted with the adaptive",
                  "Gauss-Hermite \nquadrature approximation with",
                  x$nAGQ ,"quadrature points\n\n"))
    else if(x$nAGQ <= -1)
      cat(paste("Cumulative Link Mixed Model fitted with the",
                "Gauss-Hermite \nquadrature approximation with",
                abs(x$nAGQ) ,"quadrature points"), "\n\n")
    else
        cat("Cumulative Link Mixed Model fitted with the Laplace approximation\n",
            fill=TRUE)
    if(!is.null(cl <- x$call)) {
        cat("Call:\n")
        dput(cl, control=NULL)
    }
    if(length(x$stDev)) {
        cat("\nRandom effects:\n")
        print(x$varMat, ...)
    } else {
        cat("\nNo random effects\n")
    }
### OPTION: Should the number of obs. and the number of groups be
### displayed as in clmm?
    coef <- format(round(x$coefficients, digits=digits))
    coef[,4] <- format.pval(x$coefficients[, 4])
    p <- length(x$beta); nxi <- length(x$xi)
    k <- length(x$zeta); u <- x$estimLambda
    if(p > 0) {
        cat("\nLocation coefficients:\n")
        print(coef[nxi + 1:p, , drop=FALSE],
              quote = FALSE, ...)
    } else {
        cat("\nNo location coefficients\n")
    }
    if(k > 0) {
      cat("\nScale coefficients:\n")
      print(coef[(nxi+p+1):(nxi+p+k), , drop=FALSE],
            quote = FALSE, ...)
    } else {
      cat("\nNo scale coefficients\n")
    }
    if(x$estimLambda > 0) {
        cat("\nLink coefficients:\n")
        print(coef[(nxi+p+k+1):(nxi+p+k+u), , drop=FALSE],
              quote = FALSE, ...)
    }
    if(nxi > 0) {
        cat("\nThreshold coefficients:\n")
        print(coef[1:nxi, -4, drop=FALSE], quote = FALSE, ...)
    }

    cat("\nlog-likelihood:", format(x$logLik, nsmall=2), "\n")
    cat("AIC:", format(-2*x$logLik + 2*x$edf, nsmall=2), "\n")
    cat("Condition number of Hessian:", format(x$condHess, nsmall=2), "\n")
    if(nzchar(mess <- naprint(x$na.action))) cat("(", mess, ")\n", sep="")
    if(!is.null(correl <- x$correlation)) {
        cat("\nCorrelation of Coefficients:\n")
        ll <- lower.tri(correl)
        correl[ll] <- format(round(correl[ll], digits))
        correl[!ll] <- ""
        print(correl[-1, -ncol(correl)], quote = FALSE, ...)
    }
    invisible(x)
}

## ranef.clmm2 <- function(x) {
##     x$ranef
## }

## Trace <- function(iter, stepFactor, val, maxGrad, par, first=FALSE) {
##     t1 <- sprintf(" %3d:  %-5e:    %.3f:   %1.3e:  ",
##                   iter, stepFactor, val, maxGrad)
##     t2 <- formatC(par)
##     if(first)
##         cat("iter:  step factor:     Value:     max|grad|:   Parameters:\n")
##     cat(t1, t2, "\n")
## }

gauss.hermite <- function (n)
{
    n <- as.integer(n)
    if (n < 0)
        stop("need non-negative number of nodes")
    if (n == 0)
        return(list(nodes = numeric(0), weights = numeric(0)))
    i <- 1:n
    i1 <- i[-n]
    muzero <- sqrt(pi)
    a <- rep(0, n)
    b <- sqrt(i1/2)

    A <- rep(0, n * n)
    A[(n + 1) * (i1 - 1) + 2] <- b
    A[(n + 1) * i1] <- b
    dim(A) <- c(n, n)
    vd <- eigen(A, symmetric = TRUE)
    w <- rev(as.vector(vd$vectors[1, ]))
    w <- muzero * w^2
    x <- rev(vd$values)
    list(nodes = x, weights = w)
}

profile.clmm2 <-
    function(fitted, alpha = 0.01, range, nSteps = 20, trace = 1, ...)
{
    estimStDev <- !("sdFixed" %in% names(as.list(fitted$call)))
    if(!estimStDev) ##  || is.null(fitted$Hessian))
        fitted <- update(fitted, Hess = TRUE, sdFixed = NULL)
    MLogLik <- fitted$logLik
    MLstDev <- fitted$stDev
    if(missing(range) && is.null(fitted$Hessian))
        stop("'range' should be specified or model fitted with 'Hess = TRUE'")
    if(missing(range) && !is.null(fitted$Hessian)) {
        range <- log(fitted$stDev) + qnorm(1 - alpha/2) *
            c(-1, 1) * sqrt(vcov(fitted)[fitted$edf, fitted$edf])
        range <- exp(range)
        pct <- paste(round(100*c(alpha/2, 1-alpha/2), 1), "%")
        ci <- array(NA, dim = c(1, 2),
                    dimnames = list("stDev", pct))
        ci[] <- range
    }
    stopifnot(all(range > 0))
    logLik <- numeric(nSteps)
    stDevSeq <- seq(min(range), max(range), length.out = nSteps)
    if(trace) message("Now profiling stDev with ", nSteps,
                      " steps: i =")
    if(trace) cat(1, "")
    rho <- update(fitted, Hess = FALSE, sdFixed = min(range))
    logLik[1] <- rho$logLik
    start <- as.vector(rho$coefficients)

    for(i in 2:nSteps){
        if(trace) cat(i, "")
        rho <- update(rho, sdFixed = stDevSeq[i], start = start)
        logLik[i] <- rho$logLik
        start <- as.vector(rho$coefficients)
    }
    if(trace) cat("\n")

    if(any(logLik > fitted$logLik))
        warning("Profiling found a better optimum,",
                "  so original fit had not converged")
    sgn <- 2*(stDevSeq > MLstDev) -1
    Lroot <- sgn * sqrt(2) * sqrt(-logLik + MLogLik)
    res <- data.frame("Lroot" = c(0, Lroot),
                      "stDev" = c(MLstDev, stDevSeq))
    res <- res[order(res[,1]),]
    if(!all(diff(res[,2]) > 0))
        warning("likelihood is not monotonically decreasing from maximum,\n",
                "  so profile may be unreliable for stDev")
    val <- structure(list(stDev = res), original.fit = fitted)
    if(exists("ci", inherits=FALSE)) attr(val, "WaldCI") <- ci
    class(val) <- c("profile.clmm2", "profile")
    val
}

confint.profile.clmm2 <-
    function(object, parm = seq_along(Pnames), level = 0.95, ...)
{
    Pnames <- names(object)
    confint.profile.clm2(object, parm = parm, level = level, ...)
}

plot.profile.clmm2 <-
    function(x, parm = seq_along(Pnames), level = c(0.95, 0.99),
             Log = FALSE, relative = TRUE, fig = TRUE, n = 1e3, ...,
             ylim = NULL)
{
    Pnames <- names(x)
    plot.profile.clm2(x, parm = parm, level = level, Log = Log,
                     relative = relative, fig = fig,
                     n = n, ...,  ylim = ylim)
}

update.clmm2 <-
    function(object, formula., location, scale, nominal, ...,
             evaluate = TRUE)
{
    call <- object$call
    if (is.null(call))
        stop("need an object with call component")
    extras <- match.call(expand.dots = FALSE)$...
    if (!missing(location))
        call$location <-
            update.formula(formula(attr(object$location, "terms")),
                                   location)
    if (!missing(scale))
        call$scale <-
            if(!is.null(object$scale))
                update.formula(formula(attr(object$scale, "terms")), scale)
            else
                scale

    if (!missing(nominal))
        call$nominal <-
            if(!is.null(object$nominal))
                update.formula(formula(attr(object$nominal, "terms")), nominal)
            else
                nominal

    if (length(extras)) {
        existing <- !is.na(match(names(extras), names(call)))
        for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
        if (any(!existing)) {
            call <- c(as.list(call), extras[!existing])
            call <- as.call(call)
        }
    }
    if (evaluate)
        eval(call, parent.frame())
    else call
}

Try the ordinal package in your browser

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

ordinal documentation built on Nov. 17, 2022, 1:06 a.m.