R/hett.R

In hett: Heteroscedastic t-Regression

"tlm" <- function(lform = formula(data), sform =  ~ 1, data = NULL, subset = NULL, contrasts = NULL, na.action = na.fail, start = NULL, control = tlm.control(...), obs = FALSE, estDof = FALSE, ...)
{
  dofse <- NA
  tcall <- match.call()
  sform[3] <- sform[2]
  sform[2] <- lform[2]
  fixedFrameArgs <- list(lform, na.action = na.action, contrasts = contrasts, data = data)
  dispFrameArgs <- list(sform, na.action = na.action, data = data)
  if(!missing(subset)) {
    fixedFrameArgs[["subset"]] <- subset
    dsipFrameArgs[["subset"]] <- subset
  }
  fixedFrame <- do.call("model.frame", fixedFrameArgs)
  dispFrame <- do.call("model.frame", dispFrameArgs)
  X <- model.matrix(lform, fixedFrame)
  y <- model.extract(fixedFrame, "response")
  Z <- model.matrix(sform, dispFrame)
  if(!is.null(start)) {
    if(is.null(namStart <- names(start)))
      stop("Starting values must have names for assignation")
    namList <- c("beta", "lambda", "dof", "omega")
    if(!all(pmatch(namStart, namList, nomatch = 0)))
      stop("Names of starting values must be \"beta\", \"lambda\", \"dof\" or \"omega\"")
  }
  else namStart <- ""
  if((bets <- pmatch("beta", namStart, nomatch = 0)) != 0) {
    fixedCoef <- start[[bets]]
    sqResid <- (y - crossprod(t(X), fixedCoef))^2
  }
  else {
    initFit <- lm.fit(X, y, method = "qr")
    fixedCoef <- coef(initFit)
    sqResid <- resid(initFit)^2
  }
  dispFamily <- Gamma(link = "log")
  if((lams <- pmatch("lambda", namStart, nomatch = 0)) != 0) {
    dispCoef <- start[[lams]]
    sigmaI <- dispFamily$linkinv(crossprod(t(Z), dispCoef))
  }
  else {
    dispFit <- lm.fit(Z, dispFamily$linkfun((sqResid + (sqResid == 0)/6)))
    zLambda <- dispFit$fitted.values
    dispCoef <- coef(dispFit)
    sigmaI <- dispFamily$linkinv(zLambda)
  }
  if((dofs <- pmatch("dof", namStart, nomatch = 0)) != 0){
    dof <- start[[dofs]]
    if(length(dof) != 1)
      stop(paste("initial value for degrees of freedom must be of length 1"))
  }
  else
    dof <- 3
  if((oms <- pmatch("omega", namStart, nomatch = 0)) != 0){
    randCoef <- start[[oms]]
    if(length(randCoef) != length(y))
      stop(paste("Starting values for scale random effects must be the same length as y"))
  }
  else
    randCoef <- (dof + 1)/(dof + sqResid/sigmaI)
  n <- length(y)
  ## control parameters for algorithm
  epsilon <- control$epsilon
  maxit <- control$maxit
  trace <- control$trace
  verbose <- control$verboseLev
  iter <- 0
  initTime <- proc.time()[1]
  ##
  ## iterative scheme in 4 parts : location -> fixed scale -> random scale -> degrees of freedom
  ##
  max.fit <- expression({
    const <-  - n * lgamma(1/2) - n * lgamma(dof/2) - (n/2) * log(dof) + n * lgamma((dof + 1)/2)
    logLik <- const - (sum(log(sigmaI)))/2 - sum(((dof + 1)/2) * (log(1 + (sqResid/sigmaI)/dof)))
    logLikTemp <- logLik + 1
    orthoComp <- apply(crossprod(solve(crossprod(Z)), t(Z)),1,sum)
    orthoCoef <- dispCoef - orthoComp
    while(abs(logLik - logLikTemp) > epsilon && iter < maxit) {
      ##
      ## location parameter likelihood
      ##
      if(!obs) {
        mscale <- (dof + 1)/(dof + 3)
        meanY <-  (y - X %*% fixedCoef)*(randCoef/mscale) + crossprod(t(X), fixedCoef)
        meanFit <- lm.wfit(X, meanY, rep(mscale, length(y))/sigmaI, method = "qr")
      }
      else meanFit <- lm.wfit(X, y, randCoef/sigmaI, method = "qr")
      fitm <- meanFit$fitted.values
      fixedCoef <- coef(meanFit)
      sqResid <- (as.vector(y - crossprod(t(X), fixedCoef)))^2
      if(trace) cat("\nLocation parameters :",
                    format(round(fixedCoef, 4)))
      ##
      ## scale parameter estimation
      ##
      tscale <- (1 - 3/(dof + 3))
      scale <- (1/tscale)
      dispY <- scale * (((sqResid * randCoef)/sigmaI) - rep(1, length(y))) + crossprod(t(Z), dispCoef)
      dispFit <- lm.wfit(Z, dispY, rep(tscale, length(y)), method = "qr")
      dispCoef <- coef(dispFit)
      sigmaI <- dispFamily$linkinv(dispFit$fitted.values)
      orthoI <- as.vector(sigmaI/exp(Z %*% (orthoComp*(2*(log(dof) - log(dof + 1))))))
      if(trace) cat("\nScale parameters :", format(round(
                                                              dispCoef, 4)))
      ##
      ## random component maximum likelihood
      ##
      randCoef <- (dof + 1)/(dof + sqResid/sigmaI)
      if(trace && verbose == 2)
          cat("\nRandom scale effects :", format(round(randCoef, 4)))
      logLikTemp <- logLik
      ##
      ## degrees of freedom estimation
        ##
      if(estDof) {
        opt <- nlm(dof.profile, dof, n = n, sqResid = sqResid, orthoI = orthoI, X = X, Z = Z)
        dof <- opt$estimate
        if(trace)
          cat("\ndegrees of freedom is ", dof)
        ## expected information for degrees of freedom standard error.
        expinfDof <-   n*trigamma(dof/2)/4 - n*trigamma((dof + 1)/2)/4 - n*(dof + 5)/(2 * dof * (dof + 1)*(dof + 3)) - (2/(dof*(dof + 3)*(dof + 1)^2))*(t(rep(1, length(y))) %*% Z %*% orthoComp)
        logLik <- - opt$minimum
        dofse <- 1/sqrt(expinfDof)
      }
      else
        logLik <- const - sum(log(sigmaI))/2 - sum(((dof + 1)/2) * (log(1 + (sqResid/sigmaI)/dof)))
      iter <- iter + 1
      if(trace)
        cat("\nHeteroscedastic t: - 2 x Maximum Likelihood iteration", iter, " : ", format(round(-2 * logLik, 4)), "\n")
    }})
  estMethod <-  eval(max.fit)
  ##
  ## Coverged or not?
  ##
  endTime <- proc.time()[1] - initTime
  if(trace && iter >= maxit) cat("\nMaximum", maxit,
                "iterations reached in", endTime, "\n")
  meanFit$method <- "maximum likelihood"
  random <- randCoef
  meanFit$formula <- lform
  meanFit$terms <- attr(fixedFrame, "terms")
  meanFit$iter <- iter
  meanFit$call <- tcall
  meanFit$residuals <- as.vector(y - crossprod(t(X), fixedCoef))
  dispFit$formula <- sform
  dispFit$terms <- attr(dispFrame, "terms")
  dispFit$iter <- iter
  dispFit$fitted.values <- sigmaI
  dispFit$call <- tcall
  res <- list(loc.fit = meanFit, scale.fit = dispFit, random =
              random, dof = dof, dofse = dofse, iter = iter, logLik =
              logLik, endTime = endTime)
  class(res) <- c("tlm", "glm", "lm")
  res
}

