svyciprop adjustment

DESCRIPTION

@@ -1,9 +1,9 @@
 Package: surveytable
-Title: Streamlined output from the `survey` package
+Title: Formatted Survey Estimates
 Version: 0.9
 Authors@R: 
     person("Alex", "Strashny", , "[email protected]", role = c("aut", "cre"))
-Description: Streamlined output from the `survey` package.
+Description: Formatted Survey Estimates.
 Date/Publication: 2023
 License: Apache License (>= 2)
 Encoding: UTF-8

NAMESPACE

@@ -8,6 +8,7 @@ export(show_options)
 export(show_output)
 export(show_survey)
 export(survey_subset)
+export(svyciprop_adjusted)
 export(tab)
 export(tab_cross)
 export(tab_rate)

R/svyciprop_adjusted.R

@@ -0,0 +1,79 @@
+#' Confidence intervals for proportions, adjusted for degrees of freedom
+#'
+#' A version of `survey::svyciprop` that adjusts for the degrees of freedom when `method == "beta"`.
+#'
+#' Written by Makram Talih in 2019.
+#'
+#' To use this function in tabulations, type: `options(surveytable.adjust_svyciprop = TRUE)`.
+#'
+#' @param formula see `survey::svyciprop`.
+#' @param design see `survey::svyciprop`.
+#' @param method see `survey::svyciprop`.
+#' @param level see `survey::svyciprop`.
+#' @param df_method how `df` should be calculated: "default" or "NHIS".
+#' @param ... see `survey::svyciprop`.
+#'
+#' `df_method`: for "default", `df = degf(design)`; for "NHIS", `df = nrow(design) - 1`.
+#'
+#' @return The point estimate of the proportion, with the confidence interval as an attribute.
+#' @export
+#'
+#' @examples
+#' set_survey("vars2019")
+#' options(surveytable.adjust_svyciprop = TRUE)
+#' tab("AGER")
+#' options(surveytable.adjust_svyciprop = FALSE)
+#' tab("AGER")
+svyciprop_adjusted = function(formula
+                    , design
+                    , method = c("logit", "likelihood", "asin", "beta", "mean")
+                    , level = 0.95
+                    , df_method
+                    , ...) {
+  assert_that(df_method %in% c("default", "NHIS"))
+  df = switch(df_method,
+     default = degf(design),
+     NHIS = nrow(design) - 1
+  )
+  method = match.arg(method)
+  if (method == "beta") {
+    m = eval(bquote(svymean(~as.numeric(.(formula[[2]])), design, ...)))
+    rval = coef(m)[1]
+
+    #Effective sample size
+    n.eff = coef(m) * (1 - coef(m))/stats::vcov(m)
+
+    attr(rval, "var") = stats::vcov(m)
+    alpha = 1 - level
+
+    # Degrees of freedom used only for adjusting effective sample size, below
+
+    # For NHIS, provisional guideline is to override the R default
+    # df = nrow(design) - 1  ## uncomment this row to override R default
+
+    # R-default from degf(design) uses subdomain Strata and PSUs
+
+    if (df >0) {     #Korn-Graubard df-adjustment factor
+      rat.squ = (qt(alpha/2, nrow(design) - 1)/qt(alpha/2, df))^2
+    } else {
+      rat.squ = 0 # limit case: set to zero
+    }
+
+    if (rval > 0) {  #Adjusted effective sample size
+      n.eff = min(nrow(design), n.eff*rat.squ)
+    } else {
+      n.eff = nrow(design) #limit case: set to sample size
+    }
+
+    ci = c(stats::qbeta(alpha/2, n.eff*rval, n.eff*(1-rval)+1), stats::qbeta(1-alpha/2, n.eff*rval+1, n.eff*(1 - rval)))
+  }
+  else {
+    svyciprop(formula, design, method, level, df, ...)
+  }
+  halfalpha = (1 - level)/2
+  names(ci) = paste(round(c(halfalpha, (1 - halfalpha))*100, 1), "%", sep = "")
+  names(rval) = deparse(formula[[2]])
+  attr(rval, "ci") = ci
+  class(rval) = "svyciprop"
+  rval
+}

R/tab.R

