1

我正在尝试在 R 的逻辑回归模型中计算分类变量随时间推移的预测概率值和边际效应值（使用 p 值）。基本上，我想知道 1）响应变量的预测概率（事件发生）每年 2 个类别之一的样本站点和 2）一个站点在 1 个类别中与每年另一个类别的平均边际效应。我可以使用 ggeffects 包和 margins 包中的边际效应值来获得预测的概率值，但我还没有找到一种方法来从单个包中获取两组值。

所以我的问题是 1）是否有一个包/方法来获取这两组值，以及 2）如果我从 ggeffects 获得预测概率值和从边际获得边际效应值，这些值是否兼容？或者软件包处理模型的方式是否存在差异，这意味着我不能假设一个模型的边际效应对应于另一个模型的预测概率？3）在margins包中，如何获得两个因子变量随时间相互作用的平均边际效应？4）如何让margins()处理大型数据集？

以下是一些示例数据：

### Make dataset
df <- data.frame(year = rep(2001:2010, each = 100), 
                 state = rep(c("montana", "idaho", 
                               "colorado", "wyoming", "utah"),
                             times = 10, each = 20), 
                 site_id = as.factor(rep(1:100, times = 10)),
                 cat_variable = as.factor(rep(0:1, times = 5, each = 10)),
                 ind_cont_variable = rnorm(100, mean = 20, sd = 5),
                 event_occurred = as.factor(sample(c(0, 1), 
                                                    replace = TRUE, 
                                                    size = 1000)))

### Add dummy columns for states
library(fastDummies)
df <- dummy_cols(df, 
                 select_columns = "state",
                 remove_first_dummy = TRUE)

我感兴趣的是状态和分类变量对事件发生概率的影响，以及状态和分类变量的影响如何随时间变化。这是模型：

library(lme4)
fit_state <- glmer(event_occurred ~ ind_cont_variable +
                  cat_variable*year*state +
                  (1|site_id),
                data = df, 
                family = binomial(link = "logit"),
                nAGQ = 0,
                control = glmerControl(optimizer = "nloptwrap"))

我可以使用 ggeffects 来获取每个状态和类别组合随时间推移的预测概率值：

library(ggeffects)
fit_pp_state <- data.frame(ggpredict(fit_state, 
                          terms = c("year [all]",
                                    "cat_variable",
                                    "state")))

head(fit_pp_state)
### x = year, predicted = predicted probability, group = categorical variable level, facet = state
#    x predicted std.error  conf.low conf.high group    facet
# 2001 0.2835665 0.3981910 0.1535170 0.4634655     0 colorado
# 2001 0.5911911 0.3762090 0.4089121 0.7514289     0    idaho
# 2001 0.5038673 0.3719418 0.3288209 0.6779708     0  montana
# 4 2001 0.7101610 0.3964843 0.5297327 0.8420101     0     utah
# 5 2001 0.5714579 0.3747205 0.3901606 0.7354088     0  wyoming
# 6 2001 0.6788503 0.3892568 0.4963910 0.8192719     1 colorado

这对于可视化 5 个状态下预测概率随时间的变化非常有用。但我不知道如何使用 ggeffects 从这些值估计边际效应。使用 margins 包，我可以得到分类变量随时间的边际效应，但我不确定如何一起解释两个不同包的输出，或者这是否合适（我的前两个问题）。此外，我不知道如何获得边际以给我样本站点的边际效应，随着时间的推移处于分类变量级别/状态的每个组合中（带我到我的第三个问题）：

library(margins)
fit_state_me <- summary(margins(fit_state, 
                                at = list(year = 2001:2010),
                                variables = "cat_variable"))
head(fit_state_me)
#       factor      year     AME     SE       z      p   lower
# cat_variable1 2001.0000  0.0224 0.0567  0.3953 0.6926 -0.0887
# cat_variable1 2002.0000  0.0146 0.0490  0.2978 0.7659 -0.0814
# cat_variable1 2003.0000  0.0062 0.0418  0.1478 0.8825 -0.0757
# cat_variable1 2004.0000 -0.0026 0.0359 -0.0737 0.9413 -0.0731
# cat_variable1 2005.0000 -0.0117 0.0325 -0.3604 0.7186 -0.0754
# cat_variable1 2006.0000 -0.0208 0.0325 -0.6400 0.5222 -0.0845

我使用的实际数据集相当大（原始数据的 csv 为 1.51 GB，当我将其保存为 .rds 文件时，回归模型对象为 1.29 GB）。当我尝试对我的数据使用margins()时，我收到一条错误消息：

Error: cannot allocate vector of size 369.5 Gb

有关解决此问题的任何建议，以便我可以在我的数据上使用此功能？

我将不胜感激任何提示——我应该检查的包、我在代码中犯的错误或我的概念理解等。谢谢！

4

0 回答 0