"summary.tlm" <- function(object, correlation = FALSE, ...)
{
  locOut <- object$loc.fit
  scaleOut <- object$scale.fit
  locSum <- tsum(locOut, dispersion = 1, correlation = correlation, ...)
  scaleSum <- tsum(scaleOut, dispersion = 2, correlation = correlation, ...)
  structure(list(loc.summary = locSum, scale.summary = scaleSum, iter =
                 object$iter, dof = object$dof, dofse = object$dofse,
                 logLik = object$logLik, method = locOut$method,
                 estTime = object$endTime), class = c("summary.tlm"))
}

"print.tlm" <- function(x, ...)
  {
    cat("\nCall: ")
    cat(format(x$loc.fit$call), "\n\n")
    cat("Location Coefficients:\n")
    print.default(format(x$loc.fit$coefficients, digits = options()$digits), print.gap = 2,
                  quote = FALSE)
    cat("Scale Coefficients:\n")
    print.default(format(x$scale.fit$coefficients, digits = options()$digits), print.gap = 2,
                  quote = FALSE)
    cat("Degrees of freedom:\n")
    cat(format(x$dof), "\n")
    cat("Log-Likelihood\n")
    cat(format(x$logLik), "\n")
  }

"print.summary.tlm" <- function(x, ...)
{
  cat("Location model :\n")
  print(x$loc.summary, scale = FALSE, ...)
  cat("\nScale Model :\n")
  print(x$scale.summary, scale = TRUE, ...)
  if(is.null(x$dofse))
    cat("\nFixed degrees of freedom parameter: ", format(x$dof))
  else {
    cat("\nEst. degrees of freedom parameter: ", format(x$dof))
    cat("\nStandard error for d.o.f: ", format(x$dofse))
  }
  cat("\nNo. of iterations of model :", format(x$iter), "in", format(
                                                                          round(x$estTime, options()$digits)))
  cat("\nHeteroscedastic t Likelihood :", format(x$logLik), "\n")
}

