Какими альтернативными способами можно указать биномиальные успехи/испытания в формуле?
Предположим, вы моделируете биномиальные данные, где каждый ответ представляет собой ряд успехов (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)