多部分公式解析:处理 NA 和最小化对象副本

multipart formula parsing: handling NAs and minimizing object copies

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

我希望能够构建这些对象,而无需不必要地制作额外的数据副本(例如,大 N 和大量仪器会使内存成本过高 usage/time)。我还想考虑整个数据中的 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. 同时尽量减少数据副本的数量(理想情况下只是复制到矩阵中)

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

这并不完美,但它确实有效。遗憾的是几乎没有关于如何在 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))
}

我怀疑每次都使用 valid 对数据进行子集化是相当昂贵的。但在上面的测试中,它似乎起作用了。