代码之家 › 专栏 › 技术社区 › dan

在glm中设置对比度

contrast glm effects r

0

dan · 技术社区 · 10 年前

我有二项式计数数据,来自一组过度分散的条件。为了模拟它们,我使用由 rbetabinom 的功能 emdbook R 包裹:

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 ). 我认为 glm.nb 的模型 MASS R 包允许建模:

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

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

2 回复 | 直到 10 年前

1

Ben Bolker 10 年前

两件事:(1)如果你想对比平均值,使用 contr.sum 而不是默认值 contr.treatment ; (2) 你可能不应该用负二项模型拟合β二项数据;改用β二项模型(例如,via 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)

具有 :

 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公司 (使用logit概率和分散参数)稍微容易一点:

 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公司 )为了在不同的条件下变化,我们将得到两种参数化的相同的对数可能性。

2

1

Hack-R 10 年前

影响必须相对于某个基准面进行估计。具有这三个条件中的任何一个条件的效果都将与回归中的常数相同。

由于截距是当 cond 对于两个估计水平(即。 "B" 和 "C" ),仅为参考组的平均值(即。 "A" ).

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

比较组的平均值是截距加上比较组的系数。因此,如您所知,比较组的系数为您提供了相对于参考组的比较组=1的效果(请记住,分类变量的每个级别都是一个虚拟变量,当该级别存在时,该虚拟变量=1)。

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

希望这能为您提供所需的所有信息。如果不是,那么你需要准确地问自己你想要的是什么信息。