11

假设您正在对二项式数据进行建模，其中每个响应是来自具有一些解释变量（a 和 b）的大量试验 (N) 的成功次数 (y)。有几个函数可以做这种事情，它们似乎都使用不同的方法来指定 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 一些含义，使其明确。

有没有人有更好的想法？这样做的正确方法是什么？

4

3 回答 3

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其中暗示了估计中的情况

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

或者，您可以为glm.fit或glm使用model.frame. 请参阅中的...论点?model.frame

...对于model.frame方法，混合进一步的参数，例如数据，，na.action传递subset给默认方法。到达默认方法的任何附加参数（例如offset和weights或其他命名参数）用于在模型框架中创建更多列，带有括号的名称，例如 "(offset)".

评论

之后我看到了 Ben Bolker 的评论。以上是他指出的。

于 2017-11-25T08:22:43.833 回答
0

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

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

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

于 2016-08-10T00:05:49.613 回答
-2

来自 r 在 glm 的帮助页面：“......或作为两列矩阵，其中的列给出成功和失败的数量”

所以它必须是 cbind(Y, NY)

于 2013-10-29T07:42:30.800 回答