Установить контрасты в glm

Question

У меня есть данные биномиального счета, исходящие из набора условий, которые чрезмерно перегружены. Для того, чтобы имитировать их я использую бета-биномиальное распределение, реализованный rbetabinom функции пакета emdbookR:

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 модель пакета MASSR позволяет моделировать, что:

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

Мой вопрос заключается в том, чтобы установить контрасты таким образом, что я получаю эффект каждого условия относительно средних эффектов по всем условиям, а не относительный к состоянию dummyA?

Answer 1

Две вещи: (1) если вы хотите контрастов относительно среднего значения, используйте contr.sum, а не по умолчанию contr.treatment; (2) вы, вероятно, не должны соответствовать бета-биномиальным данным с отрицательной биномиальной моделью; вместо этого используйте бета-биномиальную модель (например, через 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 парметрирование немного лучше); теоретически, если бы мы меняли разные условия, как форма1, так и форма 2 (VGAM) или lprob и ltheta (bbmle), мы получали бы те же логарифмические вероятности для обеих параметризаций.

Answer 2

Эффекты необходимо оценить относительно базового уровня. Эффект от любого из трех условий будет таким же, как константа в регрессии.

Поскольку перехват является ожидаемым среднее значение, когда cond является = 0 для обоих расчетных уровней (т.е. "B" и "C"), то есть среднее значение только для контрольной группы (т.е. "A").

Таким образом, вы в основном уже имеете эту информацию в своей модели или, по крайней мере, как можно ближе к ней.

Среднее значение группы сравнения - это перехват плюс коэффициент группы сравнения. Таким образом, коэффициенты групп сравнения, как вы знаете, дают вам эффект наличия группы сравнения = 1 (учитывая, что каждый уровень вашей категориальной переменной является фиктивной переменной, которая = 1, когда этот уровень присутствует) относительно ссылки группа.

Итак, ваши результаты дают вам средства и относительные эффекты каждого уровня. Вы можете, конечно, отключить контрольный уровень в соответствии с вашим присутствием.

Это, надеюсь, даст вам всю необходимую информацию. Если нет, тогда вам нужно спросить себя, какая именно информация вам нужна.