有哪些替代方式可以在公式中指定二项式成功/试验次数？

Question

假设您正在对具有解释变量a和b的二项数据进行建模，其中每个响应是一定数量试验次数（N）中成功的数量（y）。有一些函数可以执行此类操作，它们似乎都使用不同的方法来指定y和N。

在glm中，您需要执行glm(cbind(y,N-y)~a+b, data = d)（左侧为成功/失败矩阵）。

在inla中，您需要执行inla(y~a+b, Ntrials=d$N, data=d)（单独指定试验次数）。

在glmmBUGS中，您需要执行glmmBUGS(y+N~a+b,data=d)（将成功和试验次数作为左侧的项进行指定）。

当编写新方法时，我一直认为最好遵循glm的做法，因为这是人们通常首次遇到二项响应数据的地方。然而，我总是记不住是否应该使用cbind(y,N-y)还是cbind(y,N)，而且我通常似乎具有数据中的成功/试验次数而非成功/失败次数-您的情况可能会有所不同。

当然，也有其他方法。例如，在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，使其变得明确。

有人有更好的想法吗？做这件事的正确方式是什么？

- Spacedman

---

我喜欢使用 %of%，因为它可以使顺序更加明显。但是它不太标准... - Ari B. Friedman

glm 还允许在左侧指定比例，使用额外的 weights 参数来提供分母。 - Ben Bolker

Answer 1

您还可以让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的代码片段，它给出了估计过程中正在进行的提示。

如果通过给出两列响应来指定binomialglm模型，则由prior.weights返回的权重是案例的总数（由提供的案例权重分解），结果的组件y是成功的比例。

另外，您可以使用model.frame创建glm.fit或glm的包装器。请参阅?model.frame中的...参数。

对于model.frame方法，进一步混合其他参数，如data、na.action、subset传递给默认方法。任何额外的参数（例如offset和weights或其他命名参数），将到达默认方法用于创建进一步的模型框架列，其名称带括号，例如"(offset)"。保留html标签。

评论

我看到Ben Bolker在后面的评论。上面是他指出的内容。

- Benjamin Christoffersen

Answer 2

我喜欢glm文档中的这种方法：

对于二项式和准二项式家族，响应也可以被指定为因子（当第一级表示失败时，所有其他级别表示成功）

这与我的经验中成功和失败的方式非常契合。其中一个是万能的（例如“没有投票”），而实现另一个的方法有多种（例如“投票给A”，“投票给B”）。希望从我表述的方式清楚地看出，“成功”和“失败”由glm定义可以任意定义，以便第一级对应于“失败”，而所有其他级别都是“成功”。

- rcorty

Answer 3

从glm的帮助页面上可以看到： "...或者作为一个两列矩阵，其中列分别给出了成功和失败的次数"

所以应该使用cbind(Y, N-Y)。

- Aviad Klein

---

这对原帖的问题有什么补充？我认为他已经说过了。 - Ben Bolker

我猜这与OP的以下评论有关：“我永远记不住是cbind(y，N-y)还是cbind(y，N)...”。尽管如此，这并不是对OP问题的回答。 - Benjamin Christoffersen