小编Dyl*_*mes的帖子

最近，管道操作符的 RStudio 快捷方式%>%停止工作。这似乎是随机发生的，因为我没有更新任何内容。有谁知道发生了什么事以及我该如何解决这个问题？

> sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] dplyr_1.0.7

loaded via a namespace (and not attached):
 [1] rstudioapi_0.13  knitr_1.34       magrittr_2.0.1   tidyselect_1.1.1
 [5] R6_2.5.1         rlang_0.4.11     fastmap_1.1.0    fansi_0.5.0     
 [9] tools_4.1.0      xfun_0.25        utf8_1.2.2       cli_3.0.1       
[13] DBI_1.1.1        htmltools_0.5.2  ellipsis_0.3.2   yaml_2.2.1 …

我正在尝试在数据之上手动绘制模型估计值。我真正的问题远比这复杂得多，所以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, …

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

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 …

我正在尝试为stan_glmer.nb( rstanarm) 输出创建一个表，但是model_parameters从包中parameters抛出一个奇怪的错误，我不确定如何解决。也许这是一个错误。

sessionInfo()版本信息的缩短输出：

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

parameters_0.8.2
rstanarm_2.21.1

一个可重现的例子：

library(rstanarm)
library(parameters)

x<-rnorm(500)
dat<-data.frame(x=x,z=rep(c("A","B","C","D","E"),100), y=.2+x*.7)

mod1<-stan_glmer(y~x+(x|z),data=dat)

model_parameters(mod1, effects="all")

我将在此处为您节省输出，因为它并不重要，但该功能有效。现在负二项式模型：

dat.nb<-data.frame(x=rnorm(500),z=rep(c("A","B","C","D","E"),100),
                y=rnbinom(500,size=1,prob = .5))

mod2<-stan_glmer.nb(y~x+(x|z),data=dat.nb)

model_parameters(mod2, effects="all")

现在出现错误消息：

Error in `$<-.data.frame`(`*tmp*`, "parameter", value = c("(Intercept)",  : 
  replacement has 3 rows, data has 1

尽管在parameters0.10.1 版本中，@BenBolker 得到一个空白输出，而不是错误（见评论）。不管怎样，这个函数似乎不适用于rstanarm离散分布（见评论）。据我在帮助文档中看到，没有任何内容表明需要指定负二项式模型。此外，该功能适用​​于lme4模型：

library(lme4)
mod1<-lmer(y~x+(x|z),data=dat)

model_parameters(mod1, effects="all")

mod2<-glmer.nb(y~x+(x|z),data=dat.nb)

model_parameters(mod2, effects="all") …