How to predict and extract R Squared with .lm.fit?-CodePudding

As the title suggest, I have seen some user mentioned that .lm.fit() functions has an advantage of more speed than a regular lm(), but when i look deeper at .lm.fit() in help, it is supposed to be a fitter functions, it returns a set of list instead of a model, which makes me to think is it still possible to extract components like R squared, Adj R Squared, and lastly do a predict() out of it?

Below is sample data and executions:

test_dat <- data.frame(y = rnorm(780, 20, 10))
for(b in 1:300){
  name_var <- paste0("x",b)
  test_dat[[name_var]] <- rnorm(780, 0.01 * b, 5)
}

tic()
obj_lm <- lm(y ~ ., data = test_dat)
print(class(obj_lm))
print(summary(obj_lm)$r.squared)
print(summary(obj_lm)$adj.r.squared)
predict(obj_lm)
toc() #approximately 0.4 seconds

tic()
datm <- as.matrix(test_dat)
obj_lm_fit <- .lm.fit(cbind(1,datm[,-1]), datm[,1])
print(class(obj_lm_fit))
toc() #approximately 0.2 seconds

CodePudding user response：

Functions predict and resid are generic and since .lm.fit returns an object of class "list", all you have to do is to write methods implementing the definitions of what you want. Below are methods to compute fitted values, residuals and R^2.

set.seed(2023)    # make the results reproducible
test_dat <- data.frame(y = rnorm(780, 20, 10))
for(b in 1:300){
  name_var <- paste0("x",b)
  test_dat[[name_var]] <- rnorm(780, 0.01 * b, 5)
}

obj_lm <- lm(y ~ ., data = test_dat)

datm <- as.matrix(test_dat)
obj_lm_fit <- .lm.fit(cbind(1,datm[,-1]), datm[,1])

#------------------------------------------------------------------------
# the methods for objects of class "list"
#
fitted.list <- function(object, X) {
  X %*% object$coefficients
}
resid.list <- residuals.list <- function(object, X, y) {
  y_fitted <- fitted(object, X)
  y - y_fitted
}
rsquared <- function(x, ...) UseMethod("rsquared")
rsquared.default <- function(x, ...) {
  summary(x)$r.squared
}
rsquared.list <- function(object, X, y) {
  e <- resid.list(object, X, y)
  1 - sum(e^2)/sum( (y - mean(y))^2 )
}

rsquared(obj_lm_fit, cbind(1,datm[,-1]), datm[,1])
#> [1] 0.3948863
rsquared(obj_lm)
#> [1] 0.3948863

^{Created on 2023-01-03 with reprex v2.0.2}