我正在对需要转换的数据进行线性回归,因为我正在使用Box-Cox电源转换,然后使用反向转换来使用原始比例编写报告。我一直试图用emmeans
包来做这个,然后我按照emmeans
package vignette中描述的步骤进行,但是,我发现估计均值的汇总结果与未转换的数据完全不相似。实际上,输出根本没有变换。
以下是使用emmeans
包中的示例的可重现示例:
require(emmeans)
# Fit a model using an oddball transformation:
bctran <- make.tran("boxcox", 0.368)
warp.bc <- with(bctran,
lm(linkfun(breaks) ~ wool * tension, data = warpbreaks))
# Obtain back-transformed LS means:
emmeans(warp.bc, ~ tension | wool, type = "response")
# Fit a model without transformation:
warp <- lm(breaks ~ wool * tension, data = warpbreaks)
# Obtain LS means:
emmeans(warp, ~ tension | wool)
返回:
> emmeans(warp.bc, ~ tension | wool, type = "response")
wool = A:
tension emmean SE df lower.CL upper.CL
L 8.07 0.419 48 7.23 8.92
M 5.91 0.419 48 5.07 6.75
H 5.94 0.419 48 5.10 6.79
wool = B:
tension emmean SE df lower.CL upper.CL
L 6.45 0.419 48 5.61 7.29
M 6.53 0.419 48 5.69 7.37
H 5.22 0.419 48 4.38 6.07
Confidence level used: 0.95
> emmeans(warp, ~ tension | wool)
wool = A:
tension emmean SE df lower.CL upper.CL
L 44.6 3.65 48 37.2 51.9
M 24.0 3.65 48 16.7 31.3
H 24.6 3.65 48 17.2 31.9
wool = B:
tension emmean SE df lower.CL upper.CL
L 28.2 3.65 48 20.9 35.6
M 28.8 3.65 48 21.4 36.1
H 18.8 3.65 48 11.4 26.1
Confidence level used: 0.95
实际上估计的张力平均值:L应该是42.37,使用公式计算:
> origin + (1 + param * pmax(eta))^(1/param)
> 0 + (1 + 0.368 * pmax(8.07))^(1/0.368)
[1] 42.37179
有什么我缺少或不正确理解?
Hmmmm。我转载了这个问题。我不确定是什么问题,但到目前为止,我可以说bctran
本身是有序的:
> emm = as.data.frame(emmeans(warp.bc, ~tension|wool))
> emm
tension wool emmean SE df lower.CL upper.CL
1 L A 8.074761 0.4192815 48 7.231739 8.917783
2 M A 5.911710 0.4192815 48 5.068688 6.754732
3 H A 5.942335 0.4192815 48 5.099313 6.785357
4 L B 6.449869 0.4192815 48 5.606847 7.292891
5 M B 6.531085 0.4192815 48 5.688063 7.374107
6 H B 5.224939 0.4192815 48 4.381917 6.067961
> bctran$linkinv(emm$emmean)
[1] 42.42263 23.10060 23.32407 27.22827 27.88877 18.43951
所以这些反向转换的EMM是有序的。我将跟踪代码,看看为什么结果不会被反向转换。
几个月前我发现了一个修订版的逻辑错误,如果转换是字符(例如,"log"
),它可以正常工作,但如果它是一个列表(例如,你的bctran
),它将被忽略。我在下一个版本中修复了该错误以推送到github site(版本> = 1.3.3.0999902),修复将在下一个CRAN更新(版本> 1.3.3)中。
> emmeans(warp.bc, ~ tension | wool)
wool = A:
tension emmean SE df lower.CL upper.CL
L 8.07 0.419 48 7.23 8.92
M 5.91 0.419 48 5.07 6.75
H 5.94 0.419 48 5.10 6.79
wool = B:
tension emmean SE df lower.CL upper.CL
L 6.45 0.419 48 5.61 7.29
M 6.53 0.419 48 5.69 7.37
H 5.22 0.419 48 4.38 6.07
Results are given on the Box-Cox (lambda = 0.368) (not the response) scale.
Confidence level used: 0.95
> emmeans(warp.bc, ~ tension | wool, type = "response")
wool = A:
tension response SE df lower.CL upper.CL
L 42.4 4.48 48 34.0 52.0
M 23.1 3.05 48 17.5 29.8
H 23.3 3.07 48 17.7 30.0
wool = B:
tension response SE df lower.CL upper.CL
L 27.2 3.38 48 20.9 34.6
M 27.9 3.44 48 21.5 35.3
H 18.4 2.65 48 13.6 24.3
Confidence level used: 0.95
Intervals are back-transformed from the Box-Cox (lambda = 0.368) scale
请注意,即使没有反向转换,也会有一个注释。您的结果中根本没有注释的事实是一个提示。