使用 svyrecvar 获取“调查”R 包中统计量的方差

我正在写一个扩展surveyR 包。我的函数使用复杂的调查数据估计峰度，但我不确定如何获得这些统计数据的标准误差，例如现有函数的输出，例如svymean and svytotal.

svykurt <- function(
    x,
    design,
    na.rm = FALSE,
    excess = TRUE
) {

  if (!inherits(design, "survey.design"))
    stop("design is not a survey design")

  x <- model.frame(x, design$variables, na.action = na.pass)
  x <- as.matrix(x)

  if (ncol(x) > 1)
    stop("Only calculate kurtosis one variable at a time")

  if(na.rm){
    x <- x[!is.na(x)]
  }

  pweights <- 1/design$prob
  psum <- sum(pweights)
  mean_x <- svymean(x, design, na.rm = na.rm)
  var_x <- svyvar(x, design, na.rm = na.rm)

  m4 <- sum(pweights * (x - mean_x)^4) / psum
  kurt <- m4 / var_x^2

  if (excess) {
    kurt <- kurt - 3
  }

  class(kurt) <- "svykurt"

  return(kurt)

}

print.svykurt <- function(x) {
  m <- as.matrix(x, ncol = 1)
  rownames(m) <- names(x)
  colnames(m) <- "kurtosis"

  print(m)
}

目前，我只是打印没有标准错误的峰度。

data(api)
dclus1<-svydesign(id=~dnum, weights=~pw, data=apiclus1, fpc=~fpc)
svykurt(x = ~api00, dclus1, na.rm = TRUE)

但是，我最终想将输出分配给类svystat并打印包含这两个值的矩阵，如svymean(x = ~api00, dclus1, na.rm = TRUE).

看着survey 源代码 https://rdrr.io/cran/survey/src/R/multistage.R， 我理解了svymean and svytotal使用一个更高级别的函数，称为svrecvar, e.g.,

attr(total, "var")<-v<-svyrecvar(x/design$prob,design$cluster,
                                     design$strata, design$fpc,
                                   postStrata=design$postStrata)

我只是不确定如何将其应用到我的案例中。

如果你想要峰度，你最好与svycontrast比与svyrecvar。您可以根据原始矩写出方差和第四中心矩（在答案之后）here https://stats.stackexchange.com/questions/237558/kurtosis-expressed-in-raw-moments）然后变换

> moments<-svymean(~enroll+I(enroll^2)+I(enroll^3)+I(enroll^4), dstrat)
> moments
                  mean         SE
enroll      5.9528e+02 1.8509e+01
I(enroll^2) 5.4895e+05 4.0250e+04
I(enroll^3) 7.5137e+08 1.0066e+08
I(enroll^4) 1.3348e+12 2.7683e+11
> 
> central_moments<-svycontrast(moments, list(mu4=quote(-3*enroll^4+6*enroll^2*`I(enroll^2)`-4*enroll*`I(enroll^3)`+`I(enroll^4)`),
sigma2=quote(`I(enroll^2)`-enroll^2)))
> central_moments
            nlcon         SE
mu4    3.3618e+11 1.0458e+11
sigma2 1.9459e+05 2.3116e+04
> 
> 
> svycontrast(central_moments, quote(mu4/(sigma2*sigma2)))
          nlcon     SE
contrast 8.8783 1.5678

作为一般指南，如果您有关于矩或单元格计数或比例的明确统计公式，则更容易使用svycontrast如果您有估计方程或最大化问题，则更容易使用svyrecvar.

To use svyrecvar您需要计算或查找峰度的影响函数，这可能是可行的，但有点烦人。如果你确实计算出了影响函数，那么使用起来实际上更简单svytotal计算方差

例如，在svyVGAM:::svy_vglm.survey.design有很多代码来拟合 VGAM 模型并计算出其影响函数，接下来是

    v <- vcov(svytotal(inffuns, design))
    dimnames(v) <- list(names(coef(fit)), names(coef(fit)))

进行方差计算。