估计后分析与诊断

先查看各命令展开的示例,再按需展开其他选项及 R 实现。改写代码前,请阅读适用条件及差异说明。

brier

计算总体Brier均方概率评分

Stata

brier b p

R 包:stats

mean((data$b-data$p)^2)

与 Stata 有差异 此例仅返回总体评分;Stata的brier还报告分解分量。

dfbeta

计算删除观测对指定斜率的标准化影响

Stata

dfbeta x

R 包:stats

stats::dfbetas(fit)[,"x"]

与 Stata 有差异 返回dfbetas中x的标准化删点系数诊断;须与Stata模型统一系数顺序、缺失行和诊断名称。

estat

用学生化统计量检验与指定变量相关的异方差

Stata

estat hettest x, iid

R 包:lmtest

lmtest::bptest(fit,varformula=~x,data=data,studentize=TRUE)

与 Stata 有差异 采用学生化Koenker版Breusch-Pagan检验,对应所示iid变体。须统一辅助回归变量和估计样本;bptest默认值不同于Stata默认hettest。

另有 9 个示例及替代实现

运行经典Breusch-Pagan异方差检验

Stata

estat hettest, rhs

R 包:lmtest

lmtest::bptest(fit,studentize=FALSE)

与 Stata 有差异 studentize=FALSE选择正态理论版本;bptest默认使用Koenker版本。

计算惩罚似然模型评价准则

Stata

estat ic

R 包:stats

c(AIC=stats::AIC(fit),BIC=stats::BIC(fit))

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

运行White形式辅助回归异方差检验

Stata

estat imtest, white

R 包:lmtest

lmtest::bptest(fit,varformula=~x+z+I(x^2)+I(z^2)+x:z,data=data,studentize=TRUE)

此例仅对应White异方差分量,不包含完整的信息矩阵分解。

使用拟合值幂次运行Ramsey RESET检验

Stata

estat ovtest

R 包:lmtest

lmtest::resettest(fit,power=2:4,type="fitted")

与 Stata 有差异 按所示代码对拟合值或解释变量的二至四次幂实施RESET。因子解释变量展开和辅助项秩不足处理可能不同于Stata。

使用解释变量幂次运行RESET检验

Stata

estat ovtest, rhs

R 包:lmtest

lmtest::resettest(fit,power=2:4,type="regressor")

与 Stata 有差异 按所示代码对拟合值或解释变量的二至四次幂实施RESET。因子解释变量展开和辅助项秩不足处理可能不同于Stata。

汇总实际估计样本中的变量

Stata

estat summarize

R 包:stats

summary(stats::model.frame(fit))

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

提取回归系数协方差矩阵

Stata

estat vce

R 包:stats

stats::vcov(fit)

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

提取系数估计之间的相关系数

Stata

estat vce, correlation

R 包:stats

stats::cov2cor(stats::vcov(fit))

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

计算最小二乘模型方差膨胀因子

Stata

estat vif

R 包:car

car::vif(fit)

与 Stata 有差异 car对单自由度项返回普通VIF,对多参数项返回广义VIF。Stata的estat vif报告列层面VIF,因子项输出不能直接互换。

estimates

列出已保存模型名称

Stata

estimates dir

R 包:stats

names(models)

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

另有 6 个示例及替代实现

重新显示线性模型摘要

Stata

estimates replay

R 包:stats

summary(fit)

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

从模型集合恢复指定模型

Stata

estimates restore model1

R 包:stats

fit <- models[["model1"]]

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

保存并重新载入拟合模型对象

Stata

estimates save model1

R 包:base

{ path <- tempfile(fileext=".rds"); saveRDS(fit,path); readRDS(path) }

此例写入R临时文件;长期保存应选择明确的项目路径。

比较已保存模型的AIC

Stata

estimates stats model1 model2

R 包:stats

do.call(stats::AIC,unname(models))

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

按名称保存拟合模型

Stata

estimates store model1

R 包:stats

model1 <- fit

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

比较不同模型的系数与标准误

Stata

estimates table model1 model2, b se

R 包:stats

lapply(models,function(m) summary(m)$coefficients[,1:2,drop=FALSE])

R以拟合对象保存模型;示例中的模型名称是固定流程标签。

lincom

估计两个斜率之和及标准误

Stata

lincom x + z

R 包: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 - z

R 包: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

通过拟合值及平方项检查模型设定

Stata

linktest

R 包:stats

