#' Linear regression
#'
#' Runs an OLS regression not unlike \code{\link{lm}}
#'
#' @param y response vector (1 x n)
#' @param X covariate matrix (p x n) with no intercept
#'
#' @return A list with 4 elements: coefficients, vcov, sigma, df
#'
#' @examples
#' data(mtcars)
#' X <- as.matrix(mtcars[, c("cyl", "disp", "hp")])
#' y <- mtcars[, "mpg"]
#' linmodEst(y, X)
#'
#' @export
#'
linmodEst <- function(x, y) {
## CC: crossprod or a QR decomposition (as in the original version) are more efficient
coef <- solve(t(x) %*% x) %*% t(x) %*% y
## degrees of freedom and standard deviation of residuals
df <- nrow(x) - ncol(x)
sigma2 <- sum((y - x %*% coef) ^ 2) / df
## compute sigma^2 * (x’x)^-1
vcov <- sigma2 * solve(t(x) %*% x)
colnames(vcov) <- rownames(vcov) <- colnames(x)
list(
coefficients = coef,
vcov = vcov,
sigma = sqrt(sigma2),
df = df
)
}
linmod <- function(x, ...)
UseMethod("linmod")
linmod.default <- function(x, y, ...) {
x <- as.matrix(x)
y <- as.numeric(y)
est <- linmodEst(x, y)
est$fitted.values <- as.vector(x %*% est$coefficients)
est$residuals <- y - est$fitted.values
est$call <- match.call()
class(est) <- "linmod"
return(est)
}
print.linmod <- function(x, ...) {
cat("Call:\n")
print(x$call)
cat("\nCoefficients:\n")
print(x$coefficients)
}
linmod.formula <- function(formula, data = list(), ...) {
mf <- model.frame(formula = formula, data = data)
x <- model.matrix(attr(mf, "terms"), data = mf)
y <- model.response(mf)
est <- linmod.default(x, y, ...)
est$call <- match.call()
est$formula <- formula
return(est)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.