在 glm 中设置对比

Question

我有来自一组过度分散的条件的二项式计数数据。为了模拟它们，我使用了由 emdbook R 包的 rbetabinom 函数实现的 beta 二项分布：

library(emdbook)
set.seed(1)
df <- data.frame(p = rep(runif(3,0,1)),
                 n = as.integer(runif(30,100,200)),
                 theta = rep(runif(3,1,5)),
                 cond = rep(LETTERS[1:3],10),
                 stringsAsFactors=F)
df$k <- sapply(1:nrow(df), function(x) rbetabinom(n=1, prob=df$p[x], size=df$n[x],theta = df$theta[x], shape1=1, shape2=1))

我想找出每个条件 (cond) 对计数 (k) 的影响。我认为 MASS R 包的 glm.nb 模型允许建模：

library(MASS)
fit <- glm.nb(k ~ cond + offset(log(n)), data = df)

我的问题是如何设置对比，以便我得到每个条件的效果相对于所有条件的平均效果，而不是相对于 dummy 条件 A?

Answer 1

效果必须相对于某个基础水平进行估计。具有这 3 个条件中的任何一个的效果将与回归中的常数相同。

由于截距是当cond为两个估计水平（即"B"和"C"）=0时的预期平均值，因此它是仅供参考的平均值组（即 "A"）。

因此，您的模型中基本上已经有了这些信息，或者至少尽可能接近它。

对照组的平均值是截距加上对照组的系数。如您所知，比较组的系数因此会给您带来比较组 = 1 的效果（请记住，分类变量的每个级别都是一个虚拟变量，当该级别存在时 = 1）相对于参考团体。

因此，您的结果为您提供了每个级别的均值和相对效果。你当然可以根据自己的存在来切换参考电平。

这应该能为您提供所需的所有信息。如果不是，那么您需要问问自己，您到底在寻找什么信息。

Answer 2

两件事：（1）如果你想要相对于平均值的对比，使用contr.sum而不是默认的contr.treatment； (2) 您可能不应该用负二项式模型拟合 beta-二项式数据；改为使用 beta-binomial 模型（例如通过 VGAM 或 bbmle）！

library(emdbook)
set.seed(1)
df <- data.frame(p = rep(runif(3,0,1)),
             n = as.integer(runif(30,100,200)),
             theta = rep(runif(3,1,5)),
             cond = rep(LETTERS[1:3],10),
             stringsAsFactors=FALSE)
 ## slightly abbreviated
 df$k <- rbetabinom(n=nrow(df), prob=df$p,
                    size=df$n,theta = df$theta, shape1=1, shape2=1)

与VGAM:

 library(VGAM)
 ## note dbetabinom/rbetabinom from emdbook are masked
 options(contrasts=c("contr.sum","contr.poly"))
 vglm(cbind(k,n-k)~cond,data=df,
        family=betabinomialff(zero=2)
        ## hold shape parameter 2 constant
 )
 ## Coefficients:
 ## (Intercept):1 (Intercept):2         cond1         cond2 
 ##     0.4312181     0.5197579    -0.3121925     0.3011559 
 ## Log-likelihood: -147.7304

此处截距是跨层的平均形状参数； cond1 和 cond2 是级别 1 和级别 2 与平均值的差异（这不会给您级别 3 与平均值的差异，但根据构造它应该是 (-cond1-cond2 ) ...)

我发现使用 bbmle 的参数化（使用对数概率和离散参数）更容易一些：

 detach("package:VGAM")
 library(bbmle)
 mle2(k~dbetabinom(k, prob=plogis(lprob),
                   size=n,  theta=exp(ltheta)),
      parameters=list(lprob~cond),
      data=df,
      start=list(lprob=0,ltheta=0))
## Coefficients:
## lprob.(Intercept)       lprob.cond1       lprob.cond2            ltheta 
##       -0.09606536       -0.31615236        0.17353311        1.15201809 
## 
## Log-likelihood: -148.09

对数似然差不多（VGAM参数化好一点）；理论上，如果我们允许 shape1 和 shape2 (VGAM) 或 lprob 和 ltheta (bbmle) 在不同条件下发生变化，我们将获得两个参数化的相同对数似然。

在 glm 中设置对比

Set contrasts in glm

r

effects

glm

contrast