[R] Magnitude of trend in time series
hadley wickham
h.wickham at gmail.com
Sat Jan 13 20:59:09 CET 2007
> I am analyzing some climate time series data using the Mann Kendall package
> and was wondering if there was a way to calculate the trend using Sen's
> nonparametric estimator slope in R?
I think you can do it with the mblm package. But I vaguely remember
that the results didn't look reasonable. I've included some
simulation/testing code that I used below.
b1 <- 5
x <- 1:20 - mean(1:20)
y <- 0 + b1 * x
n1 <- function() y + rnorm(20, 0, 2)
ln1 <- function() y + log(rnorm(20, 0, 0.7))
ln2 <- function() y + log(rnorm(20, 0, 1.5))
fitlm <- function(f=n1) {
m <- lm(n1() ~ x)
ci <- t(sapply(c(0.9, 0.95, 0.99), function(x) confint(m, level=x)[2,]))
b1 > ci[,1] & b1 < ci[,2]
}
fitts <- function(f=n1) {
x <- x # because mblm mucks up the environments
y2 <- n1()
m <- mblm(y2 ~ x, repeated=FALSE)
ci <- t(sapply(c(0.9, 0.95, 0.99), function(x) confint(m, level=x)[2,]))
b1 > ci[,1] & b1 < ci[,2]
}
n1lm <- replicate(1000, fitlm(n1))
n1ts <- replicate(100, fitts(n1))
ln1lm <- replicate(1000, fitlm(ln1))
ln1ts <- replicate(100, fitts(ln1))
ln2lm <- replicate(1000, fitlm(ln2))
ln2ts <- replicate(100, fitts(ln2))
Hadley
More information about the R-help
mailing list