计算线性回归的 stderr、t 值、p 值、预测值
Calculate stderr, t-value, p-value, predict value for linear regression
我正在用 MatrixModels:::lm.fit.sparse
和 MatrixModels::glm4
(也是稀疏的)拟合线性模型。
但是,这些功能仅 return coeff
、residuals
和 fitted.values
。
get/calculate 其他值(例如 stderr
、t-value
、p-value
、predict
值的最快和最简单的方法是什么?
我使用 MatrixModels:::lm.fit.sparse
示例中的数据。
我构建了一个自定义函数 summary_sparse
来执行此模型的汇总。
所有矩阵运算都使用 Matrix
包执行。
结果与密集型模型进行比较。
注意 lm.fit.sparse
必须用 method = "chol"
进行评估才能获得正确的结果。
函数:
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
}
申请:
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
我正在用 MatrixModels:::lm.fit.sparse
和 MatrixModels::glm4
(也是稀疏的)拟合线性模型。
但是,这些功能仅 return coeff
、residuals
和 fitted.values
。
get/calculate 其他值(例如 stderr
、t-value
、p-value
、predict
值的最快和最简单的方法是什么?
我使用 MatrixModels:::lm.fit.sparse
示例中的数据。
我构建了一个自定义函数 summary_sparse
来执行此模型的汇总。
所有矩阵运算都使用 Matrix
包执行。
结果与密集型模型进行比较。
注意 lm.fit.sparse
必须用 method = "chol"
进行评估才能获得正确的结果。
函数:
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
}
申请:
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