emmeans 在计算转换后的结果变量的置信区间时是否使用 lm_robust 聚类稳健标准误差?
Does emmeans use lm_robust cluster-robust standard errors when calculating confidence intervals for transformed outcome variables?
我正在使用 emmeans
包检查两个连续预测变量之间的相互作用。我正在使用 estimatr
包中的 lm_robust()
来执行线性回归并获得聚类稳健标准误差。结果变量居中并缩放为 SD 单位方差。例如:
fit <- lm_robust(scale(Y) ~ X1 * X2 + X3 + X4, data = mydata, cluster = school, se_type = 'CR2')
然后我可以使用类似于以下的代码在 X2
的三个级别执行成对对比或可视化线条:
emmip(fit, X2 ~ X1, CIs = TRUE, at = list(X2 = c(mean(X2) - sd(X2),
mean(X2),
mean(X2) + sd(X2))))
我不希望将结果变量反向转换为其原始比例。
我的问题是 emmeans
是否使用聚类稳健标准误差来计算其报告的置信区间或 p 值,这种行为是否取决于结果变量是否在其原始尺度上或转变? estimatr
包创建者 website suggests that lm_robust
objects can be used with emmeans
, but I can't see lm_robust
listed as a supported model on the "Models supported by emmeans" vignette page 或包文档的简短示例。
我相信 lm_robust 对象是 lm 的扩展,因此它使用 lm 的 emmeans 支持。反过来,这意味着估计值是通过 coef(model) 获得的,而它们的 SE 是使用 vcov(model) 得出的。因此,如果 vcov() returns 您需要稳健的方差,emmeans 将使用它们。
对于大多数转换,它将按照转换小插图中的描述进行工作。特别是,指定 type = "response" 会导致对估计值和置信限进行反向转换,保留 P 值,并通过 delta 方法计算 SE(但 不会 在 CI 和测试中使用)。
附加信息
首先,我发现lm_robust
并没有继承自lm
;相反,estimatr 包包含了它自己对 emmeans 的支持。没有给出很多细节,但是 estimatr 的开发者必须相信所提供的必须是合适的。
scale()
转换不是内置的,因为它很复杂。只是说我们使用了 "scale"
并不像说它是 "log"
那么简单,比如说,因为要处理 scale()
结果,我们需要知道用什么来居中和划分结果.
解决方法是创建emmeans()
及其亲属需要反转转换的对象;这是 stats::make.link()
或 emmeans::make.tran()
返回的函数列表。这是一个可以达到这个目的的函数:
make.scaletran = function(y, ...) {
sy = scale(y, ...)
if(is.null(m <- attr(sy, "scaled:center")))
m = 0
if(is.null(s <- attr(sy, "scaled:scale")))
s = 1
list(
linkfun = function(mu) (mu - m) / s,
linkinv = function(eta) s * eta + m,
mu.eta = function(eta) s,
valideta = function(eta) TRUE,
name = paste0("scale(", signif(m, 3), ", ", signif(s, 3), ")")
)
}
要使用它,您需要手动指定转换,因为它不会被自动检测到。这是一个使用 R 中已有的 warpbreaks
数据的示例:
> warp.lmr = lm_robust(scale(breaks) ~ tension, cluster = wool,
+ se_type = 'CR2', data = warpbreaks)
> tran = make.scaletran(warpbreaks$breaks)
> emmeans(warp.lmr, "tension", tran = tran)
tension emmean SE df lower.CL upper.CL
L 0.624 0.619 51 -0.618 1.8666
M -0.133 0.181 51 -0.497 0.2301
H -0.491 0.219 51 -0.930 -0.0517
Results are given on the scale(28.1, 13.2) (not the response) scale.
Confidence level used: 0.95
> emmeans(warp.lmr, "tension", tran = tran, type = "response")
tension response SE df lower.CL upper.CL
L 36.4 8.17 51 20.0 52.8
M 26.4 2.39 51 21.6 31.2
H 21.7 2.89 51 15.9 27.5
Confidence level used: 0.95
Intervals are back-transformed from the scale(28.1, 13.2) scale
OP 中 emmip()
调用的代码不正确,因为它使用 emmeans()
的规范,而不是 emmip()
。
我会考虑在以后的更新中 emmeans::make.tran()
添加这个比例转换选项。
我正在使用 emmeans
包检查两个连续预测变量之间的相互作用。我正在使用 estimatr
包中的 lm_robust()
来执行线性回归并获得聚类稳健标准误差。结果变量居中并缩放为 SD 单位方差。例如:
fit <- lm_robust(scale(Y) ~ X1 * X2 + X3 + X4, data = mydata, cluster = school, se_type = 'CR2')
然后我可以使用类似于以下的代码在 X2
的三个级别执行成对对比或可视化线条:
emmip(fit, X2 ~ X1, CIs = TRUE, at = list(X2 = c(mean(X2) - sd(X2),
mean(X2),
mean(X2) + sd(X2))))
我不希望将结果变量反向转换为其原始比例。
我的问题是 emmeans
是否使用聚类稳健标准误差来计算其报告的置信区间或 p 值,这种行为是否取决于结果变量是否在其原始尺度上或转变? estimatr
包创建者 website suggests that lm_robust
objects can be used with emmeans
, but I can't see lm_robust
listed as a supported model on the "Models supported by emmeans" vignette page 或包文档的简短示例。
我相信 lm_robust 对象是 lm 的扩展,因此它使用 lm 的 emmeans 支持。反过来,这意味着估计值是通过 coef(model) 获得的,而它们的 SE 是使用 vcov(model) 得出的。因此,如果 vcov() returns 您需要稳健的方差,emmeans 将使用它们。
对于大多数转换,它将按照转换小插图中的描述进行工作。特别是,指定 type = "response" 会导致对估计值和置信限进行反向转换,保留 P 值,并通过 delta 方法计算 SE(但 不会 在 CI 和测试中使用)。
附加信息
首先,我发现lm_robust
并没有继承自lm
;相反,estimatr 包包含了它自己对 emmeans 的支持。没有给出很多细节,但是 estimatr 的开发者必须相信所提供的必须是合适的。
scale()
转换不是内置的,因为它很复杂。只是说我们使用了 "scale"
并不像说它是 "log"
那么简单,比如说,因为要处理 scale()
结果,我们需要知道用什么来居中和划分结果.
解决方法是创建emmeans()
及其亲属需要反转转换的对象;这是 stats::make.link()
或 emmeans::make.tran()
返回的函数列表。这是一个可以达到这个目的的函数:
make.scaletran = function(y, ...) {
sy = scale(y, ...)
if(is.null(m <- attr(sy, "scaled:center")))
m = 0
if(is.null(s <- attr(sy, "scaled:scale")))
s = 1
list(
linkfun = function(mu) (mu - m) / s,
linkinv = function(eta) s * eta + m,
mu.eta = function(eta) s,
valideta = function(eta) TRUE,
name = paste0("scale(", signif(m, 3), ", ", signif(s, 3), ")")
)
}
要使用它,您需要手动指定转换,因为它不会被自动检测到。这是一个使用 R 中已有的 warpbreaks
数据的示例:
> warp.lmr = lm_robust(scale(breaks) ~ tension, cluster = wool,
+ se_type = 'CR2', data = warpbreaks)
> tran = make.scaletran(warpbreaks$breaks)
> emmeans(warp.lmr, "tension", tran = tran)
tension emmean SE df lower.CL upper.CL
L 0.624 0.619 51 -0.618 1.8666
M -0.133 0.181 51 -0.497 0.2301
H -0.491 0.219 51 -0.930 -0.0517
Results are given on the scale(28.1, 13.2) (not the response) scale.
Confidence level used: 0.95
> emmeans(warp.lmr, "tension", tran = tran, type = "response")
tension response SE df lower.CL upper.CL
L 36.4 8.17 51 20.0 52.8
M 26.4 2.39 51 21.6 31.2
H 21.7 2.89 51 15.9 27.5
Confidence level used: 0.95
Intervals are back-transformed from the scale(28.1, 13.2) scale
OP 中 emmip()
调用的代码不正确,因为它使用 emmeans()
的规范,而不是 emmip()
。
我会考虑在以后的更新中 emmeans::make.tran()
添加这个比例转换选项。