tests/testthat/test-skewed_distributions.R

In DjalmaPessoa/convey: Income Concentration Analysis with Complex Survey Samples

library(survey)


test_that("functions work on weird distributions" , {
  skip_on_cran()

  no.na <-
    function(z , value = FALSE) {
      z[is.na(z)] <- value
      z
    }

  all_funs <-
    c(
		"svyarpr" ,
      "svyrmir" ,
      "svyqsr" ,
      "svyarpt" ,
      "svyarpt" ,
      "svyatk" ,
      "svyzenga" ,
      "svyfgt" ,
      "svygini" ,
      "svygpg" ,
      "svyiqalpha" ,
      "svyisq" ,
      "svypoormed" ,
      "svygei" ,
      "svyrmpg" ,
      "svylorenz" ,
      "svyjdiv" ,
      "svyrich" ,
      "svywatts"
    )
	
  for (n in c(50 , 1000)) {
    set.seed(n)

    dist_frame <-
      data.frame(
        wt_binom_two = rbinom(n , 1 , .5) ,
        wt_binom = rbinom(n , 1 , .5) ,
        wt_unif = runif(n) ,
        beta_050_050 = rbeta(n , 0.5 , 0.5) ,
        beta_500_100 = rbeta(n , 5 , 1) ,
        beta_100_300 = rbeta(n , 1 , 3) ,
        beta_200_200 = rbeta(n , 2 , 2) ,
        beta_200_500 = rbeta(n , 2 , 5) ,
        pois1 = rpois(n , 1) ,
        pois100 = rpois(n , 100) ,
        sex = ifelse(rbinom(n , 1 , 0.5) , "male" , "female") ,
        age = sample(0:99 , n , replace = TRUE)
      )

    dist_frame[-which(names(dist_frame) == 'sex')] <-
      sapply(dist_frame[-which(names(dist_frame) == 'sex')] , as.numeric)


    for (this_function in all_funs) {
      FUN <- get(this_function)

      unwtd_des <-
        convey_prep(suppressWarnings(svydesign(~ 1 , data = dist_frame)))
      binom_des <-
        convey_prep(svydesign(~ 1 , data = dist_frame , weight = ~ wt_binom))
      unif_des <-
        convey_prep(svydesign(~ 1 , data = dist_frame , weight = ~ wt_unif))

      unwtd_rep <-
        convey_prep(as.svrepdesign(unwtd_des , type = 'bootstrap'))
      binom_rep <-
        convey_prep(as.svrepdesign(binom_des , type = 'bootstrap'))
      unif_rep <-
        convey_prep(as.svrepdesign(unif_des , type = 'bootstrap'))


      for (this_prefix in c('unwtd' , 'binom' , 'unif')) {
        lin_des <- get(paste0(this_prefix , "_des"))
        rep_des <- get(paste0(this_prefix , "_rep"))

        for (this_formula in list(
          ~ beta_050_050 ,
          ~ beta_500_100 ,
          ~ beta_100_300 ,
          ~ beta_200_200 ,
          ~ beta_200_500 ,
          ~ pois1 ,
          ~ pois100
        )) {
          # function-specific parameters here

          if (any(unlist(lapply(list(svygei , svyatk , svyjdiv) , function(z)
            identical(z , FUN))))) {
            lin_des <-
              subset(lin_des , no.na(lin_des$variables[, as.character(this_formula)[2]] > 0))
            rep_des <-
              subset(rep_des , no.na(rep_des$variables[, as.character(this_formula)[2]] > 0))

          }

          lin_params_list <- list(this_formula , lin_des)
          rep_params_list <- list(this_formula , rep_des)


          if (identical(FUN , svyfgt)) {
            lin_params_list <-
              c(lin_params_list ,
                list(
                  g = 0,
                  type_thresh = "abs",
                  abs_thresh = 10000
                ))
            rep_params_list <-
              c(rep_params_list ,
                list(
                  g = 0,
                  type_thresh = "abs",
                  abs_thresh = 10000
                ))

          }

          if (identical(FUN , svygpg)) {
            lin_params_list <- c(lin_params_list , list(sex = ~ sex))
            rep_params_list <-
              c(rep_params_list , list(sex = ~ sex))

          }

          if (identical(FUN , svyiqalpha) |
              identical(FUN , svyisq)) {
            lin_params_list <- c(lin_params_list , list(alpha = 0.5))
            rep_params_list <-
              c(rep_params_list , list(alpha = 0.5))

          }

          if (identical(FUN , svyrmir)) {
            lin_params_list <- c(lin_params_list , list(age = ~ age))
            rep_params_list <-
              c(rep_params_list , list(age = ~ age))

          }

          if (identical(FUN , svyrich)) {
            lin_params_list <-
              c(
                lin_params_list ,
                list(
                  type_measure = "FGTT1",
                  g = 1,
                  type_thresh = "abs",
                  abs_thresh = 30000
                )
              )
            rep_params_list <-
              c(
                rep_params_list ,
                list(
                  type_measure = "FGTT1",
                  g = 1,
                  type_thresh = "abs",
                  abs_thresh = 30000
                )
              )

          }


          if (identical(FUN , svywatts)) {
            lin_params_list <-
              c(lin_params_list ,
                list(type_thresh = "abs", abs_thresh = 10000))
            rep_params_list <-
              c(rep_params_list ,
                list(type_thresh = "abs", abs_thresh = 10000))

          }

          if (this_prefix == 'unwtd')
            wt_vec <- seq(nrow(dist_frame))
          if (this_prefix == 'binom')
            wt_vec <- which(dist_frame$wt_binom > 0)
          if (this_prefix == 'unif')
            wt_vec <- which(dist_frame$wt_unif > 0)



          if ((identical(svyrmpg , FUN) |
               identical(svypoormed , FUN)) &
              ((identical(this_formula , ~ pois1)) |
               (coef(
                 svyarpt(this_formula , lin_des)
               ) < min(dist_frame[wt_vec , as.character(this_formula)[2]])))) {
            expect_equal(as.numeric(coef(
              do.call(FUN , lin_params_list)
            )), NA_real_)
          } else if (((identical(svyrmpg , FUN)) &
                      coef(svyarpt(this_formula , lin_des)) < min(dist_frame[wt_vec , as.character(this_formula)[2]])) |
                     ((identical(svywatts , FUN) | identical(svyqsr , FUN)) &
                      identical(this_formula , ~ pois1))) {
            expect_error(do.call(FUN , lin_params_list))

          } else {
            lin_res <- do.call(FUN , lin_params_list)
            rep_res <- do.call(FUN , rep_params_list)

            print(
              paste(
                "testing functions work on weird distributions" ,
                this_function ,
                this_prefix ,
                as.character(this_formula)[2]
              )
            )

            # result should give some non-missing number
            expect_that(coef(lin_res) , is.numeric)
            expect_that(coef(lin_res) , function(w)
              all(!is.na(w)))

            # coefficients should be equal
            expect_equal(coef(lin_res) , coef(rep_res))

            # difference between SEs for the designs should be less than 25% of the coef
            # but only for larger designs, since small ones can be unstable
            if (n > 50 &
                all(coef(lin_res) * .5 > SE(lin_res))) {
              expect_true(all(abs(SE(
                lin_res
              ) - SE(
                rep_res
              )) <= max(coef(
                lin_res
              )) * 0.3))
            }



          }

        }
      }

    }
  }


})