1. Data
1.1 Helpers
clean_prov <- function(x) {
x <- as.character(x)
x <- gsub("内蒙古自治区", "内蒙古", x)
x <- gsub("(壮族|回族|维吾尔|藏族)?自治区$|省$|市$|特别行政区$", "", x)
trimws(x)
}
file_for <- function(keyword, year) {
fs <- list.files(data_folder, pattern = keyword, full.names = TRUE, recursive = TRUE)
fs <- fs[grepl("\\.xls$|\\.xlsx$", fs)]
fs <- fs[grepl(paste0(year, "年"), fs)]
if (length(fs) == 0) NA else fs[1]
}
keep_prov <- function(out) out[nchar(out$province) >= 2 & nchar(out$province) <= 4, ]
1.2 Readers
# Marital status. 2011 has one "currently married" column; the
# "cohabiting" column (from 2016) is not part of the total.
read_marriage <- function(file_path, year) {
raw <- as.data.frame(read_excel(file_path, col_names = FALSE))
d <- raw[6:nrow(raw), ]
prov <- d[[9]]
never <- as.numeric(d[[12]])
first <- as.numeric(d[[14]])
if (year == 2011) {
rem <- 0
div <- as.numeric(d[[16]])
wid <- as.numeric(d[[18]])
} else {
rem <- as.numeric(d[[16]])
div <- as.numeric(d[[18]])
wid <- as.numeric(d[[20]])
}
out <- data.frame(year = year, province = clean_prov(prov),
never = never, first = first, rem = rem, div = div, wid = wid)
out <- out[!is.na(out$never) & !is.na(out$first), ]
out <- keep_prov(out)
out$n_marital <- out$never + out$first + out$rem + out$div + out$wid
out$never_rate <- round(out$never / out$n_marital * 100, 2)
out$first_rate <- round(out$first / out$n_marital * 100, 2)
out[, c("year", "province", "never_rate", "first_rate", "n_marital")]
}
# Sex ratio (males per 100 females)
read_gender <- function(file_path, year) {
raw <- as.data.frame(read_excel(file_path, col_names = FALSE))
d <- raw[6:nrow(raw), ]
male <- as.numeric(d[[12]])
female <- as.numeric(d[[14]])
out <- data.frame(year = year, province = clean_prov(d[[9]]),
sex_ratio = round(male / female * 100, 1),
n_gender = male + female)
out <- out[!is.na(out$sex_ratio), ]
keep_prov(out)
}
# Youth share = ages 15-29 / ages 15-59
read_age <- function(file_path, year) {
raw <- as.data.frame(read_excel(file_path, col_names = FALSE))
header <- as.character(unlist(raw[2, ]))
col_of <- function(seg) {
h <- which(grepl(seg, header) & grepl("人数", header))
if (length(h)) h[1] else NA
}
segs <- c("15-19", "20-24", "25-29", "30-34", "35-39", "40-44", "45-49", "50-54", "55-59")
idxs <- sapply(segs, col_of)
if (any(is.na(idxs))) return(NULL)
d <- raw[6:nrow(raw), ]
counts <- sapply(idxs, function(c) as.numeric(d[[c]]))
out <- data.frame(year = year, province = clean_prov(d[[9]]),
young_share = round(rowSums(counts[, 1:3]) / rowSums(counts) * 100, 1))
out <- out[!is.na(out$young_share), ]
keep_prov(out)
}
# Tertiary-education share (junior college or above)
read_edu <- function(file_path, year) {
raw <- as.data.frame(read_excel(file_path, col_names = FALSE))
header <- as.character(unlist(raw[2, ]))
find_col <- function(label) {
h <- which(grepl(label, header) & grepl("人数", header))
if (!length(h)) h <- which(grepl(label, header) & !grepl("比重", header))
if (length(h)) h[1] else NA
}
d <- raw[6:nrow(raw), ]
gc <- function(lab) {
c <- find_col(lab)
if (is.na(c)) rep(0, nrow(d)) else as.numeric(d[[c]])
}
higher <- gc("大学专科") + gc("大学本科") + gc("研究生")
hs <- gc("高中/中专")
if (all(hs == 0, na.rm = TRUE)) hs <- gc("高中") + gc("中专")
base <- gc("未上过学") + gc("小学") + gc("初中") + higher
out <- data.frame(year = year, province = clean_prov(d[[9]]),
higher_edu = round(higher / (base + hs) * 100, 1))
out <- out[!is.na(out$higher_edu), ]
keep_prov(out)
}
# Agricultural-hukou share, computed from counts (the published share
# columns are defined differently across years)
read_hukou <- function(file_path, year) {
raw <- as.data.frame(read_excel(file_path, col_names = FALSE))
header <- gsub("户口性质_", "", as.character(unlist(raw[2, ])))
d <- raw[6:nrow(raw), ]
cs <- function(name) {
cc <- which(header == name)
if (!length(cc)) rep(0, nrow(d)) else as.numeric(d[[cc[1]]])
}
agri <- cs("农业人数") + cs("农业转居民人数")
total <- cs("农业人数") + cs("非农业人数") + cs("其他户口人数") + cs("其他") +
cs("农业转居民人数") + cs("非农业转居民人数") + cs("居民人数") +
cs("其他(含无户口)人数")
out <- data.frame(year = year, province = clean_prov(d[[9]]),
agri_share = round(agri / total * 100, 1), total = total)
out <- out[!is.na(out$agri_share) & out$total > 0, ]
keep_prov(out)[, c("year", "province", "agri_share")]
}
# 2014 tables have a different layout
read_2014 <- function() {
rd <- function(keyword) {
raw <- as.data.frame(read_excel(file_for(keyword, 2014), col_names = FALSE))
nm <- as.character(unlist(raw[3, ]))
sub <- as.character(unlist(raw[4, ]))
cols <- which(!is.na(nm) & sub == "人数")
d <- raw[5:nrow(raw), ]
m <- sapply(cols, function(c) as.numeric(d[[c]]))
colnames(m) <- gsub("～", "-", nm[cols])
list(prov = clean_prov(d[[1]]), m = m)
}
mar <- rd("婚姻"); gen <- rd("性别"); age <- rd("年龄"); edu <- rd("教育"); huk <- rd("户口")
pick <- function(x, pattern) rowSums(x$m[, grepl(pattern, colnames(x$m)), drop = FALSE])
out <- data.frame(
year = 2014, province = mar$prov,
never_rate = round(mar$m[, "未婚"] / rowSums(mar$m) * 100, 2),
first_rate = round(mar$m[, "初婚"] / rowSums(mar$m) * 100, 2),
n_marital = rowSums(mar$m),
sex_ratio = round(gen$m[, "男"] / gen$m[, "女"] * 100, 1),
n_gender = rowSums(gen$m),
young_share = round(pick(age, "15-19|20-24|25-29") / rowSums(age$m) * 100, 1),
higher_edu = round(pick(edu, "专科|本科|研究生") / rowSums(edu$m) * 100, 1),
agri_share = round(pick(huk, "^农业$|^农业转居民$") / rowSums(huk$m) * 100, 1)
)
out[!(out$province %in% c("合计", "兵团")) & !is.na(out$never_rate), ]
}
1.3 Panel
stack <- function(reader, keyword) bind_rows(lapply(analysis_years, function(y) {
f <- file_for(keyword, y)
if (is.na(f)) return(NULL)
reader(f, y)
}))
marriage <- stack(read_marriage, "婚姻")
gender <- stack(read_gender, "性别")
age <- stack(read_age, "年龄")
edu <- stack(read_edu, "教育")
hukou <- stack(read_hukou, "户口")
name_map <- c(
"安徽"="Anhui","北京"="Beijing","福建"="Fujian","甘肃"="Gansu",
"广东"="Guangdong","广西"="Guangxi","贵州"="Guizhou","海南"="Hainan",
"河北"="Hebei","河南"="Henan","黑龙江"="Heilongjiang","湖北"="Hubei",
"湖南"="Hunan","吉林"="Jilin","江苏"="Jiangsu","江西"="Jiangxi",
"辽宁"="Liaoning","内蒙古"="Inner Mongolia","宁夏"="Ningxia",
"青海"="Qinghai","山东"="Shandong","山西"="Shanxi","陕西"="Shaanxi",
"上海"="Shanghai","四川"="Sichuan","天津"="Tianjin","西藏"="Tibet",
"新疆"="Xinjiang","云南"="Yunnan","浙江"="Zhejiang","重庆"="Chongqing"
)
# Inner Mongolia counted as Central, Guangxi as West
east <- c("Beijing","Tianjin","Hebei","Liaoning","Shanghai","Jiangsu",
"Zhejiang","Fujian","Shandong","Guangdong","Hainan")
central <- c("Shanxi","Jilin","Heilongjiang","Anhui","Jiangxi","Henan",
"Hubei","Hunan","Inner Mongolia")
finish <- function(p) {
p$prov_en <- unname(name_map[p$province])
p <- p[!is.na(p$prov_en), ]
p <- p[!duplicated(p[, c("year", "prov_en")]), ]
p$region <- ifelse(p$prov_en %in% east, "East",
ifelse(p$prov_en %in% central, "Central", "West"))
p
}
panel <- marriage %>%
left_join(gender, by = c("year", "province")) %>%
left_join(age, by = c("year", "province")) %>%
left_join(edu, by = c("year", "province")) %>%
left_join(hukou, by = c("year", "province")) %>%
finish()
p2014 <- finish(read_2014())
panel <- rbind(panel, p2014[, names(panel)])
panel <- panel[order(panel$year, panel$prov_en), ]
rownames(panel) <- NULL
cat("Panel rows:", nrow(panel), "| provinces:", length(unique(panel$prov_en)),
"| years:", paste(sort(unique(panel$year)), collapse = ","), "\n")
Panel rows: 248 | provinces: 31 | years: 2011,2012,2013,2014,2015,2016,2017,2018
stopifnot(nrow(panel) == 248, length(unique(panel$prov_en)) == 31,
!anyNA(panel[, c("never_rate", "young_share", "sex_ratio", "higher_edu", "agri_share")]))
vars <- c("never_rate", "young_share", "sex_ratio", "higher_edu", "agri_share")
core <- c("young_share", "sex_ratio", "higher_edu", "agri_share")
2. Data checks
cat("Table totals by year\n")
Table totals by year
print(aggregate(cbind(n_marital, n_gender) ~ year, data = panel, FUN = sum))
year n_marital n_gender
1 2011 126000 126000
2 2012 154605 154605
3 2013 194795 194795
4 2014 196937 196938
5 2015 203998 203998
6 2016 346352 440017
7 2017 348037 445392
8 2018 149311 151757
cat("\nMarital-table total by province (selected)\n")
Marital-table total by province (selected)
sel <- c("Guangdong", "Zhejiang", "Anhui", "Hebei", "Ningxia", "Tibet")
print(round(xtabs(n_marital ~ prov_en + year, data = panel[panel$prov_en %in% sel, ])))
year
prov_en 2011 2012 2013 2014 2015 2016 2017 2018
Anhui 4000 3999 5000 5000 5000 4656 4350 1892
Guangdong 10000 11992 12000 11998 15000 90253 95272 44303
Hebei 4000 3994 8000 4998 9998 2777 2585 1037
Ningxia 2000 3996 4000 3997 4000 716 622 259
Tibet 2000 3871 3999 3984 4000 412 326 231
Zhejiang 6000 9995 14000 13999 14000 55492 62874 26967
wide <- function(d, v) { w <- reshape(d[, c("prov_en", "year", v)], idvar = "prov_en",
timevar = "year", direction = "wide"); rownames(w) <- w$prov_en; w[, -1] }
cat("\n2014: cross-province correlation with neighbouring years\n")
2014: cross-province correlation with neighbouring years
chk <- t(sapply(vars, function(v) { w <- wide(panel, v)
c(r_2013_2014 = cor(w[[paste0(v, ".2013")]], w[[paste0(v, ".2014")]]),
r_2014_2015 = cor(w[[paste0(v, ".2014")]], w[[paste0(v, ".2015")]]),
r_2012_2013 = cor(w[[paste0(v, ".2012")]], w[[paste0(v, ".2013")]])) }))
print(round(chk, 2))
r_2013_2014 r_2014_2015 r_2012_2013
never_rate 0.95 0.92 0.93
young_share 0.86 0.83 0.85
sex_ratio 0.78 0.68 0.88
higher_edu 0.95 0.96 0.92
agri_share 0.97 0.98 0.96
ws <- wide(panel, "sex_ratio")
ratio14 <- ws$sex_ratio.2014 / ((ws$sex_ratio.2013 + ws$sex_ratio.2015) / 2)
names(ratio14) <- rownames(ws)
cat("2014 sex ratio / mean of 2013 and 2015:\n"); print(round(sort(ratio14), 2))
2014 sex ratio / mean of 2013 and 2015:
Shanghai Guangdong Shandong Tibet Fujian
1.00 1.09 1.11 1.12 1.15
Henan Zhejiang Jiangsu Jiangxi Hunan
1.15 1.16 1.19 1.21 1.21
Hebei Xinjiang Liaoning Shaanxi Yunnan
1.21 1.21 1.24 1.24 1.25
Hainan Inner Mongolia Sichuan Chongqing Guangxi
1.26 1.27 1.27 1.27 1.27
Anhui Gansu Heilongjiang Qinghai Jilin
1.30 1.30 1.32 1.33 1.36
Hubei Shanxi Ningxia Tianjin Beijing
1.37 1.38 1.39 1.39 1.40
Guizhou
1.46
cat("Provinces >10% above:", sum(ratio14 > 1.10), "of 31 | range:", round(range(ratio14), 2), "\n")
Provinces >10% above: 29 of 31 | range: 1 1.46
# Keep the raw panel for one robustness check; set the 2014 sex ratio to missing
panel_raw <- panel
panel$sex_ratio[panel$year == 2014] <- NA
cat("\nExtreme sex ratios (2014 excluded)\n")
Extreme sex ratios (2014 excluded)
print(panel[!is.na(panel$sex_ratio) & (panel$sex_ratio > 150 | panel$sex_ratio < 80), c("year", "prov_en", "sex_ratio")])
year prov_en sex_ratio
10 2011 Hebei 165.4
23 2011 Shandong 79.1
41 2012 Hebei 169.9
54 2012 Shandong 77.2
72 2013 Hebei 164.1
85 2013 Shandong 77.3
152 2015 Tibet 152.4
# Within-year z-scores
zpanel <- panel
for (v in vars) zpanel[[v]] <- ave(panel[[v]], panel$year,
FUN = function(x) (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE))
3. Change over time and provincial differences
trend <- aggregate(panel_raw[, c("never_rate", "first_rate", core)], by = list(year = panel_raw$year), FUN = mean)
trend[, -1] <- round(trend[, -1], 1)
print(trend)
year never_rate first_rate young_share sex_ratio higher_edu agri_share
1 2011 20.6 77.7 38.9 116.0 7.6 84.8
2 2012 21.7 75.0 38.9 117.0 8.7 84.6
3 2013 21.5 75.2 38.1 118.3 9.4 84.7
4 2014 21.7 74.5 38.7 145.6 12.8 84.7
5 2015 18.1 77.6 37.5 113.8 12.4 84.3
6 2016 16.8 79.4 33.7 109.4 12.4 84.1
7 2017 15.9 79.6 31.1 108.2 13.5 82.9
8 2018 12.0 81.7 25.1 105.9 17.3 81.3
y2018 <- panel[panel$year == 2018, c("prov_en", "never_rate")]
y2018 <- y2018[order(y2018$never_rate), ]
cat("\n2018 lowest five:\n"); print(head(y2018, 5))
2018 lowest five:
prov_en never_rate
218 Anhui 3.49
232 Inner Mongolia 5.20
230 Hubei 6.25
237 Ningxia 6.56
242 Shanxi 6.66
cat("2018 highest five:\n"); print(tail(y2018, 5))
2018 highest five:
prov_en never_rate
247 Yunnan 15.45
234 Jiangxi 15.60
223 Guangdong 18.45
229 Henan 19.03
245 Tibet 35.06
cat("2018 range:", round(diff(range(y2018$never_rate)), 1), "pp | without Tibet:",
round(diff(range(y2018$never_rate[y2018$prov_en != "Tibet"])), 1), "pp\n")
2018 range: 31.6 pp | without Tibet: 15.5 pp
w <- reshape(panel[panel$year %in% c(2011, 2018), c("prov_en", "year", "never_rate")],
idvar = "prov_en", timevar = "year", direction = "wide")
w$change <- round(w$never_rate.2018 - w$never_rate.2011, 1)
w <- w[order(w$change), c("prov_en", "change")]
cat("\nChange 2011-2018, largest falls:\n"); print(head(w, 5))
Change 2011-2018, largest falls:
prov_en change
19 Liaoning -18.8
12 Henan -16.2
4 Fujian -16.0
31 Zhejiang -13.1
10 Hebei -12.6
cat("Smallest falls / rises:\n"); print(tail(w, 5))
Smallest falls / rises:
prov_en change
7 Guangxi -5.0
1 Anhui -4.3
20 Ningxia -3.5
29 Xinjiang -3.0
28 Tibet 2.4
desc <- t(sapply(panel[vars], function(x) c(N = sum(!is.na(x)), Mean = mean(x, na.rm = TRUE),
SD = sd(x, na.rm = TRUE), Min = min(x, na.rm = TRUE), Median = median(x, na.rm = TRUE), Max = max(x, na.rm = TRUE))))
print(round(desc, 1))
N Mean SD Min Median Max
never_rate 248 18.5 6.7 3.5 18.1 42.3
young_share 248 35.2 6.3 17.9 36.0 53.3
sex_ratio 217 112.7 12.9 77.2 111.6 169.9
higher_edu 248 11.8 5.8 2.1 10.7 36.1
agri_share 248 83.9 6.7 61.5 85.7 94.4
for (v in vars) { i <- which.min(panel[[v]]); j <- which.max(panel[[v]])
cat(sprintf("%-12s min: %s %d | max: %s %d\n", v, panel$prov_en[i], panel$year[i], panel$prov_en[j], panel$year[j])) }
never_rate min: Anhui 2018 | max: Tibet 2014
young_share min: Inner Mongolia 2018 | max: Henan 2012
sex_ratio min: Shandong 2012 | max: Hebei 2012
higher_edu min: Tibet 2013 | max: Beijing 2018
agri_share min: Beijing 2016 | max: Henan 2015
par(mar = c(4.5, 4.5, 1, 1))
plot(trend$year, trend$first_rate, type = "o", col = "steelblue", lwd = 2, pch = 19,
ylim = c(0, 100), xlab = "Year", ylab = "Share of floating population (%)", xaxt = "n")
axis(1, 2011:2018)
lines(trend$year, trend$never_rate, type = "o", col = "tomato", lwd = 2, pch = 19)
text(trend$year, trend$first_rate + 5, sprintf("%.1f", trend$first_rate), cex = 0.7)
text(trend$year, trend$never_rate + 5, sprintf("%.1f", trend$never_rate), cex = 0.7)
legend("right", c("First-married", "Never-married"), col = c("steelblue", "tomato"), lwd = 2, pch = 19, bty = "n", cex = 0.8)
par(mar = c(4.5, 9, 1, 2))
bp <- barplot(y2018$never_rate, names.arg = y2018$prov_en, horiz = TRUE, las = 1, col = "steelblue",
border = NA, cex.names = 0.85, xlab = "Never-married share, 2018 (%)", xlim = c(0, 40))
text(y2018$never_rate + 1.8, bp, sprintf("%.1f", y2018$never_rate), cex = 0.75)
abline(v = mean(y2018$never_rate), lty = 2)
text(mean(y2018$never_rate), max(bp) + 1, sprintf("mean %.1f", mean(y2018$never_rate)), pos = 4, cex = 0.75, xpd = NA)
4. Correlations
cm <- cor(panel[vars], use = "pairwise.complete.obs")
print(round(cm, 2))
never_rate young_share sex_ratio higher_edu agri_share
never_rate 1.00 0.62 0.41 -0.34 0.18
young_share 0.62 1.00 0.18 -0.30 0.40
sex_ratio 0.41 0.18 1.00 -0.42 0.18
higher_edu -0.34 -0.30 -0.42 1.00 -0.62
agri_share 0.18 0.40 0.18 -0.62 1.00
cat("N per variable:", colSums(!is.na(panel[vars])), "\n")
N per variable: 248 248 217 248 248
cat("\nAfter within-year standardisation:\n")
After within-year standardisation:
print(round(cor(zpanel[vars], use = "pairwise.complete.obs")["never_rate", core], 2))
young_share sex_ratio higher_edu agri_share
0.43 0.34 -0.20 0.16
cat("R2 of youth share alone (N = 248):", round(summary(lm(never_rate ~ young_share, panel))$r.squared, 2), "\n")
R2 of youth share alone (N = 248): 0.38
cat("Province means (31 provinces):\n")
Province means (31 provinces):
prov_avg <- aggregate(panel[vars], by = list(prov_en = panel$prov_en), FUN = mean, na.rm = TRUE)
rownames(prov_avg) <- prov_avg$prov_en
print(round(cor(prov_avg[vars])["never_rate", core], 2))
young_share sex_ratio higher_edu agri_share
0.39 0.40 -0.21 0.14
partial_corr <- function(df, x, y, z) {
d <- df[complete.cases(df[, c(x, y, z)]), ]
cor(residuals(lm(d[[x]] ~ d[[z]])), residuals(lm(d[[y]] ~ d[[z]])))
}
raw_r <- sapply(core[-1], function(v) cor(panel$never_rate, panel[[v]], use = "complete.obs"))
par_r <- sapply(core[-1], function(v) partial_corr(panel, v, "never_rate", "young_share"))
cat("\nRaw vs partial (controlling for youth share):\n")
Raw vs partial (controlling for youth share):
print(round(rbind(raw = raw_r, partial = par_r), 2))
sex_ratio higher_edu agri_share
raw 0.41 -0.34 0.18
partial 0.39 -0.21 -0.09
blue <- rgb(0.27, 0.51, 0.71, 0.55)
sc <- function(x, xlab) {
par(mar = c(4.5, 4.5, 1, 1))
d <- panel_raw; used <- !(x == "sex_ratio" & d$year == 2014)
plot(d[[x]], d$never_rate, type = "n", xlab = xlab, ylab = "Never-married share (%)")
points(d[[x]][used], d$never_rate[used], pch = 21, bg = blue, col = "white")
if (any(!used)) {
points(d[[x]][!used], d$never_rate[!used], pch = 1, col = "grey55")
legend("bottomright", c("Used", "2014 (not used)"), pch = c(21, 1), pt.bg = c(blue, NA),
col = c("white", "grey55"), bty = "n", cex = 0.8)
}
abline(lm(d$never_rate[used] ~ d[[x]][used]), lwd = 2, lty = 2)
legend("topleft", bty = "n", cex = 0.9, legend = sprintf("r = %.2f (N = %d)",
cor(d[[x]][used], d$never_rate[used]), sum(used)))
}
sc("young_share", "Youth share, ages 15-29 (%)")
sc("sex_ratio", "Sex ratio (males per 100 females)")
sc("higher_edu", "Tertiary-education share (%)")
sc("agri_share", "Agricultural-hukou share (%)")
lab5 <- c("Never-married", "Youth 15-29", "Sex ratio", "Tertiary edu.", "Agri. hukou")
par(mar = c(7, 7, 1, 1))
image(1:5, 1:5, t(cm[5:1, ]), axes = FALSE, xlab = "", ylab = "",
col = colorRampPalette(c("steelblue", "white", "tomato"))(50), zlim = c(-1, 1))
axis(1, 1:5, lab5, las = 2, cex.axis = 0.75); axis(2, 1:5, rev(lab5), las = 1, cex.axis = 0.75)
for (i in 1:5) for (j in 1:5) text(j, 6 - i, sprintf("%.2f", cm[i, j]), cex = 0.8)
rp <- rbind(raw_r, par_r)
bp <- barplot(rp, beside = TRUE, ylim = c(-0.5, 0.6),
names.arg = c("Sex ratio", "Tertiary education", "Agricultural hukou"),
col = c("steelblue", "tomato"), border = NA, ylab = "Correlation with never-married share")
text(bp, rp + ifelse(rp >= 0, 0.04, -0.04), sprintf("%+.2f", rp), cex = 0.7)
legend("topright", c("Raw correlation", "Controlling for youth share"), fill = c("steelblue", "tomato"), border = NA, bty = "n", cex = 0.8)
abline(h = 0)
5. Regressions
Inference uses province-clustered standard errors (CR1, t with G − 1
df) and wild cluster bootstrap p-values (WCR, Rademacher weights, null
imposed; Cameron, Gelbach & Miller, 2008).
reg <- panel[complete.cases(panel[, vars]), c(vars, "year", "prov_en", "n_marital")]
stopifnot(nrow(reg) == 217)
cl_vcov <- function(X, e, g, w) {
n <- nrow(X); k <- ncol(X); G <- length(unique(g))
bread <- solve(crossprod(X * sqrt(w)))
S <- rowsum(X * (w * e), g)
bread %*% crossprod(S) %*% bread * (G / (G - 1)) * ((n - 1) / (n - k))
}
wild_p <- function(X, y, g, w, j, t0, B) {
Xr <- X[, -j, drop = FALSE]
br <- solve(crossprod(Xr * sqrt(w)), crossprod(Xr, w * y))
yr <- drop(Xr %*% br); er <- y - yr
A <- solve(crossprod(X * sqrt(w)), t(X * w))
gi <- as.integer(factor(g)); G <- max(gi)
tb <- numeric(B)
for (b in seq_len(B)) {
ys <- yr + er * sample(c(-1, 1), G, replace = TRUE)[gi]
bs <- A %*% ys; es <- ys - drop(X %*% bs)
tb[b] <- bs[j] / sqrt(cl_vcov(X, es, g, w)[j, j])
}
mean(abs(tb) >= abs(t0))
}
fit_table <- function(formula, dat, w = NULL, B = 9999, seed = 2026) {
if (is.null(w)) w <- rep(1, nrow(dat))
dat$.w <- w
m <- lm(formula, data = dat, weights = .w)
X <- model.matrix(m); y <- dat$never_rate; g <- dat$prov_en; G <- length(unique(g))
b <- coef(m); V <- cl_vcov(X, residuals(m), g, w); se <- sqrt(diag(V))
set.seed(seed)
out <- t(sapply(core, function(v) {
j <- which(colnames(X) == v); t0 <- b[j] / se[j]
c(Estimate = unname(b[j]), SE_cluster = unname(se[j]),
CI_low = unname(b[j] - qt(0.975, G - 1) * se[j]), CI_high = unname(b[j] + qt(0.975, G - 1) * se[j]),
p_cluster = unname(2 * pt(-abs(t0), G - 1)),
p_wild = wild_p(X, y, g, w, j, unname(t0), B),
p_ordinary = summary(m)$coefficients[v, "Pr(>|t|)"])
}))
list(model = m, table = round(out, 3), r2 = summary(m)$r.squared)
}
show_fit <- function(res, label, n) {
cat("\n", label, " | N = ", n, " | R2 = ", round(res$r2, 3), "\n", sep = ""); print(res$table)
}
5.1 Main models
f1 <- never_rate ~ young_share + sex_ratio + higher_edu + agri_share
specs <- list("Model 1: pooled OLS" = f1,
"Model 2: + year FE" = update(f1, . ~ . + factor(year)),
"Model 3: two-way FE" = update(f1, . ~ . + factor(year) + factor(prov_en)))
fits <- list()
for (nm in names(specs)) {
fits[[nm]] <- fit_table(specs[[nm]], reg)
show_fit(fits[[nm]], nm, nrow(reg))
}
Model 1: pooled OLS | N = 217 | R2 = 0.513
Estimate SE_cluster CI_low CI_high p_cluster p_wild p_ordinary
young_share 0.599 0.097 0.401 0.797 0.000 0.000 0.000
sex_ratio 0.129 0.037 0.054 0.204 0.001 0.025 0.000
higher_edu -0.222 0.128 -0.484 0.040 0.093 0.086 0.003
agri_share -0.201 0.121 -0.448 0.047 0.108 0.137 0.001
Model 2: + year FE | N = 217 | R2 = 0.564
Estimate SE_cluster CI_low CI_high p_cluster p_wild p_ordinary
young_share 0.908 0.146 0.610 1.205 0.000 0.000 0
sex_ratio 0.141 0.038 0.063 0.219 0.001 0.024 0
higher_edu -0.433 0.174 -0.789 -0.077 0.019 0.017 0
agri_share -0.370 0.138 -0.651 -0.088 0.012 0.021 0
Model 3: two-way FE | N = 217 | R2 = 0.917
Estimate SE_cluster CI_low CI_high p_cluster p_wild p_ordinary
young_share 0.814 0.117 0.575 1.052 0.000 0.000 0.000
sex_ratio 0.075 0.032 0.009 0.141 0.027 0.007 0.000
higher_edu -0.112 0.122 -0.361 0.137 0.365 0.363 0.235
agri_share -0.152 0.099 -0.354 0.049 0.134 0.200 0.055
cat("\nYear effects in Model 2 (reference 2011):\n")
Year effects in Model 2 (reference 2011):
yc <- coef(fits[["Model 2: + year FE"]]$model); print(round(yc[grepl("factor\\(year\\)", names(yc))], 2))
factor(year)2012 factor(year)2013 factor(year)2015 factor(year)2016
1.40 2.02 0.94 3.67
factor(year)2017 factor(year)2018
5.32 8.23
dm <- function(v) residuals(lm(v ~ factor(year) + factor(prov_en), data = reg))
within <- lm(dm(reg$never_rate) ~ dm(reg$young_share) + dm(reg$sex_ratio) + dm(reg$higher_edu) + dm(reg$agri_share))
cat("Within R2 of Model 3:", round(summary(within)$r.squared, 3), "\n")
Within R2 of Model 3: 0.509
5.2 Robustness checks (Model 3)
m3f <- specs[["Model 3: two-way FE"]]
wrel <- reg$n_marital / ave(reg$n_marital, reg$year) # sample size relative to its year
rob <- list(
"Weighted by relative sample size" = list(d = reg, w = wrel),
"Without 2016-2017" = list(d = reg[!reg$year %in% c(2016, 2017), ]),
"Without Tibet" = list(d = reg[reg$prov_en != "Tibet", ]),
"2011-2015 only" = list(d = reg[reg$year <= 2015, ]),
"Adding 2014 with its raw sex ratio" = list(d = panel_raw[, c(vars, "year", "prov_en", "n_marital")]))
rob_fits <- list("Main Model 3" = fits[["Model 3: two-way FE"]])
rob_n <- c("Main Model 3" = nrow(reg))
for (nm in names(rob)) {
r <- fit_table(m3f, rob[[nm]]$d, w = rob[[nm]]$w, B = 4999)
show_fit(r, nm, nrow(rob[[nm]]$d))
rob_fits[[nm]] <- r; rob_n[nm] <- nrow(rob[[nm]]$d)
}
Weighted by relative sample size | N = 217 | R2 = 0.926
Estimate SE_cluster CI_low CI_high p_cluster p_wild p_ordinary
young_share 0.772 0.081 0.606 0.938 0.000 0.000 0.000
sex_ratio 0.044 0.026 -0.010 0.097 0.106 0.066 0.017
higher_edu 0.001 0.143 -0.290 0.293 0.992 0.991 0.987
agri_share -0.043 0.142 -0.334 0.247 0.762 0.762 0.578
Without 2016-2017 | N = 155 | R2 = 0.943
Estimate SE_cluster CI_low CI_high p_cluster p_wild p_ordinary
young_share 0.643 0.149 0.338 0.948 0.000 0.000 0.000
sex_ratio 0.054 0.032 -0.012 0.119 0.104 0.046 0.019
higher_edu -0.023 0.149 -0.328 0.281 0.876 0.875 0.830
agri_share -0.066 0.162 -0.397 0.265 0.687 0.681 0.558
Without Tibet | N = 210 | R2 = 0.909
Estimate SE_cluster CI_low CI_high p_cluster p_wild p_ordinary
young_share 0.783 0.118 0.541 1.025 0.000 0.000 0.000
sex_ratio 0.076 0.034 0.006 0.145 0.033 0.008 0.000
higher_edu -0.069 0.121 -0.316 0.178 0.571 0.563 0.466
agri_share -0.139 0.102 -0.347 0.069 0.181 0.236 0.074
2011-2015 only | N = 124 | R2 = 0.964
Estimate SE_cluster CI_low CI_high p_cluster p_wild p_ordinary
young_share 0.784 0.093 0.595 0.973 0.000 0.000 0.000
sex_ratio 0.023 0.021 -0.020 0.066 0.282 0.244 0.251
higher_edu -0.171 0.193 -0.566 0.223 0.382 0.454 0.235
agri_share -0.074 0.205 -0.492 0.344 0.719 0.697 0.569
Adding 2014 with its raw sex ratio | N = 248 | R2 = 0.92
Estimate SE_cluster CI_low CI_high p_cluster p_wild p_ordinary
young_share 0.832 0.117 0.593 1.071 0.000 0.000 0.000
sex_ratio 0.068 0.029 0.009 0.127 0.025 0.012 0.000
higher_edu -0.122 0.133 -0.393 0.149 0.365 0.368 0.169
agri_share -0.147 0.109 -0.370 0.076 0.188 0.234 0.055
cat("\nSummary by variable\n")
Summary by variable
for (v in core) {
s <- t(sapply(rob_fits, function(r) r$table[v, c("Estimate", "SE_cluster", "CI_low", "CI_high", "p_cluster", "p_wild")]))
s <- data.frame(N = rob_n[rownames(s)], s, check.names = FALSE)
cat("\n", v, ":\n", sep = ""); print(s)
}
young_share:
N Estimate SE_cluster CI_low CI_high
Main Model 3 217 0.814 0.117 0.575 1.052
Weighted by relative sample size 217 0.772 0.081 0.606 0.938
Without 2016-2017 155 0.643 0.149 0.338 0.948
Without Tibet 210 0.783 0.118 0.541 1.025
2011-2015 only 124 0.784 0.093 0.595 0.973
Adding 2014 with its raw sex ratio 248 0.832 0.117 0.593 1.071
p_cluster p_wild
Main Model 3 0 0
Weighted by relative sample size 0 0
Without 2016-2017 0 0
Without Tibet 0 0
2011-2015 only 0 0
Adding 2014 with its raw sex ratio 0 0
sex_ratio:
N Estimate SE_cluster CI_low CI_high
Main Model 3 217 0.075 0.032 0.009 0.141
Weighted by relative sample size 217 0.044 0.026 -0.010 0.097
Without 2016-2017 155 0.054 0.032 -0.012 0.119
Without Tibet 210 0.076 0.034 0.006 0.145
2011-2015 only 124 0.023 0.021 -0.020 0.066
Adding 2014 with its raw sex ratio 248 0.068 0.029 0.009 0.127
p_cluster p_wild
Main Model 3 0.027 0.007
Weighted by relative sample size 0.106 0.066
Without 2016-2017 0.104 0.046
Without Tibet 0.033 0.008
2011-2015 only 0.282 0.244
Adding 2014 with its raw sex ratio 0.025 0.012
higher_edu:
N Estimate SE_cluster CI_low CI_high
Main Model 3 217 -0.112 0.122 -0.361 0.137
Weighted by relative sample size 217 0.001 0.143 -0.290 0.293
Without 2016-2017 155 -0.023 0.149 -0.328 0.281
Without Tibet 210 -0.069 0.121 -0.316 0.178
2011-2015 only 124 -0.171 0.193 -0.566 0.223
Adding 2014 with its raw sex ratio 248 -0.122 0.133 -0.393 0.149
p_cluster p_wild
Main Model 3 0.365 0.363
Weighted by relative sample size 0.992 0.991
Without 2016-2017 0.876 0.875
Without Tibet 0.571 0.563
2011-2015 only 0.382 0.454
Adding 2014 with its raw sex ratio 0.365 0.368
agri_share:
N Estimate SE_cluster CI_low CI_high
Main Model 3 217 -0.152 0.099 -0.354 0.049
Weighted by relative sample size 217 -0.043 0.142 -0.334 0.247
Without 2016-2017 155 -0.066 0.162 -0.397 0.265
Without Tibet 210 -0.139 0.102 -0.347 0.069
2011-2015 only 124 -0.074 0.205 -0.492 0.344
Adding 2014 with its raw sex ratio 248 -0.147 0.109 -0.370 0.076
p_cluster p_wild
Main Model 3 0.134 0.200
Weighted by relative sample size 0.762 0.762
Without 2016-2017 0.687 0.681
Without Tibet 0.181 0.236
2011-2015 only 0.719 0.697
Adding 2014 with its raw sex ratio 0.188 0.234
5.3 Diagnostics
vif_of <- function(d) sapply(core, function(v)
1 / (1 - summary(lm(reformulate(setdiff(core, v), response = v), data = d))$r.squared))
cat("VIF, pooled\n"); print(round(vif_of(reg), 2))
VIF, pooled
young_share sex_ratio higher_edu agri_share
1.23 1.24 1.97 1.77
dmd <- as.data.frame(lapply(reg[core], dm))
cat("VIF, after year and province fixed effects\n"); print(round(vif_of(dmd), 2))
VIF, after year and province fixed effects
young_share sex_ratio higher_edu agri_share
1.09 1.09 1.53 1.54
m3 <- fits[["Model 3: two-way FE"]]$model
big <- order(-abs(residuals(m3)))[1:5]
cat("\nFive largest residuals of Model 3\n")
Five largest residuals of Model 3
print(data.frame(reg[big, c("year", "prov_en", "never_rate")], residual = round(residuals(m3)[big], 1)))
year prov_en never_rate residual
81 2013 Liaoning 36.78 6.6
240 2018 Shandong 7.18 6.0
147 2015 Shandong 7.87 -5.1
245 2018 Tibet 35.06 4.9
225 2018 Guizhou 10.69 -4.4
par(mfrow = c(1, 2), mar = c(4.5, 4.5, 2, 1))
plot(fitted(m3), residuals(m3), pch = 21, bg = blue, col = "white",
xlab = "Fitted never-married share (%)", ylab = "Residual", main = "Residuals vs fitted")
abline(h = 0, lty = 2)
qqnorm(residuals(m3), pch = 21, bg = blue, col = "white", main = "Normal Q-Q"); qqline(residuals(m3), lty = 2)
5.4 Inference checks for the sex ratio
G <- length(unique(reg$prov_en))
tcr <- qt(0.975, G - 1)
X3 <- model.matrix(m3); e3 <- residuals(m3); gi <- as.integer(factor(reg$prov_en))
bread3 <- solve(crossprod(X3)); j_sr <- which(colnames(X3) == "sex_ratio")
S3 <- rowsum(X3 * e3, gi)
contrib <- (S3 %*% bread3[, j_sr])^2
share <- sort(setNames(drop(contrib) / sum(contrib), levels(factor(reg$prov_en))), decreasing = TRUE)
cat("Share of the clustered variance of the sex-ratio coefficient, top 5 provinces\n")
Share of the clustered variance of the sex-ratio coefficient, top 5 provinces
print(round(head(share, 5), 2))
Liaoning Hebei Shandong Anhui Inner Mongolia
0.34 0.23 0.23 0.04 0.02
loo <- sapply(levels(factor(reg$prov_en)), function(p)
coef(lm(m3f, data = reg[reg$prov_en != p, ]))["sex_ratio"])
names(loo) <- levels(factor(reg$prov_en))
cat("Coefficient when one province is dropped (range):", round(range(loo), 3), "\n")
Coefficient when one province is dropped (range): 0.056 0.094
cat("Dropping the three largest contributors:", round(loo[names(head(share, 3))], 3), "\n")
Dropping the three largest contributors: 0.056 0.094 0.057
cr3_se <- sqrt((G - 1) / G * sum((loo - mean(loo))^2))
cat("CR3 (jackknife) SE:", round(cr3_se, 3),
"| p =", round(2 * pt(-abs(coef(m3)["sex_ratio"] / cr3_se), G - 1), 3), "\n")
CR3 (jackknife) SE: 0.035 | p = 0.039
wild_dist <- function(j, B = 9999, restricted = TRUE, webb = FALSE, seed = 2026) {
set.seed(seed)
y <- reg$never_rate; b <- coef(m3)
if (restricted) { Xr <- X3[, -j]; yr <- drop(Xr %*% solve(crossprod(Xr), crossprod(Xr, y))) } else yr <- drop(X3 %*% b)
er <- y - yr; A <- bread3 %*% t(X3); wts <- if (webb) c(-sqrt(1.5), -1, -sqrt(0.5), sqrt(0.5), 1, sqrt(1.5)) else c(-1, 1)
sapply(seq_len(B), function(i) {
ys <- yr + er * sample(wts, G, replace = TRUE)[gi]; bs <- A %*% ys; es <- ys - drop(X3 %*% bs)
(bs[j] - if (restricted) 0 else b[j]) / sqrt(cl_vcov(X3, es, gi, rep(1, length(es)))[j, j])
})
}
t_sr <- coef(m3)[j_sr] / sqrt(cl_vcov(X3, e3, gi, rep(1, length(e3)))[j_sr, j_sr])
for (cfg in list(c("WCR, Rademacher", TRUE, FALSE), c("WCR, Webb", TRUE, TRUE), c("WCU, Rademacher", FALSE, FALSE))) {
tb <- wild_dist(j_sr, restricted = as.logical(cfg[2]), webb = as.logical(cfg[3]))
cat(sprintf("%-16s p = %.3f | SD of t* = %.2f | 95th percentile of |t*| = %.2f (t(30) critical value %.2f)\n",
cfg[1], mean(abs(tb) >= abs(t_sr)), sd(tb), quantile(abs(tb), 0.95), tcr))
}
WCR, Rademacher p = 0.006 | SD of t* = 0.97 | 95th percentile of |t*| = 1.83 (t(30) critical value 2.04)
WCR, Webb p = 0.006 | SD of t* = 0.98 | 95th percentile of |t*| = 1.82 (t(30) critical value 2.04)
WCU, Rademacher p = 0.041 | SD of t* = 1.10 | 95th percentile of |t*| = 2.26 (t(30) critical value 2.04)
cat("Model 3: 95% CIs (clustered SE with t(G-1); wild bootstrap-t, WCU)\n")
Model 3: 95% CIs (clustered SE with t(G-1); wild bootstrap-t, WCU)
V3 <- cl_vcov(X3, e3, gi, rep(1, length(e3)))
ci_tab <- t(sapply(core, function(v) {
j <- which(colnames(X3) == v); b <- coef(m3)[j]; s <- sqrt(V3[j, j])
q <- quantile(abs(wild_dist(j, B = 4999, restricted = FALSE)), 0.95)
c(Estimate = unname(b), CI_low = unname(b - tcr * s), CI_high = unname(b + tcr * s),
boot_low = unname(b - q * s), boot_high = unname(b + q * s))
}))
print(round(ci_tab, 3))
Estimate CI_low CI_high boot_low boot_high
young_share 0.814 0.575 1.052 0.572 1.055
sex_ratio 0.075 0.009 0.141 0.002 0.147
higher_edu -0.112 -0.361 0.137 -0.345 0.120
agri_share -0.152 -0.354 0.049 -0.328 0.023
cat("\nCorrelations: 95% CIs from a province bootstrap (2000 resamples)\n")
Correlations: 95% CIs from a province bootstrap (2000 resamples)
pc_stats <- function(d, v) {
raw <- cor(d$never_rate, d[[v]])
par <- cor(residuals(lm(d$never_rate ~ d$young_share)), residuals(lm(d[[v]] ~ d$young_share)))
c(raw = raw, partial = par, change = raw - par)
}
provs <- unique(panel$prov_en); set.seed(2026)
for (v in c("agri_share", "higher_edu", "sex_ratio")) {
d0 <- panel[complete.cases(panel[, c("never_rate", "young_share", v)]), ]
est <- pc_stats(d0, v)
bs <- t(replicate(2000, { pick <- sample(provs, length(provs), replace = TRUE)
pc_stats(do.call(rbind, lapply(pick, function(p) d0[d0$prov_en == p, ])), v) }))
ci <- apply(bs, 2, quantile, c(0.025, 0.975))
cat(sprintf("%-11s raw %.3f [%.3f, %.3f] | partial %.3f [%.3f, %.3f] | change %.3f [%.3f, %.3f]\n", v,
est[1], ci[1, 1], ci[2, 1], est[2], ci[1, 2], ci[2, 2], est[3], ci[1, 3], ci[2, 3]))
}
agri_share raw 0.183 [-0.077, 0.394] | partial -0.093 [-0.425, 0.172] | change 0.276 [0.139, 0.464]
higher_edu raw -0.340 [-0.532, -0.187] | partial -0.205 [-0.449, 0.061] | change -0.134 [-0.323, -0.033]
sex_ratio raw 0.410 [0.224, 0.560] | partial 0.390 [0.183, 0.538] | change 0.019 [-0.074, 0.135]
6. Between and within provinces
6.1 Mundlak model
mund <- reg
for (v in core) { mund[[paste0(v, "_b")]] <- ave(reg[[v]], reg$prov_en); mund[[paste0(v, "_w")]] <- reg[[v]] - mund[[paste0(v, "_b")]] }
fm <- as.formula(paste("never_rate ~", paste(c(paste0(core, "_w"), paste0(core, "_b")), collapse = " + "), "+ factor(year)"))
mm <- lm(fm, data = mund); Xm <- model.matrix(mm); Vm <- cl_vcov(Xm, residuals(mm), gi, rep(1, nrow(Xm)))
mtab <- t(sapply(core, function(v) {
jw <- which(colnames(Xm) == paste0(v, "_w")); jb <- which(colnames(Xm) == paste0(v, "_b"))
bw <- coef(mm)[jw]; bb <- coef(mm)[jb]; sw <- sqrt(Vm[jw, jw]); sb <- sqrt(Vm[jb, jb])
sd_ <- sqrt(Vm[jb, jb] + Vm[jw, jw] - 2 * Vm[jb, jw])
c(within = bw, w_low = bw - tcr * sw, w_high = bw + tcr * sw, p_within = 2 * pt(-abs(bw / sw), G - 1),
between = bb, b_low = bb - tcr * sb, b_high = bb + tcr * sb, p_between = 2 * pt(-abs(bb / sb), G - 1),
p_difference = 2 * pt(-abs((bb - bw) / sd_), G - 1))
}))
colnames(mtab) <- c("within", "w_low", "w_high", "p_within", "between", "b_low", "b_high", "p_between", "p_diff")
print(round(mtab, 3))
within w_low w_high p_within between b_low b_high p_between
young_share 0.814 0.591 1.036 0.000 0.977 0.485 1.468 0.000
sex_ratio 0.075 0.014 0.136 0.018 0.185 0.047 0.323 0.010
higher_edu -0.112 -0.345 0.120 0.332 -0.461 -0.889 -0.033 0.036
agri_share -0.152 -0.340 0.036 0.109 -0.429 -0.819 -0.038 0.033
p_diff
young_share 0.498
sex_ratio 0.128
higher_edu 0.132
agri_share 0.245
Rm <- matrix(0, 4, ncol(Xm)); for (i in 1:4) { Rm[i, which(colnames(Xm) == paste0(core[i], "_b"))] <- 1; Rm[i, which(colnames(Xm) == paste0(core[i], "_w"))] <- -1 }
rb <- Rm %*% coef(mm); Fj <- drop(t(rb) %*% solve(Rm %*% Vm %*% t(Rm)) %*% rb) / 4
cat("Joint test, between = within for all four: F =", round(Fj, 2), "| p =", round(pf(Fj, 4, G - 1, lower.tail = FALSE), 3), "\n")
Joint test, between = within for all four: F = 1.32 | p = 0.287
par(mar = c(4.5, 9, 1, 1))
lab4 <- c("Youth share", "Sex ratio", "Tertiary edu.", "Agri. hukou"); yy <- 4:1
plot(NA, xlim = range(mtab[, c("w_low", "w_high", "b_low", "b_high")], 0), ylim = c(0.5, 4.5), yaxt = "n",
ylab = "", xlab = "Coefficient (percentage points per unit), 95% CI")
axis(2, yy, lab4, las = 1); abline(v = 0, lty = 2)
segments(mtab[, "w_low"], yy + 0.12, mtab[, "w_high"], yy + 0.12, col = "steelblue", lwd = 2); points(mtab[, "within"], yy + 0.12, pch = 19, col = "steelblue")
segments(mtab[, "b_low"], yy - 0.12, mtab[, "b_high"], yy - 0.12, col = "tomato", lwd = 2); points(mtab[, "between"], yy - 0.12, pch = 17, col = "tomato")
legend("bottomright", c("Within provinces", "Between provinces"), col = c("steelblue", "tomato"), pch = c(19, 17), lwd = 2, bty = "n", cex = 0.8)
6.2 Variance decomposition
tot <- sum((reg$never_rate - mean(reg$never_rate))^2)
cat("Between-province share:", round(sum((ave(reg$never_rate, reg$prov_en) - mean(reg$never_rate))^2) / tot, 2),
"| year share:", round(sum((ave(reg$never_rate, reg$year) - mean(reg$never_rate))^2) / tot, 2), "\n")
Between-province share: 0.57 | year share: 0.26
r2 <- function(f, d) summary(lm(f, data = d))$r.squared
cat("Between provinces (31 means), R2 of single variables:",
paste(core, round(sapply(core, function(v) r2(reformulate(v, "never_rate"), prov_avg)), 3)), "\n")
Between provinces (31 means), R2 of single variables: young_share 0.149 sex_ratio 0.161 higher_edu 0.046 agri_share 0.02
full_b <- r2(reformulate(core, "never_rate"), prov_avg); cat("R2 of all four:", round(full_b, 3), "\n")
R2 of all four: 0.438
shap <- sapply(core, function(v) { others <- setdiff(core, v); tot_ <- 0
for (k in 0:3) for (S in if (k == 0) list(character(0)) else combn(others, k, simplify = FALSE)) {
w <- factorial(length(S)) * factorial(3 - length(S)) / factorial(4)
with_v <- r2(reformulate(c(S, v), "never_rate"), prov_avg)
without <- if (length(S)) r2(reformulate(S, "never_rate"), prov_avg) else 0
tot_ <- tot_ + w * (with_v - without) }
tot_ })
cat("Shapley shares of the between-province R2:", paste(core, round(shap, 3)), "\n")
Shapley shares of the between-province R2: young_share 0.201 sex_ratio 0.134 higher_edu 0.054 agri_share 0.05
cat("Within provinces, R2 youth only:",
round(summary(lm(dm(reg$never_rate) ~ dm(reg$young_share)))$r.squared, 3),
"| all four:", round(summary(within)$r.squared, 3), "\n")
Within provinces, R2 youth only: 0.451 | all four: 0.509
cat("Cross-sectional R2 of youth share by year:\n")
Cross-sectional R2 of youth share by year:
print(round(sapply(split(panel, panel$year), function(d) r2(never_rate ~ young_share, d)), 2))
2011 2012 2013 2014 2015 2016 2017 2018
0.28 0.33 0.19 0.26 0.09 0.03 0.17 0.27
6.3 Influence of single provinces
between_stats <- function(pa) {
r <- sapply(core, function(v) cor(pa$never_rate, pa[[v]]))
r2 <- sapply(core, function(v) summary(lm(reformulate(v, "never_rate"), data = pa))$r.squared)
c(setNames(r, paste0("r_", core)), setNames(r2, paste0("R2_", core)),
R2_all4 = summary(lm(reformulate(core, "never_rate"), data = pa))$r.squared)
}
pa31 <- prov_avg[, c("prov_en", vars)]
print(round(rbind(all_31 = between_stats(pa31),
without_Tibet = between_stats(pa31[pa31$prov_en != "Tibet", ]),
without_Hebei_Shandong = between_stats(pa31[!pa31$prov_en %in% c("Hebei", "Shandong"), ])), 3))
r_young_share r_sex_ratio r_higher_edu r_agri_share
all_31 0.386 0.401 -0.215 0.140
without_Tibet 0.472 0.309 -0.054 0.048
without_Hebei_Shandong 0.401 0.305 -0.198 0.169
R2_young_share R2_sex_ratio R2_higher_edu R2_agri_share
all_31 0.149 0.161 0.046 0.020
without_Tibet 0.223 0.095 0.003 0.002
without_Hebei_Shandong 0.161 0.093 0.039 0.028
R2_all4
all_31 0.438
without_Tibet 0.490
without_Hebei_Shandong 0.382
cat("\nCorrelation of province means when one province is dropped\n")
Correlation of province means when one province is dropped
for (v in core) {
rr <- sapply(pa31$prov_en, function(p) { d <- pa31[pa31$prov_en != p, ]; cor(d$never_rate, d[[v]]) })
cat(sprintf("%-11s full %.3f | min %.3f (without %s) | max %.3f (without %s)\n", v,
cor(pa31$never_rate, pa31[[v]]), min(rr), names(which.min(rr)), max(rr), names(which.max(rr))))
}
young_share full 0.386 | min 0.276 (without Henan) | max 0.472 (without Tibet)
sex_ratio full 0.401 | min 0.309 (without Tibet) | max 0.470 (without Henan)
higher_edu full -0.215 | min -0.251 (without Beijing) | max -0.054 (without Tibet)
agri_share full 0.140 | min 0.044 (without Henan) | max 0.223 (without Liaoning)
cat("\nSpearman correlation of province means\n")
Spearman correlation of province means
print(round(sapply(core, function(v) cor(pa31$never_rate, pa31[[v]], method = "spearman")), 3))
young_share sex_ratio higher_edu agri_share
0.428 0.345 -0.113 0.131
cat("\nMundlak model without Tibet\n")
Mundlak model without Tibet
mund_nt <- mund[mund$prov_en != "Tibet", ]
for (v in core) { mund_nt[[paste0(v, "_b")]] <- ave(mund_nt[[v]], mund_nt$prov_en)
mund_nt[[paste0(v, "_w")]] <- mund_nt[[v]] - mund_nt[[paste0(v, "_b")]] }
mm_nt <- lm(fm, data = mund_nt); Xn <- model.matrix(mm_nt); gn <- mund_nt$prov_en; Gn <- length(unique(gn))
Vn <- cl_vcov(Xn, residuals(mm_nt), gn, rep(1, nrow(Xn))); tn <- qt(0.975, Gn - 1)
print(round(t(sapply(core, function(v) {
jw <- which(colnames(Xn) == paste0(v, "_w")); jb <- which(colnames(Xn) == paste0(v, "_b"))
bw <- coef(mm_nt)[jw]; bb <- coef(mm_nt)[jb]; sw <- sqrt(Vn[jw, jw]); sb <- sqrt(Vn[jb, jb])
sdd <- sqrt(Vn[jb, jb] + Vn[jw, jw] - 2 * Vn[jb, jw])
c(within = unname(bw), w_low = unname(bw - tn * sw), w_high = unname(bw + tn * sw),
between = unname(bb), b_low = unname(bb - tn * sb), b_high = unname(bb + tn * sb),
p_between = unname(2 * pt(-abs(bb / sb), Gn - 1)), p_diff = unname(2 * pt(-abs((bb - bw) / sdd), Gn - 1)))
})), 3))
within w_low w_high between b_low b_high p_between p_diff
young_share 0.783 0.557 1.009 1.015 0.546 1.485 0.000 0.291
sex_ratio 0.076 0.011 0.141 0.143 0.030 0.256 0.015 0.307
higher_edu -0.069 -0.300 0.161 -0.356 -0.726 0.014 0.059 0.207
agri_share -0.139 -0.333 0.055 -0.448 -0.835 -0.060 0.025 0.195
7. How much variation can age structure explain?
7.1 Five age groups
grp_def <- list(a15_19 = "15-19", a20_24 = "20-24", a25_29 = "25-29",
a30_39 = c("30-34", "35-39"), a40_59 = c("40-44", "45-49", "50-54", "55-59"))
read_age_groups <- function(file_path, year) {
raw <- as.data.frame(read_excel(file_path, col_names = FALSE))
header <- as.character(unlist(raw[2, ])); d <- raw[6:nrow(raw), ]
num <- function(seg) { h <- which(grepl(seg, header, fixed = TRUE) & grepl("人数", header))
if (length(h)) as.numeric(d[[h[1]]]) else rep(0, nrow(d)) }
out <- data.frame(year = year, province = clean_prov(d[[9]]))
for (g in names(grp_def)) out[[g]] <- Reduce(`+`, lapply(grp_def[[g]], function(s) num(paste0(s, "岁"))))
out$a60plus <- num("60岁及以上")
out$a0_14 <- num("0-4岁") + num("5-9岁") + num("10-14岁")
out <- out[!is.na(out$a15_19), ]
keep_prov(out)
}
ag <- bind_rows(lapply(analysis_years, function(y) read_age_groups(file_for("年龄", y), y)))
raw14 <- as.data.frame(read_excel(file_for("年龄", 2014), col_names = FALSE))
nm14 <- gsub("～", "-", as.character(unlist(raw14[3, ]))); sub14 <- as.character(unlist(raw14[4, ]))
d14 <- raw14[5:nrow(raw14), ]
col14 <- function(seg) { h <- which(grepl(seg, nm14, fixed = TRUE) & sub14 == "人数"); as.numeric(d14[[h[1]]]) }
ag14 <- data.frame(year = 2014, province = clean_prov(d14[[1]]))
for (g in names(grp_def)) ag14[[g]] <- Reduce(`+`, lapply(grp_def[[g]], function(s) col14(s)))
ag14$a60plus <- 0; ag14$a0_14 <- 0
ag14 <- ag14[!(ag14$province %in% c("合计", "兵团")) & !is.na(ag14$a15_19), ]
ag <- rbind(ag, ag14)
ag$prov_en <- unname(name_map[ag$province]); ag <- ag[!is.na(ag$prov_en) & !duplicated(ag[, c("year", "prov_en")]), ]
grp <- names(grp_def)
ag$n15_59 <- rowSums(ag[, grp])
for (g in grp) ag[[paste0("s_", g)]] <- ag[[g]] / ag$n15_59
P18 <- merge(panel, ag[, c("year", "prov_en", grp, "a60plus", "a0_14", "n15_59", paste0("s_", grp))],
by = c("year", "prov_en"))
stopifnot(nrow(P18) == 248)
cat("Youth share rebuilt from age groups, max abs difference:",
round(max(abs(100 * (P18$s_a15_19 + P18$s_a20_24 + P18$s_a25_29) - P18$young_share)), 2), "\n")
Youth share rebuilt from age groups, max abs difference: 0.05
cat("\nDenominators by year\n")
Denominators by year
print(aggregate(cbind(n_marital, age_15_59 = n15_59, age_60plus = a60plus, age_0_14 = a0_14) ~ year, data = P18, FUN = sum))
year n_marital age_15_59 age_60plus age_0_14
1 2011 126000 126000 0 0
2 2012 154605 154605 0 0
3 2013 194795 194795 0 0
4 2014 196937 196938 0 0
5 2015 203998 199365 4633 0
6 2016 346352 337125 11777 91106
7 2017 348037 337710 15545 92131
8 2018 149311 146331 5426 0
7.2 Composition benchmark and Goodman regression
m_nr <- mean(P18$never_rate) / 100; m_ys <- mean(P18$young_share) / 100
cat(sprintf("Mean never-married share %.3f | mean youth share %.3f | implied r(15-29) if all 30+ married: %.3f\n",
m_nr, m_ys, m_nr / m_ys))
Mean never-married share 0.185 | mean youth share 0.352 | implied r(15-29) if all 30+ married: 0.526
# Goodman ecological regression, rates bounded in [0, 1]
bls <- function(X, y) {
f <- function(b) sum((y - X %*% b)^2); g <- function(b) -2 * drop(t(X) %*% (y - X %*% b))
optim(rep(0.3, ncol(X)), f, g, method = "L-BFGS-B", lower = 0, upper = 1)$par
}
Xs <- as.matrix(P18[, paste0("s_", grp)]); ynr <- P18$never_rate / 100
b_pool <- setNames(bls(Xs, ynr), grp)
b_ols <- setNames(coef(lm(ynr ~ 0 + Xs)), grp)
b_year <- t(sapply(sort(unique(P18$year)), function(y) { k <- P18$year == y; bls(Xs[k, ], ynr[k]) }))
dimnames(b_year) <- list(sort(unique(P18$year)), grp)
cat("\nImplied never-married rate by age group (248 province-years)\n")
Implied never-married rate by age group (248 province-years)
print(round(rbind(pooled_bounded = b_pool, pooled_unconstrained = b_ols), 3))
a15_19 a20_24 a25_29 a30_39 a40_59
pooled_bounded 0.725 1.000 0.000 0.000 0.098
pooled_unconstrained 0.269 1.469 -0.208 -0.062 0.179
cat("By year (bounded):\n"); print(round(b_year, 3))
By year (bounded):
a15_19 a20_24 a25_29 a30_39 a40_59
2011 0.997 0.859 0.000 0.000 0.082
2012 0.805 1.000 0.000 0.000 0.091
2013 1.000 0.975 0.000 0.000 0.109
2014 1.000 1.000 0.000 0.000 0.138
2015 0.952 0.831 0.000 0.000 0.125
2016 0.627 0.482 0.000 0.012 0.212
2017 0.348 1.000 0.000 0.044 0.090
2018 1.000 1.000 0.089 0.000 0.054
P18$expected <- 100 * drop(Xs %*% b_pool)
P18$expected_y <- 100 * sapply(seq_len(nrow(P18)), function(i) sum(Xs[i, ] * b_year[as.character(P18$year[i]), ]))
r2v <- function(y, f) 1 - sum((y - f)^2) / sum((y - mean(y))^2)
pm <- aggregate(cbind(never_rate, expected, expected_y) ~ prov_en, data = P18, FUN = mean)
cat(sprintf("R2 of expected vs actual: all province-years %.3f (pooled rates) / %.3f (year-specific rates)\n",
r2v(P18$never_rate, P18$expected), r2v(P18$never_rate, P18$expected_y)))
R2 of expected vs actual: all province-years 0.602 (pooled rates) / 0.648 (year-specific rates)
cat(sprintf("Between provinces (31 means): R2 = %.3f / %.3f; correlation = %.3f / %.3f\n",
r2v(pm$never_rate, pm$expected), r2v(pm$never_rate, pm$expected_y),
cor(pm$never_rate, pm$expected), cor(pm$never_rate, pm$expected_y)))
Between provinces (31 means): R2 = 0.492 / 0.511; correlation = 0.738 / 0.753
pm$gap <- round(pm$never_rate - pm$expected_y, 1); pm <- pm[order(pm$gap), ]
cat("Actual minus age-expected (year-specific rates), lowest 5:\n"); print(head(pm[, c("prov_en", "gap")], 5))
Actual minus age-expected (year-specific rates), lowest 5:
prov_en gap
1 Anhui -6.1
13 Hubei -5.1
16 Jiangsu -3.7
20 Ningxia -3.7
21 Qinghai -3.6
cat("Highest 5:\n"); print(tail(pm[, c("prov_en", "gap")], 5))
Highest 5:
prov_en gap
7 Guangxi 2.8
10 Hebei 3.3
12 Henan 5.0
19 Liaoning 6.4
28 Tibet 12.4
pma <- aggregate(P18[, c("never_rate", paste0("s_", grp))], by = list(prov_en = P18$prov_en), FUN = mean)
fa <- lm(reformulate(paste0("s_", grp[-5]), "never_rate"), data = pma)
cat(sprintf("Between provinces, OLS on four age-group shares: R2 = %.3f, adjusted R2 = %.3f (youth share alone: R2 = %.3f)\n",
summary(fa)$r.squared, summary(fa)$adj.r.squared,
summary(lm(never_rate ~ young_share, data = aggregate(P18[, c("never_rate", "young_share")], by = list(p = P18$prov_en), FUN = mean)))$r.squared))
Between provinces, OLS on four age-group shares: R2 = 0.688, adjusted R2 = 0.640 (youth share alone: R2 = 0.149)
P18$young1519 <- P18$a15_19 / (P18$a15_19 + P18$a20_24 + P18$a25_29)
cat(sprintf("Correlation of youth share with the 15-19 share among ages 15-29: pooled %.3f | province means %.3f | within %.3f\n",
cor(P18$young_share, P18$young1519),
with(aggregate(cbind(young_share, young1519) ~ prov_en, P18, mean), cor(young_share, young1519)),
cor(residuals(lm(young_share ~ factor(year) + factor(prov_en), P18)), residuals(lm(young1519 ~ factor(year) + factor(prov_en), P18)))))
Correlation of youth share with the 15-19 share among ages 15-29: pooled 0.200 | province means 0.106 | within 0.169
7.3 Variance explained by the youth share
var_parts <- function(d) {
tot <- sum((d$never_rate - mean(d$never_rate))^2)
bp <- sum((ave(d$never_rate, d$prov_en) - mean(d$never_rate))^2); by <- sum((ave(d$never_rate, d$year) - mean(d$never_rate))^2)
c(province = bp / tot, year = by / tot, remainder = 1 - (bp + by) / tot)
}
print(round(rbind(N248 = var_parts(panel), N217 = var_parts(reg)), 3))
province year remainder
N248 0.593 0.245 0.162
N217 0.570 0.261 0.169
vp <- var_parts(reg)
pm217 <- aggregate(reg[, vars], by = list(prov_en = reg$prov_en), FUN = mean)
ym217 <- aggregate(reg[, vars], by = list(year = reg$year), FUN = mean)
r2_b <- summary(lm(never_rate ~ young_share, pm217))$r.squared
r2_y <- summary(lm(never_rate ~ young_share, ym217))$r.squared
r2_w <- summary(lm(dm(reg$never_rate) ~ dm(reg$young_share)))$r.squared
cat(sprintf("Youth share, R2 within each component: between provinces %.3f | between years %.3f | within %.3f\n", r2_b, r2_y, r2_w))
Youth share, R2 within each component: between provinces 0.141 | between years 0.919 | within 0.451
cat(sprintf("As a share of total variance: between provinces %.3f | between years %.3f | within %.3f\n",
r2_b * vp["province"], r2_y * vp["year"], r2_w * vp["remainder"]))
As a share of total variance: between provinces 0.080 | between years 0.240 | within 0.076
cat("\nProvince means from the 217-observation sample\n")
Province means from the 217-observation sample
shapley <- function(d) sapply(core, function(v) { others <- setdiff(core, v); s <- 0
for (k in 0:3) for (S in if (k == 0) list(character(0)) else combn(others, k, simplify = FALSE)) {
w <- factorial(length(S)) * factorial(3 - length(S)) / factorial(4)
s <- s + w * (summary(lm(reformulate(c(S, v), "never_rate"), d))$r.squared -
if (length(S)) summary(lm(reformulate(S, "never_rate"), d))$r.squared else 0) }
s })
for (lab in c("all 31", "without Tibet")) {
d <- if (lab == "all 31") pm217 else pm217[pm217$prov_en != "Tibet", ]
cat("\n[", lab, "]\n")
print(round(rbind(r = sapply(core, function(v) cor(d$never_rate, d[[v]])),
R2 = sapply(core, function(v) summary(lm(reformulate(v, "never_rate"), d))$r.squared),
Shapley = shapley(d),
VIF = sapply(core, function(v) 1 / (1 - summary(lm(reformulate(setdiff(core, v), v), d))$r.squared))), 3))
cat("R2 of all four:", round(summary(lm(reformulate(core, "never_rate"), d))$r.squared, 3),
"| correlation tertiary vs agricultural hukou:", round(cor(d$higher_edu, d$agri_share), 3), "\n")
}
[ all 31 ]
young_share sex_ratio higher_edu agri_share
r 0.375 0.409 -0.221 0.164
R2 0.141 0.167 0.049 0.027
Shapley 0.179 0.139 0.048 0.041
VIF 2.052 1.425 3.388 3.664
R2 of all four: 0.408 | correlation tertiary vs agricultural hukou: -0.654
[ without Tibet ]
young_share sex_ratio higher_edu agri_share
r 0.471 0.318 -0.063 0.075
R2 0.221 0.101 0.004 0.006
Shapley 0.280 0.101 0.027 0.058
VIF 2.051 1.337 3.073 3.529
R2 of all four: 0.466 | correlation tertiary vs agricultural hukou: -0.637
7.4 Cross-validation with age groups
reg18 <- merge(reg, P18[, c("year", "prov_en", paste0("s_", grp))], by = c("year", "prov_en"))
reg18 <- reg18[order(match(reg18$prov_en, unique(reg$prov_en))), ]
cv2 <- list("Year means only" = never_rate ~ factor(year),
"Year + youth share" = never_rate ~ factor(year) + young_share,
"Year + youth share + sex ratio" = never_rate ~ factor(year) + young_share + sex_ratio,
"Year + four features" = never_rate ~ factor(year) + young_share + sex_ratio + higher_edu + agri_share,
"Year + age groups" = never_rate ~ factor(year) + s_a15_19 + s_a20_24 + s_a25_29 + s_a30_39,
"Year + age groups + sex ratio" = never_rate ~ factor(year) + s_a15_19 + s_a20_24 + s_a25_29 + s_a30_39 + sex_ratio)
po <- unlist(lapply(unique(reg18$prov_en), function(p) rep(p, sum(reg18$prov_en == p))))
e2 <- lapply(cv2, function(f) unlist(lapply(unique(reg18$prov_en), function(p) {
tr <- reg18[reg18$prov_en != p, ]; te <- reg18[reg18$prov_en == p, ]; te$never_rate - predict(lm(f, tr), te) })))
b1 <- mean(e2[[1]]^2); b1n <- mean(e2[[1]][po != "Tibet"]^2)
print(round(t(sapply(e2, function(e) c(RMSE = sqrt(mean(e^2)), R2_oos = 1 - mean(e^2) / b1,
RMSE_without_Tibet = sqrt(mean(e[po != "Tibet"]^2)), R2_oos_without_Tibet = 1 - mean(e[po != "Tibet"]^2) / b1n))), 3))
RMSE R2_oos RMSE_without_Tibet
Year means only 5.707 0.000 5.104
Year + youth share 5.164 0.181 4.444
Year + youth share + sex ratio 4.882 0.268 4.293
Year + four features 4.792 0.295 4.168
Year + age groups 3.888 0.536 3.178
Year + age groups + sex ratio 3.729 0.573 3.096
R2_oos_without_Tibet
Year means only 0.000
Year + youth share 0.242
Year + youth share + sex ratio 0.292
Year + four features 0.333
Year + age groups 0.612
Year + age groups + sex ratio 0.632
7.5 Other characteristics after controlling for age groups
agevars <- paste0("p_", grp[-5])
for (i in seq_along(agevars)) { reg18[[agevars[i]]] <- 100 * reg18[[paste0("s_", grp[i])]] }
others <- c("sex_ratio", "higher_edu", "agri_share")
G18 <- length(unique(reg18$prov_en)); t18 <- qt(0.975, G18 - 1)
ci_rows <- function(m, keep, g) {
X <- model.matrix(m); V <- cl_vcov(X, residuals(m), g, rep(1, nrow(X)))
t(sapply(keep, function(v) { j <- which(colnames(X) == v); b <- coef(m)[j]; s <- sqrt(V[j, j])
c(Estimate = unname(b), CI_low = unname(b - t18 * s), CI_high = unname(b + t18 * s), p = unname(2 * pt(-abs(b / s), G18 - 1))) }))
}
m3a <- lm(reformulate(c(agevars, others, "factor(year)", "factor(prov_en)"), "never_rate"), data = reg18)
cat("Model 3 with age groups, N =", nrow(reg18), "| within R2 =",
round(summary(lm(residuals(lm(never_rate ~ factor(year) + factor(prov_en), reg18)) ~
sapply(c(agevars, others), function(v) residuals(lm(reformulate(c("factor(year)", "factor(prov_en)"), v), reg18)))))$r.squared, 3), "\n")
Model 3 with age groups, N = 217 | within R2 = 0.719
print(round(ci_rows(m3a, c(agevars, others), reg18$prov_en), 3))
Estimate CI_low CI_high p
p_a15_19 1.292 0.818 1.766 0.000
p_a20_24 1.126 0.802 1.450 0.000
p_a25_29 0.190 -0.250 0.630 0.384
p_a30_39 0.224 -0.088 0.536 0.153
sex_ratio 0.073 0.017 0.129 0.012
higher_edu -0.029 -0.274 0.216 0.812
agri_share -0.034 -0.184 0.116 0.648
mu <- reg18
for (v in c(agevars, others)) { mu[[paste0(v, "_b")]] <- ave(mu[[v]], mu$prov_en); mu[[paste0(v, "_w")]] <- mu[[v]] - mu[[paste0(v, "_b")]] }
mma <- lm(reformulate(c(paste0(c(agevars, others), "_w"), paste0(c(agevars, others), "_b"), "factor(year)"), "never_rate"), data = mu)
cat("\nMundlak with age groups: within and between coefficients\n")
Mundlak with age groups: within and between coefficients
print(round(rbind(within = ci_rows(mma, paste0(others, "_w"), mu$prov_en), between = ci_rows(mma, paste0(others, "_b"), mu$prov_en)), 3))
Estimate CI_low CI_high p
sex_ratio_w 0.073 0.021 0.126 0.008
higher_edu_w -0.029 -0.259 0.202 0.801
agri_share_w -0.034 -0.175 0.107 0.627
sex_ratio_b 0.099 -0.004 0.201 0.058
higher_edu_b -0.212 -0.688 0.265 0.372
agri_share_b -0.093 -0.323 0.138 0.418
pc_age <- function(d, v) cor(residuals(lm(reformulate(agevars, "never_rate"), d)), residuals(lm(reformulate(agevars, v), d)))
pm18 <- aggregate(reg18[, c("never_rate", agevars, others)], by = list(prov_en = reg18$prov_en), FUN = mean)
cat("\nPartial correlations controlling for the age-group shares\n")
Partial correlations controlling for the age-group shares
print(round(rbind(pooled_217 = sapply(others, function(v) pc_age(reg18, v)),
province_means = sapply(others, function(v) pc_age(pm18, v)),
province_means_without_Tibet = sapply(others, function(v) pc_age(pm18[pm18$prov_en != "Tibet", ], v))), 3))
sex_ratio higher_edu agri_share
pooled_217 0.313 0.039 -0.053
province_means 0.332 -0.222 0.100
province_means_without_Tibet 0.190 0.177 -0.271
cat("Between provinces: R2 age groups =", round(summary(lm(reformulate(agevars, "never_rate"), pm18))$r.squared, 3),
"| + sex ratio =", round(summary(lm(reformulate(c(agevars, "sex_ratio"), "never_rate"), pm18))$r.squared, 3),
"| + three features =", round(summary(lm(reformulate(c(agevars, others), "never_rate"), pm18))$r.squared, 3), "\n")
Between provinces: R2 age groups = 0.703 | + sex ratio = 0.736 | + three features = 0.744
cat("Without Tibet: R2 age groups =", round(summary(lm(reformulate(agevars, "never_rate"), pm18[pm18$prov_en != "Tibet", ]))$r.squared, 3), "\n")
Without Tibet: R2 age groups = 0.781
7.6 Monotone Goodman regression
Cum <- t(apply(Xs, 1, cumsum)); colnames(Cum) <- paste0("C", 1:5)
fit_mono <- function(C, y) { # r_a = sum_{k >= a} d_k, d_k >= 0, sum(d) <= 1
f <- function(d) sum((y - C %*% d)^2); g <- function(d) -2 * drop(t(C) %*% (y - C %*% d))
d <- constrOptim(rep(0.05, 5), f, g, ui = rbind(diag(5), -rep(1, 5)), ci = c(rep(0, 5), -1) - 1e-10)$par
rev(cumsum(rev(d)))
}
fit_box <- function(X, y) bls(X, y)
fit_free <- function(X, y) unname(coef(lm(y ~ 0 + X)))
rates <- rbind(unconstrained = fit_free(Xs, ynr), bounded_0_1 = fit_box(Xs, ynr), monotone = fit_mono(Cum, ynr))
colnames(rates) <- grp; print(round(rates, 3))
a15_19 a20_24 a25_29 a30_39 a40_59
unconstrained 0.269 1.469 -0.208 -0.062 0.179
bounded_0_1 0.725 1.000 0.000 0.000 0.098
monotone 0.998 0.998 0.020 0.020 0.020
rates_y <- t(sapply(sort(unique(P18$year)), function(y) { k <- P18$year == y; fit_mono(Cum[k, ], ynr[k]) }))
dimnames(rates_y) <- list(sort(unique(P18$year)), grp)
cat("Monotone, by year:\n"); print(round(rates_y, 3))
Monotone, by year:
a15_19 a20_24 a25_29 a30_39 a40_59
2011 0.999 0.988 0.000 0.000 0.000
2012 1.000 1.000 0.014 0.014 0.014
2013 1.000 1.000 0.030 0.030 0.030
2014 1.000 1.000 0.045 0.045 0.045
2015 1.000 1.000 0.017 0.017 0.017
2016 0.695 0.661 0.068 0.068 0.068
2017 0.768 0.768 0.050 0.050 0.050
2018 1.000 1.000 0.134 0.017 0.017
fit_stats <- function(fun, Xmat, yearly) {
pred <- numeric(nrow(P18))
if (yearly) { for (y in unique(P18$year)) { k <- P18$year == y; pred[k] <- Xs[k, ] %*% fun(Xmat[k, ], ynr[k]) } }
else pred <- drop(Xs %*% fun(Xmat, ynr))
pm_ <- aggregate(cbind(a = ynr, p = pred) ~ P18$prov_en, FUN = mean)
cvp <- numeric(nrow(P18))
for (pr in unique(P18$prov_en)) {
te <- P18$prov_en == pr
if (yearly) { for (y in unique(P18$year)) { k <- P18$year == y
tr <- k & !te; tt <- k & te; cvp[tt] <- Xs[tt, , drop = FALSE] %*% fun(Xmat[tr, ], ynr[tr]) } }
else cvp[te] <- Xs[te, , drop = FALSE] %*% fun(Xmat[!te, ], ynr[!te])
}
c(R2_all = r2v(ynr, pred), R2_province_means = r2v(pm_$a, pm_$p),
RMSE_LOPO_pp = 100 * sqrt(mean((ynr - cvp)^2)))
}
fs <- rbind(
unconstrained_pooled = fit_stats(fit_free, Xs, FALSE), bounded_pooled = fit_stats(fit_box, Xs, FALSE), monotone_pooled = fit_stats(fit_mono, Cum, FALSE),
unconstrained_yearly = fit_stats(fit_free, Xs, TRUE), bounded_yearly = fit_stats(fit_box, Xs, TRUE), monotone_yearly = fit_stats(fit_mono, Cum, TRUE))
print(round(fs, 3))
R2_all R2_province_means RMSE_LOPO_pp
unconstrained_pooled 0.648 0.602 4.209
bounded_pooled 0.602 0.492 4.335
monotone_pooled 0.577 0.452 4.466
unconstrained_yearly 0.723 0.662 4.127
bounded_yearly 0.648 0.511 4.237
monotone_yearly 0.623 0.460 4.261
cat("LOPO RMSE of 'Year means only' =", round(sqrt(mean(e2[[1]]^2)), 3), "| 'Year + age groups' =", round(sqrt(mean(e2[["Year + age groups"]]^2)), 3), "\n")
LOPO RMSE of 'Year means only' = 5.707 | 'Year + age groups' = 3.888
7.7 Age composition within the youth group
P18$y2024 <- P18$a20_24 / (P18$a15_19 + P18$a20_24 + P18$a25_29)
P18$y2529 <- P18$a25_29 / (P18$a15_19 + P18$a20_24 + P18$a25_29)
dm2 <- function(v) residuals(lm(v ~ factor(year) + factor(prov_en), P18))
pmy <- aggregate(cbind(young_share, young1519, y2024, y2529, never_rate) ~ prov_en, P18, mean)
print(round(rbind(
corr_with_youth_share = c(pooled = sapply(c("young1519", "y2024", "y2529"), function(v) cor(P18$young_share, P18[[v]])),
province_means = sapply(c("young1519", "y2024", "y2529"), function(v) cor(pmy$young_share, pmy[[v]])),
within = sapply(c("young1519", "y2024", "y2529"), function(v) cor(dm2(P18$young_share), dm2(P18[[v]])))),
corr_with_never_rate = c(pooled = sapply(c("young1519", "y2024", "y2529"), function(v) cor(P18$never_rate, P18[[v]])),
province_means = sapply(c("young1519", "y2024", "y2529"), function(v) cor(pmy$never_rate, pmy[[v]])),
within = sapply(c("young1519", "y2024", "y2529"), function(v) cor(dm2(P18$never_rate), dm2(P18[[v]]))))), 3))
pooled.young1519 pooled.y2024 pooled.y2529
corr_with_youth_share 0.200 0.693 -0.591
corr_with_never_rate 0.317 0.761 -0.711
province_means.young1519 province_means.y2024
corr_with_youth_share 0.106 0.465
corr_with_never_rate 0.371 0.817
province_means.y2529 within.young1519 within.y2024
corr_with_youth_share -0.361 0.169 0.515
corr_with_never_rate -0.732 0.486 0.504
within.y2529
corr_with_youth_share -0.478
corr_with_never_rate -0.684
k15 <- P18$year <= 2015
cat(sprintf("2011-2015: mean never-married share %.3f | mean youth share %.3f | implied r(15-29) if all 30+ married %.3f\n",
mean(P18$never_rate[k15]) / 100, mean(P18$young_share[k15]) / 100, mean(P18$never_rate[k15]) / mean(P18$young_share[k15])))
2011-2015: mean never-married share 0.207 | mean youth share 0.384 | implied r(15-29) if all 30+ married 0.540
7.8 Indirect standardisation with the 2010 census
Census counts: Tabulation on the 2010 Population Census of
China, Table 5-3a (city), long form, both sexes.
census <- data.frame(
age = c("15-19", "20-24", "25-29", "30-34", "35-39", "40-44", "45-49", "50-54", "55-59"),
total = c(2728757, 4276600, 3484422, 3299594, 3713982, 3541583, 3058752, 2311229, 2158466),
never = c(2712072, 3633319, 1285204, 341651, 150538, 75994, 43230, 24912, 18518))
cen_grp <- list(a15_19 = 1, a20_24 = 2, a25_29 = 3, a30_39 = 4:5, a40_59 = 6:9)
R_cen <- sapply(cen_grp, function(i) sum(census$never[i]) / sum(census$total[i]))
cat("Never-married rate by age group (city)\n"); print(round(R_cen, 3))
Never-married rate by age group (city)
a15_19 a20_24 a25_29 a30_39 a40_59
0.994 0.850 0.369 0.070 0.015
cat("Composition bound for each age-group coefficient:\n"); print(round(R_cen[1:4] - R_cen[5], 3))
Composition bound for each age-group coefficient:
a15_19 a20_24 a25_29 a30_39
0.979 0.835 0.354 0.055
P18$exp_census <- 100 * drop(Xs %*% R_cen)
cat(sprintf("Mean actual %.2f | mean census-expected %.2f\n", mean(P18$never_rate), mean(P18$exp_census)))
Mean actual 18.55 | mean census-expected 24.65
pmc <- aggregate(cbind(never_rate, exp_census) ~ prov_en, P18, mean)
cat(sprintf("Correlation actual vs expected: all %.3f | province means %.3f (R2 %.3f) | within %.3f\n",
cor(P18$never_rate, P18$exp_census), cor(pmc$never_rate, pmc$exp_census), cor(pmc$never_rate, pmc$exp_census)^2,
cor(dm2(P18$never_rate), dm2(P18$exp_census))))
Correlation actual vs expected: all 0.703 | province means 0.545 (R2 0.297) | within 0.798
cat(sprintf("Without Tibet, province means: R2 %.3f\n", with(pmc[pmc$prov_en != "Tibet", ], cor(never_rate, exp_census)^2)))
Without Tibet, province means: R2 0.400
rc <- merge(reg, P18[, c("year", "prov_en", "exp_census")], by = c("year", "prov_en"))
gc <- rc$prov_en; Gc <- length(unique(gc)); tc <- qt(0.975, Gc - 1)
slope_ci <- function(m, v) { X <- model.matrix(m); V <- cl_vcov(X, residuals(m), gc, rep(1, nrow(X))); j <- which(colnames(X) == v)
b <- coef(m)[j]; s <- sqrt(V[j, j]); c(slope = unname(b), CI_low = unname(b - tc * s), CI_high = unname(b + tc * s)) }
rc_b <- aggregate(cbind(never_rate, exp_census) ~ prov_en, rc, mean)
cat("Slope of actual on census-expected (N = 217):\n")
Slope of actual on census-expected (N = 217):
print(round(rbind(two_way_FE = slope_ci(lm(never_rate ~ exp_census + factor(year) + factor(prov_en), rc), "exp_census"),
year_FE = slope_ci(lm(never_rate ~ exp_census + factor(year), rc), "exp_census"),
between_31_means = c(slope = unname(coef(lm(never_rate ~ exp_census, rc_b))[2]),
confint(lm(never_rate ~ exp_census, rc_b))[2, ]),
two_way_FE_plus_sex_ratio = slope_ci(lm(never_rate ~ exp_census + sex_ratio + factor(year) + factor(prov_en), rc), "exp_census")), 3))
slope CI_low CI_high
two_way_FE 1.146 0.917 1.375
year_FE 1.052 0.707 1.397
between_31_means 1.003 0.389 1.618
two_way_FE_plus_sex_ratio 1.209 1.001 1.416
cat("Sex ratio controlling for the census-expected share (two-way FE):\n")
Sex ratio controlling for the census-expected share (two-way FE):
print(round(slope_ci(lm(never_rate ~ exp_census + sex_ratio + factor(year) + factor(prov_en), rc), "sex_ratio"), 3))
slope CI_low CI_high
0.078 0.021 0.134
pmc$gap <- round(pmc$never_rate - pmc$exp_census - mean(pmc$never_rate - pmc$exp_census), 1); pmc <- pmc[order(pmc$gap), ]
cat("Actual minus census-expected, centred (province means), lowest 5:\n"); print(head(pmc[, c("prov_en", "gap")], 5))
Actual minus census-expected, centred (province means), lowest 5:
prov_en gap
1 Anhui -8.2
13 Hubei -5.8
23 Shandong -5.7
20 Ningxia -5.0
16 Jiangsu -4.5
cat("Highest 5:\n"); print(tail(pmc[, c("prov_en", "gap")], 5))
Highest 5:
prov_en gap
10 Hebei 3.7
18 Jilin 4.0
12 Henan 4.7
19 Liaoning 7.7
28 Tibet 14.5
cat("Spearman correlation of province means (217 sample):\n")
Spearman correlation of province means (217 sample):
print(round(sapply(core, function(v) cor(pm217$never_rate, pm217[[v]], method = "spearman")), 3))
young_share sex_ratio higher_edu agri_share
0.430 0.390 -0.161 0.204
fa217 <- lm(reformulate(agevars, "never_rate"), data = pm18)
cat(sprintf("Between provinces, four age-group shares: R2 = %.3f, adjusted R2 = %.3f\n",
summary(fa217)$r.squared, summary(fa217)$adj.r.squared))
Between provinces, four age-group shares: R2 = 0.703, adjusted R2 = 0.657
cat(sprintf("Five age groups, share of total variance via between provinces: %.3f\n", summary(fa217)$r.squared * vp["province"]))
Five age groups, share of total variance via between provinces: 0.401
wR2 <- function(xs) summary(lm(residuals(lm(never_rate ~ factor(year) + factor(prov_en), reg18)) ~
sapply(xs, function(v) residuals(lm(reformulate(c("factor(year)", "factor(prov_en)"), v), reg18)))))$r.squared
cat(sprintf("Within provinces R2: age groups only %.3f | age groups + sex ratio %.3f\n", wR2(agevars), wR2(c(agevars, "sex_ratio"))))
Within provinces R2: age groups only 0.677 | age groups + sex ratio 0.719
cat(sprintf("Five age groups, share of total variance via within provinces: %.3f\n", wR2(agevars) * vp["remainder"]))
Five age groups, share of total variance via within provinces: 0.114
cat("MAE of the cross-validation specifications:\n")
MAE of the cross-validation specifications:
print(round(sapply(e2, function(e) mean(abs(e))), 3))
Year means only Year + youth share
4.154 3.796
Year + youth share + sex ratio Year + four features
3.754 3.646
Year + age groups Year + age groups + sex ratio
2.816 2.735
7.9 Four census standards
Tables 5-3 (national), 5-3a (city), 5-3b (town) and 5-3c (rural).
cen <- list(
national = list(total = c(9190962, 11604416, 9278248, 8999299, 11099483, 11881111, 10159832, 7658776, 8011414),
never = c(9069370, 8698138, 2682763, 811360, 458064, 293283, 182178, 136332, 147629)),
city = list(total = census$total, never = census$never),
town = list(total = c(2021480, 2067283, 1754241, 1835721, 2315452, 2418522, 1990457, 1409108, 1434409),
never = c(2002917, 1511207, 429858, 116887, 62132, 38454, 22856, 16602, 18349)),
rural = list(total = c(4440725, 5260533, 4039585, 3863984, 5070049, 5921006, 5110623, 3938439, 4418539),
never = c(4354381, 3553612, 967701, 352822, 245394, 178835, 116092, 94818, 110762)))
chk <- with(cen, c(total = all(national$total == city$total + town$total + rural$total),
never = all(national$never == city$never + town$never + rural$never)))
cat("Check: national = city + town + rural:", chk, "\n")
Check: national = city + town + rural: TRUE TRUE
rate_of <- function(x) sapply(cen_grp, function(i) sum(x$never[i]) / sum(x$total[i]))
R_all <- t(sapply(cen, rate_of))
print(round(R_all, 3))
a15_19 a20_24 a25_29 a30_39 a40_59
national 0.987 0.750 0.289 0.063 0.020
city 0.994 0.850 0.369 0.070 0.015
town 0.991 0.731 0.245 0.043 0.013
rural 0.981 0.676 0.240 0.067 0.026
lopo_r2 <- function(f, d) {
e <- unlist(lapply(unique(d$prov_en), function(p) { tr <- d[d$prov_en != p, ]; te <- d[d$prov_en == p, ]
te$never_rate - predict(lm(f, tr), te) }))
e0 <- unlist(lapply(unique(d$prov_en), function(p) { tr <- d[d$prov_en != p, ]; te <- d[d$prov_en == p, ]
te$never_rate - predict(lm(never_rate ~ factor(year), tr), te) }))
c(R2_oos = 1 - mean(e^2) / mean(e0^2), RMSE = sqrt(mean(e^2)))
}
std_tab <- t(sapply(names(cen), function(nm) {
d <- merge(reg, P18[, c("year", "prov_en", paste0("s_", grp))], by = c("year", "prov_en"))
d$expv <- 100 * drop(as.matrix(d[, paste0("s_", grp)]) %*% R_all[nm, ])
gc <- d$prov_en
sl <- function(m) { X <- model.matrix(m); V <- cl_vcov(X, residuals(m), gc, rep(1, nrow(X))); j <- which(colnames(X) == "expv")
b <- coef(m)[j]; s <- sqrt(V[j, j]); c(unname(b), unname(b - tc * s), unname(b + tc * s)) }
w <- sl(lm(never_rate ~ expv + factor(year) + factor(prov_en), d))
ws <- sl(lm(never_rate ~ expv + sex_ratio + factor(year) + factor(prov_en), d))
db <- aggregate(cbind(never_rate, expv) ~ prov_en, d, mean); mb <- lm(never_rate ~ expv, db)
c(mean_expected = mean(d$expv), within_slope = w[1], within_low = w[2], within_high = w[3],
within_slope_sr = ws[1], within_sr_low = ws[2], within_sr_high = ws[3],
between_slope = unname(coef(mb)[2]), between_low = confint(mb)[2, 1], between_high = confint(mb)[2, 2],
R2_between = cor(db$never_rate, db$expv)^2,
R2_between_noTibet = with(db[db$prov_en != "Tibet", ], cor(never_rate, expv)^2),
R2_within = cor(residuals(lm(never_rate ~ factor(year) + factor(prov_en), d)), residuals(lm(expv ~ factor(year) + factor(prov_en), d)))^2,
lopo_offset = lopo_r2(never_rate ~ factor(year) + offset(expv), d),
lopo_slope = lopo_r2(never_rate ~ factor(year) + expv, d))
}))
cat("\nlopo_offset: slope fixed at 1; lopo_slope: one slope estimated\n")
lopo_offset: slope fixed at 1; lopo_slope: one slope estimated
print(round(t(std_tab), 3))
national city town rural
mean_expected 21.713 24.413 19.835 20.231
within_slope 1.264 1.146 1.262 1.370
within_low 1.025 0.917 1.030 1.120
within_high 1.503 1.375 1.494 1.621
within_slope_sr 1.324 1.209 1.320 1.428
within_sr_low 1.105 1.001 1.105 1.195
within_sr_high 1.543 1.416 1.534 1.660
between_slope 1.158 1.003 1.210 1.279
between_low 0.489 0.389 0.548 0.559
between_high 1.826 1.618 1.872 1.999
R2_between 0.302 0.278 0.325 0.313
R2_between_noTibet 0.423 0.394 0.448 0.437
R2_within 0.651 0.635 0.656 0.661
lopo_offset.R2_oos 0.371 0.357 0.387 0.370
lopo_offset.RMSE 4.526 4.577 4.469 4.531
lopo_slope.R2_oos 0.361 0.337 0.381 0.372
lopo_slope.RMSE 4.561 4.647 4.490 4.522
cat("Reference: LOPO R2_oos of 'Year + youth share' =", round(1 - mean(e2[["Year + youth share"]]^2) / mean(e2[[1]]^2), 3),
"| 'Year + age groups' =", round(1 - mean(e2[["Year + age groups"]]^2) / mean(e2[[1]]^2), 3), "\n")
Reference: LOPO R2_oos of 'Year + youth share' = 0.181 | 'Year + age groups' = 0.536
7.10 What age groups cannot explain
cat(sprintf("Unexplained between-province variance: %.3f | unexplained out-of-sample variance: %.3f\n",
1 - summary(fa217)$r.squared, 1 - (1 - mean(e2[["Year + age groups"]]^2) / mean(e2[[1]]^2))))
Unexplained between-province variance: 0.297 | unexplained out-of-sample variance: 0.464
ea <- e2[["Year + age groups"]]
pres <- data.frame(prov_en = unique(po), mean_error = tapply(ea, po, mean)[unique(po)],
same_sign_years = tapply(ea, po, function(e) mean(sign(e) == sign(mean(e))))[unique(po)])
pres <- pres[order(pres$mean_error), ]; pres$mean_error <- round(pres$mean_error, 1); pres$same_sign_years <- round(pres$same_sign_years, 2)
cat("Held-out prediction error (actual - predicted), lowest 5\n"); print(head(pres, 5), row.names = FALSE)
Held-out prediction error (actual - predicted), lowest 5
prov_en mean_error same_sign_years
Hubei -4.3 1.00
Zhejiang -4.0 1.00
Jiangsu -3.8 0.86
Qinghai -3.8 1.00
Anhui -3.6 1.00
cat("Highest 5\n"); print(tail(pres, 5), row.names = FALSE)
Highest 5
prov_en mean_error same_sign_years
Hebei 3.1 1.00
Guangxi 4.0 1.00
Henan 4.2 1.00
Liaoning 4.8 0.86
Tibet 12.0 1.00
cat("Provinces whose 7 yearly errors all have the same sign:", sum(pres$same_sign_years == 1), "of 31\n")
Provinces whose 7 yearly errors all have the same sign: 14 of 31
8. Province types
Input: province means of within-year z-scores of the four
characteristics (the never-married share is excluded). The number of
types is the K in 2–5 with the largest average silhouette width.
prov_z <- aggregate(zpanel[vars], by = list(prov_en = zpanel$prov_en), FUN = mean, na.rm = TRUE)
stopifnot(identical(prov_z$prov_en, prov_avg$prov_en))
z4 <- scale(prov_z[core]); z5 <- scale(prov_z[vars])
relabel <- function(cl) { o <- as.integer(names(sort(table(cl), decreasing = TRUE))); match(cl, o) } # 1 = largest
km_run <- function(z, k, seed = 42) { set.seed(seed); relabel(kmeans(z, centers = k, nstart = 100)$cluster) }
sil_width <- function(cl) mean(silhouette(cl, dist(z4))[, "sil_width"])
ksel <- t(sapply(1:6, function(k) {
set.seed(42); km <- kmeans(z4, centers = k, nstart = 100)
c(k = k, within_SS = km$tot.withinss, explained = 1 - km$tot.withinss / km$totss,
silhouette = if (k == 1) NA else sil_width(km$cluster))
}))
print(round(ksel, 3))
k within_SS explained silhouette
[1,] 1 120.000 0.000 NA
[2,] 2 75.780 0.368 0.304
[3,] 3 53.877 0.551 0.308
[4,] 4 40.045 0.666 0.315
[5,] 5 33.111 0.724 0.295
[6,] 6 28.167 0.765 0.295
K_main <- (2:5)[which.max(ksel[2:5, "silhouette"])]
cat("K =", K_main, "\n")
K = 4
par(mfrow = c(1, 2), mar = c(4.5, 4.5, 2, 1))
plot(ksel[, "k"], ksel[, "within_SS"], type = "b", pch = 19, xlab = "Number of clusters (k)",
ylab = "Within-cluster sum of squares", main = "Elbow plot")
plot(ksel[-1, "k"], ksel[-1, "silhouette"], type = "b", pch = 19, xlab = "Number of clusters (k)",
ylab = "Average silhouette width", main = "Silhouette")
ari <- function(a, b) {
tab <- table(a, b); s <- sum(choose(tab, 2)); sa <- sum(choose(rowSums(tab), 2))
sb <- sum(choose(colSums(tab), 2)); e <- sa * sb / choose(length(a), 2)
(s - e) / ((sa + sb) / 2 - e)
}
describe_partition <- function(main, k, label) {
cat("\n", label, " (K = ", k, ")\n", sep = "")
for (t in 1:k) cat("Type", t, "(", sum(main == t), "):", prov_avg$prov_en[main == t], "\n")
print(round(aggregate(prov_avg[vars], by = list(type = main), FUN = mean), 1))
alts <- list(km5 = km_run(z5, k), ward = relabel(cutree(hclust(dist(z4), method = "ward.D2"), k = k)),
pam = relabel(pam(z4, k = k)$clustering))
cat("Adjusted Rand index: K-means with outcome =", round(ari(main, alts$km5), 3),
"| Ward =", round(ari(main, alts$ward), 3), "| PAM =", round(ari(main, alts$pam), 3), "\n")
for (alt in names(alts)) {
tab <- table(main, alts[[alt]]); best <- apply(tab, 2, which.max)
moved <- prov_avg$prov_en[best[as.character(alts[[alt]])] != main]
cat("Placed differently by", alt, ":", if (length(moved)) moved else "none", "\n")
}
# Bootstrap stability (Jaccard; Hennig 2007)
set.seed(2026); Bb <- 1000; jac <- matrix(NA, Bb, k)
for (b in seq_len(Bb)) {
idx <- sample(31, 31, replace = TRUE); u <- unique(idx)
cb <- kmeans(z4[idx, ], centers = k, nstart = 20)$cluster
lab <- cb[match(u, idx)]
for (t in 1:k) {
A <- u[main[u] == t]; if (!length(A)) next
jac[b, t] <- max(sapply(unique(lab), function(h) { Bset <- u[lab == h]
length(intersect(A, Bset)) / length(union(A, Bset)) }))
}
}
cat("Bootstrap Jaccard by type:\n")
print(round(rbind(mean = colMeans(jac, na.rm = TRUE), share_gt_0.75 = colMeans(jac > 0.75, na.rm = TRUE)), 2))
invisible(alts)
}
prov_avg$type <- km_run(z4, K_main)
prov_avg$type3 <- km_run(z4, 3)
alts_main <- describe_partition(prov_avg$type, K_main, "Main partition")
Main partition (K = 4)
Type 1 ( 11 ): Anhui Gansu Guangdong Henan Hubei Hunan Jiangsu Jiangxi Ningxia Shaanxi Shandong
Type 2 ( 10 ): Guangxi Hainan Heilongjiang Inner Mongolia Jilin Liaoning Shanxi Sichuan Tianjin Xinjiang
Type 3 ( 7 ): Fujian Guizhou Hebei Qinghai Tibet Yunnan Zhejiang
Type 4 ( 3 ): Beijing Chongqing Shanghai
type never_rate young_share sex_ratio higher_edu agri_share
1 1 17.5 37.5 107.5 11.7 87.2
2 2 17.1 31.8 113.6 11.7 80.1
3 3 22.4 36.6 123.8 7.2 89.3
4 4 18.1 35.3 102.5 22.8 71.9
Adjusted Rand index: K-means with outcome = 0.822 | Ward = 1 | PAM = 0.891
Placed differently by km5 : Henan Jiangxi
Placed differently by ward : none
Placed differently by pam : Guangxi
Bootstrap Jaccard by type:
[,1] [,2] [,3] [,4]
mean 0.68 0.80 0.77 0.87
share_gt_0.75 0.41 0.59 0.55 0.79
prov_avg$km5 <- alts_main$km5
alts_k3 <- describe_partition(prov_avg$type3, 3, "Three types")
Three types (K = 3)
Type 1 ( 15 ): Fujian Gansu Guangdong Guizhou Hebei Henan Hunan Jiangsu Jiangxi Ningxia Qinghai Shaanxi Tibet Yunnan Zhejiang
Type 2 ( 13 ): Anhui Guangxi Hainan Heilongjiang Hubei Inner Mongolia Jilin Liaoning Shandong Shanxi Sichuan Tianjin Xinjiang
Type 3 ( 3 ): Beijing Chongqing Shanghai
type never_rate young_share sex_ratio higher_edu agri_share
1 1 21.2 37.6 116.5 9.5 88.2
2 2 15.5 32.5 110.6 11.8 81.7
3 3 18.1 35.3 102.5 22.8 71.9
Adjusted Rand index: K-means with outcome = 0.879 | Ward = 0.668 | PAM = 0.879
Placed differently by km5 : Ningxia
Placed differently by ward : Anhui Hubei Shandong
Placed differently by pam : Anhui
Bootstrap Jaccard by type:
[,1] [,2] [,3]
mean 0.78 0.79 0.76
share_gt_0.75 0.60 0.59 0.67
ggplot(prov_avg, aes(young_share, never_rate, colour = factor(type), shape = factor(type), label = prov_en)) +
geom_point(size = 2.6) +
geom_text_repel(size = 2.7, colour = "grey25", seed = 1, max.overlaps = Inf, min.segment.length = 0) +
scale_colour_manual(values = c("tomato", "steelblue", "seagreen", "darkorange", "purple")[1:K_main], labels = paste("Type", 1:K_main)) +
scale_shape_manual(values = c(16, 15, 17, 18, 8)[1:K_main], labels = paste("Type", 1:K_main)) +
labs(x = "Youth share, ages 15-29 (%), province mean", y = "Never-married share (%), province mean",
colour = NULL, shape = NULL) +
theme_bw(base_size = 10) + theme(legend.position = "bottom", panel.grid.minor = element_blank())
# eta-squared with permutation tests (10,000 relabellings, same group sizes)
eta2 <- function(y, g) { m <- ave(y, g); sum((m - mean(y))^2) / sum((y - mean(y))^2) }
perm_test <- function(g, R = 10000) {
y <- prov_avg$never_rate; obs <- eta2(y, g)
null <- replicate(R, eta2(y, sample(g)))
c(eta2 = obs, p_perm = mean(null >= obs), null_95 = unname(quantile(null, 0.95)),
max_gap = diff(range(tapply(y, g, mean))))
}
prov_avg$region <- ifelse(prov_avg$prov_en %in% east, 1, ifelse(prov_avg$prov_en %in% central, 2, 3))
yr_rank <- rank(-prov_avg$young_share, ties.method = "first")
prov_avg$youth3 <- cut(yr_rank, c(0, 10, 21, 31), labels = FALSE)
sz <- cumsum(table(prov_avg$type)); prov_avg$youth_sz <- cut(yr_rank, c(0, sz), labels = FALSE)
for (k in 2:5) prov_avg[[paste0("k", k)]] <- km_run(z4, k)
set.seed(2026)
print(round(t(sapply(c(types_main = "type", types_k2 = "k2", types_k3 = "k3", types_k5 = "k5",
types_main_with_outcome = "km5", regions = "region",
youth_tertiles = "youth3", youth_same_sizes_as_main = "youth_sz"),
function(v) perm_test(prov_avg[[v]]))), 3))
eta2 p_perm null_95 max_gap
types_main 0.164 0.178 0.250 5.235
types_k2 0.087 0.101 0.124 3.032
types_k3 0.279 0.009 0.193 5.702
types_k5 0.341 0.027 0.304 7.934
types_main_with_outcome 0.354 0.009 0.249 7.519
regions 0.048 0.503 0.191 2.806
youth_tertiles 0.079 0.314 0.187 3.452
youth_same_sizes_as_main 0.102 0.402 0.246 5.224
# Permutation test including the choice of K
parts <- as.matrix(prov_avg[, paste0("k", 2:5)])
yprov <- prov_avg$never_rate
obs_k <- apply(parts, 2, function(g) eta2(yprov, g))
set.seed(2026); R <- 10000
null_max <- replicate(R, { yp <- sample(yprov); max(apply(parts, 2, function(g) eta2(yp, g))) })
cat("\nIncluding the choice of K:\n")
Including the choice of K:
print(round(obs_k, 3))
k2 k3 k4 k5
0.087 0.279 0.164 0.341
cat("Largest eta-squared:", round(max(obs_k), 3), "( K =", which.max(obs_k) + 1, ")",
"| 95th percentile of null maximum:", round(quantile(null_max, 0.95), 3),
"| p =", round(mean(null_max >= max(obs_k)), 3), "\n")
Largest eta-squared: 0.341 ( K = 5 ) | 95th percentile of null maximum: 0.304 | p = 0.026
region_means <- aggregate(never_rate ~ region, data = panel, FUN = mean)
region_means$never_rate <- round(region_means$never_rate, 1)
cat("\nNever-married share by region:\n"); print(region_means)
Never-married share by region:
region never_rate
1 Central 17.0
2 East 18.6
3 West 19.8
res <- panel[, c("prov_en", "young_share", "never_rate")]
res$resid <- residuals(lm(never_rate ~ young_share, res))
pr <- aggregate(resid ~ prov_en, res, mean); pr <- pr[order(pr$resid), ]; pr$resid <- round(pr$resid, 1)
cat("\nBelow age-predicted level (top 5):\n"); print(head(pr, 5))
Below age-predicted level (top 5):
prov_en resid
1 Anhui -9.7
23 Shandong -7.5
13 Hubei -6.8
20 Ningxia -4.9
16 Jiangsu -4.9
cat("Above age-predicted level (top 5):\n"); print(tail(pr, 5))
Above age-predicted level (top 5):
prov_en resid
18 Jilin 3.8
10 Hebei 3.9
12 Henan 5.8
19 Liaoning 8.2
28 Tibet 15.1
