R/S3methods.R

In baruuum/jars: Just Another Robust Standard errors package

#' Extract variance-covariance matrix
#' 
#' @param object a \code{robustSE} class object
#' @param se logical; if TRUE, returns vector of standard errors
#' @param ... further options to be passed (not used)
#' @return returns the estimated covariance matrix
#' @export
vcov.robustSE = function(object, se = FALSE, ...) { 
    
    res = object$cov 
    
    if (se) res = sqrt(diag(res))
    
    return(res)
    
}

#' Extract regression coefficients
#' @param object a \code{robustSE} class object
#' @param ... further options (not used)
#' @return returns the estimated regression coefficients
#' @export
coef.robustSE = function(object, ...) { 
    
    return(object$coefficients)
    
}


#' Print function for \code{robustSE} object
#'
#' @param x an object of class \code{robustSE}
#' @param digits number of digits to print
#' @param stars logical; if TRUE, prints stars
#' @param ... further options to be passed to print
#' @return void
#' @export
print.robustSE = function (
          x,
          digits = 3L,
          stars = FALSE,
          ...)
{
    cat("\nCall:\n", paste(deparse(x$call), sep = "\n",
        collapse = "\n"), "\n", sep = "")


    if (ncol(x$coefficients) == 0L) {

        cat("\nNo Coefficients\n")

    } else {

        c_type = x$cov_type
        
        if (c_type[1] == "robust") {
            
            cat(
                paste0(
                    "\nRobust covariance matrix of type ", 
                    c_type[2], " is used"
                )
            )
            
        } else if (c_type[1] == "cluster") {
            
            cat(
                paste0(
                    "\nCluster-robust covariance matrix of type ",
                    c_type[2], " is used\n",
                    "(Number of clusters: ", x$m, ")"
                )
            )
            
        } else {
            
            stop("cannot find vcov type")
            
        }
        
        cat("\n\nCoefficients:\n")
        coefs = x$coefficients
        if (stars) {
            
            printCoefmat(
                coefs, 
                digits = digits, 
                signif.stars = TRUE,
                dig.tst = digits,
                na.print = "NA", 
                ...
            )
        
        } else {
            
            printCoefmat(
                coefs[, c(1, 2)], 
                digits = digits, 
                signif.stars = FALSE,
                dig.tst = digits,
                na.print = "NA", 
                ...
            )
            
        }

    }
    cat("\nN = ", x$n, ",  K = ", x$k, sep = "")
    cat(",  RMSE = ", formatC(x$rmse, digits = digits), ",  ", sep = "")
    cat("R-squared = ", formatC(x$r.squared, digits = digits), sep = "")
    if (x$cluster_warn) {
        
        cat("\nNOTE: observations used to calculate standard errors and coefficients differ due to missing values in the clustering variable.")
        
    } else cat("")
    cat("\n")
    invisible(x)
}

