context("FGT decomposition output survey.design and svyrep.design")

skip_on_cran()

library(laeken)
library(survey)

data(api)
apistrat[, sapply(apistrat, is.integer)] <-
  apply(apistrat[, sapply(apistrat, is.integer)], 2, as.numeric)
dstrat1 <- convey_prep(svydesign(id =  ~ 1, data = apistrat))
for (this_thresh in c("abs" , "relm" , "relq")) {
  for (this_g in 2:3) {
    test_that("svyfgtdec works on unweighted designs", {
      svyfgtdec(
        ~ api00,
        design = dstrat1,
        g = this_g,
        type_thresh = this_thresh,
        percent = 1,
        abs_thresh = 600 ,
        na.rm  = FALSE
      )
    })
  }
}

test_that("output svyfgtdec", {
  skip_on_cran()

  data(eusilc)
  names(eusilc) <- tolower(names(eusilc))
  eusilc[, sapply(eusilc, is.integer)] <-
    apply(eusilc[, sapply(eusilc, is.integer)], 2, as.numeric)

  des_eusilc <-
    svydesign(
      ids = ~ rb030,
      strata =  ~ db040,
      weights = ~ rb050,
      data = eusilc
    )
  des_eusilc <- convey_prep(des_eusilc)
  des_eusilc_rep <-
    as.svrepdesign(des_eusilc, type = "bootstrap" , replicates = 20)
  des_eusilc_rep <- convey_prep(des_eusilc_rep)

  # database-backed design
  library(RSQLite)
  library(DBI)
  dbfile <- tempfile()
  conn <- dbConnect(RSQLite::SQLite() , dbfile)
  dbWriteTable(conn , 'eusilc' , eusilc)

  dbd_eusilc <-
    svydesign(
      ids = ~ rb030 ,
      strata = ~ db040 ,
      weights = ~ rb050 ,
      data = "eusilc",
      dbname = dbfile,
      dbtype = "SQLite"
    )

  dbd_eusilc <- convey_prep(dbd_eusilc)

  # create a hacky database-backed svrepdesign object
  # mirroring des_eusilc_rep
  dbd_eusilc_rep <-
    svrepdesign(
      weights = ~ rb050,
      repweights = des_eusilc_rep$repweights ,
      scale = des_eusilc_rep$scale ,
      rscales = des_eusilc_rep$rscales ,
      type = "bootstrap" ,
      data = "eusilc" ,
      dbtype = "SQLite" ,
      dbname = dbfile ,
      combined.weights = FALSE
    )

  dbd_eusilc_rep <- convey_prep(dbd_eusilc_rep)


  for (this_thresh in c("abs" , "relm" , "relq")) {
    for (this_g in 2:3) {
      a1 <-
        svyfgtdec(
          ~ eqincome,
          design = des_eusilc,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE
        )
      a2 <-
        svyby(
          ~ eqincome,
          by = ~ rb090,
          design = des_eusilc,
          FUN = svyfgtdec,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE,
          deff = FALSE
        )
      b1 <-
        svyfgtdec(
          ~ eqincome,
          design = des_eusilc_rep,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE
        )
      b2 <-
        svyby(
          ~ eqincome,
          by = ~ rb090,
          design = des_eusilc_rep,
          FUN = svyfgtdec,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE,
          deff = FALSE
        )

      se_dif1 <- max(abs(SE(a1) - SE(b1)))
      se_diff2 <- max(abs(SE(a2) - SE(b2)))

      expect_is(coef(a1), "numeric")
      expect_is(coef(a2), "numeric")
      expect_is(coef(b1), "numeric")
      expect_is(coef(b2), "numeric")
      expect_equal(coef(a1), coef(b1))
      expect_equal(coef(a2), coef(b2))
      expect_is(SE(a1), "numeric")
      expect_is(SE(a2), "svyby")
      expect_is(SE(b1), "numeric")
      expect_is(SE(b2), "svyby")
      expect_lte(confint(a1)[1, 1], coef(a1)[1])
      expect_gte(confint(a1)[1, 2], coef(a1)[1])
      expect_lte(confint(a1)[2, 1], coef(a1)[2])
      expect_gte(confint(a1)[2, 2], coef(a1)[2])
      expect_lte(confint(a1)[3, 1], coef(a1)[3])
      expect_gte(confint(a1)[3, 2], coef(a1)[3])
      expect_lte(confint(a1)[4, 1], coef(a1)[4])
      expect_gte(confint(a1)[4, 2], coef(a1)[4])
      expect_lte(confint(a1)[1, 1], coef(a1)[1])
      expect_gte(confint(b1)[1, 2], coef(b1)[1])
      expect_lte(confint(b1)[2, 1], coef(b1)[2])
      expect_gte(confint(b1)[2, 2], coef(b1)[2])
      expect_lte(confint(b1)[3, 1], coef(b1)[3])
      expect_gte(confint(b1)[3, 2], coef(b1)[3])
      expect_lte(confint(b1)[4, 1], coef(b1)[4])
      expect_gte(confint(b1)[4, 2], coef(b1)[4])
      expect_equal(sum(confint(a2)[, 1] <= coef(a2)), length(coef(a2)))
      expect_equal(sum(confint(a2)[, 2] >= coef(a2)), length(coef(a2)))
      expect_equal(sum(confint(b2)[, 1] <= coef(b2)), length(coef(b2)))
      expect_equal(sum(confint(b2)[, 2] >= coef(b2)), length(coef(b2)))

      c1 <-
        svyfgtdec(
          ~ eqincome,
          design = dbd_eusilc,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE
        )
      c2 <-
        svyby(
          ~ eqincome,
          by = ~ rb090,
          design = dbd_eusilc,
          FUN = svyfgtdec,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE,
          deff = FALSE
        )

      # database svyfgtdec
      expect_equal(coef(a1), coef(c1))
      # expect_equal(rev(coef(a2)), coef(c2)) # inverted results
      expect_equal(SE(a1), SE(c1))
      # expect_equal(rev(SE(a2)), SE(c2)) # inverted results


      # compare subsetted objects to svyby objects
      sub_des <-
        svyfgtdec(
          ~ eqincome,
          design = subset(des_eusilc, rb090 == "male"),
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE
        )
      sby_des <-
        svyby(
          ~ eqincome,
          by = ~ rb090,
          design = des_eusilc,
          FUN = svyfgtdec,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE,
          deff = FALSE
        )
      sub_rep <-
        svyfgtdec(
          ~ eqincome,
          design = subset(des_eusilc_rep, rb090 == "male"),
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE
        )
      sby_rep <-
        svyby(
          ~ eqincome,
          by = ~ rb090,
          design = des_eusilc_rep,
          FUN = svyfgtdec,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE,
          deff = FALSE
        )

      # subsets equal svyby
      expect_equal(as.numeric(coef(sub_des)), as.numeric(coef(sby_des)[grepl("^male", names(coef(sby_des)))]))
      expect_equal(as.numeric(coef(sub_rep)), as.numeric(coef(sby_rep)[grepl("^male", names(coef(sby_rep)))]))
      expect_equal(as.numeric(SE(sub_des)), as.numeric(SE(sby_des)[1,]))
      expect_equal(as.numeric(SE(sub_rep)), as.numeric(SE(sby_rep)[1,]))

      # coefficients should match across svydesign & svrepdesign
      expect_equal(coef(sby_des), coef(sby_rep))

      # coefficients of variation should be within five percent
      cv_dif <- abs(cv(sby_des) - cv(sby_rep))
      expect_lte(max(unlist(cv_dif)) , .05)




      sub_dbd <-
        svyfgtdec(
          ~ eqincome,
          design = subset(dbd_eusilc, rb090 == "male"),
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE
        )
      sby_dbd <-
        svyby(
          ~ eqincome,
          by = ~ rb090,
          design = dbd_eusilc,
          FUN = svyfgtdec,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE,
          deff = FALSE
        )
      sub_dbr <-
        svyfgtdec(
          ~ eqincome,
          design = subset(dbd_eusilc_rep, rb090 == "male"),
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE
        )
      sby_dbr <-
        svyby(
          ~ eqincome,
          by = ~ rb090,
          design = dbd_eusilc_rep,
          FUN = svyfgtdec,
          g = this_g,
          type_thresh = this_thresh,
          percent = .6,
          abs_thresh = 15000 ,
          na.rm  = FALSE,
          deff = FALSE
        )

      # compare database-backed designs to non-database-backed designs
      # dbi subsets equal non-dbi subsets
      expect_equal(coef(sub_des), coef(sub_dbd))
      expect_equal(coef(sub_rep), coef(sub_dbr))
      expect_equal(SE(sub_des), SE(sub_dbd))
      expect_equal(SE(sub_rep), SE(sub_dbr))



      # compare database-backed subsetted objects to database-backed svyby objects
      # dbi subsets equal dbi svyby
      expect_equal(as.numeric(coef(sub_dbd)), as.numeric(coef(sby_dbd[2,]))) # inverted results!
      expect_equal(as.numeric(coef(sub_dbr)), as.numeric(coef(sby_dbr[2,]))) # inverted results!
      expect_equal(as.numeric(SE(sub_dbd)), as.numeric(SE(sby_dbd[2,]))) # inverted results!
      expect_equal(as.numeric(SE(sub_dbr)), as.numeric(SE(sby_dbr[2,]))) # inverted results!


    }
  }

  dbRemoveTable(conn , 'eusilc')
  dbDisconnect(conn)

})