Create list of matrices
First, I load the data and create a list of matrices:
t <- data.table(read_dta(paste0(path, "data/tmatrices_cz.dta")))
t <- t[complete.cases(t)] # remove missing data
nrow(t)
[1] 729
# loop to create matrices
mm <- list()
for (i in 1:nrow(t)) {
tm <- matrix(NA, 5, 5)
for (h in 1:5) {
for (j in 1:5) {
tm[h, j] <- t[i, get(paste0("prob_p", h, "_k", j))]
}
}
mm[[i]] <- tm
}
I defined functions to get upper and lower mobility differences.
# functions to get lower and upper mobility
ld <- function(mat) {
am <- sum(diag(mat)) / 5
s <- apply(upper.triangle(mat), 1, sum)
t <- diag(mat)
t[1:3] <- s[1:3]
return(sum(t)/5 - am)
}
ud <- function(mat) {
am <- sum(diag(mat)) / 5
s <- apply(upper.triangle(mat), 1, sum)
t <- diag(mat)
t[3:5] <- s[3:5]
return(sum(t)/5 - am)
}
Now I can compute figures and plot them. Most of the difference is observed for lower changes (as expected).
lmd <- sapply(mm, ld)
umd <- sapply(mm, ud)
# get some plots
hist(lmd, main = "lower mobility")

hist(umd, main = "upper mobility")

Explore patterns
The correlation between scores is pretty high 0.89. This is not good if we want to distinguish the effect of lower and upper mobility.
plot(lmd, umd, xlab = "lower mobility", ylab = "upper mobility")

I can use just the extreme quantiles (remove the third):
# functions to get lower and upper mobility
ld <- function(mat) {
am <- sum(diag(mat)) / 5
s <- apply(upper.triangle(mat), 1, sum)
t <- diag(mat)
t[1:2] <- s[1:2]
return(sum(t)/5 - am)
}
ud <- function(mat) {
am <- sum(diag(mat)) / 5
s <- apply(upper.triangle(mat), 1, sum)
t <- diag(mat)
t[4:5] <- s[4:5]
return(sum(t)/5 - am)
}
Get new plots:
lmdb <- sapply(mm, ld)
umdb <- sapply(mm, ud)
# get some plots
hist(lmdb, main = "lower mobility")

hist(umdb, main = "upper mobility")

Now the correlation is lower 0.69:
plot(lmdb, umdb, xlab = "lower mobility", ylab = "upper mobility")

