理论上
首先，您可以拟合具有多个 LHS 的线性模型 https://stackoverflow.com/q/39262534/4891738.
其次，显式数据分割并不是分组回归的唯一方法（或推荐方法）。看R回归分析：分析特定种族的数据 https://stackoverflow.com/q/51407989/4891738 and R：为每个类别建立单独的模型 https://stackoverflow.com/q/43189497/4891738。所以将你的模型构建为cbind(x1, x2, x3, x4) ~ day * subject where subject是因子变量。
最后，由于您有许多因子级别并使用大型数据集，lm是不可行的。考虑使用speedglm::speedlm with sparse = TRUE, or MatrixModels::glm4 with sparse = TRUE.
事实上
Neither speedlm nor glm4正在积极开发中。 （在我看来）它们的功能很原始。
Neither speedlm nor glm4支持多个 LHS 作为lm。所以你需要做 4 个独立的模型x1 ~ day * subject to x4 ~ day * subject反而。
这两个包背后的逻辑不同sparse = TRUE.
-
speedlm首先使用标准构建密集设计矩阵model.matrix.default，然后使用is.sparse调查是否稀疏。如果为 TRUE，则后续计算可以使用稀疏方法。
-
glm4 uses model.Matrix构造设计矩阵，可以直接构造稀疏矩阵。
所以，这并不奇怪speedlm和一样糟糕lm在这个稀疏性问题中，并且glm4是我们真正想要使用的。
glm4没有一套完整、有用的通用函数来分析拟合模型。您可以通过以下方式提取系数、拟合值和残差coef, fitted and residuals，但必须自己计算所有统计数据（标准误差、t 统计量、F 统计量等）。这对于熟悉回归理论的人来说不是什么大事，但仍然相当不方便。
glm4仍然期望你使用最好的模型公式，以便可以构造最稀疏的矩阵。常规的~ day * subject确实不太好。我可能应该稍后就这个问题设立一个问答。基本上，如果你的公式有截距并且因子是对比的，你就会失去稀疏性。这是我们应该使用的：~ 0 + subject + day:subject.
测试与glm4
## use chinsoon12's data in his answer
set.seed(0L)
nSubj <- 200e3
nr <- 1e6
DF <- data.frame(subject = gl(nSubj, 5),
day = 3:7,
y1 = runif(nr),
y2 = rpois(nr, 3),
y3 = rnorm(nr),
y4 = rnorm(nr, 1, 5))
library(MatrixModels)
fit <- glm4(y1 ~ 0 + subject + day:subject, data = DF, sparse = TRUE)
在我的 1.1GHz Sandy Bridge 笔记本电脑上大约需要 6 ~ 7 秒。让我们提取它的系数：
b <- coef(fit)
head(b)
# subject1 subject2 subject3 subject4 subject5 subject6
# 0.4378952 0.3582956 -0.2597528 0.8141229 1.3337102 -0.2168463
tail(b)
#subject199995:day subject199996:day subject199997:day subject199998:day
# -0.09916175 -0.15653402 -0.05435883 -0.02553316
#subject199999:day subject200000:day
# 0.02322640 -0.09451542
你可以做B <- matrix(b, ncol = 2)，因此第一列是截距，第二列是斜率。
我的想法：我们可能需要更好的包来进行大回归
Using glm4这里并没有提供有吸引力的优势青松12的data.table解决方案 https://stackoverflow.com/a/51940464/4891738因为它基本上也只是告诉你回归系数。它也比慢一点data.table方法，因为它计算拟合值和残差。
简单回归分析不需要适当的模型拟合例程。我有一些关于如何在这种回归上做一些奇特的事情的答案，比如数据框中所有变量之间的快速成对简单线性回归 https://stackoverflow.com/q/51953709/4891738其中还给出了如何计算所有统计数据的详细信息。但当我写这个答案时，我正在思考一些关于大型回归问题的一般问题。我们可能需要更好的包，否则就没有尽头进行案例到案例的编码。
回复OP
speedglm::speedlm(x1 ~ 0 + subject + day:subject, data = df, sparse = TRUE) gives 错误：无法分配大小为 74.5 Gb 的向量
是的，因为它有一个坏处sparse logic.
MatrixModels::glm4(x1 ~ day * subject, data = df, sparse = TRUE)给出错误Cholesky(crossprod(from), LDL = FALSE)：internal_chm_factor：Cholesky 分解失败
这是因为某些数据只有一个数据subject。您至少需要两个数据才能拟合一条线。这是一个示例（在密集设置中）：
dat <- data.frame(t = c(1:5, 1:9, 1),
f = rep(gl(3,1,labels = letters[1:3]), c(5, 9, 1)),
y = rnorm(15))
级别“c”f只有一个数据/行。
X <- model.matrix(~ 0 + f + t:f, dat)
XtX <- crossprod(X)
chol(XtX)
#Error in chol.default(XtX) :
# the leading minor of order 6 is not positive definite
Cholesky 分解无法解决秩缺陷模型。如果我们使用lm的 QR 分解，我们会看到NA系数。
lm(y ~ 0 + f + t:f, dat)
#Coefficients:
# fa fb fc fa:t fb:t fc:t
# 0.49893 0.52066 -1.90779 -0.09415 -0.03512 NA
我们只能估计水平“c”的截距，而不能估计斜率。
请注意，如果您使用data.table解决方案，你最终会得到0 / 0计算该级别的斜率时，最终结果为NaN.
更新：现已提供快速解决方案
查看快速分组简单线性回归和一般配对简单线性回归 https://stackoverflow.com/q/51996021/4891738.
