1 Overview

This document contains the complete, reproducible analysis for the thesis “Inter-Provincial Variation in the Marital Composition of China’s Floating Population, 2011-2018.” It runs from the raw Excel files all the way to every figure and table in the thesis.

Data source. Province-level floating-population data from the National Earth System Science Data Center, National Science & Technology Infrastructure of China (http://www.geodata.cn), based on the China Migrants Dynamic Survey (CMDS).

Key design choices (explained where they occur):

  • Rates are computed from person counts, not from the raw “proportion” column (which is each province’s share of the national total, not the within-province composition).
  • 2010 is dropped (no marital-status data). 2014 is dropped entirely (anomalous sex ratio and a different table format). Final period: 7 years (2011-2013, 2015-2018).
  • Youth share uses the working-age (15-59) denominator, so all years are comparable (early years lack the 0-14 and 60+ bands).
  • Province names are cleaned so Inner Mongolia is preserved. The Xinjiang Production and Construction Corps and “Total” rows are removed.
  • Final sample: 31 provinces × 7 years = 217 province-years.

2 Setup

library(readxl)   # read Excel files
library(dplyr)    # data manipulation

data_folder <- "/Users/yvainee.yanganker.com/Desktop/论文"

analysis_years <- c(2011, 2012, 2013, 2015, 2016, 2017, 2018)
# Clean province names; preserve Inner Mongolia (bugfix)
clean_prov <- function(x) {
  x <- as.character(x)
  x <- gsub("内蒙古自治区", "内蒙古", x)
  x <- gsub("(壮族|回族|维吾尔|藏族)?自治区$|省$|市$|特别行政区$", "", x)
  trimws(x)
}

# Find a file by keyword + year
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(as.character(year), fs)]
  if (length(fs) == 0) NA else fs[1]
}

3 Step 1 — Marital status (outcome)

Reads each year’s marriage file and computes the within-province never-married and first-married rates from person counts. Column positions are fixed (marriage tables share the same layout, except 2011 which labels first marriage differently).

read_marriage <- function(file_path, year) {
  raw <- as.data.frame(read_excel(file_path, col_names = FALSE))
  if (year == 2011) {
    d <- raw[6:nrow(raw), ]; prov <- d[[9]]
    never <- as.numeric(d[[12]]); first <- as.numeric(d[[14]])
    rem <- 0; div <- as.numeric(d[[16]]); wid <- as.numeric(d[[18]])
  } else {
    d <- raw[6:nrow(raw), ]; prov <- d[[9]]
    never <- as.numeric(d[[12]]); first <- as.numeric(d[[14]])
    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 <- out[nchar(out$province) >= 2 & nchar(out$province) <= 4, ]
  total <- out$never + out$first + out$rem + out$div + out$wid
  out$never_rate <- round(out$never / total * 100, 2)
  out$first_rate <- round(out$first / total * 100, 2)
  out[, c("year","province","never_rate","first_rate")]
}

marriage <- bind_rows(lapply(analysis_years, function(y){
  f <- file_for("婚姻", y); if (is.na(f)) return(NULL); read_marriage(f, y)
}))
cat("Marriage rows:", nrow(marriage),
    "| provinces:", length(unique(marriage$province)), "\n")
## Marriage rows: 218 | provinces: 32

4 Step 2 — English names and region labels

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"
)
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")
marriage$prov_en <- name_map[marriage$province]
marriage$region  <- ifelse(marriage$prov_en %in% east, "East",
                    ifelse(marriage$prov_en %in% central, "Central", "West"))

5 Step 3 — Explanatory variables

5.1 Sex ratio (fixed columns)

