从 glm 模型中提取参考类别的最佳方法？答案

【问题标题】：Best way to extract a reference category from a glm model?从 glm 模型中提取参考类别的最佳方法？
【发布时间】：2014-01-12 05:24:14
【问题描述】：

我正在编写一个函数，该函数接受 full 和 reduced glm 对象来汇总感兴趣变量 varofint 和交互变量 interaction_var 的交互结果（通过执行 lrtest并在full 对象上使用svycontrast 来提取varofint 的每个级别interaction_var 的结果）。样本数据：

x <- data.frame(outcome=rbinom(100,1,.3),varofint=rnorm(100), interaction_var=sample(letters[1:3],100,replace=TRUE))

reduced <- glm(outcome~varofint+interaction_var,data=x)
full <- glm(outcome~varofint*interaction_var,data=x)

我想知道为所述 (full) glm 模型提取参考类别的最佳方法。我显然可以做类似的事情

levels(full$data$interaction_var)[1]

但是在给定contrasts 参数的输入的情况下，这将是一种提取参考类别的“安全”方法吗？看起来，如果可以选择 SAS 对比度，此方法可能会产生interactionv_var 的级别，这不是模型中用作参考类别的级别。以下会更安全吗？

mf <- model.frame(full)
setdiff(rownames(contrasts(mf[, "interaction_var"])), colnames(contrasts(mf[, "interaction_var"])))

或类似

names(which(apply(contrasts(mf[, "interaction_var"]),1,function(.v){all(.v==0)})))

我是否缺少一种更简单的方法来提取参考类别？

【问题讨论】：

如果没有参考类别怎么办？参考类别仅适用于治疗对比。
好的，这是一个好的开始。因此，如果相应的交互变量没有处理对比，该函数应该返回错误...

标签： r lm

【解决方案1】：

这里是这个任务的一个函数：

refCat <- function(model, var) {
  cs <- attr(model.matrix(model), "contrasts")[[var]]
  if (is.character(cs)) {
    if (cs == "contr.treatment")
      ref <- 1
    else stop("No treatment contrast")
  }  
  else {
    zeroes <- !cs
    ones <- cs == 1
    stopifnot(all(zeroes | ones))
    cos <- colSums(ones)
    stopifnot(all(cos == 1))
    ros <- rowSums(ones)
    stopifnot(sum(!ros) == 1 && sum(ros) != ncol(cs))
    ref <- which(!ros)
  }
  return(levels(model$data[[var]])[ref])  
}

如果变量var 未表示为治疗对比，该函数将停止。

例子：

refCat(reduced, "interaction_var")
# [1] "a"
refCat(full, "interaction_var")
# [1] "a"

【讨论】：

谢谢，这似乎有效。在什么情况下is.character(cs) 是FALSE？
@Michael cs 可以是带有对比函数名称的字符串，例如contr.treatment，也可以是数值对比矩阵。
如果使用了 relevel() 就不行了，对吧？
@Helen 我没试过。

【解决方案2】：

有点晚了，但dummy.coef() 可以工作...其输出的每个变量元素中的第一个值是参考类别。

# R 4.0.0 data.frame() does not produce factors
x <- data.frame(
  outcome = rbinom(100, 1, .3),
  varofint = rnorm(100),
  interaction_var = sample(letters[1:3], 100, replace = TRUE),
  stringsAsFactors = TRUE
)
reduced <- glm(outcome ~ varofint + interaction_var, data = x)
full <- glm(outcome ~ varofint * interaction_var, data = x)

d <- dummy.coef(full)
d
# Full coefficients are 
#                                                                 
# (Intercept):                    0.310136                        
# varofint:                    -0.07247677                        
# interaction_var:                       a           b           c
#                               0.00000000  0.07017833 -0.05891015
# varofint:interaction_var:              a           b           c
#                               0.00000000 -0.14824179 -0.04123618

d$interaction_var
#           a           b           c 
#  0.00000000  0.07017833 -0.05891015 
d$interaction_var[1]
# a 
# 0 
names(d$interaction_var[1])
# [1] "a"

【讨论】：