0

我试图了解如何使用公式对象。假设我想制作自己的 2SLS 函数，并且想将我正在使用的对象分为 4 个主要组：y = response；X = 外生变量；E = 内生变量；Z = 仪器。

我希望能够构建这些对象，而无需不必要地制作额外的数据副本（例如，大 N 和大量仪器会使内存使用/时间成本过高）。我还想考虑数据中的 NA。

让我们使用类似于 felm 的公式语法（我尝试查看那里的解析代码，但无法遵循它）。

frml = y ~ x1 + x2 + x3*x4 | (e1 | e2 ~ z1 + z2)

library(Formula)
N = 12 # be divisible by 6
data = data.frame(y=rnorm(N), x1=rnorm(N), x2=rnorm(N), x3=rnorm(N),
                  x4=factor(rep(1:2, N/2)), e1=rnorm(N), e2=rnorm(N),
                  z1=rnorm(N), z2=factor(rep(1:3, N/3)))
data[2,'y'] = data[3,'x1'] = data[4,'e1'] = data[5,'z2'] = NA

parse_frml = function(frml, data, subset=NULL) {
    frml = as.Formula(frml)
    # does not take into account NAs at all
    y = model.part(frml, data=data, subset=subset, lhs=1)
    # does not take into account NAs in other variables (y, Z, E)
    X = model.matrix(frml, data=data, subset=subset, lhs=0, rhs=1)
    Z = model.matrix(frml, data=data, subset=subset, lhs=0, rhs=2)
    #E =  # I can't figure this out at all
    return(list(y=y, X=X, E=E, Z=Z))
}

现在，我可以做类似的事情

mf = model.frame(frml, data=data, subset=subset, lhs=1, rhs=1)

这将考虑 y 和 X 中的 NA，但忽略 E 和 Z。此外，这会将数据复制到 mf 中，然后再次复制到 y 和 X 中。

所以，我有 2 个问题和 1 个约束

  1. 我如何获得E？（第二个等式的 LHS 矩阵）
  2. 如何在所有矩阵中考虑 frml 使用的数据中的 NA？
  3. 同时最小化数据副本的数量（理想情况下只是复制到矩阵中）

更一般地说，什么是理解公式、公式、术语等的好资源？我还没有发现，例如公式库包文档非常有用。

4

1 回答 1

0

这并不完美，但它有效。遗憾的是，几乎没有关于如何在 R 代码中实际处理和操作公式的信息。我的解决方案取决于 formula.tools

library(formula.tools)
parse_frml = function(frml, data, subset=NULL) {
    frml = as.Formula(frml)
    vars = all.vars(frml)
    other_vars = c(all.vars(formula(frml, lhs=1, rhs=1)),
                   rhs.vars(formula(frml, lhs=0, rhs=2)))
    e_vars = setdiff(vars, other_vars)
    valid = which(complete.cases(data[, vars]))
    if (!is.null(subset)) {
        if (class(subset) == 'logical') {
            subset = which(subset)
        }
        valid = intersect(valid, subset)
    }
    y = model.part(frml, data=data[valid,], lhs=1)
    X = model.matrix(frml, data=data[valid,], lhs=0, rhs=1)
    Z = model.matrix(frml, data=data[valid,], lhs=0, rhs=2)
    E = data.matrix(data[valid, e_vars])
    return(list(y=y, X=X, E=E, Z=Z))
}

我怀疑每次都使用有效的子集数据是相当昂贵的。但是在上面的测试中，它似乎起作用了。

于 2015-12-04T05:13:13.497 回答