请围观：关于Mann-Kendall检验结果的分析(附R分析代码)

zhouyuanshen

我按照魏凤英老师书上介绍做了M-K分析，但是我的结果好像不好解释，麻烦各位帮我支支招。
左边是原始数据序列，右边是M-K分析结果。
我的分析结果显示，UF和UB大部分时间都在置信区间外，而且大很多，这正常吗？

魏老师说UK或者UB大于0有上升趋势，但是我的数据大部分时间都是一个大于0，一个小于0，这又怎么解释呢？

UF和UB的交点大部分都在置信区间外，是不是说明我的数据序列可能就不存在突变点？

从我的M-K分析结果中，能看出来哪个时期变化显著吗？哪个时期变化平缓吗？

谢谢啦

R代码如下：

1. Mann_Kendall <- function(timeserial,start){
   
2.   Mann_Kendall_sub <- function(timeserial){
   
3.     r <- c()
   
4.     s <- c()
   
5.     U <- c()
   
6.    
   
7.     for(i in 2:length(timeserial))
   
8.     {r[i] <- 0
   
9.      
   
10.      for(j in 1:i)
    
11.      {
    
12.        if(timeserial[i]>timeserial[j]){r[i]=r[i]+1}
    
13.      }
    
14.      
    
15.      
    
16.      s[i] <- 0
    
17.      for (ii in 2:i){
    
18.        s[i] <- s[i]+r[ii]
    
19.      }
    
20.      
    
21.      
    
22.      U[i] <- 0
    
23.      U[i] <- (s[i]- (i*(i-1)/4))/sqrt(i*(i-1)*(2*i+5)/72)
    
24.      
    
25.     }
    
26.    
    
27.     r[1] <- 0
    
28.     s[1] <- 0
    
29.     U[1] <- 0
    
30.    
    
31.     LST <- list(r = r, s = s, U = U)
    
32.    
    
33.     return (LST)
    
34.   }
    
35.   timeserial_rev <- rev(timeserial)
    
36.   
    
37.   y1 <- Mann_Kendall_sub(timeserial)
    
38.   y2 <- Mann_Kendall_sub(timeserial_rev)
    
39.   
    
40.   y2$U <- -(rev(y2$U))
    
41.   
    
42.   LST <- list(UF=y1,UB=y2)
    
43.   return(LST)  
    
44. }
    
45. 
    
46. 
    
47. 
    
48. #这里是你要修改的地方
    
49. setwd("d:/")
    
50. od <- read.table("1.txt", header=T)
    
51. Variable <- c("Jan","Feb","Mar","Apr","May","Jun")
    
52. #修改上面三段代码
    
53. 
    
54. 
    
55. 
    
56. 
    
57. startyear <- as.numeric(od[1,1])
    
58. years <- od$Year
    
59. 
    
60. tiff("filename.tif", width=14.6, height=16, units="cm", res=300, family = 'serif')
    
61. par(mfrow=c(6,2),oma=c(3,0,0,0), mar=c(0,2,0,0),cex=0.7)
    
62. 
    
63. for(i in 1:length(Variable)){
    
64.   
    
65.   name <- paste("od$",Variable[i], sep="")
    
66.   value <- eval(parse(text=name))
    
67.   
    
68.   plot(value,type="l", ylab=Variable[i], cex.axis=0.6,xaxt="n",mgp=c(1,0.1,0),tck=-0.02)
    
69.   if(i==length(Variable)){
    
70.     axis(side=1, at=years, tck=-0.04, hadj=0.4, labels=years,mgp=c(1,0.4,0), cex.axis=1) # add x-axis to the last figure
    
71.     axis(side=1, at=1:length(od$Year), tck=-0.01, hadj=0.4, labels=NA, cex.axis=1) # add month labels to the x-axis
    
72.     mtext("年份",side=1,line=1.5)
    
73.   }
    
74.   
    
75.   
    
76.   d<-Mann_Kendall(value,startyear)
    
77.   yUF <- as.data.frame(d$UF[3])$U
    
78.   yUB <- as.data.frame(d$UB[3])$U
    
79.   
    
80.   plot(x=c(1:length(od$Year)),y=yUF, type="l", ylim=c(min(yUF,yUB,-1.96),max(yUF,yUB,1.96)),lwd=1, lty=5, ylab="", cex=0.5,xaxt="n",mgp=c(1,0.1,0),tck=-0.02)
    
81.   points(x=c(1:length(od$Year)),y=yUB,type="l",lty=3,col=6,lwd=1)
    
82.   abline(h=1.969,lty=4,lwd=0.5)
    
83.   abline(h=-1.96,lty=4,lwd=0.5)
    
84.   abline(h=0,col="gray",lwd=0.5)
    
85. }
    
