Merge pull request #88 from nutriverse:dev

create function for patchy test; fix #85

DESCRIPTION

@@ -30,9 +30,11 @@ Depends:
     R (>= 4.1.0)
 Imports: 
     cli,
+    methods,
     parallel,
     parallelly,
-    squeacr
+    squeacr,
+    stats
 Suggests: 
     covr,
     knitr,

NAMESPACE

@@ -2,6 +2,7 @@
 
 S3method(plot,lqasSim)
 S3method(print,lqasClass)
+export(check_coverage_homogeneity)
 export(estimate_coverage)
 export(estimate_coverage_overall)
 export(get_n_cases)
@@ -20,12 +21,16 @@ export(lqas_simulate_run)
 export(lqas_simulate_runs)
 export(lqas_simulate_test)
 importFrom(cli,cli_abort)
+importFrom(cli,cli_alert_info)
+importFrom(cli,cli_alert_warning)
 importFrom(graphics,abline)
 importFrom(graphics,legend)
 importFrom(graphics,lines)
 importFrom(graphics,plot)
 importFrom(graphics,points)
+importFrom(methods,is)
 importFrom(parallel,mclapply)
 importFrom(parallelly,availableCores)
 importFrom(stats,lowess)
+importFrom(stats,pchisq)
 importFrom(stats,phyper)

R/sleacr.R → R/00-sleacr.R

@@ -27,7 +27,9 @@
 #' @importFrom stats lowess phyper
 #' @importFrom parallel mclapply
 #' @importFrom parallelly availableCores
-#' @importFrom cli cli_abort
+#' @importFrom cli cli_abort cli_alert_warning cli_alert_info
+#' @importFrom methods is
+#' @importFrom stats pchisq
 #'
 
 "_PACKAGE"

R/06-chi-square.R

@@ -0,0 +1,123 @@
+#'
+#' Check coverage distribution
+#' 
+#' @inheritParams estimate_coverage_overall
+#' @param p Minimum p-value to test statistic. Default is 0.05.
+#' 
+#' @returns A named list of 2 lists: one for case-finding effectiveness (*cf*)
+#'   and the second for treatment coverage (*tc*). For each list, the following
+#'   values are provided:
+#'   * **statistic** - calculated chi-square statistic
+#'   * **df** - degrees of freedom
+#'   * **p** - p-value of chi-square statistic
+#' 
+#' @examples
+#' check_coverage_homogeneity(survey_data)
+#' 
+#' @export
+#' 
+
+check_coverage_homogeneity <- function(cov_df, k = 3, p = 0.05) {
+  ## Check coverage data ----
+  check_coverage_data(cov_df)
+
+  ## Check p ----
+  check_p(p)
+
+  ## Calculate chi-square statistic ----
+  x2 <- calculate_x2_stat(
+    cases_in = cov_df$cases_in, cases_out = cov_df$cases_out, 
+    rec_in = cov_df$rec_in, k = k
+  )
+
+  ## Sum of chi-square statistics ----
+  x2_cf <- sum(x2[[1]])
+  x2_tc <- sum(x2[[2]])
+
+  deg_free <- nrow(cov_df) - 1
+
+  ## Get critical value ----
+  #crit_value <- get_critical_value(df = deg_free, p = as.character(p))
+
+  p_cf <- stats::pchisq(x2_cf, df = deg_free, lower.tail = FALSE)
+  p_tc <- stats::pchisq(x2_tc, df = deg_free, lower.tail = FALSE)
+
+  ## Concatenate results ----
+  x2_results <- list(
+    cf = list(
+      statistic = x2_cf,
+      df = deg_free,
+      p = p_cf
+    ),
+    tc = list(
+      statistic = x2_tc,
+      df = deg_free,
+      p = p_tc
+    )
+  )
+
+  ## Create messages ----
+  if (p_cf < p) {
+    cli::cli_alert_warning(
+      "{.strong Case-finding effectiveness} across {nrow(cov_df)} surveys is {.strong patchy}."
+    )
+  } else {
+    cli::cli_alert_info(
+      "{.strong Case-finding effectiveness} across {nrow(cov_df)} surveys is {.strong not patchy}."
+    )
+  }
+
+  if (p_tc < p) {
+    cli::cli_alert_warning(
+      "{.strong Treatment coverage} across {nrow(cov_df)} surveys is {.strong patchy}."
+    )
+  } else {
+    cli::cli_alert_info(
+      "{.strong Treatment coverage} across {nrow(cov_df)} surveys is {.strong not patchy}."
+    )
+  }
+
+  ## Return x2_results ----
+  x2_results
+}
+
+
+#'
+#' @keywords internal
+#' 
+
+calculate_x2_stat <- function(cases_in, cases_out, rec_in, k = 3) {
+  ## Get observed ----
+  observed_cf <- cases_in
+  observed_tc <- cases_in + rec_in
+
+  ## Calculate rec_out ----
+  rec_out <- squeacr::calculate_rout(
+    cin = cases_in, cout = cases_out, rin = rec_in, k = k
+  )
+
+  ## Calculate variables needed ----
+  cases_cf <- cases_in + cases_out
+  cases_tc <- cases_in + cases_out + rec_in + rec_out
+
+  total_cases_cf <- sum(cases_cf)
+  total_cases_tc <- sum(cases_tc)
+
+  total_cases_in_cf <- sum(cases_in)
+  total_cases_in_tc <- sum(cases_in + rec_in)
+
+  ## Calculate expected ----
+  expected_cf <- cases_cf * (total_cases_in_cf / total_cases_cf)
+  expected_tc <- cases_tc * (total_cases_in_tc / total_cases_tc) 
+
+  ## Calculate chi-square ----
+  x2_cf <- ((observed_cf - expected_cf) ^ 2) / expected_cf
+  x2_tc <- ((observed_tc - expected_tc) ^ 2) / expected_tc
+
+  ## Concatenate ----
+  x2 <- list(cf = x2_cf, tc = x2_tc)
+
+  ## Return chi-square ----
+  x2
+}
+

