Как установить формулу R программно

Question

Я пытаюсь использовать пакет 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, ...)

Нужна помощь мастера, пожалуйста. Какую форму хочет этот знаменатель?

Answer 1

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. Я надеюсь, что вы найдете их полезными!