In this vignette we will explain how some functions of the package
are used in order to estimate a contingency table. We will work on the
eusilc dataset of the laeken package. All the
functions presented in the following are explained in the proposed
manuscript by Raphaël Jauslin and Yves Tillé (2021) doi:10.1016/j.jspi.2022.12.003.
We will estimate the contingency table when the factor variable which
represents the economic status pl030 is crossed with a
discretized version of the equivalized household income
eqIncome. In order to discretize the equivalized income, we
calculate percentiles (0.15,0.30,0.45,0.60,0.75,0.90) of the variable
and define the category as intervals between the values.
library(laeken)
library(sampling)
library(StratifiedSampling)
#> Loading required package: Matrix
data("eusilc")
eusilc <- na.omit(eusilc)
N <- nrow(eusilc)
# Xm are the matching variables and id are identity of the units
Xm <- eusilc[,c("hsize","db040","age","rb090","pb220a")]
Xmcat <- do.call(cbind,apply(Xm[,c(2,4,5)],MARGIN = 2,FUN = disjunctive))
Xm <- cbind(Xmcat,Xm[,-c(2,4,5)])
id <- eusilc$rb030
# categorial income splitted by the percentile
c_income <- eusilc$eqIncome
q <- quantile(eusilc$eqIncome, probs = seq(0, 1, 0.15))
c_income[which(eusilc$eqIncome <= q[2])] <- "(0,15]"
c_income[which(q[2] < eusilc$eqIncome & eusilc$eqIncome <= q[3])] <- "(15,30]"
c_income[which(q[3] < eusilc$eqIncome & eusilc$eqIncome <= q[4])] <- "(30,45]"
c_income[which(q[4] < eusilc$eqIncome & eusilc$eqIncome <= q[5])] <- "(45,60]"
c_income[which(q[5] < eusilc$eqIncome & eusilc$eqIncome <= q[6])] <- "(60,75]"
c_income[which(q[6] < eusilc$eqIncome & eusilc$eqIncome <= q[7])] <- "(75,90]"
c_income[which( eusilc$eqIncome > q[7] )] <- "(90,100]"
# variable of interests
Y <- data.frame(ecostat = eusilc$pl030)
Z <- data.frame(c_income = c_income)
# put same rownames
rownames(Xm) <- rownames(Y) <- rownames(Z)<- id
YZ <- table(cbind(Y,Z))
addmargins(YZ)
#> c_income
#> ecostat (0,15] (15,30] (30,45] (45,60] (60,75] (75,90] (90,100] Sum
#> 1 409 616 722 807 935 1025 648 5162
#> 2 189 181 205 184 165 154 82 1160
#> 3 137 90 72 75 59 52 33 518
#> 4 210 159 103 95 74 49 46 736
#> 5 470 462 492 477 459 435 351 3146
#> 6 57 25 28 30 17 11 10 178
#> 7 344 283 194 149 106 91 40 1207
#> Sum 1816 1816 1816 1817 1815 1817 1210 12107Here we set up the sampling designs and define all the quantities we will need for the rest of the vignette. The sample are selected with simple random sampling without replacement and the weights are equal to the inverse of the inclusion probabilities.
# size of sample
n1 <- 1000
n2 <- 500
# samples
s1 <- srswor(n1,N)
s2 <- srswor(n2,N)
# extract matching units
X1 <- Xm[s1 == 1,]
X2 <- Xm[s2 == 1,]
# extract variable of interest
Y1 <- data.frame(Y[s1 == 1,])
colnames(Y1) <- colnames(Y)
Z2 <- as.data.frame(Z[s2 == 1,])
colnames(Z2) <- colnames(Z)
# extract correct identities
id1 <- id[s1 == 1]
id2 <- id[s2 == 1]
# put correct rownames
rownames(Y1) <- id1
rownames(Z2) <- id2
# here weights are inverse of inclusion probabilities
d1 <- rep(N/n1,n1)
d2 <- rep(N/n2,n2)
# disjunctive form
Y_dis <- sampling::disjunctive(as.matrix(Y))
Z_dis <- sampling::disjunctive(as.matrix(Z))
Y1_dis <- Y_dis[s1 ==1,]
Z2_dis <- Z_dis[s2 ==1,]Then the harmonization step must be performed. The
harmonize function returns the harmonized weights. If by
chance the true population totals are known, it is possible to use these
instead of the estimate made within the function.
re <- harmonize(X1,d1,id1,X2,d2,id2)
# if we want to use the population totals to harmonize we can use
re <- harmonize(X1,d1,id1,X2,d2,id2,totals = c(N,colSums(Xm)))
w1 <- re$w1
w2 <- re$w2
colSums(Xm)
#> 1 2 3 4 5 6 7 8 9 10 11
#> 476 887 2340 763 1880 1021 2244 1938 558 6263 5844
#> 12 13 14 hsize age
#> 11073 283 751 36380 559915
colSums(w1*X1)
#> 1 2 3 4 5 6 7 8 9 10 11
#> 476 887 2340 763 1880 1021 2244 1938 558 6263 5844
#> 12 13 14 hsize age
#> 11073 283 751 36380 559915
colSums(w2*X2)
#> 1 2 3 4 5 6 7 8 9 10 11
#> 476 887 2340 763 1880 1021 2244 1938 558 6263 5844
#> 12 13 14 hsize age
#> 11073 283 751 36380 559915The statistical matching is done by using the otmatch
function. The estimation of the contingency table is calculated by
extracting the id1 units (respectively id2
units) and by using the function tapply with the correct
weights.
# Optimal transport matching
object <- otmatch(X1,id1,X2,id2,w1,w2)
head(object[,1:3])
#> id1 id2 weight
#> 101 101 45701 3.656193
#> 101.1 101 441101 8.373407
#> 202 202 375502 9.887632
#> 202.1 202 379001 3.185493
#> 401 401 284805 14.184304
#> 902 902 95801 2.964985
Y1_ot <- cbind(X1[as.character(object$id1),],y = Y1[as.character(object$id1),])
Z2_ot <- cbind(X2[as.character(object$id2),],z = Z2[as.character(object$id2),])
YZ_ot <- tapply(object$weight,list(Y1_ot$y,Z2_ot$z),sum)
# transform NA into 0
YZ_ot[is.na(YZ_ot)] <- 0
# result
round(addmargins(YZ_ot),3)
#> (0,15] (15,30] (30,45] (45,60] (60,75] (75,90] (90,100] Sum
#> 1 619.411 722.561 640.446 866.536 966.009 981.995 440.512 5237.469
#> 2 76.517 137.649 144.173 187.683 196.843 176.594 113.992 1033.452
#> 3 123.355 87.854 29.249 73.792 89.980 62.119 33.451 499.801
#> 4 195.170 67.552 128.444 41.511 104.097 110.730 66.849 714.353
#> 5 444.303 471.796 357.499 569.230 375.301 476.928 502.776 3197.832
#> 6 15.064 0.000 12.864 34.636 37.165 9.846 23.208 132.784
#> 7 214.936 158.993 150.540 162.010 259.129 248.027 97.672 1291.309
#> Sum 1688.755 1646.406 1463.216 1935.399 2028.524 2066.240 1278.461 12107.000As you can see from the previous section, the optimal transport
results generally do not have a one-to-one match, meaning that for every
unit in sample 1, we have more than one unit with weights not equal to 0
in sample 2. The bsmatch function creates a one-to-one
match by selecting a balanced stratified sampling to obtain a data.frame
where each unit in sample 1 has only one imputed unit from sample 2.
# Balanced Sampling
BS <- bsmatch(object)
head(BS$object[,1:3])
#> id1 id2 weight
#> 101 101 45701 3.656193
#> 202 202 375502 9.887632
#> 401 401 284805 14.184304
#> 902.1 902 322401 9.563587
#> 1 1102 1102 13.648046
#> 1201 1201 565201 8.662422
Y1_bs <- cbind(X1[as.character(BS$object$id1),],y = Y1[as.character(BS$object$id1),])
Z2_bs <- cbind(X2[as.character(BS$object$id2),],z = Z2[as.character(BS$object$id2),])
YZ_bs <- tapply(BS$object$weight/BS$q,list(Y1_bs$y,Z2_bs$z),sum)
YZ_bs[is.na(YZ_bs)] <- 0
round(addmargins(YZ_bs),3)
#> (0,15] (15,30] (30,45] (45,60] (60,75] (75,90] (90,100] Sum
#> 1 583.257 737.127 633.194 851.819 1085.673 926.664 419.736 5237.469
#> 2 77.547 160.997 124.874 224.505 190.605 149.724 105.199 1033.452
#> 3 120.232 81.493 26.789 71.281 98.796 63.447 37.763 499.801
#> 4 184.501 58.070 163.903 47.595 101.280 92.179 66.826 714.353
#> 5 435.544 513.348 357.722 530.307 357.484 456.425 547.002 3197.832
#> 6 15.064 0.000 12.864 34.636 35.340 11.671 23.208 132.784
#> 7 199.983 201.288 154.913 134.780 270.865 243.615 85.865 1291.309
#> Sum 1616.127 1752.324 1474.259 1894.923 2140.043 1943.725 1285.599 12107.000
# With Z2 as auxiliary information for stratified balanced sampling.
BS <- bsmatch(object,Z2)
Y1_bs <- cbind(X1[as.character(BS$object$id1),],y = Y1[as.character(BS$object$id1),])
Z2_bs <- cbind(X2[as.character(BS$object$id2),],z = Z2[as.character(BS$object$id2),])
YZ_bs <- tapply(BS$object$weight/BS$q,list(Y1_bs$y,Z2_bs$z),sum)
YZ_bs[is.na(YZ_bs)] <- 0
round(addmargins(YZ_bs),3)
#> (0,15] (15,30] (30,45] (45,60] (60,75] (75,90] (90,100] Sum
#> 1 633.498 725.401 656.807 821.013 976.200 1018.169 406.382 5237.469
#> 2 85.795 149.993 113.627 167.989 204.812 183.719 127.515 1033.452
#> 3 96.315 81.493 26.789 84.450 98.796 74.194 37.763 499.801
#> 4 197.660 65.862 131.390 35.677 115.094 114.100 54.571 714.353
#> 5 435.053 473.904 353.856 588.683 345.965 445.646 554.725 3197.832
#> 6 15.064 0.000 12.864 34.636 35.340 11.671 23.208 132.784
#> 7 224.715 153.933 166.966 195.771 247.355 219.853 82.716 1291.309
#> Sum 1688.099 1650.587 1462.299 1928.219 2023.562 2067.353 1286.880 12107.000
# split the weight by id1
q_l <- split(object$weight,f = object$id1)
# normalize in each id1
q_l <- lapply(q_l, function(x){x/sum(x)})
q <- as.numeric(do.call(c,q_l))
Z_pred <- t(q*disjunctive(object$id1))%*%disjunctive(Z2[as.character(object$id2),])
colnames(Z_pred) <- levels(factor(Z2$c_income))
head(Z_pred)
#> (0,15] (15,30] (30,45] (45,60] (60,75] (75,90] (90,100]
#> [1,] 0.303933 0.0000000 0.6960670 0 0.0000000 0.0000000 0
#> [2,] 0.000000 0.2436673 0.7563327 0 0.0000000 0.0000000 0
#> [3,] 0.000000 0.0000000 1.0000000 0 0.0000000 0.0000000 0
#> [4,] 0.000000 0.0000000 0.0000000 0 0.7633422 0.2366578 0
#> [5,] 0.000000 0.0000000 0.0000000 0 1.0000000 0.0000000 0
#> [6,] 0.000000 0.0000000 0.0000000 0 0.0000000 0.0000000 1