我正在写一个函数,它接受一个 full
和a reduced
glm
对象,用于总结感兴趣的变量的交互结果。varofint
和互动变量 interaction_var
(通过执行 lrtest
并使用 svycontrast
在...上 full
对象提取结果 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)})))
我是不是缺少一个更简单的方法来提取参考类别?
这是一个用于这个任务的函数。
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"
有点晚,但 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
# a b c
# 0.00000000 0.07017833 -0.05891015
d$interaction_var[1]
# a
# 0
names(d$interaction_var[1])
# [1] "a"