read_gender <- function(file_path, year) {
  raw <- as.data.frame(read_excel(file_path, col_names = FALSE))
  d <- raw[6:nrow(raw), ]; prov <- d[[9]]
  male <- as.numeric(d[[12]]); female <- as.numeric(d[[14]])
  out <- data.frame(year = year, province = clean_prov(prov),
                    sex_ratio = round(male / female * 100, 1))
  out <- out[!is.na(out$sex_ratio), ]
  out[nchar(out$province) >= 2 & nchar(out$province) <= 4, ]
}
gender <- bind_rows(lapply(analysis_years, function(y){
  f <- file_for("性别", y); if (is.na(f)) return(NULL); read_gender(f, y)
}))

5.2 Age structure — youth 15-29 over working-age 15-59 (by column name)

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), ]; prov <- d[[9]]
  counts <- sapply(idxs, function(c) as.numeric(d[[c]]))
  young <- rowSums(counts[, 1:3]); denom <- rowSums(counts)   # 15-29 / 15-59
  out <- data.frame(year = year, province = clean_prov(prov),
                    young_share = round(young / denom * 100, 1))
  out <- out[!is.na(out$young_share), ]
  out[nchar(out$province) >= 2 & nchar(out$province) <= 4, ]
}
age <- bind_rows(lapply(analysis_years, function(y){
  f <- file_for("年龄", y); if (is.na(f)) return(NULL); read_age(f, y)
}))

5.3 Education — tertiary share (by column name; handles 2016 merge)

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), ]; prov <- d[[9]]
  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)) hs <- gc("高中") + gc("中专")
  base <- gc("未上过学") + gc("小学") + gc("初中") + gc("大学专科") + gc("大学本科") + gc("研究生")
  out <- data.frame(year = year, province = clean_prov(prov),
                    higher_edu = round(higher / (base + hs) * 100, 1))
  out <- out[!is.na(out$higher_edu), ]
  out[nchar(out$province) >= 2 & nchar(out$province) <= 4, ]
}
edu <- bind_rows(lapply(analysis_years, function(y){
  f <- file_for("教育", y); if (is.na(f)) return(NULL); read_edu(f, y)
}))

5.4 Hukou — agricultural share (by exact column name)

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), ]; prov <- d[[9]]
  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("其他(含无户口)人数")
  out <- data.frame(year = year, province = clean_prov(prov),
                    agri_share = round(agri / total * 100, 1))
  out <- out[!is.na(out$agri_share) & total > 0, ]
  out[nchar(out$province) >= 2 & nchar(out$province) <= 4, ]
}
hukou <- bind_rows(lapply(analysis_years, function(y){
  f <- file_for("户口", y); if (is.na(f)) return(NULL); read_hukou(f, y)
}))

6 Step 4 — Merge into the province-year panel

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"))

# Drop non-standard / unmapped province rows and duplicate province-years
panel <- panel[!is.na(panel$prov_en), ]
panel <- panel[!duplicated(panel[, c("year","prov_en")]), ]

cat("Panel rows:", nrow(panel),
    "| provinces:", length(unique(panel$prov_en)),
    "| years:", paste(sort(unique(panel$year)), collapse=","), "\n")
## Panel rows: 217 | provinces: 31 | years: 2011,2012,2013,2015,2016,2017,2018

7 Step 5 — RQ1: Trend, and RQ2: Inter-provincial variation

trend <- aggregate(cbind(never_rate, first_rate) ~ year, data = panel, FUN = mean)
trend[,2:3] <- round(trend[,2:3], 1)
trend
##   year never_rate first_rate
## 1 2011       20.6       77.7
## 2 2012       21.7       75.0
## 3 2013       21.5       75.2
## 4 2015       18.1       77.6
## 5 2016       16.8       79.4
## 6 2017       15.9       79.6
## 7 2018       12.0       81.7
plot(trend$year, trend$first_rate, type = "o", pch = 1, col = "blue", lwd = 2,
     ylim = c(0, 100), xlab = "Year", ylab = "Share of floating population (%)",
     main = "Marital Composition of China's Floating Population\n2011-2018 (excluding 2014)")
