From: hadley wickham <h.wickham_at_gmail.com>

Date: Sat 13 Jan 2007 - 19:59:09 GMT

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Jan 14 07:04:34 2007

Date: Sat 13 Jan 2007 - 19:59:09 GMT

> 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

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun Jan 14 07:04:34 2007

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Sat 13 Jan 2007 - 20:30:38 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*