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
#> 401 401 254002 11.7841084
#> 702 702 556001 12.3098893
#> 801 801 261602 4.4896375
#> 801.1 801 389901 1.1400850
#> 801.2 801 428102 6.2008025
#> 2505 2505 197306 0.6448038
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 695.162 873.416 691.997 847.267 672.098 1008.925 443.593 5232.458
#> 2 213.368 135.464 204.948 166.648 179.516 164.385 166.979 1231.308
#> 3 69.073 48.938 65.550 56.270 53.162 114.630 64.956 472.578
#> 4 109.057 223.515 65.026 92.955 141.701 67.636 47.791 747.680
#> 5 286.413 322.544 549.285 483.227 636.903 377.311 338.946 2994.629
#> 6 29.383 4.337 13.874 31.431 0.000 18.434 26.351 123.809
#> 7 213.705 193.648 203.613 176.199 252.436 172.934 92.004 1304.539
#> Sum 1616.160 1801.862 1794.292 1853.996 1935.816 1924.254 1180.620 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
#> 401 401 254002 11.784108
#> 702 702 556001 12.309889
#> 801 801 261602 4.489638
#> 2505.1 2505 430305 9.655457
#> 5004 5004 16801 2.706425
#> 5201 5201 122301 4.029041
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 705.200 859.073 673.762 833.488 679.826 1043.420 437.689 5232.458
#> 2 200.524 138.610 216.627 165.455 153.125 206.500 150.467 1231.308
#> 3 68.552 58.721 57.686 58.450 49.960 97.175 82.033 472.578
#> 4 99.298 219.134 68.234 92.667 144.394 89.428 34.527 747.680
#> 5 304.198 294.865 550.979 514.525 611.235 359.656 359.171 2994.629
#> 6 21.665 12.055 13.874 21.445 0.000 28.174 26.597 123.809
#> 7 215.973 192.366 212.651 176.608 248.622 178.094 80.224 1304.539
#> Sum 1615.410 1774.824 1793.813 1862.637 1887.162 2002.446 1170.708 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 729.172 846.503 656.953 850.944 677.396 1038.468 433.021 5232.458
#> 2 245.039 140.102 217.025 176.646 163.699 138.538 150.260 1231.308
#> 3 68.552 46.539 57.686 71.356 47.032 111.067 70.345 472.578
#> 4 108.285 248.407 68.234 56.522 129.017 68.535 68.681 747.680
#> 5 253.677 310.272 535.667 506.093 689.101 383.731 316.088 2994.629
#> 6 21.665 12.055 13.874 21.445 0.000 28.174 26.597 123.809
#> 7 194.520 201.590 237.265 163.766 237.141 163.445 106.812 1304.539
#> Sum 1620.909 1805.468 1786.704 1846.772 1943.386 1931.957 1171.803 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.0000000 0.00000000 1.00000000 0.0000000 0.0000000 0.00000000 0.0000000
#> [2,] 0.0000000 1.00000000 0.00000000 0.0000000 0.0000000 0.00000000 0.0000000
#> [3,] 0.0000000 0.37949605 0.00000000 0.0000000 0.0000000 0.09636808 0.5241359
#> [4,] 0.9373993 0.00000000 0.06260073 0.0000000 0.0000000 0.00000000 0.0000000
#> [5,] 0.0000000 0.04272186 0.00000000 0.2379195 0.7193586 0.00000000 0.0000000
#> [6,] 0.0000000 0.00000000 0.00000000 0.3543750 0.6456250 0.00000000 0.0000000