时间序列多模型批处理及可视化分析是指通过并行运行多种时间序列模型（如ARIMA、ARCH、VAR等）对大规模数据进行批量处理，并结合可视化技术呈现分析结果的方法。其公共卫生意义主要体现在以下方面：
一、核心技术与流程
1.多模型批处理
并行建模：同时应用ARIMA（处理趋势与季节性）、ARCH（分析波动性）、VAR/VEC（多变量关联）等模型进行批量计算。
智能优化：通过AIC准则自动筛选最优模型参数（如p/d/q阶数），提升预测效率。
批量处理能力：支持对多组数据（如不同病种发病率）同步建模，减少人工干预。
2.可视化关键技术
多维呈现：采用交互式线图展示趋势，堆叠柱形图/玫瑰图对比多指标周期规律（如季度门诊量）。
认知优化：通过降噪设计（如移动平均平滑）和分层着色，降低多数据并行的解读负担。
动态仪表盘：整合预测结果与历史数据，实时监控关键指标（如床位占用率）。
二、公共卫生应用价值
1.疫情预警与资源调配
预测疾病暴发周期（如流感季住院率高峰），指导疫苗、床位及医护人员配置。
通过多变量模型分析环境因素（如气温）与传染病传播的关联性。
2.医疗效率优化
识别门诊量/手术量的长期趋势与季节性波动，优化排班与设备采购计划。
监测药品消耗规律，避免短缺或过量库存。
3.政策效果评估
可视化对比政策实施前后关键指标（如疫苗接种率），量化干预措施效果。
通过残差分析检测异常事件（如突发公共卫生事件）的影响。
三、实施挑战与建议
数据质量：需规范时间序列声明（tsset命令）并处理缺失值/断点。
模型适配：非平稳数据需差分处理，多变量场景优先选择VAR/VEC模型。
可视化设计：避免过度叠加图表，优先采用交互式动态展示提升可读性。
案例：某医院通过ARIMA预测住院需求，结合玫瑰图可视化季节性峰值，将床位准备效率提升32%。
时间序列分析模型用于预测基于时间顺序的数据模式，广泛应用于经济、医疗和环境等领域。以下分两大类介绍常用模型：
一、传统统计模型
自回归模型 (AR)：通过历史值预测当前值，适用于平稳序列，如经济指标预测。
移动平均模型 (MA)：基于随机误差项建模当前值，常用于噪声数据处理。
ARMA模型：结合AR和MA特性，优化平稳序列的预测精度。
ARIMA模型：扩展ARMA模型，通过差分处理非平稳序列，支持季节性调整（如SARIMA）。
向量自回归模型 (VAR/VEC)：处理多变量序列（如气温与疾病关联），捕捉变量间动态关系。
二、现代机器学习与深度学习模型
LSTM：通过门控机制处理长期依赖，适用于复杂时序模式（如传染病传播预测）。
Informer：优化注意力机制，提升长序列预测效率（如资源需求分析）。
Transformer-based模型：如Time-LLM和TimesNet，利用全局信息建模非线性趋势（如药品库存预警）。
模型选择需结合数据特性：传统模型适用于小规模平稳数据，而现代模型更擅长处理高维非平稳序列。
# 设置工作目录和清理环境
rm(list = ls())
dev.off()
# 加载必要的包
if (!require(pacman)) install.packages("pacman")
pacman::p_load(tidyverse, lubridate, openxlsx, ggmap, sf, maps,
mapdata, officer, flextable, ggpubr, pyramid, forecast,
prophet, bsts, e1071, mem, mgcv, dlnm, splines,
gridExtra, plotly, leaflet, tmap, htmlwidgets, webshot, rmarkdown, knitr)
# 设置工作目录为桌面
setwd("C:/Users/L/Desktop/")
# 创建结果文件夹
if (!dir.exists("Results3")) dir.create("Results3")
# 1. 读取数据
df <- read.csv("医疗卡片数据库.csv", stringsAsFactors = FALSE)
df$发病日期 <- as.Date(df$发病日期, format = "%Y/%m/%d")
# 2. 生成时间序列数据 - 日、周、月尺度
# 按日汇总
daily_ts <- df %>%
group_by(日期 = as.Date(发病日期)) %>%
summarise(病例数 = n(), .groups = "drop") %>%
complete(日期 = seq.Date(min(日期), max(日期), by = "day"),
fill = list(病例数 = 0))
# 按周汇总
weekly_ts <- df %>%
mutate(周 = floor_date(发病日期, "week")) %>%
group_by(周) %>%
summarise(病例数 = n(), .groups = "drop") %>%
complete(周 = seq.Date(min(周), max(周), by = "week"),
fill = list(病例数 = 0))
# 按月汇总
monthly_ts <- df %>%
mutate(月 = floor_date(发病日期, "month")) %>%
group_by(月) %>%
summarise(病例数 = n(), .groups = "drop") %>%
complete(月 = seq.Date(min(月), max(月), by = "month"),
fill = list(病例数 = 0))
# 保存时间序列数据
library(writexl)
write_xlsx(list(日数据 = daily_ts, 周数据 = weekly_ts, 月数据 = monthly_ts),
"Results3/时间序列数据.xlsx")
# 3. 时间序列分析函数（修复版）
run_time_series_analysis <- function(data, time_unit) {
results <- list()
# 创建时间序列对象
ts_data <- ts(data$病例数, frequency = switch(time_unit, "day" = 365, "week" = 52, "month" = 12))
# 移动平均/中位数
results$移动平均 <- tryCatch({
forecast::ma(ts_data, order = switch(time_unit, "day" = 7, "week" = 4, "month" = 3))
}, error = function(e) {
message("移动平均计算失败: ", e$message)
NULL
})
results$移动中位数 <- tryCatch({
k <- switch(time_unit, "day" = 7, "week" = 4, "month" = 3)
# 确保k是奇数
if (k %% 2 == 0) k <- k + 1
runmed(ts_data, k = k)
}, error = function(e) {
message("移动中位数计算失败: ", e$message)
NULL
})
# ARIMA模型
results$ARIMA <- tryCatch({
auto.arima(ts_data)
}, error = function(e) {
message("ARIMA模型拟合失败: ", e$message)
NULL
})
# SARIMA模型
results$SARIMA <- tryCatch({
auto.arima(ts_data, seasonal = TRUE)
}, error = function(e) {
message("SARIMA模型拟合失败: ", e$message)
NULL
})
# 支持向量回归 - 简化模型
results$SVR <- tryCatch({
svr_data <- data.frame(y = ts_data, t = 1:length(ts_data))
svm(y ~ t, data = svr_data)
}, error = function(e) {
message("SVR模型拟合失败: ", e$message)
NULL
})
# 移动流行区间法 (MEM) - 仅周尺度
if(time_unit == "week") {
results$MEM <- tryCatch({
memmodel(i.data = matrix(ts_data, ncol = 1), i.seasons = 3)
}, error = function(e) {
message("MEM模型拟合失败: ", e$message)
NULL
})
}
# GAM模型 - 简化实现
results$GAM <- tryCatch({
data_gam <- data.frame(cases = ts_data, time = 1:length(ts_data))
gam(cases ~ s(time, k = 12), data = data_gam, family = poisson)
}, error = function(e) {
message("GAM模型拟合失败: ", e$message)
NULL
})
return(results)
}
# 4. 执行不同尺度的时间序列分析
analysis_results <- list()
# 日尺度分析
analysis_results$日尺度 <- run_time_series_analysis(daily_ts, "day")
# 周尺度分析
analysis_results$周尺度 <- run_time_series_analysis(weekly_ts, "week")
# 月尺度分析
analysis_results$月尺度 <- run_time_series_analysis(monthly_ts, "month")
# 保存分析结果
saveRDS(analysis_results, "Results3/时间序列分析结果.rds")
# 5. 结果可视化函数
visualize_results <- function(data, results, time_unit, prefix) {
plots <- list()
# 基础时间序列
base_plot <- ggplot(data, aes(x = .data[[names(data)[1]]], y = .data[[names(data)[2]]])) +
geom_line(color = "steelblue") +
labs(title = paste("原始", time_unit, "尺度时间序列"),
x = "时间", y = "病例数") +
theme_minimal()
plots[[paste0(prefix, "_原始序列")]] <- base_plot
# 移动平均
if (!is.null(results$移动平均) && length(results$移动平均) > 0) {
ma_data <- data
ma_data$移动平均 <- as.numeric(results$移动平均)
p_ma <- ggplot(ma_data, aes(x = .data[[names(ma_data)[1]]])) +
geom_line(aes(y = 移动平均), color = "red", linewidth = 1) +
geom_line(aes(y = .data[[names(ma_data)[2]]]), color = "steelblue", alpha = 0.5) +
labs(title = paste("移动平均 -", time_unit, "尺度"),
x = "时间", y = "病例数") +
theme_minimal()
plots[[paste0(prefix, "_移动平均")]] <- p_ma
}
# ARIMA拟合
if (!is.null(results$ARIMA)) {
arima_data <- data
arima_data$ARIMA拟合 <- as.numeric(fitted(results$ARIMA))
p_arima <- ggplot(arima_data, aes(x = .data[[names(arima_data)[1]]])) +
geom_line(aes(y = ARIMA拟合), color = "darkgreen", linewidth = 1) +
geom_line(aes(y = .data[[names(arima_data)[2]]]), color = "steelblue", alpha = 0.5) +
labs(title = paste("ARIMA模型拟合 -", time_unit, "尺度"),
x = "时间", y = "病例数") +
theme_minimal()
plots[[paste0(prefix, "_ARIMA")]] <- p_arima
}
# SARIMA拟合
if (!is.null(results$SARIMA)) {
sarima_data <- data
sarima_data$SARIMA拟合 <- as.numeric(fitted(results$SARIMA))
p_sarima <- ggplot(sarima_data, aes(x = .data[[names(sarima_data)[1]]])) +
geom_line(aes(y = SARIMA拟合), color = "purple", linewidth = 1) +
geom_line(aes(y = .data[[names(sarima_data)[2]]]), color = "steelblue", alpha = 0.5) +
labs(title = paste("SARIMA模型拟合 -", time_unit, "尺度"),
x = "时间", y = "病例数") +
theme_minimal()
plots[[paste0(prefix, "_SARIMA")]] <- p_sarima
}
# SVR拟合
if (!is.null(results$SVR)) {
svr_data <- data
svr_data$SVR拟合 <- predict(results$SVR)
p_svr <- ggplot(svr_data, aes(x = .data[[names(svr_data)[1]]])) +
geom_line(aes(y = SVR拟合), color = "brown", linewidth = 1) +
geom_line(aes(y = .data[[names(svr_data)[2]]]), color = "steelblue", alpha = 0.5) +
labs(title = paste("支持向量回归拟合 -", time_unit, "尺度"),
x = "时间", y = "病例数") +
theme_minimal()
plots[[paste0(prefix, "_SVR")]] <- p_svr
}
# GAM拟合
if (!is.null(results$GAM)) {
gam_data <- data
gam_data$GAM拟合 <- predict(results$GAM, type = "response")
p_gam <- ggplot(gam_data, aes(x = .data[[names(gam_data)[1]]])) +
geom_line(aes(y = GAM拟合), color = "darkblue", linewidth = 1) +
geom_line(aes(y = .data[[names(gam_data)[2]]]), color = "steelblue", alpha = 0.5) +
labs(title = paste("广义加性模型拟合 -", time_unit, "尺度"),
x = "时间", y = "病例数") +
theme_minimal()
plots[[paste0(prefix, "_GAM")]] <- p_gam
}
# MEM结果 (仅周尺度)
if(time_unit == "周" && !is.null(results$MEM)) {
tryCatch({
jpeg("Results3/MEM_temp.jpg", width = 800, height = 600)
plot(results$MEM, main = paste("移动流行区间法 -", time_unit, "尺度"))
dev.off()
plots[[paste0(prefix, "_MEM")]] <- "Results3/MEM_temp.jpg"
}, error = function(e) {
message("MEM可视化失败: ", e$message)
})
}
return(plots)
}
# 6. 生成可视化结果
all_plots <- list()
# 日尺度可视化
all_plots$日尺度 <- visualize_results(daily_ts, analysis_results$日尺度, "日", "日尺度")
# 周尺度可视化
all_plots$周尺度 <- visualize_results(weekly_ts, analysis_results$周尺度, "周", "周尺度")
# 月尺度可视化
all_plots$月尺度 <- visualize_results(monthly_ts, analysis_results$月尺度, "月", "月尺度")
# 7. 保存可视化结果
save_visualizations <- function(plots, time_unit) {
prefix <- paste0("Results3/", time_unit, "_")
for (name in names(plots)) {
if (is.ggplot(plots[[name]])) {
# 保存静态图像
ggsave(paste0(prefix, name, ".jpg"), plots[[name]], width = 10, height = 6, dpi = 300)
ggsave(paste0(prefix, name, ".pdf"), plots[[name]], width = 10, height = 6, device = cairo_pdf)
# 保存ggplot对象
saveRDS(plots[[name]], paste0(prefix, name, ".rds"))
# 保存交互式HTML
interactive_plot <- ggplotly(plots[[name]])
saveWidget(interactive_plot, paste0(prefix, name, ".html"), selfcontained = TRUE)
} elseif (is.character(plots[[name]])) {
# 处理基础图形文件
file.copy(plots[[name]], paste0(prefix, name, ".jpg"), overwrite = TRUE)
}
}
}
# 保存所有可视化
if (!is.null(all_plots$日尺度)) save_visualizations(all_plots$日尺度, "日尺度")
if (!is.null(all_plots$周尺度)) save_visualizations(all_plots$周尺度, "周尺度")
if (!is.null(all_plots$月尺度)) save_visualizations(all_plots$月尺度, "月尺度")
# 8. 生成Word报告
doc <- officer::read_docx()
# 添加报告标题
doc <- doc %>%
body_add_par("传染病时间序列分析报告", style = "heading 1") %>%
body_add_par(paste("报告生成日期:", Sys.Date()), style = "Normal") %>%
body_add_par(" ", style = "Normal") %>%
body_add_par("一、分析概述", style = "heading 2") %>%
body_add_par("本报告对传染病发病数据进行了多尺度时间序列分析，包括：", style = "Normal") %>%
body_add_par("- 日、周、月尺度时间序列分析", style = "Normal") %>%
body_add_par("- 移动平均/中位数", style = "Normal") %>%
body_add_par("- ARIMA/SARIMA模型", style = "Normal") %>%
body_add_par("- 支持向量回归(SVR)", style = "Normal") %>%
body_add_par("- 移动流行区间法(MEM)", style = "Normal") %>%
body_add_par("- 分布滞后非线性模型(DLNM)", style = "Normal") %>%
body_add_par("- 广义加性模型(GAM)", style = "Normal")
# 添加日尺度分析结果
doc <- doc %>%
body_add_par("二、日尺度分析结果", style = "heading 2")
img_files <- c(
"Results3/日尺度_日尺度_原始序列.jpg",
"Results3/日尺度_日尺度_移动平均.jpg",
"Results3/日尺度_日尺度_ARIMA.jpg"
)
for (i in seq_along(img_files)) {
if (file.exists(img_files[i])) {
doc <- doc %>%
body_add_par(paste0("图", i, ": ", c("日尺度原始序列", "日尺度移动平均", "日尺度ARIMA拟合")[i]), style = "heading 3") %>%
body_add_img(src = img_files[i], width = 6, height = 4)
}
}
# 添加周尺度分析结果
doc <- doc %>%
body_add_par("三、周尺度分析结果", style = "heading 2")
img_files <- c(
"Results3/周尺度_周尺度_原始序列.jpg",
"Results3/周尺度_周尺度_SARIMA.jpg",
"Results3/周尺度_周尺度_MEM.jpg"
)
for (i in seq_along(img_files)) {
if (file.exists(img_files[i])) {
doc <- doc %>%
body_add_par(paste0("图", i+3, ": ", c("周尺度原始序列", "周尺度SARIMA拟合", "周尺度MEM分析")[i]), style = "heading 3") %>%
body_add_img(src = img_files[i], width = 6, height = 4)
}
}
# 添加月尺度分析结果
doc <- doc %>%
body_add_par("四、月尺度分析结果", style = "heading 2")
img_files <- c(
"Results3/月尺度_月尺度_原始序列.jpg",
"Results3/月尺度_月尺度_GAM.jpg",
"Results3/月尺度_月尺度_DLNM1.jpg"
)
for (i in seq_along(img_files)) {
if (file.exists(img_files[i])) {
doc <- doc %>%
body_add_par(paste0("图", i+6, ": ", c("月尺度原始序列", "月尺度GAM拟合", "月尺度DLNM分析")[i]), style = "heading 3") %>%
body_add_img(src = img_files[i], width = 6, height = 4)
}
}
# 添加结论与建议
doc <- doc %>%
body_add_par("五、结论与建议", style = "heading 2") %>%
body_add_par("1. 时间趋势分析表明...", style = "Normal") %>%
body_add_par("2. 不同模型拟合效果比较...", style = "Normal") %>%
body_add_par("3. 建议防控措施...", style = "Normal") %>%
body_add_par("六、附录", style = "heading 2") %>%
body_add_par("完整分析结果和数据可在Results3文件夹中查看。", style = "Normal")
# 保存Word报告
print(doc, target = "Results3/传染病时间序列分析报告.docx")
# 9. 生成HTML报告
rmd_content <- c(
"---",
"title: \"传染病时间序列交互式分析报告\"",
"output: html_document",
"---",
"",
"```{r setup, include=FALSE}",
"knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)",
"library(plotly)",
"library(ggplot2)",
"library(htmltools)",
"```",
"",
"# 时间序列分析报告",
"",
"## 日尺度分析",
"",
"### 原始序列",
"```{r day1}",
"if (file.exists('Results3/日尺度_日尺度_原始序列.rds')) {",
" daily_plot <- readRDS('Results3/日尺度_日尺度_原始序列.rds')",
" ggplotly(daily_plot)",
"} else {",
" cat('图表不可用')",
"}",
"```",
"",
"### ARIMA模型拟合",
"```{r day2}",
"if (file.exists('Results3/日尺度_日尺度_ARIMA.rds')) {",
" arima_plot <- readRDS('Results3/日尺度_日尺度_ARIMA.rds')",
" ggplotly(arima_plot)",
"} else {",
" cat('图表不可用')",
"}",
"```",
"",
"## 周尺度分析",
"",
"### SARIMA模型拟合",
"```{r week1}",
"if (file.exists('Results3/周尺度_周尺度_SARIMA.rds')) {",
" sarima_plot <- readRDS('Results3/周尺度_周尺度_SARIMA.rds')",
" ggplotly(sarima_plot)",
"} else {",
" cat('图表不可用')",
"}",
"```",
"",
"## 月尺度分析",
"",
"### GAM模型拟合",
"```{r month1}",
"if (file.exists('Results3/月尺度_月尺度_GAM.rds')) {",
" gam_plot <- readRDS('Results3/月尺度_月尺度_GAM.rds')",
" ggplotly(gam_plot)",
"} else {",
" cat('图表不可用')",
"}",
"```"
)
# 保存Rmd文件
writeLines(rmd_content, "时间序列分析报告.Rmd")
# 渲染HTML报告
rmarkdown::render("时间序列分析报告.Rmd",
output_file = "Results3/传染病时间序列交互式分析报告.html")
cat("分析完成！所有结果已保存到Results3文件夹:\n")
cat("- 时间序列数据.xlsx\n")
cat("- 时间序列分析结果.rds\n")
cat("- 各种分析图表（JPG/PDF/HTML格式）\n")
cat("- 传染病时间序列分析报告.docx\n")
cat("- 传染病时间序列交互式分析报告.html\n")
医学统计数据分析分享交流SPSS、R语言、Python、ArcGis、Geoda、GraphPad、数据分析图表制作等心得。承接数据分析，论文返修，医学统计，机器学习，生存分析，空间分析，问卷分析业务。若有投稿和数据分析代做需求，可以直接联系我，谢谢！
公众号右下角-联系作者，
可加我微信，邀请入粉丝群！
【临床】有临床流行病学数据分析如（t检验、方差分析、χ2检验、logistic回归）、（重复测量方差分析与配对T检验、ROC曲线）、（非参数检验、生存分析、样本含量估计）、（筛检试验：灵敏度、特异度、约登指数等计算）、（绘制柱状图、散点图、小提琴图、列线图等）、机器学习、深度学习、生存分析等需求的同仁们，加入【临床】粉丝群。
【公卫】疾控，公卫岗位的同仁，可以加一下【公卫】粉丝群，分享生态学研究、空间分析、时间序列、监测数据分析、时空面板技巧等工作科研自动化内容。
【生信】有实验室数据分析需求的同仁们，可以加入【生信】粉丝群，交流NCBI（基因序列）、UniProt（蛋白质）、KEGG（通路）、GEO（公共数据集）等公共数据库、基因组学转录组学蛋白组学代谢组学表型组学等数据分析和可视化内容。
往期推荐：【监测预警自动化】系列教程
往期推荐：样本含量估计（样本量计算与功效分析）
往期推荐：SPSS、R语言、Python等临床数据分析专题
往期推荐：科研图表绘制专题
往期推荐：重复测量数据分析专题
往期推荐：生信分析、基因测序数据、实验室数据专题
往期推荐：生存分析及机器学习
往期推荐：二分类因变量机器学习及相关评价可视化
技术分享|如何综合评价临床预测模型？手把手教你学会9种机器学习模型×5种模型评价曲线的综合评价方法（附Python批处理代码）
一文读懂十模型lasso、贝叶斯、KNN、Logistic、决策树、随机森林、SVM、神经网络、XGBoost、lightGBM
一文读懂十模型lasso、贝叶斯、KNN、Logistic、决策树、随机森林、SVM、神经网络、XGBoost、lightGBM
【Python机器学习】十分钟读懂Logistic、决策树、随机森林、SVM、神经网络、XGBoost、lightGBM七种模型
【R语言机器学习】十分钟读懂Logistic回归、决策树、随机森林、SVM、神经网络、XGBoost、lightGBM七种模型
往期推荐：时间序列分析
往期推荐：地统计分析-GIS、地图、相关、聚类、回归
往期推荐：科研自动化探究
往期推荐：趣味阅读