#' New summary.lm function
#'
#' Adding the options \code{robust} and \code{cluster} to the summary.lm function. Loading the package will mask the summary.lm function.
#'
#' @param object an object of class "lm", usually, a result of a call to lm.
#' @param correlation logical; if TRUE, the correlation matrix of the estimated parameters is returned and printed.
#' @param symbolic.cor logical. If TRUE, print the correlations in a symbolic form (see symnum) rather than as numbers.
#' @param robust logical; if TRUE, heteroskedasticity robust standard errors are used in the hypothesis tests
#' @param cluster integer vector indicating cluster membership
#' @param type type of robust standard errors to calculate
#' @param ... further options to pass to summary
#' @details detailed descriptions regarding the difference between the standard errors can be found by \code{?summary_robust}. The summary.lm function should work as usuall when calling without specifying any of the options \code{robust} or \code{cluster}. If called with these options, on the other hand, a \code{robustSE} class object will be returned.
#' @export
summary.lm = function (object,
                       correlation = FALSE,
                       symbolic.cor = FALSE,
                       robust = FALSE,
                       cluster = NULL,
                       type = "HC1",
                       ...)
{
    
    # add conditional statement to lead to summary_robust
    if (robust || !is.null(cluster)) {

        res = .summary_robust(object, robust, cluster, type)
        return(res)

    }

    # original summary.lm function definition below ---------

    z <- object
    p <- z$rank
    rdf <- z$df.residual
    if (p == 0) {
        r <- z$residuals
        n <- length(r)
        w <- z$weights
        if (is.null(w)) {
            rss <- sum(r^2)
        }
        else {
            rss <- sum(w * r^2)
            r <- sqrt(w) * r
        }
        resvar <- rss/rdf
        ans <- z[c("call", "terms", if (!is.null(z$weights)) "weights")]
        class(ans) <- "summary.lm"
        ans$aliased <- is.na(coef(object))
        ans$residuals <- r
        ans$df <- c(0L, n, length(ans$aliased))
        ans$coefficients <- matrix(NA_real_, 0L, 4L, dimnames = list(NULL,
            c("Estimate", "Std. Error", "t value",
                "Pr(>|t|)")))
        ans$sigma <- sqrt(resvar)
        ans$r.squared <- ans$adj.r.squared <- 0
        ans$cov.unscaled <- matrix(NA_real_, 0L, 0L)
        if (correlation)
            ans$correlation <- ans$cov.unscaled
        return(ans)
    }
    if (is.null(z$terms))
        stop("invalid 'lm' object:  no 'terms' component")
    if (!inherits(object, "lm"))
        warning("calling summary.lm(<fake-lm-object>) ...")
    
    ## Small changes here since stats:::qr.lm throws a warning on CRAN.
    ## Literally copying stats:::qr.lm into this part ...
    if (is.null(Qr <- object$qr)) 
        stop("lm object does not have a proper 'qr' component.\n Rank zero or should not have used lm(.., qr=FALSE).")
    
    n <- NROW(Qr$qr)
    if (is.na(z$df.residual) || n - p != z$df.residual)
        warning("residual degrees of freedom in object suggest this is not an \"lm\" fit")
    r <- z$residuals
    f <- z$fitted.values
    w <- z$weights
    if (is.null(w)) {
        mss <- if (attr(z$terms, "intercept"))
            sum((f - mean(f))^2)
        else sum(f^2)
        rss <- sum(r^2)
    }
    else {
        mss <- if (attr(z$terms, "intercept")) {
            m <- sum(w * f/sum(w))
            sum(w * (f - m)^2)
        }
        else sum(w * f^2)
        rss <- sum(w * r^2)
        r <- sqrt(w) * r
    }
    resvar <- rss/rdf
    if (is.finite(resvar) && resvar < (mean(f)^2 + var(c(f))) *
        1e-30)
        warning("essentially perfect fit: summary may be unreliable")
    p1 <- 1L:p
    R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
    se <- sqrt(diag(R) * resvar)
    est <- z$coefficients[Qr$pivot[p1]]
    tval <- est/se
    ans <- z[c("call", "terms", if (!is.null(z$weights)) "weights")]
    ans$residuals <- r
    ans$coefficients <- cbind(Estimate = est, `Std. Error` = se,
        `t value` = tval, `Pr(>|t|)` = 2 * pt(abs(tval),
            rdf, lower.tail = FALSE))
    ans$aliased <- is.na(z$coefficients)
    ans$sigma <- sqrt(resvar)
    ans$df <- c(p, rdf, NCOL(Qr$qr))
    if (p != attr(z$terms, "intercept")) {
        df.int <- if (attr(z$terms, "intercept"))
            1L
        else 0L
        ans$r.squared <- mss/(mss + rss)
        ans$adj.r.squared <- 1 - (1 - ans$r.squared) * ((n -
            df.int)/rdf)
        ans$fstatistic <- c(value = (mss/(p - df.int))/resvar,
            numdf = p - df.int, dendf = rdf)
    }
    else ans$r.squared <- ans$adj.r.squared <- 0
    ans$cov.unscaled <- R
    dimnames(ans$cov.unscaled) <- dimnames(ans$coefficients)[c(1,
        1)]
    if (correlation) {
        ans$correlation <- (R * resvar)/outer(se, se)
        dimnames(ans$correlation) <- dimnames(ans$cov.unscaled)
        ans$symbolic.cor <- symbolic.cor
    }
    if (!is.null(z$na.action))
        ans$na.action <- z$na.action
    class(ans) <- "summary.lm"
    ans
}