#' 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
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.