如何在R中绘制含Theil-Sen斜率与最小二乘回归的时间序列图

架构老林

2026-6-15

问题描述

我有一个包含POSIXct格式日期（Date列）和结果值（Result列）的数据框，结构如下：

structure(list(Date = structure(c(1704067200, 1704153600, 1704240000, 
1704326400, 1704412800, 1704499200, 1704585600, 1704672000, 1704758400, 
1704844800, 1704931200, 1705017600, 1705104000, 1705190400, 1705276800, 
1705363200, 1705449600, 1705536000, 1705622400, 1705708800, 1705795200, 
1705881600, 1705968000, 1706054400, 1706140800, 1706227200, 1706313600, 
1706400000, 1706486400, 1706572800, 1706659200, 1706745600, 1706832000, 
1706918400, 1707004800, 1707091200, 1707177600), class = c("POSIXct", 
"POSIXt"), tzone = "UTC"), Result = c(0.895103686169796, 0.0556219135717945, 
0.807153205141622, 0.461672184316161, 0.547758011882109, 0.927777651259266, 
0.308058477041002, 0.155755538233442, 0.373906913054041, 0.337253446353939, 
0.7097912469153, 0.135931197003286, 0.933045789684956, 0.990237264379879, 
0.935909893802418, 0.0584744261368051, 0.142195186247945, 0.0778947823206183, 
0.730092948719125, 0.990958934884833, 0.430560076395454, 0.416003595506252, 
0.329657988797774, 0.923642326109456, 0.768435704916231, 0.199460363021509, 
0.0447731691906169, 0.754509722206595, 0.213058588668422, 0.176568550116404, 
0.903394659596492, 0.229033023417391, 0.410922519487293, 0.327905759041092, 
0.722040727161857, 0.178261568258274, 0.659604652931138)), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -37L))

我需要用R绘制散点图，同时添加最小二乘回归线和基于Theil-Sen斜率的回归线，并且保留x轴日期的mm/dd/yyyy格式。但直接使用POSIXct类型日期计算Theil-Sen斜率时出现问题，转为数值型后回归线会偏离绘图范围。现有脚本如下：

library(trend)
library(ggplot2)
library(zyp)

#make junk1 a time series object for mk.test
xy<-as.ts(junk1)

#perform mann-kendall test, get pval and estimates
mktest<-mk.test(xy)
pVal<-round(mktest$pvalg,4)
estimates<-mktest$estimates

#get sens slope
sSlope<-zyp.sen(junk1$Result~junk1$Date)

#fabricate regression line with sen's slope and the data from junk1
#y = mX + b
mkr<-(sSlope$coefficients[2]*junk1$Date)+sSlope$coefficients[1]

#plot 1) the data in junk 1, 2)least squares regression with lm function, 3) regression with sen's slope
ggplot()+
  geom_point(aes(x = junk1$Date,y = junk1$Result,color = "junk1"))+
  geom_smooth(method = lm, fullrange = TRUE, se = FALSE,aes(x=junk1$Date, y=junk1$Result, color = "LR"))+
  geom_line(aes(x = junk1$Date, y = mkr, color ="MKR"))+
  scale_color_manual(name=" ", values = setNames(c("black", "red", "blue"), c("junk1", "LR", "MKR")))

解决方案

问题出在手动计算Theil-Sen回归线时，POSIXct类型的日期以秒为单位存储为数值，导致截距数值异常，进而让绘制的线偏离范围。可以直接利用ggplot2的stat_smooth结合zyp.sen方法自动拟合回归线，同时设置x轴日期格式：

library(trend)
library(ggplot2)
library(zyp)
library(scales)

# 加载数据
junk1 <- structure(list(Date = structure(c(1704067200, 1704153600, 1704240000, 
1704326400, 1704412800, 1704499200, 1704585600, 1704672000, 1704758400, 
1704844800, 1704931200, 1705017600, 1705104000, 1705190400, 1705276800, 
1705363200, 1705449600, 1705536000, 1705622400, 1705708800, 1705795200, 
1705881600, 1705968000, 1706054400, 1706140800, 1706227200, 1706313600, 
1706400000, 1706486400, 1706572800, 1706659200, 1706745600, 1706832000, 
1706918400, 1707004800, 1707091200, 1707177600), class = c("POSIXct", 
"POSIXt"), tzone = "UTC"), Result = c(0.895103686169796, 0.0556219135717945, 
0.807153205141622, 0.461672184316161, 0.547758011882109, 0.927777651259266, 
0.308058477041002, 0.155755538233442, 0.373906913054041, 0.337253446353939, 
0.7097912469153, 0.135931197003286, 0.933045789684956, 0.990237264379879, 
0.935909893802418, 0.0584744261368051, 0.142195186247945, 0.0778947823206183, 
0.730092948719125, 0.990958934884833, 0.430560076395454, 0.416003595506252, 
0.329657988797774, 0.923642326109456, 0.768435704916231, 0.199460363021509, 
0.0447731691906169, 0.754509722206595, 0.213058588668422, 0.176568550116404, 
0.903394659596492, 0.229033023417391, 0.410922519487293, 0.327905759041092, 
0.722040727161857, 0.178261568258274, 0.659604652931138)), class = c("tbl_df", 
"tbl", "data.frame"), row.names = c(NA, -37L))

# 执行Mann-Kendall测试
xy <- as.ts(junk1)
mktest <- mk.test(xy)
pVal <- round(mktest$pvalg, 4)
estimates <- mktest$estimates

# 绘制图形
ggplot(junk1, aes(x = Date, y = Result)) +
  geom_point(color = "black", aes(label = "junk1")) +
  # 添加最小二乘回归线
  stat_smooth(method = lm, fullrange = TRUE, se = FALSE, color = "red", aes(label = "LR")) +
  # 添加Theil-Sen回归线
  stat_smooth(method = zyp.sen, fullrange = TRUE, se = FALSE, color = "blue", aes(label = "MKR")) +
  # 设置x轴日期格式为mm/dd/yyyy
  scale_x_datetime(labels = date_format("%m/%d/%Y")) +
  # 设置图例
  scale_color_manual(name = "", values = c("black", "red", "blue"), breaks = c("junk1", "LR", "MKR")) +
  labs(x = "日期", y = "结果值") +
  theme_bw()

关键说明

避免手动计算回归线：直接使用stat_smooth(method = zyp.sen)，让ggplot2自动处理POSIXct日期的拟合与绘图，无需手动转换数值，避免截距异常导致的线偏离问题。
x轴日期格式设置：通过scale_x_datetime(labels = date_format("%m/%d/%Y"))指定日期显示格式，需要加载scales包。
数据映射优化：将数据框直接传给ggplot的data参数，避免在aes中直接引用全局变量，让代码更规范。

内容的提问来源于stack exchange，提问作者Brian