将“聚类函数”应用于一系列线性模型答案

【问题标题】：Applying 'clustering functions' to a series of linear models将“聚类函数”应用于一系列线性模型
【发布时间】：2019-05-03 17:23:06
【问题描述】：

我想遍历线性模型列表，并使用 vcovCL 函数将“聚集”标准误差应用于每个模型。我的目标是尽可能高效地执行此操作（我正在跨数据框的许多列运行线性模型）。我的问题是尝试在匿名函数中指定其他参数。下面我模拟一些假数据。区域代表我的横截面尺寸；月代表我的时间维度（在 4 个月内观察到 5 个单位）。变量int 是干预发生时的虚拟变量。

df <- data.frame(
  precinct = c( rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4), rep(5, 4) ),
  month = rep(1:4, 5),
  crime = rnorm(20, 10, 5),
  int = c(c(0, 1, 1, 0), rep(0, 4), rep(0, 4), c(1, 1, 1, 0), rep(0, 4))
  )

df[1:10, ]

outcome <- df[3]
est <- lapply(outcome, FUN = function(x) { lm(x ~ as.factor(precinct) + as.factor(month) + int, data = df) })

se <- lapply(est, function(x) { sqrt(diag(vcovCL(x, cluster = ~ precinct + month))) })

在vcovCL 函数内添加cluster 参数时收到以下错误消息。

Error in eval(expr, envir, enclos) : object 'x' not found

据我估计，唯一的解决方法是索引数据框，即df$，然后指定“聚类”变量。这可以通过在函数调用中为df 指定一个附加参数来实现吗？ 这段代码效率高吗？

我想，也许以公式的方式指定模型方程是一种更好的方法。

任何想法/cmets 总是有帮助的 :)

【问题讨论】：

我很抱歉没有提供更好的例子。我将在数据框中使用多个列。我只使用了一个用于说明目的。
是的，它非常有帮助。你的意思是“三角形”还是“复选标记”？还是习惯了一切。
“向上”三角形表示赞成票。您会为您认为特别有用的帖子投票。帖子被赞的人将获得 10 点声望。如果人们用朝下的三角形对您的帖子投反对票，您也可以松开它们。复选标记用于接受答案。

标签： r regression lapply coding-efficiency covariance-matrix

【解决方案1】：

这是一种检索多个模型的聚集标准错误的方法：

library(sandwich)

# I am going to use the same model three times to get the "sequence" of linear models. 
mod <- lm(crime ~ as.factor(precinct) + as.factor(month) + int, data = df)

# define function to retrieve standard errors:
robust_se <- function(mod) {sqrt(diag(vcovCL(mod, cluster = list(df$precinct, df$month))))}

# apply function to all models:
se <- lapply(list(mod, mod, mod), robust_se)

如果您想调整整个输出，以下可能会有所帮助：

library(lmtest)
adj_stats <- function(mod) {coeftest(mod, vcovCL(mod, cluster = list(df$precinct, df$month)))}

adjusted_models <- lapply(list(mod, mod, mod), adj_stats)

解决多列问题：

如果您难以在多列上运行线性模型，以下内容可能会有所帮助。以上所有内容都将保持不变，除了您将模型列表传递给lapply。

首先，让我们在这里使用这个数据框：

df <- data.frame(
  precinct = c( rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4), rep(5, 4) ),
  month = rep(1:4, 5),
  crime = rnorm(20, 10, 5),
  crime2 = rnorm(20, 10, 5),
  crime3 = rnorm(20, 10, 5),
  int = c(c(0, 1, 1, 0), rep(0, 4), rep(0, 4), c(1, 1, 1, 0), rep(0, 4))
)

让我们定义结果列：

outcome_columns <- c("crime", "crime2", "crime3")

现在，让我们对每个结果进行回归：

models <- lapply(outcome_columns, 
         function(outcome) lm( eval(parse(text = paste0(outcome, " ~ as.factor(precinct) + as.factor(month) + int"))), data = df) )

然后你就打电话

adjusted_models <- lapply(models, adj_stats)

关于效率：

上面的代码很有效，因为它很容易调整并且可以快速编写。对于大多数用例，它会非常好。为了计算效率，请注意您的设计矩阵在所有情况下都是相同的，即通过预先计算公共元素（例如inv(X'X)*X'），您可以节省一些计算。但是，您会失去许多内置函数的便利性。

【讨论】：