如何从R中有分类处理的回归模型中得到均值比的标准误差？
如何从R中有分类处理的回归模型中得到均值比的标准误差？
提问于 2021-08-04 11:25:30
我有一个间作试验的一些例子数据，其中包括A和B两种作物的单作产量，并相互间作。在这种情况下，单作作物表现为与其自身的间作，即作物A与作物A=单作；作物B与作物B=单作；作物A与作物B=间作；作物B与作物B=间作。
我想要获得土地等值比率的平均值和SE，其定义为：
LER =(作物A与B间作)/(单作作物A的产量)+(作物A与B的产量)/(单作作物A的产量)
建议我对SE使用增量方法，因此我想使用来自msm的msm函数(但不太确定)
以下是数据(更多物种的大型实验的一个子集)：
structure(list(site = c("Two", "Two", "Two", "Two", "Two", "Two",
"Two", "Two", "Two", "Two", "Two", "Two", "One", "One", "One",
"One", "One", "One", "One", "One", "One", "One", "One", "One",
"Three", "Three", "Three", "Three", "Three", "Three", "Three",
"Three", "Three", "Three", "Three", "Three"), siteblock = c("Two_2",
"Two_1", "Two_3", "Two_1", "Two_3", "Two_2", "Two_1", "Two_3",
"Two_2", "Two_2", "Two_3", "Two_1", "One_2", "One_1", "One_3",
"One_3", "One_1", "One_2", "One_3", "One_1", "One_2", "One_2",
"One_3", "One_1", "Three_1", "Three_2", "Three_3", "Three_1",
"Three_2", "Three_3", "Three_1", "Three_2", "Three_3", "Three_1",
"Three_2", "Three_3"), species = c("A", "A", "A", "A", "A", "A",
"B", "B", "B", "B", "B", "B", "A", "A", "A", "A", "A", "A", "B",
"B", "B", "B", "B", "B", "A", "A", "A", "A", "A", "A", "B", "B",
"B", "B", "B", "B"), pair = c("A", "A", "A", "B", "B", "B", "A",
"A", "A", "B", "B", "B", "A", "A", "A", "B", "B", "B", "A", "A",
"A", "B", "B", "B", "A", "A", "A", "B", "B", "B", "A", "A", "A",
"B", "B", "B"), yield.pp = c(118.6648257, 80.29858998, 90.53153963,
48.64188012, 140.3934164, 103.0518444, 18.23102298, 7.104490919,
22.81524354, 21.07202039, 7.529486987, 18.63650567, 242.2567602,
202.1431331, 185.5609192, 283.144789, 241.1690115, 241.5258056,
23.78876862, 35.87524173, 41.1028137, 18.55380809, 22.46060419,
18.46323056, 242.9551749, 231.387521, 455.9878777, 288.2237713,
156.3390735, 207.4286019, 7.167311238, 34.66607289, 22.41394604,
42.22510313, 38.70415176, 57.86653817)), class = "data.frame", row.names = c(NA,
-36L))该实验在三个地点的三个街区内进行了复制，但为了简单起见，我们暂时忽略它。我创建了一个模型来估计每个组合中每种物种的产量：
mod<-lm(yield.pp~species*pair,data=df)
summary(mod)然后就可以很简单地得到每种作物在每种组合中的估计方法：
est_grid<-data.frame(expand.grid(species=levels(as.factor(df$species)),pair=levels(as.factor(df$pair))))
est_means<-predict(mod,newdata=est_grid)
est_means
## 1 2 3 4
##205.53182 23.68499 189.99091 27.27905 并计算平均土地当量比率：
(est_means[3]/est_means[1])+(est_means[2]/est_means[4])
## 3
##1.792635因此，我的问题是，如何获得符合平均比率的SE。我目前被deltamethod难住了，因为它似乎只提供估计模型参数的比率，而不是估计平均值的比率。
回答 2
Stack Overflow用户
回答已采纳
发布于 2021-08-04 14:54:54
使用emmeans使用ref_grid创建引用网格，在这种情况下，它将返回一个emmGrid对象，并且有一个vcov方法，这是从car包运行deltaMethod所需要的。
library(car)
library(emmeans)
r <- ref_grid(mod); r
## 'emmGrid' object with variables:
## species = A, B
## pair = A, B
s <- summary(r); s
## species pair prediction SE df
## A A 205.5 23.7 32
## B A 23.7 23.7 32
## A B 190.0 23.7 32
## B B 27.3 23.7 32
m <- with(s, setNames(prediction, paste0(species, pair)))
## AA BA AB BB
## 205.53182 23.68499 189.99091 27.27905
deltaMethod(m, "AB / AA + BA / BB", vcov(r))
## Estimate SE 2.5 % 97.5 %
## AB/AA + BA/BB 1.79264 1.15910 -0.47916 4.0644Stack Overflow用户
发布于 2021-08-08 03:09:48
如果您只需要比率，就意味着直接提供包。关键是将预测重新划分到日志规模：
> rlog <- ref_grid(mod, transform = "log")或者使用另一个答案的引用网格：rlog <- regrid(r, "log")。无论哪种方式，我们都转换成引用网格，就像应用了日志转换一样。使用type = "response"撤消日志转换：
> confint(rlog, type = "response")
species pair response SE df lower.CL upper.CL
A A 205.5 23.7 32 162.58 260
B A 23.7 23.7 32 3.10 181
A B 190.0 23.7 32 147.43 245
B B 27.3 23.7 32 4.66 160
Confidence level used: 0.95
Intervals are back-transformed from the log scale 这些正是@G.Grotandieck的答案中所得到的结果。现在进行两两比较；由于日志的差异是比率的日志，反向转换过程会产生比率：
> pairs(rlog, type = "response", infer = c(TRUE, TRUE), adjust = "none")
contrast ratio SE df lower.CL upper.CL null t.ratio p.value
A A / B A 8.678 8.725 32 1.1194 67.268 1 2.149 0.0393
A A / A B 1.082 0.183 32 0.7659 1.528 1 0.464 0.6460
A A / B B 7.534 6.591 32 1.2682 44.763 1 2.309 0.0276
B A / A B 0.125 0.125 32 0.0160 0.969 1 -2.069 0.0467
B A / B B 0.868 1.148 32 0.0587 12.846 1 -0.107 0.9156
A B / B B 6.965 6.102 32 1.1692 41.487 1 2.215 0.0340
Confidence level used: 0.95
Intervals are back-transformed from the log scale
Tests are performed on the log scale(我添加了infer参数以获得CIs和test。)
通过一些花哨的操作，甚至有可能回答原来的问题：
> con = contrast(rlog, list(`AB/AA` = c(-1,0,1,0), `BA/BB` = c(0,1,0,-1)))
> confint(con, type = "r")
contrast ratio SE df lower.CL upper.CL
AB/AA 0.924 0.157 32 0.6544 1.31
BA/BB 0.868 1.148 32 0.0587 12.85
Confidence level used: 0.95
Intervals are back-transformed from the log scale
> confint(contrast(regrid(con), list(`AB/AA + BA/BB` = c(1,1))))
contrast estimate SE df lower.CL upper.CL
AB/AA + BA/BB 1.79 1.16 32 -0.568 4.15
Confidence level used: 0.95置信限略有不同，因为我们使用32 d.f。而deltaMethod给出了渐近结果(无限d.f.)
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/68650393
复制相关文章
相似问题
腾讯云开发者
