tests/testthat/test-extract_svyjdivdec-pos.R

In convey: Income Concentration Analysis with Complex Survey Samples

skip_on_cran()


library(laeken)
library(survey)

data(api)
dstrat1<-convey_prep(svydesign(id=~1,data=apistrat))
test_that("svyjdivdec works on unweighted designs",{
  svygei(~api00, design=dstrat1)
})

# build designs
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_save <- des_eusilc_rep <-as.svrepdesign(des_eusilc, type= "bootstrap")
des_eusilc_rep <- convey_prep(des_eusilc_rep)
des_eusilc <- subset( des_eusilc , eqincome > 0 )
des_eusilc_rep <- subset( des_eusilc_rep , eqincome > 0 )

# point estimates
a1 <- svyjdivdec(~eqincome, subgroup = ~db040, design = des_eusilc )
a2 <- svyby(formula = ~eqincome, by =~db040 , design = des_eusilc , subgroup = ~rb090, svyjdivdec )
b1 <- svyjdivdec(~eqincome, subgroup = ~db040, design = des_eusilc_rep )
b2 <- svyby(formula = ~eqincome, by =~db040 , design = des_eusilc_rep , subgroup = ~rb090, svyjdivdec )


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

test_that("output svyjdivdec",{
  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( max( confint(a2[1,])[,1] - coef(a2[1,]) ) , 0)
  expect_gte( min( confint(a2[1,])[,2] - coef(a2[1,]) ) , 0)
  expect_lte( max( confint(a2[2,])[,1] - coef(a2[2,]) ) , 0)
  expect_gte( min( confint(a2[2,])[,2] - coef(a2[2,]) ) , 0)
  expect_lte( max( confint(a2[3,])[,1] - coef(a2[3,]) ) , 0)
  expect_gte( min( confint(a2[3,])[,2] - coef(a2[3,]) ) , 0)
  expect_lte( max( confint(a2[4,])[,1] - coef(a2[4,]) ) , 0)
  expect_gte( min( confint(a2[4,])[,2] - coef(a2[4,]) ) , 0)
  expect_lte( max( confint(a2[5,])[,1] - coef(a2[5,]) ) , 0)
  expect_gte( min( confint(a2[5,])[,2] - coef(a2[5,]) ) , 0)
  expect_lte( max( confint(a2[6,])[,1] - coef(a2[6,]) ) , 0)
  expect_gte( min( confint(a2[6,])[,2] - coef(a2[6,]) ) , 0)
  expect_lte( max( confint(a2[7,])[,1] - coef(a2[7,]) ) , 0)
  expect_gte( min( confint(a2[7,])[,2] - coef(a2[7,]) ) , 0)
  expect_lte( max( confint(b2[1,])[,1] - coef(b2[1,]) ) , 0)
  expect_gte( min( confint(b2[1,])[,2] - coef(b2[1,]) ) , 0)
  expect_lte( max( confint(b2[2,])[,1] - coef(b2[2,]) ) , 0)
  expect_gte( min( confint(b2[2,])[,2] - coef(b2[2,]) ) , 0)
  expect_lte( max( confint(b2[3,])[,1] - coef(b2[3,]) ) , 0)
  expect_gte( min( confint(b2[3,])[,2] - coef(b2[3,]) ) , 0)
  expect_lte( max( confint(b2[4,])[,1] - coef(b2[4,]) ) , 0)
  expect_gte( min( confint(b2[4,])[,2] - coef(b2[4,]) ) , 0)
  expect_lte( max( confint(b2[5,])[,1] - coef(b2[5,]) ) , 0)
  expect_gte( min( confint(b2[5,])[,2] - coef(b2[5,]) ) , 0)
  expect_lte( max( confint(b2[6,])[,1] - coef(b2[6,]) ) , 0)
  expect_gte( min( confint(b2[6,])[,2] - coef(b2[6,]) ) , 0)
  expect_lte( max( confint(b2[7,])[,1] - coef(b2[7,]) ) , 0)
  expect_gte( min( confint(b2[7,])[,2] - coef(b2[7,]) ) , 0)
})


# compare subsetted objects to svyby objects
sub_des <- svyjdivdec( ~eqincome , subgroup = ~db040 , design = subset( des_eusilc , rb090 == "male" ) )
sby_des <- svyby( ~eqincome, by = ~rb090, design = des_eusilc, FUN = svyjdivdec , subgroup = ~db040 )
sub_rep <- svyjdivdec( ~eqincome , subgroup = ~db040 , design = subset( des_eusilc_rep , rb090 == "male" ) )
sby_rep <- svyby( ~eqincome, by = ~rb090, design = des_eusilc_rep , FUN = svyjdivdec , subgroup = ~db040 )

test_that("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(coef(sub_des), coef(sub_rep))
  expect_equal(coef(sby_des), coef(sby_rep))

  # coefficients of variation should be within five percent
  cv_dif <- max(unlist(abs(cv(sby_des)-cv(sby_rep))))
  expect_lte( cv_dif, 5 ) # add warning for unreliable cv
})



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

dbd_eusilc <- subset( dbd_eusilc , eqincome > 0 )


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

dbd_eusilc_rep <- convey_prep( dbd_eusilc_rep )

dbd_eusilc_rep <- subset( dbd_eusilc_rep , eqincome > 0 )

c1 <- svyjdivdec( ~ eqincome , subgroup = ~db040 , design = dbd_eusilc )
c2 <- svyby(formula = ~eqincome, by =~db040 , design = dbd_eusilc , subgroup = ~rb090, svyjdivdec )
d1 <- svyjdivdec( ~ eqincome , subgroup = ~db040 , design = dbd_eusilc_rep )
d2 <- svyby(formula = ~eqincome, by =~db040 , design = dbd_eusilc_rep , subgroup = ~rb090, svyjdivdec )

