GBD 分析方法与复现教程

library(dplyr) calc_eapc <- function(df_rate) { stopifnot(all(c("location", "year", "val") %in% names(df_rate))) df_rate %>% filter(!is.na(val), val > 0) %>% group_by(location) %>% summarise( n_year = n, beta = coef(lm(log(val) ~ year))[2], se = summary(lm(log(val) ~ year))$coefficients[2, 2], .groups = "drop" ) %>% mutate( eapc = (exp(beta) - 1) * 100, lci = (exp(beta - 1.96 * se) - 1) * 100, uci = (exp(beta + 1.96 * se) - 1) * 100 ) %>% mutate(across(c(eapc, lci, uci), ~round(.x, 3)))
} # 示例
input_df <- read.csv("examples/eapc_input_toy.csv")
out_eapc <- calc_eapc(input_df)
write.csv(out_eapc, "examples/output_eapc_toy.csv", row.names = FALSE)

library(dplyr) prepare_joinpoint_input <- function(df, measure_name) { df %>% filter(age == "Age-standardized", metric == "Rate", measure == measure_name) %>% mutate(SE = (upper - lower) / (2 * 1.96)) %>% select(location_name = location, year, val, SE) %>% arrange(location_name, year)
} jp_in <- read.csv("examples/joinpoint_input_toy.csv")
jp_export <- prepare_joinpoint_input(jp_in, "Incidence")
write.csv(jp_export, "examples/output_joinpoint_export.csv", row.names = FALSE) # 外部 Joinpoint 运行后，将结果保存为 examples/joinpoint_aapc_result.txt
aapc_raw <- read.table("examples/joinpoint_aapc_result_toy.txt", header = TRUE)
aapc_tbl <- aapc_raw %>% transmute( location_name, val = round(AAPC, 3), lower = round(LowerCI, 3), upper = round(UpperCI, 3), AAPC = sprintf("%.3f (%.3f, %.3f)", val, lower, upper) )
write.csv(aapc_tbl, "examples/output_aapc_table.csv", row.names = FALSE)

library(dplyr)
library(ggplot2)
library(maps) map_country_names <- function(df) { name_map <- c( "United States of America" = "USA", "Russian Federation" = "Russia", "Republic of Korea" = "South Korea", "Democratic People's Republic of Korea" = "North Korea", "Viet Nam" = "Vietnam", "Cote d'Ivoire" = "Ivory Coast" ) df %>% mutate(location = recode(location, !!!name_map))
} map_df <- read.csv("examples/map_input_toy.csv")
map_df <- map_country_names(map_df)
world_data <- map_data("world")
merged_map <- left_join(world_data, map_df, by = c("region" = "location")) unmatched <- setdiff(unique(map_df$location), unique(world_data$region))
write.csv(data.frame(location = unmatched), "examples/output_unmatched_countries.csv", row.names = FALSE) ggplot(merged_map, aes(long, lat, group = group, fill = val)) + geom_polygon(color = "grey70", linewidth = 0.1) + scale_fill_viridis_c(option = "C", na.value = "grey95") + coord_fixed(1.3) + theme_void

library(dplyr)
library(cluster)
library(factoextra) clu_df <- read.csv("examples/cluster_input_toy.csv")
row.names(clu_df) <- clu_df$location
x <- clu_df %>% select(-location) %>% scale hc <- hclust(dist(x, method = "euclidean"), method = "complete") sil_scores <- sapply(2:6, function(k) { grp <- cutree(hc, k = k) mean(silhouette(grp, dist(x))[, 3])
}) best_k <- which.max(sil_scores) + 1
cluster_out <- data.frame(location = row.names(x), cluster = cutree(hc, k = best_k))
write.csv(cluster_out, "examples/output_cluster_labels.csv", row.names = FALSE) fviz_dend(hc, k = best_k, cex = 0.6, rect = TRUE, horiz = TRUE)

pkgs <- c("dplyr", "ggplot2", "maps", "cluster", "factoextra")
print(R.version.string)
print(sapply(pkgs, packageVersion))
sessionInfo

### Statistical analysis
We extracted [measure] rates for [disease/exposure] from [start_year] to [end_year] at the [age_group] level. Temporal trend was quantified using estimated annual percentage change (EAPC) by fitting a log-linear model:
ln(rate) ~ year.
EAPC was calculated as (exp(β)-1)×100, with 95% confidence interval (CI) derived from β ± 1.96×SE. For segmented temporal patterns, joinpoint regression was applied and average annual percentage change (AAPC) with 95% CI was reported.
Spatial distribution was visualized using country-level choropleth maps after harmonizing country-name mappings.
Unsupervised hierarchical clustering (Euclidean distance, complete linkage) was used to identify pattern groups across countries/regions.

### Temporal trend
From [start_year] to [end_year], the age-standardized [measure] rate of [disease] in [population/region] changed from [value_start] to [value_end] per 100,000.
The overall trend showed an EAPC of [eapc]% (95% CI: [lci] to [uci]), indicating a [decreasing/increasing/stable] pattern. Joinpoint analysis identified [n] significant turning point(s), and the corresponding AAPC was [aapc]% (95% CI: [aapc_lci] to [aapc_uci]).
Spatially, higher burdens were observed in [high_regions], whereas lower values were concentrated in [low_regions]. Cluster analysis grouped countries/regions into [k] clusters with distinct epidemiological profiles, suggesting heterogeneous burden trajectories across development contexts.

**Figure [X]. Global pattern and temporal trend of [disease/measure].**
(A) Choropleth map of age-standardized [measure] rate in [year], per 100,000 population.
(B) Temporal trend of age-standardized rate from [start_year] to [end_year], with estimated annual percentage change (EAPC) and 95% CI.
(C) Joinpoint-derived average annual percentage change (AAPC) across selected regions.
(D) Hierarchical clustering dendrogram showing pattern groups based on trend and burden indicators. Abbreviations: EAPC, estimated annual percentage change; AAPC, average annual percentage change; ASR, age-standardized rate; CI, confidence interval.