使用R语言gformula包执行干预分析时遇非函数应用错误及NA警告

问题场景 #

我用R语言的gformula包，基于含缺失值的非平衡面板数据，估计自然治疗史对因变量newbusiness的影响，运行代码后出现错误及警告，查阅包文献和help文档仍无法解决。

运行代码 #

gformula(obs_data = dt, 
         outcome_type = 'continuous_eof', 
         id= 'country_iso3', 
         time_name = 'recode_year', 
         seed = 1234, 
         outcome_name = 'newbusiness', 
         covnames = c('totalnumber', 'nkillwound', 'GDPppp'), 
         covtypes = c('zero-inflated normal', 'zero-inflated normal' ,"normal"), 
         intvars = c('totalnumber', 'nkillwound'), 
         interventions = c('natural', 'natural'), 
         histvars = list(c('totalnumber', 'nkillwound', 'GDPppp'), c('totalnumber', 'nkillwound', 'GDPppp')), 
         histories = c(lagged, cumavg), 
         covparams = list(covmodels = c(totalnumber ~ lag1_totalnumber + cumavg_totalnumber, 
                                        nkillwound ~ lag1_nkillwound + cumavg_totalnumber, 
                                        GDPppp ~ lag1_GDPppp + lag1_totalnumber + lag1_nkillwound + cumavg_totalnumber + cumavg_nkillwound)), 
         ymodel = newbusiness ~ nkillwound + GDPppp + lag1_totalnumber + lag1_nkillwound)

报错信息 #

Error in intervention[[i]][[1]](newdf, pool, intvar[i], intvals = NULL, :
attempt to apply non-function
In addition: There were 19 warnings (use warnings() to see them)

其中17条警告为：

In stats::rnorm(x, mean, est_sd) : NAs produced

修复方案 #

1. 修正histories参数格式 #

histories需要传入字符串形式的历史函数名称，原代码中lagged和cumavg未加引号，R会将其视为未定义变量，触发"非函数"错误。修改为：

histories = c("lagged", "cumavg")

2. 预处理衍生变量 #

gformula不会自动生成滞后项、累积平均项，需提前用dplyr或data.table处理数据，示例代码：

library(dplyr)
dt <- dt %>%
  group_by(country_iso3) %>%
  mutate(
    lag1_totalnumber = lag(totalnumber),
    lag1_nkillwound = lag(nkillwound),
    lag1_GDPppp = lag(GDPppp),
    cumavg_totalnumber = cummean(totalnumber),
    cumavg_nkillwound = cummean(nkillwound)
  ) %>%
  ungroup()

3. 处理缺失值 #

NA警告源于数据中的缺失值，非平衡面板的缺失会导致模拟阶段生成无效值，可选择：

• 用na.omit(dt)删除含缺失值的行（注意会丢失样本）
• 用mice等包对缺失值进行插补

4. 完善零膨胀协变量模型 #

对于zero-inflated normal类型的协变量，需同时指定计数部分和零膨胀部分的模型，修改covparams为列表形式：

covparams = list(
  covmodels = list(
    totalnumber = list(
      count = totalnumber ~ lag1_totalnumber + cumavg_totalnumber,
      zero = totalnumber ~ lag1_totalnumber + cumavg_totalnumber
    ),
    nkillwound = list(
      count = nkillwound ~ lag1_nkillwound + cumavg_totalnumber,
      zero = nkillwound ~ lag1_nkillwound + cumavg_totalnumber
    ),
    GDPppp = GDPppp ~ lag1_GDPppp + lag1_totalnumber + lag1_nkillwound + cumavg_totalnumber + cumavg_nkillwound
  )
)

内容的提问来源于stack exchange，提问作者flxflks