警告:二项式glm中的非整数#successes!(调查包)

Question

我正在使用该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上是否会更好,但总的来说我想我会先在这里尝试一下.

Answer 1

没有什么不对,glm在指定二项式(和泊松)模型时只是挑剔.如果检测到没有,它会发出警告.试验或成功是非完整的,但它仍然适用于模型.如果您想要取消警告(并且您确定它不是问题),请family=quasibinomial改用.

Answer 2

计算上没有任何错误,但从统计学角度来说,你可能没有做出有意义的事情.在这种情况下,最好使用稳健的回归方法,如果您的数据包含精确为1或精确为0的单位,这通常是比例响应数据的好主意.

Answer 3

就像@HoongOoi说，glm.fit与binomial家人预计整数计数，否则，会抛出一个警告; 如果您需要非整数计数，请使用quasi-binomial。我的其余答案将这些进行比较。

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

源代码比较

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

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

以及以下移除：

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

只能simfun有所作为，但是的源代码未glm.fit显示该功能的使用，与stats::binomial诸如mu.eta和所返回的对象中的其他字段不同link。

最小的工作例子

在这个最小的工作示例中，使用quasibinomial或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 \ k），尽管可以使用Gamma函数。这对于概率分布的定义可能是个问题，但与估计无关，因为项（n选择k）不取决于参数且不属于最佳化。

对数似然估计器是：

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

其中常数与参数theta无关，因此不符合优化要求。

标准误差比较

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

GLM的色散未在拟合过程中使用，但需要使用它来查找标准误差。如果dispersion未提供或NULL，则将色散1与binomial和Poisson族一样，并用剩余卡方统计量（根据权重非零的情况计算）除以剩余自由度来估算。

...

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