{ data$hat <- stats::fitted(fit); stats::lm(y ~ hat+I(hat^2),data=data) }

此例用于最小二乘模型;其他模型的连接检验需重新拟合对应分布族。

lrtest

比较嵌套的最大似然模型

Stata

lrtest reduced full

R 包:lmtest

lmtest::lrtest(reduced,full)

与 Stata 有差异 固定流程。

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) atmeans

R 包: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, cooksd

R 包:stats

data$result <- stats::cooks.distance(fit)

与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。

另有 12 个示例及替代实现

计算二项模型偏差残差

Stata

predict result, deviance

R 包:stats

data$result <- stats::residuals(fit,type="deviance")

与 Stata 有差异 在二项响应广义线性模型拟合后使用;预测选项因模型而异。

计算DFFITS影响统计量

Stata

predict result, dfits

R 包:stats

data$result <- stats::dffits(fit)

与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。

计算杠杆值

Stata

predict result, hat

R 包:stats

data$result <- stats::hatvalues(fit)

与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。

计算二项模型Pearson残差

Stata

predict result, pearson

R 包:stats

data$result <- stats::residuals(fit,type="pearson")

与 Stata 有差异 在二项响应广义线性模型拟合后使用;预测选项因模型而异。

计算逻辑回归预测概率

Stata

predict result, pr

R 包:stats

data$result <- stats::predict(fit,type="response")

与 Stata 有差异 在二项响应广义线性模型拟合后使用;预测选项因模型而异。

计算最小二乘原始残差

Stata

predict result, residuals

R 包:stats

data$result <- stats::residuals(fit)

与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。

计算内部标准化残差

Stata

predict result, rstandard

R 包:stats

data$result <- stats::rstandard(fit)

与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。

计算外部学生化残差

Stata

predict result, rstudent

R 包:stats

data$result <- stats::rstudent(fit)

与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。

计算包含残差方差的预测标准误

Stata

predict result, stdf

R 包:stats

data$result <- sqrt(stats::predict(fit,se.fit=TRUE)$se.fit^2+summary(fit)$sigma^2)

与 Stata 有差异 将残差方差加到OLS拟合均值标准误的平方上。这是新个体预测的标准误,不是拟合均值标准误;假定同方差误差且模型样本一致。

计算拟合均值标准误

Stata

predict result, stdp

R 包:stats

data$result <- stats::predict(fit,se.fit=TRUE)$se.fit

与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。

计算线性拟合值

Stata

predict result, xb

R 包:stats

data$result <- stats::predict(fit)

与 Stata 有差异 先运行fit <- lm(y ~ x + z, data=data)。

计算逻辑回归线性指数

Stata

predict result, xb

R 包:stats

data$result <- stats::predict(fit,type="link")

与 Stata 有差异 在二项响应广义线性模型拟合后使用;预测选项因模型而异。

roccomp

使用DeLong推断比较配对标记的ROC面积

Stata

roccomp b x z

R 包: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 x

R 包:pROC

pROC::roc(data$b,data$x,levels=c(0,1),direction="<",quiet=TRUE)

与 Stata 有差异 显式指定较大标记值表示病例;固定方向可避免按数据自动反向。

列出不同ROC阈值的敏感度与特异度

Stata

roctab b x, detail

R 包: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 = 0

R 包:car

car::linearHypothesis(fit,"x = 0")

与 Stata 有差异 对已有模型检验所示系数约束。car根据模型类选择参考统计量;GLM的卡方检验不同于Stata regress后估计的OLS F检验。

另有 2 个示例及替代实现

检验两个回归斜率相等

Stata

test x = z

R 包:car

car::linearHypothesis(fit,"x = z")

与 Stata 有差异 对已有模型检验所示系数约束。car根据模型类选择参考统计量;GLM的卡方检验不同于Stata regress后估计的OLS F检验。

联合检验两个回归系数

Stata

test x z

R 包:car

car::linearHypothesis(fit,c("x = 0","z = 0"))

与 Stata 有差异 对已有模型检验所示系数约束。car根据模型类选择参考统计量;GLM的卡方检验不同于Stata regress后估计的OLS F检验。

testnl

用非线性Wald检验检验指数化系数

Stata

testnl exp(_b[x]) = 1

R 包: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.g

R 包:car

car::linearHypothesis(fit,c("g6 = 0","g8 = 0"))

与 Stata 有差异 约束显式指定已拟合的因子系数。参考水平或对比编码不同时须修改系数名;不自动展开任意Stata因子变量表达式。