lines(trend$year, trend$never_rate, type = "o", pch = 1, col = "red", lwd = 2)
text(trend$year, trend$first_rate, labels = trend$first_rate, pos = 3, col = "blue", cex = 0.8)
text(trend$year, trend$never_rate, labels = trend$never_rate, pos = 3, col = "red", cex = 0.8)
legend("right", legend = c("Never-married","First-married"),
       col = c("red","blue"), lwd = 2, pch = 1, bg = "white")
Figure 1. Marital composition trend, 2011-2018 (excluding 2014)

Figure 1. Marital composition trend, 2011-2018 (excluding 2014) 

y2018 <- panel[panel$year == 2018, c("prov_en","never_rate")]
y2018 <- y2018[order(y2018$never_rate), ]
cat("2018 range:", round(max(y2018$never_rate) - min(y2018$never_rate), 1), "pp\n")
## 2018 range: 31.6 pp
barplot(y2018$never_rate, names.arg = y2018$prov_en, horiz = TRUE, las = 1,
        col = "darkred", xlab = "Never-married rate (%)",
        main = "Never-married Rate by Province, 2018", cex.names = 0.7)
Figure 2. Never-married rate by province, 2018

Figure 2. Never-married rate by province, 2018

8 Step 6 — RQ3: Correlation and partial correlation

vars <- c("never_rate","young_share","sex_ratio","higher_edu","agri_share")
round(cor(panel[vars], use = "complete.obs"), 2)
##             never_rate young_share sex_ratio higher_edu agri_share
## never_rate        1.00        0.63      0.41      -0.39       0.20
## young_share       0.63        1.00      0.18      -0.36       0.41
## sex_ratio         0.41        0.18      1.00      -0.42       0.17
## higher_edu       -0.39       -0.36     -0.42       1.00      -0.62
## agri_share        0.20        0.41      0.17      -0.62       1.00
age_model <- lm(never_rate ~ young_share, data = panel)

cat("R-squared for age-only model:",
    round(summary(age_model)$r.squared, 2), "\n")
## R-squared for age-only model: 0.39
partial_corr <- function(df, x, y, z) {
  dd <- df[, c(x,y,z)]; dd <- dd[complete.cases(dd), ]
  rx <- dd[[x]] - predict(lm(dd[[x]] ~ dd[[z]]))
  ry <- dd[[y]] - predict(lm(dd[[y]] ~ dd[[z]]))
  cor(rx, ry)
}
for (v in c("sex_ratio","higher_edu","agri_share"))
  cat(v, ": raw =", round(cor(panel$never_rate, panel[[v]], use="complete.obs"), 2),
      " | controlling age =", round(partial_corr(panel, v, "never_rate","young_share"), 2), "\n")
## sex_ratio : raw = 0.41  | controlling age = 0.39 
## higher_edu : raw = -0.39  | controlling age = -0.23 
## agri_share : raw = 0.2  | controlling age = -0.08

8.1 Scatter plots (Figures 3-6)

sc <- function(x, xlab, col, rr) {
  plot(panel[[x]], panel$never_rate, pch = 21, bg = col, col = "gray30",
       xlab = xlab, ylab = "Never-married rate (%)",
       main = paste0(xlab, " vs Never-married Rate (2011-2018)"))
  abline(lm(panel$never_rate ~ panel[[x]]), col = "red", lwd = 3, lty = 2)
  legend("topright", legend = paste0("r = ", rr), col = "red", lwd = 3, lty = 2)
}
sc("young_share", "Youth share 15-29 (%)", rgb(0.2,0.6,0.3,0.5),
   round(cor(panel$young_share, panel$never_rate, use="complete.obs"),2))

sc("sex_ratio", "Sex ratio (M/F)", rgb(0.6,0.4,0.7,0.5),
   round(cor(panel$sex_ratio, panel$never_rate, use="complete.obs"),2))

