tests/testthat/test-extract_svyrich.R

In convey: Income Concentration Analysis with Complex Survey Samples

context("svyrich output survey.design and svyrep.design")

skip_on_cran()

library(laeken)
library(survey)


data(api)
dstrat1<-convey_prep(svydesign(id=~1,data=apistrat))
for ( this_measure in c( "Cha" , "FGTT1" , "FGTT2" ) ){
  for ( this_thresh in c( "abs" , "relm" , "relq" ) ){
    for ( this_g in c( 0 , 1 ) ) {

      if (this_measure == "Cha" & this_g == 0 ) this_g = 1/2
      if (this_measure == "FGTT2" ) this_g = this_g + 2

      test_that( paste0(this_measure , "-" , this_g ,"-svyrich works on unweighted designs"),{
        svyrich(~api00, design=dstrat1, type_measure = this_measure , g=this_g, type_thresh=this_thresh, abs_thresh=500)
      })
    }
  }
}



test_that( "output svyrich" , {
  skip_on_cran()

  data(eusilc) ; names( eusilc ) <- tolower( names( eusilc ) )

  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" )

  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_measure in c( "Cha" , "FGTT1" , "FGTT2" ) ){
    for ( this_thresh in c( "abs" , "relm" , "relq" ) ){
      for ( this_g in c( 0 , 1 ) ) {

        if (this_measure == "Cha" & this_g == 0 ) this_g = 1/2
        if (this_measure == "FGTT2" ) this_g = this_g + 2

        a1 <- svyrich(~eqincome, design = des_eusilc, type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000 )

        a2 <- svyby(~eqincome, by = ~hsize, design = des_eusilc, FUN = svyrich, type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000, deff = FALSE)

        b1 <- svyrich(~eqincome, design = des_eusilc_rep, type_measure = this_measure , g=this_g, type_thresh= this_thresh, abs_thresh=10000)

        b2 <- svyby(~eqincome, by = ~hsize, design = des_eusilc_rep, FUN = svyrich, type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000, deff = FALSE)


        se_dif1 <- 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_lte(se_dif1, coef(a1) * 0.05 ) # the difference between CVs should be less than 5% of the coefficient, otherwise manually set it
        expect_lte(se_diff2, max( coef(a2) ) * 0.1 ) # the difference between CVs should be less than 10% of the maximum coefficient, otherwise manually set it
        expect_is(SE(a1),"matrix")
        expect_is(SE(a2), "numeric")
        expect_is(SE(b1),"numeric")
        expect_is(SE(b2),"numeric")
        expect_lte(confint(a1)[1], coef(a1))
        expect_gte(confint(a1)[2],coef(a1))
        expect_lte(confint(b1)[,1], coef(b1))
        expect_gte(confint(b1)[2], coef(b1))
        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 <- svyrich(~eqincome, design = dbd_eusilc, type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000)
        c2 <- svyby(~eqincome, by = ~hsize, design = dbd_eusilc, FUN = svyrich, type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000, deff = FALSE)


        # database svyrich
        expect_equal(coef(a1), coef(c1))
        expect_equal(coef(a2), coef(c2))
        expect_equal(SE(a1), SE(c1))
        expect_equal(SE(a2), SE(c2))


        # compare subsetted objects to svyby objects
        sub_des <- svyrich( ~eqincome , design = subset( des_eusilc , hsize == 1), type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000 )
        sby_des <- svyby( ~eqincome, by = ~hsize, design = des_eusilc, FUN = svyrich , type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000)
        sub_rep <- svyrich( ~eqincome , design = subset( des_eusilc_rep , hsize == 1), type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000 )
        sby_rep <- svyby( ~eqincome, by = ~hsize, design = des_eusilc_rep, FUN = svyrich, type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000)

        # subsets equal svyby
        expect_equal(as.numeric(coef(sub_des)), as.numeric(coef(sby_des))[1])
        expect_equal(as.numeric(coef(sub_rep)), as.numeric(coef(sby_rep))[1])
        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(as.numeric(coef(sub_des)), as.numeric(coef(sby_rep))[1])

        # coefficients of variation should be within five percent
        cv_dif <- abs(cv(sub_des)-cv(sby_rep)[1])

        expect_lte(cv_dif,5)


        sub_dbd <- svyrich( ~eqincome , design = subset( dbd_eusilc , hsize == 1), type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000 )
        sby_dbd <- svyby( ~eqincome, by = ~hsize, design = dbd_eusilc, FUN = svyrich , type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000)
        sub_dbr <- svyrich( ~eqincome , design = subset( dbd_eusilc_rep , hsize == 1), type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000 )
        sby_dbr <- svyby( ~eqincome, by = ~hsize, design = dbd_eusilc_rep, FUN = svyrich , type_measure = this_measure, g=this_g, type_thresh= this_thresh, abs_thresh=10000)



        # 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))[1])
        expect_equal(as.numeric(coef(sub_dbr)), as.numeric(coef(sby_dbr))[1])
        expect_equal(as.numeric(SE(sub_dbd)), as.numeric(SE(sby_dbd))[1])
        expect_equal(as.numeric(SE(sub_dbr)), as.numeric(SE(sby_dbr))[1])



      }
    }
  }

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

})