估计后分析与诊断
先查看各命令展开的示例,再按需展开其他选项及 R 实现。改写代码前,请阅读适用条件及差异说明。
brier
计算总体Brier均方概率评分
Stata
brier b pR 包:stats
mean((data$b-data$p)^2)与 Stata 有差异 此例仅返回总体评分;Stata的brier还报告分解分量。
dfbeta
计算删除观测对指定斜率的标准化影响
Stata
dfbeta xR 包:stats
stats::dfbetas(fit)[,"x"]与 Stata 有差异 返回dfbetas中x的标准化删点系数诊断;须与Stata模型统一系数顺序、缺失行和诊断名称。
estat
用学生化统计量检验与指定变量相关的异方差
Stata
estat hettest x, iidR 包:lmtest
lmtest::bptest(fit,varformula=~x,data=data,studentize=TRUE)与 Stata 有差异 采用学生化Koenker版Breusch-Pagan检验,对应所示iid变体。须统一辅助回归变量和估计样本;bptest默认值不同于Stata默认hettest。
另有 9 个示例及替代实现
运行经典Breusch-Pagan异方差检验
Stata
estat hettest, rhsR 包:lmtest
lmtest::bptest(fit,studentize=FALSE)与 Stata 有差异 studentize=FALSE选择正态理论版本;bptest默认使用Koenker版本。
计算惩罚似然模型评价准则
Stata
estat icR 包:stats
c(AIC=stats::AIC(fit),BIC=stats::BIC(fit))R以拟合对象保存模型;示例中的模型名称是固定流程标签。
运行White形式辅助回归异方差检验
Stata
estat imtest, whiteR 包:lmtest
lmtest::bptest(fit,varformula=~x+z+I(x^2)+I(z^2)+x:z,data=data,studentize=TRUE)此例仅对应White异方差分量,不包含完整的信息矩阵分解。
使用拟合值幂次运行Ramsey RESET检验
Stata
estat ovtestR 包:lmtest
lmtest::resettest(fit,power=2:4,type="fitted")与 Stata 有差异 按所示代码对拟合值或解释变量的二至四次幂实施RESET。因子解释变量展开和辅助项秩不足处理可能不同于Stata。
使用解释变量幂次运行RESET检验
Stata
estat ovtest, rhsR 包:lmtest
lmtest::resettest(fit,power=2:4,type="regressor")与 Stata 有差异 按所示代码对拟合值或解释变量的二至四次幂实施RESET。因子解释变量展开和辅助项秩不足处理可能不同于Stata。
汇总实际估计样本中的变量
Stata
estat summarizeR 包:stats
summary(stats::model.frame(fit))R以拟合对象保存模型;示例中的模型名称是固定流程标签。
提取系数估计之间的相关系数
Stata
estat vce, correlationR 包:stats
stats::cov2cor(stats::vcov(fit))R以拟合对象保存模型;示例中的模型名称是固定流程标签。
计算最小二乘模型方差膨胀因子
Stata
estat vifR 包:car
car::vif(fit)与 Stata 有差异 car对单自由度项返回普通VIF,对多参数项返回广义VIF。Stata的estat vif报告列层面VIF,因子项输出不能直接互换。
estimates
另有 6 个示例及替代实现
重新显示线性模型摘要
Stata
estimates replayR 包:stats
summary(fit)R以拟合对象保存模型;示例中的模型名称是固定流程标签。
从模型集合恢复指定模型
Stata
estimates restore model1R 包:stats
fit <- models[["model1"]]R以拟合对象保存模型;示例中的模型名称是固定流程标签。
保存并重新载入拟合模型对象
Stata
estimates save model1R 包:base
{ path <- tempfile(fileext=".rds"); saveRDS(fit,path); readRDS(path) }此例写入R临时文件;长期保存应选择明确的项目路径。
比较已保存模型的AIC
Stata
estimates stats model1 model2R 包:stats
do.call(stats::AIC,unname(models))R以拟合对象保存模型;示例中的模型名称是固定流程标签。
按名称保存拟合模型
Stata
estimates store model1R 包:stats
model1 <- fitR以拟合对象保存模型;示例中的模型名称是固定流程标签。
lincom
估计两个斜率之和及标准误
Stata
lincom x + zR 包:stats
{ L <- c(0,1,1); c(estimate=sum(L*stats::coef(fit)),se=sqrt(drop(t(L)%*%stats::vcov(fit)%*%L))) }与 Stata 有差异 对比向量假定系数顺序为截距、x、z。仅返回估计及按协方差计算的标准误;Stata还打印检验和置信限。
另有 1 个示例及替代实现
估计两个斜率之差及标准误
Stata
lincom x - zR 包:stats
{ L <- c(0,1,-1); c(estimate=sum(L*stats::coef(fit)),se=sqrt(drop(t(L)%*%stats::vcov(fit)%*%L))) }与 Stata 有差异 对比向量假定系数顺序为截距、x、z。仅返回估计及按协方差计算的标准误;Stata还打印检验和置信限。
linktest
lrtest
margins
在指定协变量值处计算预测
Stata
margins, at(x=(2 3 4))R 包:stats
vapply(c(2,3,4),function(value) { newdata <- data; newdata$x <- value; mean(stats::predict(fit,newdata=newdata,type="response")) },numeric(1))与 Stata 有差异 在每个x取值处,按z的实际观测分布对预测值取平均。
另有 3 个示例及替代实现
计算连续解释变量的平均边际效应
Stata
margins, dydx(x)R 包:stats
{ X <- stats::model.matrix(fit); beta <- stats::coef(fit); p <- stats::fitted(fit); gradient <- beta["x"]*colMeans(X*as.vector(p*(1-p)*(1-2*p))); gradient["x"] <- gradient["x"]+mean(p*(1-p)); c(estimate=unname(beta["x"]*mean(p*(1-p))),se=sqrt(drop(t(gradient)%*%stats::vcov(fit)%*%gradient))) }与 Stata 有差异 使用已拟合GLM及其估计数据;需明确选择参照网格与平均规则。
在协变量均值处计算边际效应
Stata
margins, dydx(x) atmeansR 包:stats
{ X <- colMeans(stats::model.matrix(fit)); beta <- stats::coef(fit); p <- stats::plogis(sum(X*beta)); c(estimate=unname(beta["x"]*p*(1-p))) }与 Stata 有差异 返回协变量均值处的导数;标准误需另外计算Delta方法梯度。
计算线性指数尺度的平均预测
Stata
margins, expression(predict(xb))R 包:stats
{ L <- colMeans(stats::model.matrix(fit)); c(estimate=sum(L*stats::coef(fit)),se=sqrt(drop(t(L)%*%stats::vcov(fit)%*%L))) }与 Stata 有差异 使用已拟合GLM及其估计数据;需明确选择参照网格与平均规则。
nlcom
使用Delta方法估计斜率比
Stata
nlcom _b[x]/_b[z]R 包:car
car::deltaMethod(fit,"x / z")与 Stata 有差异 对具名已拟合系数采用一阶delta方法。返回变换估计和标准误;Stata还打印参考检验及置信限。
另有 1 个示例及替代实现
使用Delta方法估计指数化斜率
Stata
nlcom exp(_b[x])R 包:car
car::deltaMethod(fit,"exp(x)")与 Stata 有差异 对具名已拟合系数采用一阶delta方法。返回变换估计和标准误;Stata还打印参考检验及置信限。
predict
计算Cook影响距离
Stata
predict result, cooksdR 包:stats
data$result <- stats::cooks.distance(fit)与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。
另有 12 个示例及替代实现
计算二项模型偏差残差
Stata
predict result, devianceR 包:stats
data$result <- stats::residuals(fit,type="deviance")与 Stata 有差异 在二项响应广义线性模型拟合后使用;预测选项因模型而异。
计算DFFITS影响统计量
Stata
predict result, dfitsR 包:stats
data$result <- stats::dffits(fit)与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。
计算杠杆值
Stata
predict result, hatR 包:stats
data$result <- stats::hatvalues(fit)与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。
计算二项模型Pearson残差
Stata
predict result, pearsonR 包:stats
data$result <- stats::residuals(fit,type="pearson")与 Stata 有差异 在二项响应广义线性模型拟合后使用;预测选项因模型而异。
计算逻辑回归预测概率
Stata
predict result, prR 包:stats
data$result <- stats::predict(fit,type="response")与 Stata 有差异 在二项响应广义线性模型拟合后使用;预测选项因模型而异。
计算最小二乘原始残差
Stata
predict result, residualsR 包:stats
data$result <- stats::residuals(fit)与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。
计算内部标准化残差
Stata
predict result, rstandardR 包:stats
data$result <- stats::rstandard(fit)与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。
计算外部学生化残差
Stata
predict result, rstudentR 包:stats
data$result <- stats::rstudent(fit)与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。
计算包含残差方差的预测标准误
Stata
predict result, stdfR 包:stats
data$result <- sqrt(stats::predict(fit,se.fit=TRUE)$se.fit^2+summary(fit)$sigma^2)与 Stata 有差异 将残差方差加到OLS拟合均值标准误的平方上。这是新个体预测的标准误,不是拟合均值标准误;假定同方差误差且模型样本一致。
计算拟合均值标准误
Stata
predict result, stdpR 包:stats
data$result <- stats::predict(fit,se.fit=TRUE)$se.fit与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。
计算线性拟合值
Stata
predict result, xbR 包:stats
data$result <- stats::predict(fit)与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。
计算逻辑回归线性指数
Stata
predict result, xbR 包:stats
data$result <- stats::predict(fit,type="link")与 Stata 有差异 在二项响应广义线性模型拟合后使用;预测选项因模型而异。
roccomp
使用DeLong推断比较配对标记的ROC面积
Stata
roccomp b x zR 包:pROC
pROC::roc.test(pROC::roc(data$b,data$x,levels=c(0,1),direction="<",quiet=TRUE),pROC::roc(data$b,data$z,levels=c(0,1),direction="<",quiet=TRUE),paired=TRUE,method="delong")与 Stata 有差异 配对推断要求两个标记在相同完整样本受试者上测量。
roctab
估计诊断标记的非参数ROC曲线下面积
Stata
roctab b xR 包:pROC
pROC::roc(data$b,data$x,levels=c(0,1),direction="<",quiet=TRUE)与 Stata 有差异 显式指定较大标记值表示病例;固定方向可避免按数据自动反向。
列出不同ROC阈值的敏感度与特异度
Stata
roctab b x, detailR 包:pROC
pROC::coords(pROC::roc(data$b,data$x,levels=c(0,1),direction="<",quiet=TRUE),x="all",ret=c("threshold","sensitivity","specificity"))与 Stata 有差异 在显式事件方向下返回pROC全部阈值、敏感度和特异度。阈值位置及端点约定不同于Stata的roctab detail切点。
test
检验单个回归系数
Stata
test x = 0R 包:car
car::linearHypothesis(fit,"x = 0")与 Stata 有差异 对已有模型检验所示系数约束。car根据模型类选择参考统计量;GLM的卡方检验不同于Stata regress后估计的OLS F检验。
另有 2 个示例及替代实现
检验两个回归斜率相等
Stata
test x = zR 包:car
car::linearHypothesis(fit,"x = z")与 Stata 有差异 对已有模型检验所示系数约束。car根据模型类选择参考统计量;GLM的卡方检验不同于Stata regress后估计的OLS F检验。
联合检验两个回归系数
Stata
test x zR 包:car
car::linearHypothesis(fit,c("x = 0","z = 0"))与 Stata 有差异 对已有模型检验所示系数约束。car根据模型类选择参考统计量;GLM的卡方检验不同于Stata regress后估计的OLS F检验。
testnl
用非线性Wald检验检验指数化系数
Stata
testnl exp(_b[x]) = 1R 包:car
{ d <- car::deltaMethod(fit,"exp(x)"); w <- ((d$Estimate-1)/d$SE)^2; c(chisq=w,p.value=stats::pchisq(w,df=1,lower.tail=FALSE)) }与 Stata 有差异 用一阶delta标准误构造单自由度Wald卡方。R向量不包含Stata的系数及置信区间输出。
testparm
联合检验分类解释变量
Stata
testparm i.gR 包:car
car::linearHypothesis(fit,c("g6 = 0","g8 = 0"))与 Stata 有差异 约束显式指定已拟合的因子系数。参考水平或对比编码不同时须修改系数名;不自动展开任意Stata因子变量表达式。