|
| 1 | +# Test if vardom and vardomh return the same result with the same information |
| 2 | +# where: |
| 3 | +# - vardom use the data at level 1 (for example households) |
| 4 | +# - vardomh use the data at level 2 (for example persons) |
| 5 | +# |
| 6 | +# test_fun - test for data without period |
| 7 | +# test_fun2 - test for data with period |
| 8 | + |
| 9 | +# library(testthat) |
| 10 | +# devtools::load_all(path = "vardpoor") |
| 11 | +# library(data.table) |
| 12 | +# library(vardpoor) |
| 13 | + |
| 14 | +#### test without period ---- |
| 15 | + |
| 16 | +test_fun <- function(n1 = 3000) { |
| 17 | + # data generation |
| 18 | + # g-matrix |
| 19 | + g_dat <- data.table(id_m = 1:n1, |
| 20 | + x0 = 1L, |
| 21 | + x1 = sample(0:2, n1, replace = T), |
| 22 | + x2 = sample(0:3, n1, replace = T), |
| 23 | + x3 = sample(0:3, n1, replace = T), |
| 24 | + x4 = sample(0:3, n1, replace = T), |
| 25 | + x5 = sample(0:3, n1, replace = T), |
| 26 | + g = rnorm(n1, 1, 0.1), |
| 27 | + q = runif(n1)) |
| 28 | + |
| 29 | + # vardom data |
| 30 | + nn <- round(n1 / 300, 0) # number of PSUs |
| 31 | + # dat_x <- data.table() |
| 32 | + strata <- sample(1:4, n1, replace = T,) |
| 33 | + dat_x <- data.table( |
| 34 | + id_m = 1:n1, |
| 35 | + b06 = sample(1:11, n1, replace = T, |
| 36 | + prob = c(0.1922021, 0.1302067, 0.0641481, 0.0471424, |
| 37 | + 0.01797438, 0.007695619, 0.002179529, 0.00131848, |
| 38 | + 0.0003498009, 0.000107631, 2.690776e-05)), |
| 39 | + wd = sample(1:5, n1, replace = T), |
| 40 | + strata = strata, |
| 41 | + survey = rep(1:4, each = n1/4), |
| 42 | + iec = sample(1:nn, n1, replace = T) |
| 43 | + ) |
| 44 | + |
| 45 | + # vardomh data |
| 46 | + dat_y <- data.table() |
| 47 | + l <- 1 |
| 48 | + for (i in 1:n1) { |
| 49 | + if (dat_x[i, b06] == 1) { |
| 50 | + rinda <- dat_x[i] |
| 51 | + rinda <- rinda[, id_p := l] |
| 52 | + dat_y <- rbind(dat_y, rinda) |
| 53 | + l <- l + 1 |
| 54 | + } else { |
| 55 | + j <- dat_x[i, b06] |
| 56 | + for (k in 1:j) { |
| 57 | + rinda <- dat_x[i] |
| 58 | + rinda <- rinda[, b06 := k] |
| 59 | + rinda <- rinda[, id_p := l] |
| 60 | + dat_y <- rbind(dat_y, rinda) |
| 61 | + l <- l + 1 |
| 62 | + } |
| 63 | + } |
| 64 | + } |
| 65 | + rm(i, j, k, l, rinda) |
| 66 | + |
| 67 | + n2 <- nrow(dat_y) |
| 68 | + if ((nrow(dat_y) - sum(dat_x$b06)) != 0) stop("wrong number of records") |
| 69 | + |
| 70 | + dat_y[, y := sample(0:1, n2, replace = T)] # generate y variable |
| 71 | + gg <- g_dat[, c("id_m", "g")] |
| 72 | + |
| 73 | + ycol <- dat_y[, lapply(.SD, sum), |
| 74 | + keyby = id_m, |
| 75 | + .SDcols = "y"] |
| 76 | + dat_yy <- merge(dat_y, gg, by = "id_m") |
| 77 | + dat_x$y <- ycol$y # add to vardom data |
| 78 | + dat_xx <- merge(dat_x, g_dat, by = c("id_m")) |
| 79 | + dat_xx[, wc := wd * g] |
| 80 | + dat_yy[, wc := wd * g] |
| 81 | + rm(dat_x, dat_y, gg, ycol) |
| 82 | + |
| 83 | + # calculations ---- |
| 84 | + |
| 85 | + n_h <- data.table(strata = 1:4, pop = sample(500:1500, 4)) |
| 86 | + res_1 <- vardom(Y = "y", |
| 87 | + H = "strata", |
| 88 | + PSU = "iec", |
| 89 | + w_final = "wc", |
| 90 | + fh_zero = TRUE, |
| 91 | + N_h = n_h, |
| 92 | + X = paste0("x", 0:5), |
| 93 | + g = "g", |
| 94 | + q = "q", |
| 95 | + dataset = dat_xx)$all_result |
| 96 | + |
| 97 | + res_2 <- vardomh(Y = "y", |
| 98 | + H = "strata", |
| 99 | + PSU = "iec", |
| 100 | + w_final = "wc", |
| 101 | + ID_level1 = "id_m", |
| 102 | + ID_level2 = "id_p", |
| 103 | + N_h = n_h, |
| 104 | + fh_zero = TRUE, |
| 105 | + X = paste0("x", 0:5), |
| 106 | + X_ID_level1 = "id_m", |
| 107 | + g = "g", |
| 108 | + q = "q", |
| 109 | + dataset = dat_yy, |
| 110 | + datasetX = g_dat)$all_result |
| 111 | + |
| 112 | + # Those variables differ if calculated from different levels |
| 113 | + res_1 <- res_1[, -c("respondent_count", "n_nonzero","pop_size")] |
| 114 | + res_2 <- res_2[, -c("respondent_count", "n_nonzero","pop_size")] |
| 115 | + |
| 116 | + names_a <- names(res_2) |
| 117 | + return(list(res_2, res_1[, .SD, .SDcols = names_a])) |
| 118 | +} |
| 119 | + |
| 120 | +test_that("test equal for one period, diferent levels", { |
| 121 | + results <- test_fun() |
| 122 | + expect_equal(results[[1]], results[[2]]) |
| 123 | +}) |
| 124 | + |
| 125 | + |
| 126 | +#### 2 periods, with calibration, without Z ---- |
| 127 | + |
| 128 | +test_fun2 <- function(n1){ |
| 129 | + # g-matrix |
| 130 | + g_dat <- data.table(id_m = 1:n1, |
| 131 | + x0 = 1L, |
| 132 | + x1 = sample(0:2, n1, replace = T), |
| 133 | + x2 = sample(0:3, n1, replace = T), |
| 134 | + x3 = sample(0:3, n1, replace = T), |
| 135 | + x4 = sample(0:3, n1, replace = T), |
| 136 | + x5 = sample(0:3, n1, replace = T), |
| 137 | + g = rnorm(n1, 1, 0.1), |
| 138 | + q = runif(n1), |
| 139 | + period = rep(1:2, each = n1 / 2)) |
| 140 | + |
| 141 | + # vardom data |
| 142 | + nn <- round(n1 / 300, 0) # number of PSUs |
| 143 | + # dat_x <- data.table() |
| 144 | + strata <- sample(1:4, n1, replace = T,) |
| 145 | + dat_x <- data.table( |
| 146 | + id_m = 1:n1, |
| 147 | + b06 = sample(1:11, n1, replace = T, |
| 148 | + prob = c(0.1922021, 0.1302067, 0.0641481, 0.0471424, |
| 149 | + 0.01797438, 0.007695619, 0.002179529, 0.00131848, |
| 150 | + 0.0003498009, 0.000107631, 2.690776e-05)), |
| 151 | + wd = sample(1:5, n1, replace = T), |
| 152 | + strata = strata, |
| 153 | + survey = rep(1:4, each = n1/4), |
| 154 | + iec = sample(1:nn, n1, replace = T) |
| 155 | + ) |
| 156 | + |
| 157 | + # vardomh data |
| 158 | + dat_y <- data.table() |
| 159 | + l <- 1 |
| 160 | + for (i in 1:n1) { |
| 161 | + if (dat_x[i, b06] == 1) { |
| 162 | + rinda <- dat_x[i] |
| 163 | + rinda <- rinda[, id_p := l] |
| 164 | + dat_y <- rbind(dat_y, rinda) |
| 165 | + l <- l + 1 |
| 166 | + } else { |
| 167 | + j <- dat_x[i, b06] |
| 168 | + for (k in 1:j) { |
| 169 | + rinda <- dat_x[i] |
| 170 | + rinda <- rinda[, b06 := k] |
| 171 | + rinda <- rinda[, id_p := l] |
| 172 | + dat_y <- rbind(dat_y, rinda) |
| 173 | + l <- l + 1 |
| 174 | + } |
| 175 | + } |
| 176 | + } |
| 177 | + rm(i, j, k, l, rinda) |
| 178 | + |
| 179 | + n2 <- nrow(dat_y) |
| 180 | + if ((nrow(dat_y) - sum(dat_x$b06)) != 0) stop("wrong number of records") |
| 181 | + |
| 182 | + dat_y[, y := sample(0:1, n2, replace = T)] # add y variable |
| 183 | + gg <- g_dat[, c("id_m", "g", "period")] |
| 184 | + |
| 185 | + ycol <- dat_y[, lapply(.SD, sum), |
| 186 | + keyby = id_m, |
| 187 | + .SDcols = "y"] |
| 188 | + dat_yy <- merge(dat_y, gg, by = "id_m") |
| 189 | + dat_x$y <- ycol$y # add to vardom data |
| 190 | + dat_xx <- merge(dat_x, g_dat, by = c("id_m")) |
| 191 | + dat_xx[, wc := wd * g] |
| 192 | + dat_yy[, wc := wd * g] |
| 193 | + |
| 194 | + n_h <- data.table(period = rep(1:2, each = 4), |
| 195 | + strata = rep(1:4, 2), |
| 196 | + pop = rep(sample(500:1500, 4),2)) |
| 197 | + |
| 198 | + res_1 <- vardom(Y = "y", |
| 199 | + H = "strata", |
| 200 | + PSU = "iec", |
| 201 | + w_final = "wc", |
| 202 | + fh_zero = TRUE, |
| 203 | + N_h = n_h, |
| 204 | + X = paste0("x", 0:5), |
| 205 | + g = "g", |
| 206 | + q = "q", |
| 207 | + dataset = dat_xx, |
| 208 | + period = "period")$all_result |
| 209 | + |
| 210 | + res_2 <- vardomh(Y = "y", |
| 211 | + H = "strata", |
| 212 | + PSU = "iec", |
| 213 | + w_final = "wc", |
| 214 | + ID_level1 = "id_m", |
| 215 | + ID_level2 = "id_p", |
| 216 | + N_h = n_h, |
| 217 | + fh_zero = TRUE, |
| 218 | + X = paste0("x", 0:5), |
| 219 | + X_ID_level1 = "id_m", |
| 220 | + g = "g", |
| 221 | + q = "q", |
| 222 | + dataset = dat_yy, |
| 223 | + datasetX = g_dat, |
| 224 | + period = "period", |
| 225 | + periodX = "period")$all_result |
| 226 | + |
| 227 | + # Those variables differ if calculated from different levels |
| 228 | + res_1 <- res_1[, -c("respondent_count", "n_nonzero","pop_size")] |
| 229 | + res_2 <- res_2[, -c("respondent_count", "n_nonzero","pop_size")] |
| 230 | + |
| 231 | + names_a <- names(res_2) |
| 232 | + return(list(res_2, res_1[, .SD, .SDcols = names_a])) |
| 233 | +} |
| 234 | + |
| 235 | +test_that("test equal for one period, diferent levels", { |
| 236 | + results <- test_fun2(6000) |
| 237 | + expect_equal(results[[1]], results[[2]]) |
| 238 | +}) |
0 commit comments