test_that("output database svyjdivdec",{
  expect_is(coef(c1),"numeric")
  expect_is(coef(c2), "numeric")
  expect_is(coef(d1),"numeric")
  expect_is(coef(d2),"numeric")
  expect_equal(coef(c1), coef(d1))
  expect_equal(coef(c2), coef(d2))
  expect_is(SE(c1),"numeric")
  expect_is(SE(c2), "svyby" )
  expect_is(SE(d1),"numeric")
  expect_is(SE(d2),"svyby")
  expect_lte(confint(c1)[1,1], coef(c1)[1])
  expect_gte(confint(c1)[1,2], coef(c1)[1])
  expect_lte(confint(c1)[2,1], coef(c1)[2])
  expect_gte(confint(c1)[2,2], coef(c1)[2])
  expect_lte(confint(c1)[3,1], coef(c1)[3])
  expect_gte(confint(c1)[3,2], coef(c1)[3])
  expect_lte( max( confint(c2[1,])[,1] - coef(c2[1,]) ) , 0)
  expect_gte( min( confint(c2[1,])[,2] - coef(c2[1,]) ) , 0)
  expect_lte( max( confint(c2[2,])[,1] - coef(c2[2,]) ) , 0)
  expect_gte( min( confint(c2[2,])[,2] - coef(c2[2,]) ) , 0)
  expect_lte( max( confint(c2[3,])[,1] - coef(c2[3,]) ) , 0)
  expect_gte( min( confint(c2[3,])[,2] - coef(c2[3,]) ) , 0)
  expect_lte( max( confint(c2[4,])[,1] - coef(c2[4,]) ) , 0)
  expect_gte( min( confint(c2[4,])[,2] - coef(c2[4,]) ) , 0)
  expect_lte( max( confint(c2[5,])[,1] - coef(c2[5,]) ) , 0)
  expect_gte( min( confint(c2[5,])[,2] - coef(c2[5,]) ) , 0)
  expect_lte( max( confint(c2[6,])[,1] - coef(c2[6,]) ) , 0)
  expect_gte( min( confint(c2[6,])[,2] - coef(c2[6,]) ) , 0)
  expect_lte( max( confint(c2[7,])[,1] - coef(c2[7,]) ) , 0)
  expect_gte( min( confint(c2[7,])[,2] - coef(c2[7,]) ) , 0)
  expect_lte( max( confint(d2[1,])[,1] - coef(d2[1,]) ) , 0)
  expect_gte( min( confint(d2[1,])[,2] - coef(d2[1,]) ) , 0)
  expect_lte( max( confint(d2[2,])[,1] - coef(d2[2,]) ) , 0)
  expect_gte( min( confint(d2[2,])[,2] - coef(d2[2,]) ) , 0)
  expect_lte( max( confint(d2[3,])[,1] - coef(d2[3,]) ) , 0)
  expect_gte( min( confint(d2[3,])[,2] - coef(d2[3,]) ) , 0)
  expect_lte( max( confint(d2[4,])[,1] - coef(d2[4,]) ) , 0)
  expect_gte( min( confint(d2[4,])[,2] - coef(d2[4,]) ) , 0)
  expect_lte( max( confint(d2[5,])[,1] - coef(d2[5,]) ) , 0)
  expect_gte( min( confint(d2[5,])[,2] - coef(d2[5,]) ) , 0)
  expect_lte( max( confint(d2[6,])[,1] - coef(d2[6,]) ) , 0)
  expect_gte( min( confint(d2[6,])[,2] - coef(d2[6,]) ) , 0)
  expect_lte( max( confint(d2[7,])[,1] - coef(d2[7,]) ) , 0)
  expect_gte( min( confint(d2[7,])[,2] - coef(d2[7,]) ) , 0)
})

test_that("check dbi and non-dbi svyjdivdec output",{
  expect_equal(coef(c1), coef(a1))
  expect_equal(coef(c2), coef(a2))
  expect_equal(coef(d1), coef(b1))
  expect_equal(coef(d2), coef(b2))
  expect_equal(SE(c1), SE(a1))
  expect_equal(SE(c2), SE(a2))
  expect_equal(SE(d1), SE(b1))
  expect_equal(SE(d2), SE(b2))
  expect_equal(cv(c1), cv(a1))
  expect_equal(cv(c2), cv(a2))
  expect_equal(cv(d1), cv(b1))
  expect_equal(cv(d2), cv(b2))
  expect_equal(confint(c1), confint(a1))
  expect_equal(confint(c2), confint(a2))
  expect_equal(confint(d1), confint(d1))
  expect_equal(confint(d2), confint(d2))
})

sub_dbd <- svyjdivdec( ~eqincome , subgroup = ~db040, design = subset( dbd_eusilc , rb090 == "male" ) )
sub_dbr <- svyjdivdec( ~eqincome , subgroup = ~db040, design = subset( dbd_eusilc_rep , rb090 == "male" ) )
sby_dbd <- svyby(formula = ~eqincome, by =~rb090 , design = dbd_eusilc , subgroup = ~db040, svyjdivdec )
sby_dbr <- svyby(formula = ~eqincome, by =~rb090 , design = dbd_eusilc_rep , subgroup = ~db040, svyjdivdec )

# compare database-backed designs to non-database-backed designs
test_that("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))
  expect_equal(confint(sub_des), confint(sub_dbd))
  expect_equal(confint(sub_rep), confint(sub_dbr))
})

# compare database-backed subsetted objects to database-backed svyby objects
test_that("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 )