r - 在 R 中使用 logitmfx 计算稳健集群标准误差时出错

Question

我希望运行一个 logit 回归来预测家庭规模和年龄平均值的边际效应，以及二元指标(个人是否是移民、是否有健康保险或吸烟)对患高血压的预测概率的影响.
该数据来自聚类调查，我希望在输出中包含稳健的聚类标准误差。
但是，当我添加代码以包含强大的集群 SE 时，我收到一个错误，即不再找到回归中的变量，我不确定为什么。任何建议都会很棒!谢谢。

AGE       IMMIGRANT     FAMSIZE     HLTH_INS    HYPERTEN   SMOKE    PSU
<int>       <dbl>         <int>       <dbl>       <dbl>     <dbl>  <int>
40           0              2          1            0         0      2
23           0              2          1            0         0      1
24           0              2          1            0         0      2
18           0              3          1            1         0      2
30           0              2          1            0         0      2
33           1              6          0            0         0      1

#or if this is an easier output to reproduce:
structure(list(AGE = c(40L, 23L, 24L, 18L, 30L, 33L, 32L, 63L, 
22L, 24L), IMMIGRANT = c(0, 0, 0, 0, 0, 1, 0, 0, 0, 1), FAMSIZE = c(2L, 
2L, 2L, 3L, 2L, 6L, 2L, 1L, 2L, 1L), HLTH_INS = c(1, 1, 1, 1, 
1, 0, 1, 1, 1, 0), HYPERTEN = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0), 
    SMOKE = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1), PSU = c(2L, 1L, 
    2L, 2L, 2L, 1L, 2L, 2L, 1L, 2L)), row.names = c(NA, -10L), class = "data.frame")


#The regression works without adjusting for clustered SE 
logit<-logitmfx(HYPERTEN~scale(AGE)+IMMIGRANT+scale(FAMSIZE)+HLTH_INS+
                 SMOKE,data=sample,
                atmean=TRUE,robust=T)


#However, when I add in the code to cluster SE I receive the error: "Error in scale(AGE) : object 'AGE' not found" 
logit<-logitmfx(HYPERTEN~scale(AGE)+IMMIGRANT+scale(FAMSIZE)+HLTH_INS+
                 SMOKE,data=sample,
                atmean=TRUE,robust=T,clustervar1="PSU", clustervar2=NULL,!is.null("PSU")) 

Answer 1

尝试使用源代码复制函数的步骤，Steffen Moritz 的解决方案确实应该有效。问题出现了，因为第一次logitmfx立即调用另一个函数 logitmfxest .
此函数具有相同的参数，但它还有以下代码:

if(!is.null(clustervar1)){
    if(is.null(clustervar2)){
      if(!(clustervar1 %in% names(data))){
        stop("clustervar1 not in data.frame object")
      }    
      data = data.frame(model.frame(formula, data, na.action=NULL),data[,clustervar1])
      names(data)[dim(data)[2]] = clustervar1
      data=na.omit(data)
    }
    if(!is.null(clustervar2)){
      if(!(clustervar1 %in% names(data))){
        stop("clustervar1 not in data.frame object")
      }    
      if(!(clustervar2 %in% names(data))){
        stop("clustervar2 not in data.frame object")
      }    
      data = data.frame(model.frame(formula, data, na.action=NULL),
                        data[,c(clustervar1,clustervar2)])
      names(data)[c(dim(data)[2]-1):dim(data)[2]] = c(clustervar1,clustervar2)
      data=na.omit(data)
    }
  }

在您的情况下，以下代码将被激活:

if(!is.null(clustervar1)){
    if(is.null(clustervar2)){  
      data = data.frame(model.frame(formula, data, na.action=NULL),data[,clustervar1])
      names(data)[dim(data)[2]] = clustervar1
      data=na.omit(data)
    }
  }

这将“数据”重新定义为在 model.frame 上构建的 data.frame。但是模型框架使用您公式中的名称，所以突然第 2 列被称为 规模。年龄。 第 3 列称为 规模。FAMSIZE。 .
这是个大问题因为该函数随后调用了广义线性模型:

fit = glm(formula, data=data, family = binomial(link = "logit"), x=T, 
        start = start, control = control) 

它使用包含 scale(AGE) 和 scale(FAMSIZE) 的原始公式，但使用带有重命名列的新数据框。
所以在输入之前缩放应该可以工作。事实上，正如 Steffen 提到的，任何其他函数都会导致同样的错误，因为它们会在 model.frame 时产生类似的列重命名。叫做。