"tsum" <- function (object, dispersion = NULL, correlation = FALSE, symbolic.cor = FALSE,
            ...)
  ## Function edited from summary.glm()
  ## Sets up the summary information for the location and scale model for the t
  ## Produces the correct standard errors and correlations between coefficients.
{
  Qr <- object$qr
  est.disp <- FALSE
  df.r <- object$df.residual
  if (dispersion == 1)
    est.disp <- TRUE
  p <- object$rank
  p1 <- 1:p
  coef.p <- object$coefficients[Qr$pivot[p1]]
  covmat.unscaled <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
  dimnames(covmat.unscaled) <- list(names(coef.p), names(coef.p))
  covmat <- dispersion * covmat.unscaled
  var.cf <- diag(covmat)
  s.err <- sqrt(var.cf)
  tvalue <- coef.p/s.err
  dn <- c("Estimate", "Std. Error")
  if (!est.disp) {
    pvalue <- 2 * pnorm(-abs(tvalue))
    coef.table <- cbind(coef.p, s.err, tvalue, pvalue)
    dimnames(coef.table) <- list(names(coef.p), c(dn, "z value",
                                                  "Pr(>|z|)"))
  }
  else if (df.r > 0) {
    pvalue <- 2 * pt(-abs(tvalue), df.r)
    coef.table <- cbind(coef.p, s.err, tvalue, pvalue)
    dimnames(coef.table) <- list(names(coef.p), c(dn, "t value",
                                                  "Pr(>|t|)"))
  }
  else {
    coef.table <- cbind(coef.p, Inf)
    dimnames(coef.table) <- list(names(coef.p), dn)
  }
  ans <- c(object[c("call", "terms", "df.residual", "iter")], list(resid = residuals(object), coefficients = coef.table, dispersion = dispersion, df = c(object$rank, df.r), cov.unscaled = covmat.unscaled, cov.scaled = covmat))
  if (correlation) {
    dd <- sqrt(diag(covmat.unscaled))
    ans$correlation <- covmat.unscaled/outer(dd, dd)
    ans$symbolic.cor <- symbolic.cor
  }
  class(ans) <- "tsum"
  return(ans)
}

"print.tsum" <- function (x, digits = max(3, getOption("digits") - 3), symbolic.cor = x$symbolic.cor,
            signif.stars = getOption("show.signif.stars"), scale = TRUE, ...)
{
  cat("\nCall:\n")
  cat(paste(deparse(x$call), sep = "\n", collapse = "\n"),
      "\n\n", sep = "")
  cat("Residuals: \n")
  if (x$df.residual > 5) {
    x$resid <- quantile(x$resid, na.rm = TRUE)
    names(x$resid) <- c("Min", "1Q", "Median", "3Q",
                        "Max")
  }
  print.default(x$resid, digits = digits, na.print = "",
                print.gap = 2)
  cat("\nCoefficients:\n")
  printCoefmat(x$coef, digits = digits, signif.stars = signif.stars,
                ...)
  if(scale)
    cat("\n(Scale parameter taken to be ",
        format(x$dispersion),")\n")
  else
    cat("\n(Scale parameter(s) as estimated below)\n")
  correl <- x$correlation
  if (!is.null(correl)) {
    p <- NCOL(correl)
    if (p > 1) {
      cat("\nCorrelation of Coefficients:\n")
      if (is.logical(symbolic.cor) && symbolic.cor) {
        print(symnum(correl, abbr.colnames = NULL))
      }
      else {
        correl <- format(round(correl, 2), nsmall = 2,
                         digits = digits)
        correl[!lower.tri(correl)] <- ""
        print(correl[-1, -p, drop = FALSE], quote = FALSE)
      }
    }
  }
  cat("\n")
  invisible(x)
}

"tlm.control" <- function(epsilon = 1e-07, maxit = 50, trace = FALSE, verboseLev = 1)
{
  ## tdispersion control parameters.
  if(epsilon <= 0) {
    warning("the value of epsilon supplied is zero or negative;\nthe default value of 1e-7 was used instead"
            )
    epsilon <- 1e-07
  }
  if(maxit < 1) {
    warning("the value of maxit supplied is zero or negative;\nthe default value of 50 was used instead"
            )
    maxit <- 50
  }
  if(verboseLev < 0 || verboseLev > 2)
    warning("the value of the verbose level is not in the expected range. Use \"verboseLev = 1\" for simple output and \"verboseLev = 2\" for advanced output"
            )
  list(epsilon = epsilon, maxit = maxit, trace = trace,
       verboseLev = verboseLev)
}

