我正在使用 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上发布，但总的来说我想先在这里尝试一下。
没有问题，glm只是在指定二项式（和泊松）模型时十分挑剔。如果它检测到试验数或成功次数不是整数，它会发出警告，但还是会拟合模型。如果你想要抑制警告（并且你确定这不是问题），可以使用family=quasibinomial。
model$aic 会返回 NA。 - joel.wilsonlibrary('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无关，因此它从优化中排除。
stats::summary.glm计算的标准误差对于binomial和quasibinomial系列使用不同的离散度值，如stats::summary.glm所述：
GLM的离散度在拟合过程中没有使用，但需要用于找到标准误差。如果未提供
dispersion或为NULL，则对于binomial和Poisson系列，离散度被视为1，否则通过剩余自由度除以剩余卡方统计量（根据具有非零权重的案例计算）来估计离散度。...
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
stats::summary.glm函数，就可以看到这种差异。对于二项式家族（和泊松分布），离散度被硬编码为1。对于准二项式家族，离散度以“通常”的方式计算。这导致了不同的缩放协方差矩阵，从而导致了不同的标准误差。 - nograpes计算上没有任何问题，但从统计学的角度来看，你可能在做一些不太合理的事情。在这种情况下，最好使用鲁棒回归方法，如果你的数据包括恰好为1或恰好为0的单位，则通常是比例响应数据的一个好选择。
抱歉，但从更健壮的角度来看，如果底层机制是过分离散的二项模型，那么在估计标准误差时，过分离散的二项式将会加以考虑。因此，即使点估计相同，您也将获得更好的覆盖率。
glm.fit函数中，我看不出标准误差可能会有所不同（请参见我的答案）。也许survey在二项式家族中使用了simfun的附加功能。 - miguelmorin
m30变量的类型是什么？ - James