在 glm 中设置对比度答案

【问题标题】：Set contrasts in glm在 glm 中设置对比度
【发布时间】：2016-08-06 17:05:01
【问题描述】：

我有来自一组过度分散的条件的二项式计数数据。为了模拟它们，我使用了emdbookR 包的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？

【问题讨论】：

标签： r effects glm contrast

【解决方案1】：

两件事：（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

这里 intercept 是各个级别的平均形状参数； cond1 和 cond2 是 1 级和 2 级与平均值的差异（这并没有给你 3 级与平均值的差异，但通过构造它应该是 (-cond1-cond2) ...)

我发现使用bbmle（带有logit-probability 和色散参数）的参数化更容易一些：

 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) 在不同条件下变化，我们将得到两个参数化的相同对数似然。

【讨论】：

有没有办法获取std. cond3的错误？或者它只是标准的总和。 cond1 和 cond2 的错误？
不完全是，你也需要协方差。应该是sqrt(var1+var2-2*cov(1,2))；不幸的是，大多数便利包和方法（预测 std 错误、lsmeans、效果包）在这种情况下都不起作用......
非常感谢。一个相关的问题。在我的实验中，n 是细胞计数，k 是细胞计数，其中我感染的病毒整合到基因组中（因此给出了可检测的信号），并且条件是不同的细胞类型。假设我的部分分散是由于感染的技术变异性（而不是生物变异性），对此我有每种细胞类型的经验控制数据（同样，k's 和 n's）。如何从我拟合真实数据的 beta 二项式模型的离散度中减去技术可变性？
没有什么明显的想法。尝试在stats.stackexchange.com 上提问...？

【解决方案2】：

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

因为当cond 对两个估计水平（即"B" 和"C"）都为0 时，截距是预期平均值，所以它只是参考组的平均值（即"A"）。

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

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

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

这有望为您提供所需的所有信息。如果不是，那么您需要准确地问自己，您所追求的信息是什么。

【讨论】：