Какими альтернативными способами можно указать биномиальные успехи/испытания в формуле?

Предположим, вы моделируете биномиальные данные, где каждый ответ представляет собой ряд успехов (y) из ряда испытаний (N) с некоторыми объясняющими переменными (a и b). Там есть несколько функций, которые делают такие вещи, и все они, кажется, используют разные методы для указания y и N.

В glm вы выполняете glm(cbind(y,N-y)~a+b, data = d) (матрица успеха/сбой на LHS)

В inla вы выполняете inla(y~a+b, Ntrials=d$N, data=d) (укажите количество испытаний отдельно)

В glmmBUGS вы выполняете glmmBUGS(y+N~a+b,data=d) (укажите успешные + испытания как термины на LHS)

При программировании новых методов я всегда считал, что лучше всего следить за тем, что делает glm, так как обычно люди сталкиваются с данными биномиального ответа. Тем не менее, я никогда не запомню, если его cbind(y,N-y) или cbind(y,N) - и у меня обычно, кажется, есть успех/количество проб в моих данных, а не успех/число сбоев - YMMV.

Конечно, возможны и другие подходы. Например, используя функцию в RHS, чтобы отметить, является ли переменная числом проб или числом сбоев:

myblm( y ~ a + b + Ntrials(N), data=d)
myblm( y ~ a + b + Nfails(M), data=d)  # if your dataset has succ/fail variables

или определяя оператор, чтобы просто сделать cbind, чтобы вы могли:

myblm( y %of% N ~ a + b, data=d)

таким образом, придавая значение LHS, что делает его явным.

Есть ли у кого-нибудь лучшие идеи? Какой правильный способ сделать это?

Ответ 1

Мне нравится этот метод из документации glm:

Для биномиальных и квазибиномиальных семейств ответ также может быть указанным как фактор (когда первый уровень обозначает отказ и все другие успехи)

Это хорошо сочетается с тем, что в моем опыте часто возникают успехи и неудачи. Один из них является "пойманным" (например, "не голосовал" ), и существует множество способов достижения другого (например, "проголосовали за A", "проголосовали за B" ). Я надеюсь, что это ясно из того, как я формулирую это, что "успех" и "неудача", как определено glm, могут быть определены произвольно, чтобы первый уровень соответствовал "сбою", а все остальные уровни - "успех" ".

Ответ 2

Вы также можете позволить y быть фракцией, в этом случае вам нужно поставить weights. Это не в аргументе formula, а почти равном количестве нажатий клавиш, как если бы он находился в formula. Вот пример

> set.seed(73574836)
> x <- runif(10)
> n <- sample.int(10, 2)
> y <- sapply(mapply(rbinom, size = 1, n, (1 + exp(1 - x))^-1), function(x) 
+   sum(x == 1))
> df <- data.frame(y = y, frac = y / n, x = x, weights = n)
> df
   y  frac      x weights
1  2 1.000 0.9051       2
2  5 0.625 0.3999       8
3  1 0.500 0.4649       2
4  4 0.500 0.5558       8
5  0 0.000 0.8932       2
6  3 0.375 0.1825       8
7  1 0.500 0.1879       2
8  4 0.500 0.5041       8
9  0 0.000 0.5070       2
10 3 0.375 0.3379       8
> 
> # the following two fits are identical
> summary(glm(cbind(y, weights - y) ~ x, binomial(), df))

Call:
glm(formula = cbind(y, weights - y) ~ x, family = binomial(), 
    data = df)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.731  -0.374   0.114   0.204   1.596  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.416      0.722   -0.58     0.56
x              0.588      1.522    0.39     0.70

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9.5135  on 9  degrees of freedom
Residual deviance: 9.3639  on 8  degrees of freedom
AIC: 28.93

Number of Fisher Scoring iterations: 3

> summary(glm(frac ~ x, binomial(), df, weights = weights))

Call:
glm(formula = frac ~ x, family = binomial(), data = df, weights = weights)

Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.731  -0.374   0.114   0.204   1.596  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.416      0.722   -0.58     0.56
x              0.588      1.522    0.39     0.70

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 9.5135  on 9  degrees of freedom
Residual deviance: 9.3639  on 8  degrees of freedom
AIC: 28.93

Number of Fisher Scoring iterations: 3

Причина того, что приведенные выше работы сводятся к тому, что glm действительно делает для биномиальных исходов. Он вычисляет фракцию для каждого наблюдения и вес, связанный с наблюдением, независимо от того, как вы определяете результат. Вот фрагмент от ?glm, который дает намек на то, что происходит в оценке

Если была задана модель binomial glm, давая двухколоночный ответ, веса, возвращаемые prior.weights, - это общее число случаев (с учетом прилагаемых весов корпуса) и компонента yрезультатом является доля успехов.

В качестве альтернативы вы можете сделать обертку для glm.fit или glm с помощью model.frame. См. Аргумент ... в ?model.frame

... для методов model.frame, сочетание дополнительных аргументов, таких как data, na.action, subset, чтобы перейти к методу по умолчанию. Любые дополнительные аргументы (например, offset и weights или другие имена аргументы), которые достигают метода по умолчанию, используются для создания столбцы в рамке модели с именами в скобках, такими как: "(offset)".

Комментарий

Я видел комментарий Бена Болкера после этого. Вышеупомянутое указывает на то, что он указывает.

Ответ 3

Из справочной страницы r на glm: "... или как двухколоночная матрица с столбцами, дающими число успехов и отказов"

Так что это должно быть cbind (Y, N-Y)