R/utils.R

@@ -154,5 +154,19 @@ check_pop_data <- function(pop_df) {
 }
 
 
+#'
+#' Check p value
+#' 
+#' @keywords internal
+#' 
+
+check_p <- function(p) {
+  if (!is(p, "numeric"))
+    cli::cli_abort("{.arg p} should be numeric.")
+
+  if (p >= 1 | p <= 0)
+    cli::cli_abort("{.arg p} should be greater than 0 and less than 1")
+}
+
 
 

README.Rmd

@@ -125,7 +125,7 @@ Using this dataset, per district coverage classifications can be calculated as f
 ```{r classify-coverage, eval = FALSE}
 with(
   survey_data, 
-  lqas_classify_coverage(
+  lqas_classify(
     cases_in = cases_in, cases_out = cases_out, rec_in = rec_in
   )
 )
@@ -136,7 +136,7 @@ which outputs the following results:
 ```{r classify-coverage-show, echo = FALSE}
 with(
   survey_data, 
-  lqas_classify_coverage(
+  lqas_classify(
     cases_in = cases_in, cases_out = cases_out, rec_in = rec_in
   )
 )

README.md

@@ -249,29 +249,29 @@ calculated as follows:
 ``` r
 with(
   survey_data, 
-  lqas_classify_coverage(
+  lqas_classify(
     cases_in = cases_in, cases_out = cases_out, rec_in = rec_in
   )
 )
 ```
 
 which outputs the following results:
 
-    #>     cf       tc
-    #> 1  Low Moderate
-    #> 2  Low      Low
-    #> 3  Low      Low
-    #> 4  Low      Low
-    #> 5  Low      Low
-    #> 6  Low Moderate
-    #> 7  Low      Low
-    #> 8  Low Moderate
-    #> 9  Low Moderate
-    #> 10 Low Moderate
-    #> 11 Low      Low
-    #> 12 Low      Low
-    #> 13 Low      Low
-    #> 14 Low      Low
+    #>    cf tc
+    #> 1   0  1
+    #> 2   0  0
+    #> 3   0  0
+    #> 4   0  0
+    #> 5   0  0
+    #> 6   0  1
+    #> 7   0  0
+    #> 8   1  1
+    #> 9   1  1
+    #> 10  1  1
+    #> 11  0  0
+    #> 12  0  0
+    #> 13  0  0
+    #> 14  0  0
 
 The function provides estimates for *case-finding effectiveness* and for
 *treatment coverage* as a `data.frame` object.
@@ -294,10 +294,10 @@ lqas_sim_pop <- lqas_simulate_test(
 
 ## Get classification probabilities ----
 lqas_get_class_prob(lqas_sim_pop)
-#>                     Low : 0.9549
-#>                Moderate : 0.8308
-#>                    High : 0.8317
-#>                 Overall : 0.9054
+#>                     Low : 0.9552
+#>                Moderate : 0.8305
+#>                    High : 0.8395
+#>                 Overall : 0.9062
 #> Gross misclassification : 0
 ```
 

data-raw/chi_square.R

@@ -0,0 +1,12 @@
+# Retrieve chi-square statistic table ------------------------------------------
+
+session <- rvest::session(
+  "https://people.richland.edu/james/lecture/m170/tbl-chi.html"
+)
+
+chi_square_table <- rvest::html_table(session) |>
+  (\(x) x[[1]])()
+
+usethis::use_data(
+  chi_square_table, internal = TRUE, overwrite = TRUE, compress = "xz"
+)

inst/WORDLIST

@@ -98,6 +98,7 @@ Yalenga
 Yarlenga
 Yengema
 cf
+df
 frac
 lceil
 lfloor

pkgdown/_pkgdown.yml

@@ -2,8 +2,8 @@ title: sleacr
 
 url: https://nutriverse.io/sleacr/
 
-development:
-  mode: auto
+# development:
+#   mode: auto
 
 template:
   bootstrap: 5
@@ -73,6 +73,10 @@ reference:
     contents:
     - starts_with("estimate_coverage")
 
+  - title: Test for coverage homogeneity
+    contents:
+    - check_coverage_homogeneity
+
   - title: Datasets
     contents:
     - village_list