调查包中的剩余标准误差

Residual standard error in survey package

我正在尝试使用调查包计算线性回归模型的残差标准误差。我正在处理一个复杂的设计,复杂设计的采样权重由下面代码中的 "weight" 给出。

fitM1 <- lm(med~x1+x2,data=pop_sample,weight=weight)  
fitM2 <- svyglm(med~x1+x2,data=pop_sample,design=design)

首先,如果我调用 "summary(fitM1)",我会得到以下信息:

Call: lm(formula=med~x1+x2,data=pop_sample,weights=weight)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.001787   0.042194   0.042    0.966    
x1           0.382709   0.061574   6.215 1.92e-09 ***
x2           0.958675   0.048483  19.773  < 2e-16 ***

Residual standard error: 9.231 on 272 degrees of freedom
Multiple R-squared:  0.8958,    Adjusted R-squared:  0.8931 
F-statistic: 334.1 on 7 and 272 DF,  p-value: < 2.2e-16

接下来,如果我调用 "summary(fitM2)",我会得到以下信息:

summary(fitM2)

Call: svyglm(formula=med~x1+x2,data=pop_sample,design=design)

Survey design: svydesign(id=~id_cluster,strat=~id_stratum,weight=weight,data=pop_sample)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.001787   0.043388   0.041 0.967878    
x1           0.382709   0.074755   5.120 0.000334 ***
x2           0.958675   0.041803  22.933 1.23e-10 ***

当使用"lm"时,我可以通过调用提取残差标准误差:

fitMvariance <- summary(fitM1)$sigma^2

但是,我在调查包的任何地方都找不到 "svyglm" 的类似函数。比较这两种方法时,点估计是相同的,但系数的标准误差(以及模型的残差标准误差)是不同的。

调查分析

使用 r 中的库调查来执行调查分析,它提供了广泛的函数来计算统计数据,例如百分比、下限 CI、上限 CI、人口和 RSE。

RSE

我们可以使用 survey 包中的svyby函数来获取所有统计数据,包括平方根误差

library("survey")
Survey design: svydesign(id=~id_cluster,strat=~id_stratum,weight=weight,data=pop_sample)
  
svyby(~med, ~x1+x2, design, svytotal, deff=TRUE, verbose=TRUE,vartype=c("se","cv","cvpct","var"))

cvpct 将给出平方根误差

参考了解更多信息svyby

因为 svyglm 建立在 glm 而不是 lm 的基础上,方差估计被称为 $dispersion 而不是 $sigma

> data(api)
> dstrat<-svydesign(id = ~1, strata = ~stype, weights = ~pw, data = apistrat, 
+     fpc = ~fpc)
> model<-svyglm(api00~ell+meals+mobility, design=dstrat)
> summary(model)$dispersion
     variance     SE
[1,]     5172 492.28

这是$\sigma^2$的估计值,也就是总体残差。在这个例子中我们实际上有整个人口,所以我们可以比较

> popmodel<-lm(api00~ell+meals+mobility, data=apipop)
> summary(popmodel)$sigma
[1] 70.58365
> sqrt(5172)
[1] 71.91662