"dof.profile" <- function(dof, n, sqResid, orthoI, X, Z)
{
  n <-  length(sqResid)
  sigmaI <- as.vector(orthoI*exp(Z %*% crossprod(solve(crossprod(Z)), t(Z)) %*% rep(2*(log(dof) - log(dof + 1)), n)))
  - (- n * lgamma(1/2) - n * lgamma(dof/2) - (n/2) * log(dof) + n * lgamma((dof + 1)/2) -
     (sum(log(sigmaI)))/2 - sum(((dof + 1)/2) * (log(1 + (sqResid/sigmaI)/dof))))
}

"tscore" <- function(..., data = NULL, scale = FALSE)
{
  ## General score test for the t-dispersion model for both scale and location
  listm <- list(...)
  if(is.null(data))
    stop("Score function requires the data to reproduce some model elements")
#  assign("data", data, envir = .GlobalEnv)
  dofse <- unlist(lapply(listm, function(el) el$dofse))
  if(any(is.na(dofse))){
    if(!all(is.na(dofse)))
      stop("All models must have estimated or fixed degrees of freedom\n")
  }
  fullm <- do.call("model.frame", list(listm[[length(listm)]]$loc.fit$formula, data = data))
  y <- model.extract(fullm, "response")
  k <- length(listm)
  ss <- pv <- df <- vector()
  if(scale)
    Xl <- lapply(listm, function(el, data) model.matrix(el$scale.fit$formula,  data = data), data = data)
  else
    Xl <- lapply(listm, function(el, data) model.matrix(el$loc.fit$formula, data = data), data = data)
  k <- length(listm)
  nams <- vector()
  j <- 1
  while(k > 1){
    X2 <- as.matrix(Xl[[k]][, !pmatch(dimnames(Xl[[k]])[[2]], dimnames(Xl[[k - 1]])[[2]], nomatch = 0)])
    random <- listm[[k - 1]]$random
    ssig <- (listm[[k - 1]]$scale.fit$fitted.values)^(1/2)
    Om <- random/(ssig^2)
    if(scale){
      Q <- qr.Q(qr(Xl[[k - 1]]))
      P <- - crossprod(t(Q))
      diag(P) <- 1 + diag(P)
      ss[j] <- ((listm[[k - 1]]$dof + 3)/(2 * listm[[k - 1]]$dof))*(t((resid(listm[[k - 1]]$loc.fit)^2)*Om - 1) %*% X2 %*% solve(t(X2) %*% P %*% X2)
                                                                    %*% t(X2) %*% (Om*(resid(listm[[k - 1]]$loc.fit)^2) - 1))
      df[j] <- listm[[k - 1]]$scale.fit$df.residual - listm[[k]]$scale.fit$df.residual
      nams[j] <- paste(paste(deparse(listm[[k - 1]]$scale.fit$formula[[3]]), "vs", sep = " "), deparse(listm[[k]]$scale.fit$formula[[3]]), sep = " ")
    }
    else {
      Q <- qr.Q(qr(ssig*Xl[[k - 1]]))
      P <- diag(1/ssig^2) - diag(1/ssig) %*% crossprod(t(Q)) %*% diag(1/ssig)
      ss[j] <-  ((listm[[k - 1]]$dof + 3)/(listm[[k - 1]]$dof + 1))*(t(resid(listm[[k - 1]]$loc.fit)*Om) %*% X2 %*% solve(t(X2) %*% P %*% X2)
                                                                     %*% t(X2) %*% (Om*resid(listm[[k - 1]]$loc.fit)))
      df[j] <- listm[[k - 1]]$loc.fit$df.residual - listm[[k]]$loc.fit$df.residual
      nams[j] <- paste(paste(deparse(listm[[k - 1]]$loc.fit$formula[[3]]), "vs", sep = " "), deparse(listm[[k]]$loc.fit$formula[[3]]), sep = " ")
    }
    pv[j] <- 1 - pchisq(ss[j], df[j])
    k <- k - 1
    j <- j + 1
  }
  tab <- cbind(df, ss, pv)
  colnames(tab) <- c("df", "Score Stat.", "P-value")
  rownames(tab) <- nams
  cat("Heteroscedastic t Regression:\n")
  if(scale)
    cat("\nScale score statistic tests:\n\n")
  else
    cat("\nLocation score statistic tests:\n\n")
  printCoefmat(tab)
}