sc("higher_edu", "Tertiary-education share (%)", rgb(0.4,0.5,0.7,0.5),
   round(cor(panel$higher_edu, panel$never_rate, use="complete.obs"),2))

sc("agri_share", "Agricultural-hukou share (%)", rgb(0.4,0.5,0.7,0.5),
   round(cor(panel$agri_share, panel$never_rate, use="complete.obs"),2))

8.2 Correlation matrix and partial-correlation comparison (Figures 7-8)

cm <- cor(panel[vars], use = "complete.obs")
# simple heatmap
image(1:5, 1:5, t(cm[5:1,]), axes = FALSE, xlab="", ylab="",
      col = colorRampPalette(c("steelblue","white","firebrick"))(50), zlim=c(-1,1),
      main = "Correlation Matrix of Provincial Indicators")
axis(1, 1:5, colnames(cm), las=2, cex.axis=0.8)
axis(2, 1:5, rev(colnames(cm)), las=1, cex.axis=0.8)
for (i in 1:5) for (j in 1:5) text(j, 6-i, round(cm[i,j],2), cex=0.8)

raw <- sapply(c("sex_ratio","higher_edu","agri_share"),
              function(v) cor(panel$never_rate, panel[[v]], use="complete.obs"))
ctl <- sapply(c("sex_ratio","higher_edu","agri_share"),
              function(v) partial_corr(panel, v, "never_rate","young_share"))
m <- rbind(raw, ctl)
barplot(m, beside = TRUE, names.arg = c("Sex ratio","Tertiary edu","Agri hukou"),
        col = c("steelblue","orange"), ylab = "Correlation with never-married rate",
        main = "Correlation Before vs After Controlling for Age",
        legend.text = c("Raw","Controlling age"), args.legend = list(x="topright"))
abline(h = 0)

9 Step 7 — RQ3: Regression (OLS, +year FE, two-way FE)

reg <- na.omit(panel[, c(vars, "year", "prov_en")])
cat("Regression N =", nrow(reg), "( = provinces x 7 )\n\n")
## Regression N = 217 ( = provinces x 7 )
m1 <- lm(never_rate ~ young_share + sex_ratio + higher_edu + agri_share, reg)
m2 <- lm(never_rate ~ young_share + sex_ratio + higher_edu + agri_share + factor(year), reg)
m3 <- lm(never_rate ~ young_share + sex_ratio + higher_edu + agri_share +
           factor(year) + factor(prov_en), reg)

core <- c("young_share","sex_ratio","higher_edu","agri_share")
cat("--- Model 1: pooled OLS ---\n"); print(round(summary(m1)$coefficients[core,], 3))
## --- Model 1: pooled OLS ---
##             Estimate Std. Error t value Pr(>|t|)
## young_share    0.600      0.054  11.206    0.000
## sex_ratio      0.129      0.027   4.825    0.000
## higher_edu    -0.224      0.075  -2.999    0.003
## agri_share    -0.202      0.061  -3.323    0.001
cat("R2 =", round(summary(m1)$r.squared, 3), "\n\n")
## R2 = 0.514
cat("--- Model 2: + year FE ---\n"); print(round(summary(m2)$coefficients[core,], 3))
## --- Model 2: + year FE ---
##             Estimate Std. Error t value Pr(>|t|)
## young_share    0.909      0.087  10.456        0
## sex_ratio      0.140      0.026   5.366        0
## higher_edu    -0.435      0.090  -4.805        0
## agri_share    -0.371      0.072  -5.189        0
cat("R2 =", round(summary(m2)$r.squared, 3), "\n\n")
## R2 = 0.564
cat("--- Model 3: two-way FE ---\n"); print(round(summary(m3)$coefficients[core,], 3))
## --- Model 3: two-way FE ---
##             Estimate Std. Error t value Pr(>|t|)
## young_share    0.816      0.061  13.459    0.000
## sex_ratio      0.074      0.019   4.008    0.000
## higher_edu    -0.118      0.094  -1.256    0.211
## agri_share    -0.162      0.078  -2.072    0.040
cat("R2 =", round(summary(m3)$r.squared, 3), "\n")
## R2 = 0.917

