这里有一个关于如何在 Stata 中为比例数据运行 GLM 的示例
IV 是在学校接受免费或减价膳食的学生比例。stata 模型如下所示：
glm meals yr_rnd parented api99, link(logit) family(binomial) robust nolog
我有兴趣学习如何在 R 中复制这个结果（最好使用相同的稳健方法）。让我们想象一下，我有关于获得免费膳食的学生人数（成功）和其余学生（失败）的数据。我猜 R 中的模型可能看起来像这样：
fitglm <- glm(cbind(Successes,Failures) ~ yr_rnd + parented + api99, family=binomial)
此外，有人在其他地方（Penguin_Knight）向我指出错误消息"meals has non-integer values"可能是错误的。我对这个错误一无所知...
使用 R 包sandwich，您可以复制这样的结果（我假设您已经下载了数据集）：
#-----------------------------------------------------------------------------
# Load the required packages
#-----------------------------------------------------------------------------
require(foreign)
require(sandwich)
#-----------------------------------------------------------------------------
# Load the data
#-----------------------------------------------------------------------------
dat <- read.dta("MyPath/proportion.dta")
#-----------------------------------------------------------------------------
# Inspect dataset
#-----------------------------------------------------------------------------
str(dat)
#-----------------------------------------------------------------------------
# Fit the glm
#-----------------------------------------------------------------------------
fitglm <- glm(meals ~ yr_rnd + parented + api99, family = binomial(logit), data = dat)
#-----------------------------------------------------------------------------
# Output of the model
#-----------------------------------------------------------------------------
summary(fitglm)
#-----------------------------------------------------------------------------
# Calculate robust standard errors
#-----------------------------------------------------------------------------
cov.m1 <- vcovHC(fitglm, type = "HC0")
std.err <- sqrt(diag(cov.m1))
q.val <- qnorm(0.975)
r.est <- cbind(
Estimate = coef(fitglm)
, "Robust SE" = std.err
, z = (coef(fitglm)/std.err)
, "Pr(>|z|) "= 2 * pnorm(abs(coef(fitglm)/std.err), lower.tail = FALSE)
, LL = coef(fitglm) - q.val * std.err
, UL = coef(fitglm) + q.val * std.err
)
r.est
使用稳健标准误差的模型输出为：
Estimate Robust SE z Pr(>|z|) LL UL
(Intercept) 6.801682703 0.072368970 93.98618 0.000000e+00 6.659842129 6.943523277
yr_rndYes 0.048252657 0.032167588 1.50004 1.336041e-01 -0.014794657 0.111299970
parented -0.766259824 0.039066917 -19.61403 1.173462e-85 -0.842829574 -0.689690073
api99 -0.007304603 0.000215534 -33.89072 9.127821e-252 -0.007727042 -0.006882164
估计和标准误差与使用 Stata 计算的非常相似。我不知道为什么拦截是不同的。Stata输出是：
------------------------------------------------------------------------------
| Robust
meals | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
yr_rnd | .0482527 .0321714 1.50 0.134 -.0148021 .1113074
parented | -.7662598 .0390715 -19.61 0.000 -.8428386 -.6896811
api99 | -.0073046 .0002156 -33.89 0.000 -.0077271 -.0068821
_cons | 6.75343 .0896767 75.31 0.000 6.577667 6.929193
------------------------------------------------------------------------------
该函数有多种方法可用vcovHC。有关详细信息，请参阅的帮助文件vcovHC。
请注意，如果您使用 option family = quasibinomial(logit)，则不会出现错误消息（请参阅此处）。
您可以复制加州大学洛杉矶分校关于比例的常见问题解答（以百分比作为因变量），如下所示：
require(foreign);require(lmtest);require(sandwich)
meals <- read.dta("http://www.ats.ucla.edu/stat/stata/faq/proportion.dta")
fitperc <- glm(meals ~ yr_rnd + parented + api99, family = binomial, data=meals)
## Warning message:
## In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
我不知道上面的警告是否是这里的问题。由于某种原因，截距在 R 和 Stata 中不匹配，但由于我们通常不会在 logit/probit 中解释它，所以它应该无关紧要。
summary(fitperc)
##
## Call:
## glm(formula = meals ~ yr_rnd + parented + api99, family = binomial,
## data = meals, na.action = na.exclude)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.77722 -0.18995 -0.01649 0.18692 1.60959
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.801683 0.231914 29.329 <2e-16 ***
## yr_rndYes 0.048253 0.104210 0.463 0.643
## parented -0.766260 0.090733 -8.445 <2e-16 ***
## api99 -0.007305 0.000506 -14.435 <2e-16 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1953.94 on 4256 degrees of freedom
## Residual deviance: 395.81 on 4253 degrees of freedom
## (164 observations deleted due to missingness)
## AIC: 2936.7
##
## Number of Fisher Scoring iterations: 5
在 R 中，使用的小样本校正与 Stata 中的不同，但稳健的 SE 非常相似：
coeftest(fitperc, function(x) vcovHC(x, type = "HC1"))
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.80168270 0.07240299 93.9420 <2e-16 ***
## yr_rndYes 0.04825266 0.03218271 1.4993 0.1338
## parented -0.76625982 0.03908528 -19.6048 <2e-16 ***
## api99 -0.00730460 0.00021564 -33.8748 <2e-16 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
要使用完全相同的小样本校正，您需要关注这篇文章：
sandwich1 <- function(object, ...) sandwich(object) * nobs(object) / (nobs(object) - 1)
coeftest(fitperc, vcov = sandwich1)
##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.80168270 0.07237747 93.9751 <2e-16 ***
## yr_rndYes 0.04825266 0.03217137 1.4999 0.1336
## parented -0.76625982 0.03907151 -19.6117 <2e-16 ***
## api99 -0.00730460 0.00021556 -33.8867 <2e-16 ***
## ---
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
对数似然和置信区间（略有不同，因为估计过程似乎不同）：
logLik(fitperc)
## 'log Lik.' -1464.363 (df=4)
confint(fitperc)
## Waiting for profiling to be done...
## 2.5 % 97.5 %
## (Intercept) 6.352788748 7.262067304
## yr_rndYes -0.155529338 0.253123151
## parented -0.944775733 -0.588903012
## api99 -0.008303668 -0.006319185
要获得预测：
meals_pred <- data.frame(api99=rep(c(500,600,700), 2),
yr_rnd=rep(c("No", "Yes"), times=1, each=3),
parented=rep(2.5, 6))
cbind(meals_pred, pred=predict(fitperc, meals_pred, "response"))
## api99 yr_rnd parented pred
## 1 500 No 2.5 0.7744710
## 2 600 No 2.5 0.6232278
## 3 700 No 2.5 0.4434458
## 4 500 Yes 2.5 0.7827873
## 5 600 Yes 2.5 0.6344891
## 6 700 Yes 2.5 0.4553849
有关相关讨论，请参阅此问题：