@@ -147,7 +147,7 @@ tab = function(...
 	##
 	counts = svyby(frm, frm, design, unwtd.count)$counts
 	assert_that(length(counts) == nlv)
-	if (getOption("surveytable.do_present")) {
+	if (getOption("surveytable.check_present")) {
 	  pro = getOption("surveytable.present_restricted") %>% do.call(list(counts))
 	} else {
 		pro = list(flags = rep("", length(counts)), has.flag = c())
@@ -170,7 +170,7 @@ tab = function(...
   mmcr$ll = exp(mmcr$lnx - mmcr$k)
   mmcr$ul = exp(mmcr$lnx + mmcr$k)
 
-	if (getOption("surveytable.do_present")) {
+	if (getOption("surveytable.check_present")) {
 	  pco = getOption("surveytable.present_count") %>% do.call(list(mmcr))
 	} else {
 		pco = list(flags = rep("", nrow(mmcr)), has.flag = c())
@@ -188,7 +188,12 @@ tab = function(...
 		design$variables$.tmp = NULL
 		design$variables$.tmp = (design$variables[,vr] == lv)
 		# Korn and Graubard, 1998
-		xp = svyciprop(~ .tmp, design, method="beta", level = 0.95)
+		xp = if ( getOption("surveytable.adjust_svyciprop") ) {
+		  svyciprop_adjusted(~ .tmp, design, method="beta", level = 0.95
+		      , df_method = getOption("surveytable.adjust_svyciprop.df_method"))
+		} else {
+		  svyciprop(~ .tmp, design, method="beta", level = 0.95)
+		}
 		ret1 = data.frame(Proportion = xp %>% as.numeric
 		                  , SE = attr(xp, "var") %>% as.numeric %>% sqrt)
 
@@ -204,7 +209,7 @@ tab = function(...
 	}
 	ret$degf = df1
 
-	if (getOption("surveytable.do_present")) {
+	if (getOption("surveytable.check_present")) {
 	  ppo = getOption("surveytable.present_prop") %>% do.call(list(ret))
 	} else {
 		nlvs = design$variables[, vr] %>% nlevels
@@ -238,8 +243,10 @@ tab = function(...
 }
 
 .add_flags = function(df1, has.flag) {
-	if (is.null(has.flag)) {
-	  attr(df1, "footer") = "(No flags.)"
+  if (!getOption("surveytable.check_present")) {
+    attr(df1, "footer") = NULL
+  } else if (is.null(has.flag)) {
+	  attr(df1, "footer") = "(Checked presentation standards. Nothing to report.)"
 	} else {
 		v1 = c()
 		for (ff in has.flag) {

R/total.R

@@ -26,7 +26,7 @@ total = function(screen = getOption("surveytable.screen")
 
   ##
   counts = nrow(design$variables)
-  if (getOption("surveytable.do_present")) {
+  if (getOption("surveytable.check_present")) {
     pro = getOption("surveytable.present_restricted") %>% do.call(list(counts))
   } else {
     pro = list(flags = rep("", length(counts)), has.flag = c())
@@ -49,7 +49,7 @@ total = function(screen = getOption("surveytable.screen")
   mmcr$ll = exp(mmcr$lnx - mmcr$k)
   mmcr$ul = exp(mmcr$lnx + mmcr$k)
 
-  if (getOption("surveytable.do_present")) {
+  if (getOption("surveytable.check_present")) {
     pco = getOption("surveytable.present_count") %>% do.call(list(mmcr))
   } else {
     pco = list(flags = rep("", nrow(mmcr)), has.flag = c())

R/zzz.R

@@ -32,7 +32,7 @@
     surveytable.survey = ""
     , surveytable.survey_label = ""
 
-    , surveytable.do_present = TRUE
+    , surveytable.check_present = TRUE
     , surveytable.present_restricted = ".present_restricted"
     , surveytable.present_count = ".present_count"
     , surveytable.present_prop = ".present_prop"
@@ -50,6 +50,9 @@
 
     , surveytable.rate_per = 100
     , surveytable.tx_rate = ".tx_rate"
+
+    , surveytable.adjust_svyciprop = FALSE
+    , surveytable.adjust_svyciprop.df_method = "NHIS"
   )
   set_count_1k()
 #  set_output(csv = "", screen = TRUE, max_levels = 20)

README.Rmd

@@ -13,7 +13,7 @@ knitr::opts_chunk$set(
 )
 ```
 
-# Survey table: streamlined output from the `survey` package
+# Survey table: formatted survey estimates
 
 <!-- badges: start -->
 <!-- badges: end -->
@@ -29,11 +29,8 @@ The `surveytable` package is easier to use than using `survey` directly. With fe
 You can install `surveytable` like so:
 
 ``` r
-install.packages("remotes")
-remotes::install_git(
-	url = "https://github.com/CDCgov/surveytable"
-	, upgrade = "never")
-
+install.packages(c("remotes", "git2r"))
+remotes::install_github("CDCgov/surveytable", upgrade = "never")
 ```
 
 ## Documentation

README.md

@@ -1,7 +1,7 @@
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
-# Survey table: streamlined output from the `survey` package
+# Survey table: formatted survey estimates
 
 <!-- badges: start -->
 <!-- badges: end -->
@@ -36,10 +36,8 @@ package also performs checking for presentation standards.
 You can install `surveytable` like so:
 
 ``` r
-install.packages("remotes")
-remotes::install_git(
-    url = "https://github.com/CDCgov/surveytable"
-    , upgrade = "never")
+install.packages(c("remotes", "git2r"))
+remotes::install_github("CDCgov/surveytable", upgrade = "never")
 ```
 
 ## Documentation
@@ -290,7 +288,7 @@ Under 15 years
 </tr>
 <tr>
 <td colspan="9" style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0.8pt 0pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">
-(No flags.)
+(Checked presentation standards. Nothing to report.)
 </td>
 </tr>
 </table>

_pkgdown.yml

@@ -35,3 +35,4 @@ reference:
   - show_options
   - surveytable-options
   - survey_subset
+  - svyciprop_adjusted