Reading the models. Age structure and sex ratio stay significant and robust across all three models. Education becomes not significant under two-way fixed effects (p ≈ 0.21). Hukou shows a sign reversal relative to its simple correlation (+0.20 in raw correlation, about −0.20 in the regression), because it is confounded with age and education.

10 Step 8 — RQ4: Clustering and residual analysis

prov_avg <- aggregate(panel[vars], by = list(prov_en = panel$prov_en), FUN = mean)
rownames(prov_avg) <- prov_avg$prov_en
set.seed(42)
km <- kmeans(scale(prov_avg[vars]), centers = 3, nstart = 25)
prov_avg$cluster <- km$cluster

for (c in 1:3) { cat("Cluster", c, ":", prov_avg$prov_en[prov_avg$cluster==c], "\n\n") }
## Cluster 1 : Fujian Guangdong Guizhou Hebei Henan Hunan Jiangsu Jiangxi Qinghai Shaanxi Tibet Yunnan Zhejiang 
## 
## Cluster 2 : Anhui Gansu Guangxi Hainan Heilongjiang Hubei Inner Mongolia Jilin Liaoning Ningxia Shandong Shanxi Sichuan Tianjin Xinjiang 
## 
## Cluster 3 : Beijing Chongqing Shanghai
cluster_means <- aggregate(prov_avg[vars], by = list(cluster = prov_avg$cluster), FUN = mean)
cluster_means[vars] <- round(cluster_means[vars], 1)
cluster_means
##   cluster never_rate young_share sex_ratio higher_edu agri_share
## 1       1       21.5        37.4     117.2        9.4       88.4
## 2       2       15.3        32.5     110.8       11.4       82.3
## 3       3       17.5        34.7     102.5       22.4       71.8
plot(prov_avg$young_share, prov_avg$never_rate, col = prov_avg$cluster, pch = 19, cex = 1.5,
     xlab = "Youth share (%)", ylab = "Never-married rate (%)",
     main = "Provincial Clusters by Floating-Population Profile")
text(prov_avg$young_share, prov_avg$never_rate, prov_avg$prov_en, pos = 3, cex = 0.6)

cat("Never-married rate by region:\n")
## Never-married rate by region:
region_means <- aggregate(never_rate ~ region, data = panel, FUN = mean)
region_means$never_rate <- round(region_means$never_rate, 1)
print(region_means)
##    region never_rate
## 1 Central       16.6
## 2    East       18.2
## 3    West       19.2
cluster_gap <- max(cluster_means$never_rate) - min(cluster_means$never_rate)
region_gap  <- max(region_means$never_rate) - min(region_means$never_rate)

cat("\nCluster gap:",
    round(cluster_gap, 1), "percentage points\n")
## 
## Cluster gap: 6.2 percentage points
cat("Region gap:",
    round(region_gap, 1), "percentage points\n")
## Region gap: 2.6 percentage points
res <- na.omit(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("Over-married (actual < age-predicted), top 4:\n")
## Over-married (actual < age-predicted), top 4:
print(head(pr, 4))
##     prov_en resid
## 1     Anhui  -9.4
## 23 Shandong  -6.8
## 13    Hubei  -6.5
## 20  Ningxia  -4.8
cat("\nOver-single (actual > age-predicted), top 4:\n")
## 
## Over-single (actual > age-predicted), top 4:
print(tail(pr, 4))
##     prov_en resid
## 10    Hebei   4.2
## 12    Henan   5.6
## 19 Liaoning   7.5
## 28    Tibet  14.4

11 Data source statement

The province-level floating-population data used in this study were obtained from the National Earth System Science Data Center, National Science & Technology Infrastructure of China (http://www.geodata.cn).