我使用1.3.0包中的函数ggpredict()和函数来计算混合效应模型的平均估计值和置信区间（以下简称：CI）。这些功能依赖于ggemmeans()ggeffectspredict()emmeans()使其输出对 ggplot 友好。这两个函数预测/估计的值的平均值和 CI 均不同。为什么？

以下可重现的示例基于数据集 RIKZ（Janssen e Mulder 2005；Zuur 等人 2007），该数据集着眼于与平均潮汐水位（NAP，以米为单位）相比，物种丰富度（物种数量）如何随采样站高度变化）和暴露水平（具有三个级别的因素：低、中、高）：

rm(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() …

我正在尝试在数据之上手动绘制模型估计值。我真正的问题远比这复杂得多，所以predict如果可以的话，我想避免使用，并且更愿意了解如何计算这些预测而不是依赖某个包。

（底部可重现示例的数据。）

所以我首先运行一个模型，并获取模型估计和标准误差：

library(glmmTMB)
glmmLep<-glmmTMB(Lepidoptera ~ DayL50, 
                 data=Dat, family=nbinom2(link="log") )
dB_est<-(summary(glmmLep)$coeff$cond[2,1])
dB_SE<-(summary(glmmLep)$coeff$cond[2,2])
Int<-(summary(glmmLep)$coeff$cond[1,1])
Int_SE<-(summary(glmmLep)$coeff$cond[1,2])

然后，我创建了一系列 x 值来预测

x<-seq(from=min(Dat$DayL50),to=max(Dat$DayL50),length.out = length(Dat$DayL50))

然后我用两种不同的方法预测 y 值（使用predict和编写应该做同样事情的方程）

ypred<-exp(dB_est*x+Int)
y<-predict(glmmLep,list(DayL50=x),type="response",se.fit = T)

我们绘制两条预测线（一条是顶部较小的红线）：

ggplot(aes(x=DayL50,y=Lepidoptera),data=Dat)+
  geom_point(size=2)+
  geom_line(aes(y=y$fit,x=x),size=2)+
  geom_ribbon(aes(ymax=y$fit+1.96*y$se.fit,ymin=y$fit-1.96*y$se.fit,x=x),alpha=0.2)+
  geom_line(aes(y=ypred,x=x),size=1,color="red")+
  # geom_ribbon(aes(ymax=ymax,ymin=ymin,x=x),alpha=0.2,color="red")+
  coord_cartesian(ylim=c(0,1000))

我们看到我写的方程与predict函数的作用相同。都好。但是，当我在该行周围添加 SE / 95% CI 功能区时，我在尝试重新创建它时遇到了问题（这里我保留为 SE，因为 95%CI 会导致更多笨拙的情节）。我用许多不同的方式玩过这个公式，但似乎无法理解。出于某种原因，我似乎找不到任何关于它的帖子，但也许我没有使用正确的搜索词。任何人都可以向我解释我在这里缺少什么。似乎我在错误功能区（以红色标出）中遗漏了相当多的复杂性。

ymin<-exp((dB_est-dB_SE)*x+(Int))
ymax<-exp((dB_est+dB_SE)*x+(Int))


ggplot(aes(x=DayL50,y=Lepidoptera),data=Dat)+
  geom_point(size=2)+
  geom_line(aes(y=y$fit,x=x),size=2)+
  geom_ribbon(aes(ymax=y$fit+1.96*y$se.fit,ymin=y$fit-1.96*y$se.fit,x=x),alpha=0.2)+
  geom_line(aes(y=ypred,x=x),size=1,color="red")+
  geom_ribbon(aes(ymax=ymax,ymin=ymin,x=x),alpha=0.2,color="red")+
  coord_cartesian(ylim=c(0,1000))

或者使用 95% CI，就像我的predict色带一样，它甚至更远：

ymin<-exp((dB_est-1.96*dB_SE)*x+(Int))
ymax<-exp((dB_est+1.96*dB_SE)*x+(Int))

ggplot(aes(x=DayL50,y=Lepidoptera),data=Dat)+
  geom_point(size=2)+
  geom_line(aes(y=y$fit,x=x),size=2)+
  geom_ribbon(aes(ymax=y$fit+1.96*y$se.fit,ymin=y$fit-1.96*y$se.fit,x=x),alpha=0.2)+
  geom_line(aes(y=ypred,x=x),size=1,color="red")+
  geom_ribbon(aes(ymax=ymax,ymin=ymin,x=x),alpha=0.2,color="red")+
  coord_cartesian(ylim=c(0,1000))

Dat<-structure(list(Lepidoptera = c(0L, 0L, 1L, 0L, …

我正在运行 R 4.02 和 RStudio 1.3.1073。使用 glmmTMB 运行模型时出现错误。我应该更新还是恢复到其他版本？

.Call("FreeADFunObject", ptr, PACKAGE = DLL) 中的错误："FreeADFunObject" 不适用于包 "glmmTMB" 的 .Call()

我在安装时也收到此错误：

.Call("FreeADFunObject", ptr, PACKAGE = DLL) 中的错误：
“FreeADFunObject”不适用于包“glmmTMB”的 .Call() 尝试 URL 'https://cran.rstudio.com/bin/macosx/contrib/ 4.0/glmmTMB_1.0.2.1.tgz' 内容类型 'application/x-gzip' 长度 11403989 字节 (10.9 MB)

我有半连续数据（许多精确的零和连续的正结果），我正在尝试建模。我从 Zuur 和 Ieno 的《R 中零膨胀模型初学者指南》中学到了关于大量零质量的建模数据的知识，该指南区分了零膨胀伽玛模型和他们所描述的“零改变”伽玛模型作为障碍模型，结合了零点的二项式分量和正连续结果的伽玛分量。我一直在探索包ziGamma中选项的使用glmmTMB，并将所得系数与我按照 Zuur 书中的说明（第 128-129 页）构建的障碍模型进行比较，但它们并不相符。我无法理解为什么不这样做，因为我知道伽玛分布不能呈现零值，所以我认为每个零膨胀伽玛模型在技术上都是一个障碍模型。谁能为我阐明这一点？请参阅代码下方有关模型的更多注释。

library(tidyverse)
library(boot)
library(glmmTMB)
library(parameters)

### DATA

id <- rep(1:75000)
age <- sample(18:88, 75000, replace = TRUE)
gender <- sample(0:1, 75000, replace = TRUE)
cost <- c(rep(0, 30000), rgamma(n = 37500, shape = 5000, rate = 1), 
          sample(1:1000000, 7500, replace = TRUE))
disease <- sample(0:1, 75000, replace = TRUE)
time <- sample(30:3287, 75000, replace = TRUE)

df <- data.frame(cbind(id, disease, age, gender, cost, time))

# …

我想知道是否有一种简单的方法可以在不重新运行大型模型的情况下改变截距中的值，也许是数学上的。举个例子：

mtcars$cyl<-as.factor(mtcars$cyl)

summary(
  lm(mpg~cyl+hp,data=mtcars)
)

输出：

    Call:
lm(formula = mpg ~ cyl + hp, data = mtcars)

Residuals:
   Min     1Q Median     3Q    Max 
-4.818 -1.959  0.080  1.627  6.812 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 28.65012    1.58779  18.044  < 2e-16 ***
cyl6        -5.96766    1.63928  -3.640  0.00109 ** 
cyl8        -8.52085    2.32607  -3.663  0.00103 ** 
hp          -0.02404    0.01541  -1.560  0.12995    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.146 on 28 …

我正在使用 glmmTMB 库在 R 降价文档中运行混合模型。我运行的任何模型，都会收到以下警告：

'giveCsparse' has been deprecated; setting 'repr = "T"' for you'giveCsparse' has been deprecated; setting 'repr = "T"' for you'giveCsparse' has been deprecated; setting 'repr = "T"' for you

然后，如果我在控制台中运行代码，则会收到以下警告：

Warning messages:
1: In Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j = integer(0),  :
  'giveCsparse' has been deprecated; setting 'repr = "T"' for you
2: In Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j = integer(0),  :
  'giveCsparse' has been deprecated; setting 'repr = "T"' …

我有以下数据，并使用 R 中的包 glmmTMB 创建了一个模型，用于植物直径〜植物密度（植物数量），并具有随机绘图效果：

d <- data.frame (diameter = c(17,16,15,13,11, 19,17,15,11,11, 19,15,14,11,8),
                      plant_density = c(1000,2000,3000,4000,5000, 1000,2000,3000,4000,5000, 1000,2000,3000,4000,5000),
                      plot = c(1,1,1,1,1, 2,2,2,2,2, 3,3,3,3,3))

glmm.model <- glmmTMB(diameter ~  plant_density + (1|plot),
                        data = d,
                        na.action = na.omit,
                        family="gaussian",
                        ziformula = ~ 0)

我的目的是创建一个图，其中包含不同植物密度的预测直径数据，并包含随机图效果。所以我尝试预测数据：

new.dat <- data.frame(diameter= d$diameter,
                      plant_density = d$plant_density,
                      plot= d$plot) 

new.dat$prediction <- predict(glmm.model, new.data = new.dat, 
                              type = "response", re.form = NA)

不幸的是，我得到了每个地块的输出，但想要对直径〜植物密度进行广义预测。

我的目标是创建一个像这里一样的图，但使用 glmmTMB 的回归模型来考虑随机效应。

感谢你的帮助！

例子

library(glmmTMB)
library(ggeffects)

## Zero-inflated negative binomial model
(m <- glmmTMB(count ~ spp + mined + (1|site),
              ziformula=~spp + mined, 
              family=nbinom2, 
              data=Salamanders, 
              na.action = "na.fail"))
summary(m)

ggemmeans(m, terms="spp")
spp   | Predicted |       95% CI
--------------------------------
GP    |      1.11 | [0.66, 1.86]
PR    |      0.42 | [0.11, 1.59]
DM    |      1.32 | [0.81, 2.13]
EC-A  |      0.75 | [0.37, 1.53]
EC-L  |      1.81 | [1.09, 3.00]
DES-L |      2.00 | [1.25, 3.21]
DF    |      0.99 | [0.61, 1.62]

ggeffects::ggeffect(m, terms="spp")
spp …

我正在使用 glmmTMB 包运行混合模型，并使用预测函数使用以下代码计算预测平均值：

运行模型

model_1 <- glmmTMB(Step.rate ~ Treatment*Week + 
    (1|Treatment.Group/Lamb.ID) +  (1|Plot),
     data = data.df, family = nbinom1) 

创建新的数据框

new.dat <- data.frame(Treatment = data.df$Treatment,
                      Week = data.df$Week, Plot = data.df$Plot, 
                      Treatment.Group = data.df$Treatment.Group,
                      Lamb.ID = data.df$Lamb.ID) 

预测平均值

new.dat$prediction <- predict(model_1, new.data = new.dat, 
       type = "response", re.form = NA) 

这段代码工作正常，但是当我添加Interval =“confidence”来计算置信区间时，它似乎不起作用。R 忽略代码的最后部分，仅计算预测平均值。

new.dat$prediction <- predict(model_1, new.data = new.dat, 
     type = "response", re.form = NA, intervals = "confidence")

为什么间隔=“置信度”不起作用？这可能是与 glmmTMB 包相关的问题吗？

我将 pscl 中的 glmmTMB 和 Zeroinfl 应用到同一数据集。我获得了条件部分的相同系数，但二进制部分的系数有些不同。您对造成差异的潜在因素有什么想法吗？

谢谢！

这是 Zeroinfl 的结果：

# > summary(pscl.res)\n# \n# Call:\n#   zeroinfl(formula = outside_treatment ~ group + baseline.risk + Age.group + \n#              offset(log(Follow.up)) | group, data = final, dist = "negbin", \n#            link = "logit", trace = TRUE)\n# \n# Pearson residuals:\n#     Min      1Q  Median      3Q     Max \n# -1.0332 -0.7003 -0.4041  0.1675 12.5728 \n\n# Count model coefficients (negbin with log link):\n#                Estimate Std. Error z value Pr(>|z|)    \n# …

我正在运行一个零膨胀glmmTMB模型。我有兴趣在条件成分和零通胀成分的不同因子水平之间进行成对比较。有条件的部分，我可以用通常的emmeans方法轻松完成。我一直在尝试使用（相对）新创建的glmmTMB:::emm_basis.glmmTMB，但无法弄清楚该函数采用的一些参数，也找不到示例......

这是我目前所处位置的一个玩具示例。我专门poly()向模型添加了一个组件 - 我的完整模型同时具有poly()和ns()，因此需要弄清楚它们在这里是如何工作的。

所以问题如下：1）我trms提供的论点是否正确？2）函数需要什么参数xlev和grid参数emm_basis.glmmTMB？

library(glmmTMB)

data(Salamanders)
mod <- glmmTMB(count ~ spp + mined + poly(cover, 2) + (1|site), zi=~spp + mined, Salamanders, 
   family=nbinom2)

tt <- y ~ spp + mined + poly(cover, 2)
tt <- delete.response(terms(tt))

glmmTMB:::emm_basis.glmmTMB(mod,  trms = tt)

非常感谢您的任何想法！