Calculate stderr, t-value, p-value, predict value for linear regression

I use the data from MatrixModels:::lm.fit.sparse example. I built a custom function summary_sparse to perform a summary for this model. All matrix operations are performed with Matrix package. Results are compared with dense type model.

Note lm.fit.sparse have to be evaluated with method = "chol" to get proper results.

Functions:

summary_sparse <- function(l, X) {
  
  XXinv <- Matrix::chol2inv(Matrix::chol(Matrix::crossprod(X)))
  se <- sqrt(Matrix::diag(XXinv*sum(l$residuals**2)/(nrow(X)-ncol(X))))
  ts <- l$coef/se
  pvals <- 2*c(1 - pnorm(abs(ts)))
  
  list(coef = l$coef, se = se, t = ts, p = pvals)
  
}

predict_sparse <- function(X, coef) {
  X %*% coef
}

Application:

dd <- expand.grid(a = as.factor(1:3),
                  b = as.factor(1:4),
                  c = as.factor(1:2),
                  d= as.factor(1:8))
n <- nrow(dd <- dd[rep(seq_len(nrow(dd)), each = 10), ])

set.seed(17)

dM <- cbind(dd, x = round(rnorm(n), 1))
## randomly drop some
n <- nrow(dM <- dM[- sample(n, 50),])
dM <- within(dM, { A <- c(2,5,10)[a]
                   B <- c(-10,-1, 3:4)[b]
                   C <- c(-8,8)[c]
                   D <- c(10*(-5:-2), 20*c(0, 3:5))[d]
   Y <- A + B + A*B + C + D + A*D + C*x + rnorm(n)/10
   wts <- sample(1:10, n, replace=TRUE)
   rm(A,B,C,D)
})

X <- Matrix::sparse.model.matrix( ~ (a+b+c+d)^2 + c*x, data = dM)

Xd <- as(X,"matrix")

fmDense <- lm(dM[,"Y"]~Xd-1)

ss <- summary(fmDense)

r1 <- MatrixModels:::lm.fit.sparse(X, y = dM[,"Y"], method = "chol")

f <- summary_sparse(r1, X)

all.equal(do.call(cbind, f), ss$coefficients, check.attributes = F)
#TRUE

all.equal(predict_sparse(X, r1$coef)@x, predict(fmDense), check.attributes = F, check.names=F)
#TRUE

CLICK HERE to find out more related problems solutions.

Leave a Comment

Your email address will not be published.

Scroll to Top