sens.slope () против mblm () в R - PullRequest
0 голосов
/ 09 мая 2020

Я работал над оценкой Тейла Сена в R и наткнулся на две функции, которые давали бы чувствительный наклон, однако они дают разные значения. Кто-нибудь знает, почему это так? Первый метод, который я использовал из sens.slope(), и второй метод, который я использовал из mblm(). В обоих методах использовался один и тот же фрейм данных

tm14D<-data.frame(date=as.Date(date),temperature=temperature,stringsAsFactors=FALSE)
tm14D$date<-as.numeric(tm14D$date) #date must be changed to numerical values to avoid difftime error
tm14D.zs<-sens.slope(as.numeric(tm14D$temperature)) #method 1
tm14D.ts_fit<-mblm(temperature~date,data=tm14D)     #method 2

И вот результаты:
Метод 1
tm14D.zs$estimates
Наклон Сена
0,0471944

Метод 2
ts_fittm14D$coefficients[2]
дата
0,0001055271

Какой правильный, а какой один не прав? Или мне чего-то не хватает?

1 Ответ

0 голосов
/ 09 мая 2020

Также есть zyp.sen из пакета zyp. zyp.sen и sens.slope похоже согласны.

library(trends)
library(zyp)
library(mblm)
set.seed(1)
df <- data.frame(x = 1:100,y = 1:100-runif(100,-20,20))

coef(zyp.sen(y~x,df))["x"]
       x 
1.003835 

coef(mblm(y ~ x, df),repeated = FALSE)['x'] 
       x 
1.031631 

sens.slope(ts(df)[,-1])$estimate
Sen's slope 
   1.003835

Если мы рассмотрим источник mblm, мы обнаружим, что есть ошибка и что он возвращает результаты повторного метода независимо от аргумента repeated =.

coef(mblm(y ~ x, df),repeated = TRUE)['x'] 
       x 
1.031631 

Вот отрывок из кода mblm:

    x = df$x
    y = df$y
#    if (length(term) > 2) {
#        stop("Only linear models are accepted")
#    }
    xx = sort(x)
    yy = y[order(x)]
    n = length(xx)
    slopes = c()
    intercepts = c()
    smedians = c()
    imedians = c()
#    if (repeated) {
#        <snip>
#    }
#    else {
        for (i in 1:(n - 1)) {
            for (j in i:n) {
                if (xx[j] != xx[i]) {
                  slopes = c(slopes, (yy[j] - yy[i])/(xx[j] - 
                    xx[i]))
                }
            }
        }
        slope = median(slopes)
        intercepts = yy - slope * xx
        intercept = median(intercepts)
#    }
slope
[1] 1.003835

И, таким образом, все согласны.

...