Any ideas or suggestions? To use weights? What about the direction of mobility?
LS0tCnRpdGxlOiAiQ2hldHR5J3MgTW9iaWxpdHkgTWF0cmljZXMiCm91dHB1dDogaHRtbF9ub3RlYm9vawphdXRob3I6IFNlYmFzdGlhbiBEYXphCi0tLQoKCiMgQ3JlYXRlIGxpc3Qgb2YgbWF0cmljZXMgCgpgYGB7ciwgaW5jbHVkZT1GQUxTRX0Kcm0obGlzdD1scyhhbGw9VFJVRSkpCmxpYnJhcnkoaGF2ZW4pCmxpYnJhcnkoc2RhemFyKQpsaWJyYXJ5KG1hdHJpeGNhbGMpCmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeShtY2x1c3QpCnBhdGggPC0gIi9Vc2Vycy9zZGF6YS9Hb29nbGUgRHJpdmUvMDBEaXNzZXJ0YXRpb24vQ2hhcHRlcnMvMDIvIgpgYGAKCkZpcnN0LCBJIGxvYWQgdGhlIGRhdGEgYW5kIGNyZWF0ZSBhIGxpc3Qgb2YgbWF0cmljZXM6IAoKYGBge3J9CnQgPC0gZGF0YS50YWJsZShyZWFkX2R0YShwYXN0ZTAocGF0aCwgImRhdGEvdG1hdHJpY2VzX2N6LmR0YSIpKSkKdCA8LSB0W2NvbXBsZXRlLmNhc2VzKHQpXSAjIHJlbW92ZSBtaXNzaW5nIGRhdGEKbnJvdyh0KQoKIyBsb29wIHRvIGNyZWF0ZSBtYXRyaWNlcyAKbW0gPC0gbGlzdCgpCmZvciAoaSBpbiAxOm5yb3codCkpIHsKICAgICAgdG0gPC0gbWF0cml4KE5BLCA1LCA1KQogICAgICAgICAgICBmb3IgKGggaW4gMTo1KSB7CiAgICAgICAgICAgICAgICAgIGZvciAoaiBpbiAxOjUpIHsKICAgICAgICAgICAgICAgICAgICAgICAgdG1baCwgal0gPC0gdFtpLCBnZXQocGFzdGUwKCJwcm9iX3AiLCBoLCAiX2siLCBqKSldCiAgICAgICAgICAgICAgICAgIH0KICAgICAgICAgICAgfQogICAgICBtbVtbaV1dIDwtIHRtCn0KYGBgCgpJIGRlZmluZWQgZnVuY3Rpb25zIHRvIGdldCAgdXBwZXIgYW5kIGxvd2VyIG1vYmlsaXR5IGRpZmZlcmVuY2VzLgoKYGBge3J9CiMgZnVuY3Rpb25zIHRvIGdldCBsb3dlciBhbmQgdXBwZXIgbW9iaWxpdHkKbGQgPC0gZnVuY3Rpb24obWF0KSB7CiAgICAgIGFtIDwtIHN1bShkaWFnKG1hdCkpIC8gNSAKICAgICAgcyA8LSBhcHBseSh1cHBlci50cmlhbmdsZShtYXQpLCAxLCBzdW0pCiAgICAgIHQgPC0gZGlhZyhtYXQpCiAgICAgIHRbMTozXSA8LSBzWzE6M10KICAgICAgcmV0dXJuKHN1bSh0KS81IC0gIGFtKQp9Cgp1ZCA8LSBmdW5jdGlvbihtYXQpIHsKICAgICAgYW0gPC0gc3VtKGRpYWcobWF0KSkgLyA1IAogICAgICBzIDwtIGFwcGx5KHVwcGVyLnRyaWFuZ2xlKG1hdCksIDEsIHN1bSkKICAgICAgdCA8LSBkaWFnKG1hdCkKICAgICAgdFszOjVdIDwtIHNbMzo1XQogICAgICByZXR1cm4oc3VtKHQpLzUgLSAgYW0pCn0KCmBgYAoKTm93IEkgY2FuIGNvbXB1dGUgZmlndXJlcyBhbmQgcGxvdCB0aGVtLiBNb3N0IG9mIHRoZSBkaWZmZXJlbmNlIGlzIG9ic2VydmVkIGZvciBsb3dlciBjaGFuZ2VzIChhcyBleHBlY3RlZCkuIAoKYGBge3J9CmxtZCA8LSBzYXBwbHkobW0sIGxkKQp1bWQgPC0gc2FwcGx5KG1tLCB1ZCkKIyBnZXQgc29tZSBwbG90cwpoaXN0KGxtZCwgbWFpbiA9ICJsb3dlciBtb2JpbGl0eSIpCmhpc3QodW1kLCBtYWluID0gInVwcGVyIG1vYmlsaXR5IikKYGBgCgojIEV4cGxvcmUgcGF0dGVybnMKClRoZSBjb3JyZWxhdGlvbiBiZXR3ZWVuIHNjb3JlcyBpcyBwcmV0dHkgaGlnaCBgciByb3VuZChjb3IodW1kLCBsbWQpLCAyKWAuIFRoaXMgaXMgbm90IGdvb2QgaWYgd2Ugd2FudCB0byBkaXN0aW5ndWlzaCB0aGUgKmVmZmVjdCogb2YgbG93ZXIgYW5kIHVwcGVyIG1vYmlsaXR5LiAKCmBgYHtyfQpwbG90KGxtZCwgdW1kLCB4bGFiID0gImxvd2VyIG1vYmlsaXR5IiwgeWxhYiA9ICJ1cHBlciBtb2JpbGl0eSIpCmBgYAoKSSBjYW4gdXNlIGp1c3QgdGhlIGV4dHJlbWUgcXVhbnRpbGVzIChyZW1vdmUgdGhlIHRoaXJkKTogCgoKYGBge3J9CiMgZnVuY3Rpb25zIHRvIGdldCBsb3dlciBhbmQgdXBwZXIgbW9iaWxpdHkKbGQgPC0gZnVuY3Rpb24obWF0KSB7CiAgICAgIGFtIDwtIHN1bShkaWFnKG1hdCkpIC8gNSAKICAgICAgcyA8LSBhcHBseSh1cHBlci50cmlhbmdsZShtYXQpLCAxLCBzdW0pCiAgICAgIHQgPC0gZGlhZyhtYXQpCiAgICAgIHRbMToyXSA8LSBzWzE6Ml0gCiAgICAgIHJldHVybihzdW0odCkvNSAtICBhbSkKfQoKdWQgPC0gZnVuY3Rpb24obWF0KSB7CiAgICAgIGFtIDwtIHN1bShkaWFnKG1hdCkpIC8gNSAKICAgICAgcyA8LSBhcHBseSh1cHBlci50cmlhbmdsZShtYXQpLCAxLCBzdW0pCiAgICAgIHQgPC0gZGlhZyhtYXQpCiAgICAgIHRbNDo1XSA8LSBzWzQ6NV0gCiAgICAgIHJldHVybihzdW0odCkvNSAtICBhbSkKfQpgYGAKCkdldCBuZXcgcGxvdHM6ICAKCmBgYHtyfQpsbWRiIDwtIHNhcHBseShtbSwgbGQpCnVtZGIgPC0gc2FwcGx5KG1tLCB1ZCkKIyBnZXQgc29tZSBwbG90cwpoaXN0KGxtZGIsIG1haW4gPSAibG93ZXIgbW9iaWxpdHkiKQpoaXN0KHVtZGIsIG1haW4gPSAidXBwZXIgbW9iaWxpdHkiKQpgYGAKCk5vdyB0aGUgY29ycmVsYXRpb24gaXMgbG93ZXIgYHIgcm91bmQoY29yKHVtZGIsIGxtZGIpLCAyKWA6IAoKYGBge3J9CnBsb3QobG1kYiwgdW1kYiwgeGxhYiA9ICJsb3dlciBtb2JpbGl0eSIsIHlsYWIgPSAidXBwZXIgbW9iaWxpdHkiKQpgYGAKCkFueSBpZGVhcyBvciBzdWdnZXN0aW9ucz8gVG8gdXNlIHdlaWdodHM/IFdoYXQgYWJvdXQgdGhlIGRpcmVjdGlvbiBvZiBtb2JpbGl0eT8KCg==