Add davidfunctions

R/direct_estimation.R

@@ -104,6 +104,7 @@ direct <- function(y,
                    design = NULL,
                    threshold = NULL,
                    var = FALSE, 
+                   HT = FALSE, 
                    boot_type = "naive", 
                    B = 50,
                    seed = 123,
@@ -112,7 +113,7 @@ direct <- function(y,
                    custom_indicator = NULL,
                    na.rm = FALSE){
 
-
+  smp_data <- as.data.frame(smp_data)
 
   direct_check(y = y, smp_data = smp_data, smp_domains = smp_domains, 
                 weights = weights, design = design, threshold = threshold, 
@@ -161,6 +162,7 @@ direct <- function(y,
                            smp_data = framework$smp_data,  
                            smp_domains = framework$smp_domains_vec,
                            design = design,
+                           HT = HT, 
                            bootType = boot_type, 
                            B = B,
                            seed = seed,

R/direct_variance.R

@@ -5,7 +5,8 @@ direct_variance <- function(direct_estimator,
                             smp_data,
                             smp_domains,
                             design,
-                            indicator, 
+                            indicator,
+                            HT, 
                             bootType, 
                             B = B,
                             seed,
@@ -50,6 +51,8 @@ direct_variance <- function(direct_estimator,
     B <- as.integer(B[1])
   }
 
+
+
   bootType #<- match.arg(bootType)
 
   # if calibrate bootstrap is selected
@@ -73,6 +76,38 @@ direct_variance <- function(direct_estimator,
   smp_data$weight <- weights
   smp_data$Domain <- smp_domains
 
+  # separate code for HT estimation 
+  if (HT == TRUE) {
+
+
+
+    domain_var <- function(df) {
+      tosum <-df[,2]*(df[,2]-1)*df[,1]^2 
+      nrow(df)^-2*sum(tosum)
+    }
+
+
+
+  if (indicator_name == "Mean") {
+    smp_data$indicator <- smp_data$y
+  }
+  else if (indicator_name == "Head_Count") {
+  smp_data$indicator <- as.integer(smp_data$y<threshold)
+}
+  else {
+    smp_data$indicator <- NA 
+  }
+
+    var <- as.vector(by(data=smp_data[c("indicator","weight")],INDICES=smp_data$Domain,FUN=domain_var))
+
+
+    varByDomain <- data.frame(Domain = rs, var = var)
+    indicator$varMethod <- "HT"
+  }
+
+  else {
+
+
   # set seed for bootstrap
   if (!is.null(seed)) {
     set.seed(seed)
@@ -97,7 +132,6 @@ direct_variance <- function(direct_estimator,
                      envir = envir,
                      indicator_name = indicator_name)
                      #, ...)
-
   # if variance is calculated by domain
   if (byDomain) {
     var <- apply(b$t, 2, var, na.rm = TRUE)
@@ -106,9 +140,11 @@ direct_variance <- function(direct_estimator,
   } else {
     var <- var(b$t[, 1], na.rm = TRUE)
   }
+
 
   # preparation of return
   indicator$varMethod <- "bootstrap"
+  } # close HT else loop 
   indicator$var <- var
   if (byDomain) {
     indicator$varByDomain <- varByDomain