86. 
    
87. axis(side=1, at=years, tck=-0.04, hadj=0.4, labels=years,mgp=c(1,0.4,0), cex.axis=1) # add x-axis to the last figure
    
88. axis(side=1, at=1:length(od$Year), tck=-0.01, hadj=0.4, labels=NA, cex.axis=1) # add month labels to the x-axis
    
89. mtext("年份",side=1,line=1.5)
    
90. dev.off()

复制代码

zhouyuanshen

愚遥 发表于 2015-7-24 13:46
代码有问题啊。
第23行
应该是这样才对

谢谢您的提醒，我对着魏老师的书看了一遍，好像是输错了。在这里给大家道歉了，希望对大家的影响不大。
我试了一下，用魏老师的数据和我自己的数据重新运行了，错误代码的结果似乎比较保守。。。。。下面两张图，左边是经过您提醒后修改后的代码，右边是原来的错误代码

魏凤英老师数据

自己的数据

再次谢谢您

way_liu

电视机 发表于 2016-12-12 16:26
为什么你的和魏老师的图不一样呢？

我用Matlab程序测试了你提供的魏老师的数据，做出来的结果和魏老师的一致。现在把这个matlab程序传上来，供大家分享吧。

romeo19850806

zhouyuanshen 发表于 2014-5-26 23:23
我是用R做的，有需要的我可以共享代码

您能帮我看看这个数据怎么做MK检验吗？ 1.xlsx (18.85 KB, 下载次数: 30)

2014-5-27 08:06 上传

点击文件名下载附件

高地不保

zhouyuanshen 发表于 2015-5-23 18:58
论坛的好友都在回复里面了，您可以浏览一下。大部分在第一页

可是看的不是很明白，我想知道你最后是不是用MK做的突变分析，还是用其他的方法做的，求回复

278803532

不在区间内的，楼主要检验是否是突变点，因为有可能是，也有可能不是，我记得有一片文献上有说，具体我忘记了。按照你变化趋势图和MK中的UF曲线看结果是没有错的，具体原因说不清楚。根据魏老师的书，你主要是看UF线，超出检验的是变化显著的。看上升或下降趋势，你要看一条线，而不是两条线，UF>0,UB<o是上升趋势，反之是下降的，这个可以根据你趋势图看，结果是对应的。

Kongxiangye

有些文献里说，其实这个突变检验只是个方法而已，不一定准~~~可以用其他突变检验方法试试~

zhouyuanshen

278803532 发表于 2014-5-25 07:39
不在区间内的，楼主要检验是否是突变点，因为有可能是，也有可能不是，我记得有一片文献上有说，具体我忘记 ...

虽然有点拗口，但是还是非常明白你的意思，QQQ

zhouyuanshen

Kongxiangye 发表于 2014-5-25 09:39
有些文献里说，其实这个突变检验只是个方法而已，不一定准~~~可以用其他突变检验方法试试~

恩，好像文献是这么说的，我也只是作为一个辅助，不是重点研究突变。

romeo19850806

用什么软件做的？？？？？？

zhouyuanshen

romeo19850806 发表于 2014-5-26 17:14
用什么软件做的？？？？？？

我是用R做的，有需要的我可以共享代码

zhouyuanshen

romeo19850806 发表于 2014-5-27 08:06
您能帮我看看这个数据怎么做MK检验吗？

我只画了一月到6月的图，你对着代码修改修改就能画剩下的了
只能帮你到这啦，

木林森(lynn)

过来学习一下