分解分析 (Decomposition Analysis)

library(data.table)
library(dplyr) # 读取数据：需要两个时间点的数据（早年和近年）
# 格式：location, year, age, sex, measure, metric, value # 示例：2010年和2021年的中国数据
df_2010 <- fread("examples/decomp_2010.csv")
df_2021 <- fread("examples/decomp_2021.csv") # 合并两个时间点
df <- rbind(df_2010[, .(location, age, sex, measure, value_2010 = value)], df_2021[, .(location, age, sex, measure, value_2021 = value)], by = c("location", "age", "sex", "measure")) # 填充缺失值
df[is.na(value_2010), value_2010 := 0]
df[is.na(value_2021), value_2021 := 0] # 计算总变化
df[, delta := value_2021 - value_2010]

# 读取人口数据（需要各年龄组的人口数）
pop_2010 <- fread("examples/population_2010.csv")
pop_2021 <- fread("examples/population_2021.csv") # 合并疾病数据和人口数据
df <- merge(df, pop_2010, by = c("location", "age", "sex"), suffixes = c("", "_pop"))
df <- merge(df, pop_2021, by = c("location", "age", "sex"), suffixes = c("", "_2021")) # 计算疾病率 (rate = cases / population)
df[, rate_2010 := value_2010 / pop_2010]
df[, rate_2021 := value_2021 / pop_2021] # 分解计算（使用经典的 Das Gupta 方法简化版）
# 1. 人口增长效应 = (pop_2021 - pop_2010) * rate_2010
df[, effect_growth := (pop_2021 - pop_2010) * rate_2010] # 2. 老龄化效应 = pop_2021 * (rate_2010_age - rate_2010_total)
# （需要按年龄加权计算，简化版略过） # 3. 疾病率变化效应 = pop_2021 * (rate_2021 - rate_2010)
df[, effect_rate := pop_2021 * (rate_2021 - rate_2010)] # 汇总各效应
results <- df[, .( growth = sum(effect_growth, na.rm = TRUE), rate_change = sum(effect_rate, na.rm = TRUE), total_change = sum(delta, na.rm = TRUE)
), by = .(location, measure)] # 计算贡献百分比
results[, pct_growth := growth / total_change * 100]
results[, pct_rate := rate_change / total_change * 100]

library(ggplot2) # 绘制堆叠条形图
ggplot(results, aes(x = measure, y = total_change, fill = factor(1))) + geom_bar(stat = "identity", position = "stack") + geom_bar(aes(y = growth), fill = "steelblue", stat = "identity") + geom_bar(aes(y = rate_change), fill = "#e31a1c", stat = "identity") + labs(x = "Measure", y = "Change in Cases", fill = "Effect") + theme_minimal # 绘制贡献百分比
results_long <- melt(results, id.vars = c("location", "measure"), measure.vars = c("pct_growth", "pct_rate"), variable.name = "effect_type", value.name = "percentage") ggplot(results_long, aes(x = measure, y = percentage, fill = effect_type)) + geom_bar(stat = "identity", position = "stack") + scale_fill_manual(values = c("steelblue", "#e31a1c"), labels = c("Population Growth", "Rate Change")) + labs(x = "Measure", y = "Contribution (%)", fill = "Effect Type") + theme_minimal + geom_hline(yintercept = 0) # 保存结果
fwrite(results, "examples/output_decomposition.csv")