R –第四部分:趋势处理入门

嘿,

welcome to part 4 of our short introduction to 极值分析 using the extRemes package in R.

在最近的帖子中 块最大值法阈值超额法,我们仅假设满足了所有极值分析的假设。但是,在处理环境变量时,情况极有可能并非如此。特别是在许多情况下,可能会违反平稳性的假设。在全球气候变化的背景下,气象或其他环境变量的时间序列可能有相当大的趋势。当然,这种趋势必须纳入分析中,因为最终的回报水平会随着时间而变化。

The following code shows a short practical example of fitting a generalized pareto distribution to a time series of precipitation 数据 using the extRemes package in R. The sample 数据 set features precipitation 数据 of Bärnkopf (Lower Austria) from 1971 through 2013 and is provided by the hydrographic services of Austria via via H.

# load required packages
library(extRemes)
library(xts)
library(ggplot2)

# get 数据 from H
ehyd_url <- "http://ehyd.gv.at/eHYD/MessstellenExtraData/nlv?id=107540&file=2"
precipitation_xts <- read_ehyd(ehyd_url)

# derive AMS for maximum precipitation
ams <- apply.yearly(precipitation_xts, max)

# check stationarity within AMS:
# Mann-Kendall 趋势 test
Kendall::MannKendall(ams)

# simple linear model
ams_df <- fortify(ams)
colnames(ams_df)[1] <- "date"
summary(lm(ams~date, 数据=ams_df))

p <- ggplot(ams_df, aes(x = date, y = ams)) + geom_point() + geom_line()
p + stat_smooth(method = "lm", formula = y ~ x, size = 1)

趋势

拟合线性模型的结果和该图给出的图形印象都表明年最大降水量呈上升趋势。 Mann-Kendall趋势检验返回非常小的p值,从而确认了该趋势。因此,必须应用趋势校正以解决随时间变化的回报水平。

# maximum likelihood estimation
mle_trend <- fevd(x = ams, 数据 = ams_df, location.fun=~date, method = "MLE", type="GEV")
rl_trend <- return.level(mle_trend, conf = 0.05, return.period= c(2,5,10,20,50,100))

# return level plot
plot(mle_trend, type="rl", main="Return Level Plot for Bärnkopf w/ MLE")

return_levels_trend
与以前帖子的返回水平图(无趋势)相比,此返回水平图看起来有所不同。它显示5年和100年回报水平随时间的变化。

参考:Coles,Stuart(2001)。极值统计建模简介。统计中的Springer系列。伦敦:施普林格。

关于作者

马蒂亚斯在维也纳自然资源与生命科学大学学习了环境信息管理,并获得了环境统计博士学位。他的论文的重点是罕见(极端)事件的统计建模,作为对关键基础设施进行漏洞评估的基础。他目前在奥地利国家气象和地球物理服务局(ZAMG)和BOKU大学山区风险工程研究所工作。他目前专注于(统计)不良天气事件和自然灾害以及减少灾害风险的评估。他的主要兴趣是环境现象的统计建模以及用于数据科学,地理信息和遥感的开源工具。

9条留言

您可以在这篇文章中发表评论。


  • 非常好,非常有帮助!

    您如何获得返回水平x轴的标签“Year”?即使我运行您的代码,我也会得到“index.”
    提前致谢!

    麦克风 5年前 回复


    • Try adding xlab="年" inside the plot function.

      马丁 5年前 回复


      • That makes sense, but something in fevd doesn’t like that.

        plot(mle.trend, type="rl", ylim=c(0,200),
        +      main="Return Level Plot for Bärnkopf w/ MLE",xlab="年")
        Error in plot.default(y, type = "l", xlab = "index", ylab = ylb, ...) : 
          formal argument "xlab" matched by multiple actual arguments
        

        同样,您完整的示例(数据,代码和图形)对像我这样的初学者很有帮助。

        麦克风 5年前 回复


        • As a matter of fact, I am using customized plotting functions for fevd objects. In this particular case, it’s not that easy to modify the plots, as the fevd function is somehow bloated, and so are its plotting functions. So, even though the code for the plots is quite simple, you have to find the correct position where you have to make your modifications. In your code, the xlab-definition gets passed on to a plot()-function where the argument xlab is already defined as xlab = "index". If you want to customize your plot, you have to modify the function plot.fevd.mle. Around line 550 you will find the actual plot function that gets called in this specific case. You can change type = "b", xlab = "年", pch = 16 to get a similar plot to the one I used in this post. In addition, you can dig into this plotting function if you want to change colors, return levels or plot multiple fitting methods in one plot (e.g. comparison of MLE and L-moments in one plot).

          马蒂亚斯 5年前 回复


  • I’使其起作用,并了解其重要性,但是我’我对如何消除趋势有些迷惑。换句话说,究竟是什么从数据中删除了趋势,还是只是随时间将趋势应用于收益水平。

    卡梅伦 3年前 回复


  • 获得的时间为$ t $的趋势校正估计值$ z $为$ \ hat {z} _ {t} = y_ {t}-\ hat {y} _ {t} + \ hat {y} _ {l } $,其中$ y_ {t} $是时间$ t $的度量,$ \ hat {y} _ {t} $是时间$ t $从线性趋势模型$ \ hat {y}获得的趋势_ {t} = \ beta_ {0} + \ beta_ {1} t $,截距$ \ beta_ {0} $和斜率$ \ beta_ {1} $和$ \ hat {y} _ {l} $是最近一年的趋势估计。

    这意味着您(1)进行每次测量,(2)减去通过线性回归获得的该测量值的估计值,以及(3)添加通过线性回归获得的时间序列中最后一年的估计值。

    马蒂亚斯 3年前 回复


  • 嗨,您好
    I’我的命令有问题“read_ehyd(ehyd_url)”.
    您能解释一下该命令吗? ðŸ™,
    提前致谢。

    托马斯·汉森(Tomas Hansson) 3年前 回复


    • 嗨,托马斯,

      这是一个自定义功能,可以从奥地利的数字水文档案中读取数据。
      在帖子中有解释 使用R从eHYD读取数据.

      问候,
      马蒂亚斯

      马蒂亚斯 3年前 回复



      • 谢谢,马蒂亚斯。

        问候,
        托马斯

        托马斯·汉森(Tomas Hansson) 3年前 回复


发表回复

*