这是运行单个固定效应方法的三种不同方法，它们给出或多或少相同的结果（见下文）。我的主要问题是如何使用第二个模型 ( model_plm) 或第三个模型 ( model_felm) 获得预测概率或平均边际效应。我知道如何使用第一个模型 ( model_lm) 来做到这一点，并使用下面的示例来展示ggeffects，但这仅在我有一个小样本时才有效。
由于我有超过一百万人，我的模型只能使用model_plm和来工作model_felm。如果我使用model_lm，则需要花费大量时间来运行一百万个人，因为它们是在模型中受到控制的。我还收到以下错误：Error: vector memory exhausted (limit reached?)。我检查了 StackOverflow 上的许多线程来解决该错误，但似乎没有任何解决方案。
我想知道是否有有效的方法来解决这个问题。我的主要兴趣是提取交互的预测概率residence*union。我通常使用以下软件包之一提取预测概率或平均边际效应：ggeffects、emmeans或margins。
library(lfe)
library(plm)
library(ggeffects)
data("Males")
model_lm = lm(wage ~ exper + residence+health + residence*union +factor(nr)-1, data=Males)
model_plm = plm(wage ~ exper + residence + health + residence*union,model = "within", index=c("nr", "year"), data=Males)
model_felm = felm(wage ~ exper + residence + health …我使用1.3.0包中的函数ggpredict()和函数来计算混合效应模型的平均估计值和置信区间（以下简称：CI）。这些功能依赖于ggemmeans()ggeffectspredict()emmeans()使其输出对 ggplot 友好。这两个函数预测/估计的值的平均值和 CI 均不同。为什么？
以下可重现的示例基于数据集 RIKZ（Janssen e Mulder 2005；Zuur 等人 2007），该数据集着眼于与平均潮汐水位（NAP，以米为单位）相比，物种丰富度（物种数量）如何随采样站高度变化）和暴露水平（具有三个级别的因素：低、中、高）：
\nrm(list=ls())\nif (!require(pacman)) install.packages(\'pacman\'); library(pacman)\np_load(emmeans)\np_load(ggplot2)\np_load(ggpubr)\np_load(ggeffects)\np_load(lme4, lmerTest, glmmTMB)\np_load(RCurl)\n# get data:\nRIKZ <- read.csv(text = RCurl::getURL(\n"https://raw.githubusercontent.com/marcoplebani85/datasets/master/RIKZ.csv"))\nstr(RIKZ)\n# "Exposure" is a factor:\nRIKZ$Exposure <- as.factor(RIKZ$Exposure)\n在这里，我使用泊松分布残差将广义混合效应模型拟合到数据中glmmTMB()：
mem1 <- glmmTMB(Richness ~ NAP+Exposure + (1 | Beach),\n family="poisson",\n data = RIKZ, REML=T)\n模型预测和 CI 根据ggeffects::ggpredict()，不考虑随机效应的不确定性（参见本页了解为何考虑或不考虑）：
richness.predicted <- ggpredict(mem1, \nterms=c("NAP", "Exposure"), type="fixed")\n同一模型的预测和 CI根据ggeffects::ggemmeans() …
我正在使用该emmeans软件包，但无法理解为什么它会起作用
library(emmeans)\nmod <- lm(mpg ~ disp + hp, data = mtcars)\nl1 <- emmeans(mod, list("disp","hp"))\n\nfor (x in l1) {print(data.frame(x))}\n\n disp emmean SE df lower.CL upper.CL\n1 230.7219 20.09062 0.5527102 29 18.96021 21.22104\n hp emmean SE df lower.CL upper.CL\n1 146.6875 20.09062 0.5527102 29 18.96021 21.22104\n但这不是吗？
\nlapply(l1, function(x) data.frame(x))\n\n Error in as.data.frame.default(x[[i]], optional = TRUE) : \ncannot coerce class \xe2\x80\x98"function"\xe2\x80\x99 to a data.frame\n使用普通列表，我得到预期的输出（列表中的数据帧）：
\nl2 <- list(a=matrix(c(1,2,3,4)), b=matrix(c(5,6,7,8)))\nlapply(l2, function(x) data.frame(x))\n$a\n x\n1 1\n2 2\n3 1\n4 3\n\n$b\n x\n1 …我正在使用我的响应变量作为百分比（0-100）进行 GAM。我使用反正弦变换来改进模型拟合（asin(sqrt(myvariable/100))）。我现在想评估原始规模上解释因子变量水平之间的对比。我一直在尝试使用 emmeans 并按照转换和链接函数小插图中的步骤来设置我的模型，并以 emmeans 可以读取的格式进行转换。但是，当我运行 emmeans 函数时，出现以下错误： link$mu.eta(object@bhat[estble]) 中的错误：尝试应用非函数。
我这样设置转换对象：
tran <- make.tran("asin.sqrt", 100)
我相信 bit 可以工作，因为当我用 emmeans 在线性模型上尝试它时，它工作了：
warp.t <- with(tran, lm(linkfun(breaks)~wool*tension, warpbreaks))
emmeans(warp.t, ~wool|tension, type="response")
tension = L:
wool response SE df lower.CL upper.CL
A 44.2 4.00 48 36.3 52.3
B 27.7 3.61 48 20.8 35.3
tension = M:
wool response SE df lower.CL upper.CL
A 23.5 3.41 48 17.0 30.7
B 28.4 3.63 48 21.4 35.9
tension = H:
wool response SE df lower.CL …我不知道如何提取emmean 和SE从列emmGrid 的emmeans R包。下面给出了 MWE。
library(emmeans)
warp.lm <- lm(breaks ~ wool * tension, data = warpbreaks)
Test <- emmeans(warp.lm, specs = "wool")
Test
wool emmean SE df lower.CL upper.CL
A 31.03704 2.105459 48 26.80373 35.27035
B 25.25926 2.105459 48 21.02595 29.49257
Results are averaged over the levels of: tension
Confidence level used: 0.95
class(Test)
[1] "emmGrid"
attr(,"package")
[1] "emmeans"
我对计算混合模型的 SE 很感兴趣。为此，首先我在一个更简单的模型中使用了包中包含的数据集之一。
\n\npigs$percent <- as.factor(pigs$percent)\nDoc_lm_1 <- lm(conc~percent, pigs) \nsummary(Doc_lm_1)\nemmeans(Doc_lm_1, pairwise~percent)$emmeans\n输出：
\n\npercent emmean SE df lower.CL upper.CL\n9 32.7 2.92 25 26.7 38.7\n12 38.0 2.76 25 32.3 43.7\n15 40.1 3.12 25 33.7 46.6\n18 39.9 3.70 25 32.3 47.6\n当我尝试使用平衡数据集时，所有组的 SE 都是相同的，并且与手工制作的 SE 不匹配。我想在这种情况下不会考虑任何因素，但它仍然应该与手工制作的SE相匹配
\n\n难道SE就是参数的SE吗？正如我们在表中看到的，当数据不平衡时，组间的 SE 会有所不同。我的假设基于以下事实：该包的 cran 项目网站表明（https://cran.r-project.org/web/packages/emmeans/vignettes/basics.html#backstory）：
\n\n估计边际均值基于模型 \xe2\x80\x93，而不是直接基于数据”
\n\n所以我问我，SE是如何计算的？添加随机因子将如何改变这个计算？提前致谢。
\n我可以通过下面的代码得到两两比较的意义
m <- lmer(angle ~ recipe*temp + (1|replicate), data=cake)
emtrends(m, pairwise~recipe, var="temp")
$emtrends
recipe temp.trend SE df lower.CL upper.CL
A 0.1537143 0.02981898 250 0.09498586 0.2124427
B 0.1645714 0.02981898 250 0.10584300 0.2232999
C 0.1558095 0.02981898 250 0.09708110 0.2145379
$contrasts
contrast estimate SE df t.ratio p.value
A - B -0.010857143 0.0421704 250 -0.257 0.9641
A - C -0.002095238 0.0421704 250 -0.050 0.9986
B - C 0.008761905 0.0421704 250 0.208 0.9765
但是，如果我对每个趋势recipe本身是否重要感兴趣怎么办？我怎样才能得到 的意义$emtrends？
我正在分析一个数据集，其中 y 值为蛋白质水平范围为 0-3.6nmol/L。我对线性模型进行了 sqrt 变换，并尝试进行事后测试，但我的 emmean 值和置信区间没有意义，特别是因为其中一个平均值以某种方式结果为负。
我这样初始化我的模型：
fasting.sqrt <- lm(sqrt(FCP) ~ Type * BMI_Percentile, data = Database_Mexican_children_with_different_types_of_DM)
然后制作我的数据框以供事后专门查看类型：
fasting.sqrt.emm = emmeans(fasting.sqrt, specs="Type") %>% as.data.frame()
返回这个：
# Type emmean SE df lower.CL upper.CL
# 1 0.732 0.0980 129 0.538 0.926
# 2 -0.391 0.6511 129 -1.679 0.897
# 3 0.708 0.0728 129 0.564 0.852
# 4 0.260 0.8955 129 -1.512 2.032
并且，使用此图形代码：
ggplot() +
geom_point(data=Database_Mexican_children_with_different_types_of_DM, aes(x=Type, y=FCP), alpha=0.3) +
geom_errorbar(data=fasting.sqrt.emm, aes(x=Type, ymin=lower.CL*abs(lower.CL), ymax=(upper.CL)^2), width=0.1) +
geom_point(data=fasting.sqrt.emm, aes(x=Type, y=emmean*abs(emmean)),
color="red", …其他关键字： 最佳线性无偏估计量 (BLUE)、调整均值、混合模型、固定效应、线性组合、对比度、R
mmer()使用sommer包拟合模型后- 是否可以从对象获得估计边际均值 ( emmeans()) /最小二乘均值 (LS-means)mmer？也许类似于predict()ASReml-R v3 的功能？
实际上，我想要很多东西，也许分开要求会更清楚：
predict(..., vcov=T)）emmeans(mod, pairwise ~ effect, adjust="Tukey")predict(..., sed=T)）所以是的，基本上 和 的混合predict()将emmeans()是这里的目标。
提前致谢！
我正在使用该包检查两个连续预测变量之间的相互作用emmeans。我使用包lm_robust()中的estimatr内容来执行线性回归并获得集群稳健的标准误差。结果变量以 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包创建者网站上的一个简短示例表明lm_robust对象可以与 一起使用，但我在“emmeans 支持的模型”小插图页面或包文档中看emmeans不到将其列为受支持的模型。lm_robust