没问题！我来帮你用统计模型复现你用aggregate算出的分组均值结果。先从构造示例数据开始，方便我们一起验证：

# 构造你提供的数据集
dat <- data.frame(
  EXIST = c(0, 0, 0, 1, 0),
  DATE = c("10/2015", "01/2016", "01/2014", "05/2015", "11/2015"),
  VAR1 = c(6, 6, 5, 5, 6),
  VAR2 = c(4, 4, 4, 4, 4)
)

首先，先运行你原来的aggregate代码，确认基准结果：

# 你的原分组均值计算代码
ag <- data.frame(as.matrix(aggregate(EXIST ~ VAR1 + VAR2, data = dat, function(x) c(mean = mean(x)))))
print(ag)
# 输出结果：
#   VAR1 VAR2 EXIST.mean
# 1    5    4        0.5
# 2    6    4        0.0

接下来，我们用**线性模型（lm）**来复现这个结果。因为你要的是VAR1和VAR2组合分组下EXIST的均值，而EXIST是0/1二分类变量，分组均值本质上就是每个组的“事件发生概率”，线性模型可以直接拟合出这个均值。

方法1：使用交叉项去掉全局截距 #

我们可以把VAR1和VAR2转为因子，拟合包含它们交叉项且无全局截距的模型，这样每个系数就对应一个分组的均值：

# 拟合线性模型
model <- lm(EXIST ~ factor(VAR1)*factor(VAR2) - 1, data = dat)
# 提取系数并整理成数据框
model_results <- data.frame(
  VAR1 = as.integer(sub("factor\\(VAR1\\)", "", names(coef(model)))),
  VAR2 = as.integer(sub(".*factor\\(VAR2\\)", "", names(coef(model)))),
  EXIST_mean = coef(model)
)
print(model_results)
# 输出结果和aggregate完全一致：
#   VAR1 VAR2 EXIST_mean
# 1    5    4        0.5
# 2    6    4        0.0

方法2：使用interaction函数 #

如果你觉得交叉项写法有点繁琐，也可以用interaction()函数直接生成分组变量，同样去掉全局截距：

# 用interaction构造分组变量
model2 <- lm(EXIST ~ interaction(VAR1, VAR2) - 1, data = dat)
# 整理结果
model2_results <- data.frame(
  group = names(coef(model2)),
  EXIST_mean = coef(model2)
)
# 拆分group列得到VAR1和VAR2
model2_results$VAR1 <- as.integer(sapply(strsplit(as.character(model2_results$group), "\\."), "[", 1))
model2_results$VAR2 <- as.integer(sapply(strsplit(as.character(model2_results$group), "\\."), "[", 2))
# 调整列顺序
model2_results <- model2_results[, c("VAR1", "VAR2", "EXIST_mean")]
print(model2_results)
# 同样得到和aggregate一致的结果

补充：二分类变量的广义线性模型（glm） #

如果你的数据量更大或者想考虑二分类变量的分布特性，也可以用glm的identity链接来拟合，结果和lm完全一致：

glm_model <- glm(EXIST ~ factor(VAR1)*factor(VAR2) - 1, data = dat, family = binomial(link = "identity"))
print(coef(glm_model))
# 系数同样是0.5和0.0，和之前的结果匹配

这样就完美复现了你用aggregate得到的分组均值结果啦！

内容的提问来源于stack exchange，提问作者Maximilian