时间戳

首页 > 解决方案 > 如何以编程方式设置 R 公式

r - 如何以编程方式设置 R 公式

问题描述

我正在尝试使用 R 包“ipw”进行逆概率加权。我有一些名为“covar.1”、“covar.2”、“covar.3”的列......所以我想为它们制定一个公式。从上一个问题中，我可以使用glm,matchit和其他功能。但是使用ipw，它不起作用。如果我手动复制并粘贴print(f1)分母处的输出，它会起作用，所以我尝试不使用它，as.formula但它仍然不起作用。要重现，运行

library(ipw)

betaz <- c(0.75, -0.5,  0.25)
betay <- c(0.5, 1.0, -1.5)

X <- matrix(rnorm(3 * 250), 250)
ps <- pnorm(X %*% betaz)
Z <- rbinom(250, 1, ps)
epsilon <- rnorm(250, 0.0, 0.5)

Y0 <- X %*% betay + epsilon
Y1 <- X %*% betay + 0.5 + epsilon
Y <- Y0 * (1 - Z) + Y1 * Z
df <- data.frame(id = seq(250), covar = X, group = Z, metric = Y)
print(df[1:10,])

cols <- colnames(df)
covars <- cols[grep("covar", colnames(df))]
f <- as.formula(paste('group','~', paste(covars, collapse="+")))
psmodel <- glm(f, family = binomial(), data=df)
pscore <- psmodel$fitted.values

f1 <- as.formula(paste('~', paste(covars, collapse="+")))
print(f1)
weightmodel <- ipwpoint(
  exposure = group, family = "binomial", link = "logit", 
  denominator = f1,
  data = df, trunc = .01
)

随着as.formula，它抱怨object 'groupf1' not found。不知道为什么要进行这种连接。基本上我需要一种使用变量动态设置 f1 的方法。

从回溯我看到源代码

glm(formula = eval(
  parse(
    text = paste(
      deparse(tempcall$exposure, width.cutoff = 500), 
      deparse(tempcall$denominator, width.cutoff = 500), sep = ""))), 
  family = lf, data = data, na.action = na.fail, ...)

需要R大师的帮助。这个分母想要什么形式？

标签： rformula

解决方案

ipw以这样一种方式编写，很难动态输入公式。这是我必须编写WeightIt包的动机之一，它具有相同的功能（除少数情况外）。此外，在我的cobalt包中，有一个函数f.build()可以从其输入中创建一个公式。

您可以将代码的最后几行替换为以下内容：

f1 <- f.build("group", covars)
w.out <- weightit(f1, data = df, estimand = "ATE")
w.out2 <- trim(w.out, .01, lower = TRUE)

这f1是您的公式，由f.build. 这样，您可以在第一个参数中循环遍历多个处理变量。第二个参数可以是协变量名称的向量，也可以是协变量本身的 data.frame。w.out是weightit包含由 估计的权重的对象weightit()。默认为逻辑回归，但可以更改。（我注意到真正的治疗倾向是使用概率模型生成的，可以在 中请求weightit。link = "probit"）

似乎您想截断第一个和第 99 个百分位的权重，这就是这样trim做的。默认情况下，它只修剪最高的权重，所以我设置lower = TRUE为也修剪较低的权重。一般来说，您应该在修剪前检查协变量平衡和权重的可变性，以防未修剪的权重足够。cobalt旨在评估平衡并与WeightIt. 以下是如何评估weightit对象的平衡：

bal.tab(w.out, un = TRUE)

您也可以比较修剪和未修剪的重量：

bal.tab(f1, data = df, un = TRUE, 
         weights = list(untrimmed = w.out$weights,
                         trimmed = w.out2$weights))

当您准备好估计您的治疗效果时，您可以从weightit对象中提取权重。我使用jtools包来获得强大的标准误差，这是 PS 加权所必需的：

w1 <- w.out$weights
jtools::summ(lm(metric ~ group, data = df, weights = w1),
             robust = TRUE, confint = TRUE)

有大量关于WeightIt和的文档cobalt。我希望你觉得它们有用！

推荐阅读