56

我正在使用该twang包来创建倾向得分，这些得分在二项式 glm 中用作权重，使用survey::svyglm. 代码看起来像这样：

pscore <- ps(ppci ~ var1+var2+.........., data=dt....)

dt$w <- get.weights(pscore, stop.method="es.mean")

design.ps <- svydesign(ids=~1, weights=~w, data=dt,)

glm1 <- svyglm(m30 ~ ppci, design=design.ps,family=binomial)

这会产生以下警告：

Warning message:
   In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!

有谁知道我可能做错了什么？

我不确定这条消息在 stats.SE 上是否会更好，但总的来说，我想我会先在这里尝试。

4

4 回答 4

77

没有错，glm只是在指定二项式（和泊松）模型时很挑剔。如果检测到否，它会发出警告。试验或成功的次数是不完整的，但它会继续前进并且无论如何都适合模型。如果您想抑制警告（并且您确定这不是问题），请family=quasibinomial改用。

于 2012-10-18T12:00:53.100 回答
18

就像@HoongOoi 所说，glm.fit家庭binomial期望整数计数，否则会发出警告；如果您想要非整数计数，请使用quasi-binomial. 我的其余答案比较了这些。

R for 中的准二项式与系数估计glm.fit完全相同binomial（如@HongOoi 的评论中所述），但与标准误差不同（如@nograpes 的评论中所述）。

源码对比

源代码的差异stats::binomial并stats::quasibinomial显示以下更改：

  • 文本“二项式”变为“准二项式”
  • aic 函数返回 NA 而不是计算的 AIC

以及以下删除：

  • 当权重 = 0 时将结果设置为 0
  • 检查权重的完整性
  • simfun模拟数据的功能

onlysimfun可以有所作为，但是 的源代码glm.fit显示没有使用该函数，这与返回的对象中的其他字段不同，stats::binomial例如mu.etaand link。

最小的工作示例

对于这个最小工作示例中的系数，使用quasibinomialor的结果是相同的：binomial

library('MASS')
library('stats')

gen_data <- function(n=100, p=3) {

  set.seed(1)  
  weights <- stats::rgamma(n=n, shape=rep(1, n), rate=rep(1, n))
  y <- stats::rbinom(n=n, size=1, prob=0.5)
  theta <- stats::rnorm(n=p, mean=0, sd=1)
  means <- colMeans(as.matrix(y) %*% theta)
  x <- MASS::mvrnorm(n=n, means, diag(1, p, p))

  return(list(x=x, y=y, weights=weights, theta=theta))  
}

fit_glm <- function(family) {
  data <- gen_data()
  fit <- stats::glm.fit(x = data$x,
                        y = data$y,
                        weights = data$weights,
                        family = family)
  return(fit)
}

fit1 <- fit_glm(family=stats::binomial(link = "logit"))
fit2 <- fit_glm(family=stats::quasibinomial(link = "logit"))

all(fit1$coefficients == fit2$coefficients)

与拟二项式概率分布的比较

该线程表明准二项式分布与具有附加参数的二项式分布不同phi。但它们在统计和R.

首先，源代码中没有quasibinomial提到该附加phi参数的地方。

其次，准概率类似于概率，但不是正确的概率。在这种情况下，当数字是非整数时，无法计算项 (n \choose k)，尽管可以使用 Gamma 函数。这对于概率分布的定义可能是一个问题，但与估计无关，因为项 (n 选择 k) 不依赖于参数并且不属于优化。

对数似然估计是：

对数似然估计器

作为二项式族的 theta 函数的对数似然是：

二项式家庭的对数似然

其中常数与参数 theta 无关，因此它不在优化范围内。

标准误比较

如stats::summary.glm中所述，通过对和系列stats::summary.glm使用不同的离散值计算的标准误差：binomialquasibinomial

拟合过程中不使用 GLM 的离散度，但需要找到标准误差。如果dispersion未提供 或NULL，则将离散视为和族，否则通过残差卡方统计量（根据非零权重的情况计算）除以残差自由度来估计1。binomialPoisson

...

cov.unscaled: 估计系数的未缩放 ( dispersion = 1) 估计协方差矩阵。

cov.scaled: 同上，按比例缩放dispersion。

使用上面的最小工作示例：

summary1 <- stats::summary.glm(fit1)
summary2 <- stats::summary.glm(fit2)

print("Equality of unscaled variance-covariance-matrix:")
all(summary1$cov.unscaled == summary2$cov.unscaled)

print("Equality of variance-covariance matrix scaled by `dispersion`:")
all(summary1$cov.scaled == summary2$cov.scaled)

print(summary1$coefficients)
print(summary2$coefficients)

显示相同的系数、相同的未缩放方差-协方差矩阵和不同的缩放方差-协方差矩阵：

[1] "Equality of unscaled variance-covariance-matrix:"
[1] TRUE
[1] "Equality of variance-covariance matrix scaled by `dispersion`:"
[1] FALSE
       Estimate Std. Error   z value   Pr(>|z|)
[1,] -0.3726848  0.1959110 -1.902317 0.05712978
[2,]  0.5887384  0.2721666  2.163155 0.03052930
[3,]  0.3161643  0.2352180  1.344133 0.17890528
       Estimate Std. Error   t value   Pr(>|t|)
[1,] -0.3726848  0.1886017 -1.976042 0.05099072
[2,]  0.5887384  0.2620122  2.246988 0.02690735
[3,]  0.3161643  0.2264421  1.396226 0.16583365
于 2019-03-21T16:03:29.230 回答
8

计算上没有错，但从统计上看，你可能没有做一些有意义的事情。在这种情况下，最好使用稳健的回归方法，如果您的数据包含正好为 1 或正好为 0 的单位，这对于比例响应数据通常是一个好主意。

于 2015-07-23T18:42:23.070 回答
0

抱歉，但从某种意义上说，如果底层机制是过度分散的二项式模型，则在估计标准误差时，过度分散的二项式将考虑到它。因此，即使点估计值相同，您也会获得更好的覆盖率。

于 2017-06-01T19:09:32.257 回答