Skip to content

Instantly share code, notes, and snippets.

@dholstius
Created Sep 23, 2014
Embed
What would you like to do?
Theil-Sen regression
# Ported by David Holstius <dholstius@baaqmd.gov>
# from http://skip.ucsc.edu/leslie_MOUSE/programs/plotting/tsreg.r
# Module tsreg
# Author: E. A. Houseman
# Last update July 2004
# AR(q) time series regression assuming regular intervals
# Support for cholesky residuals [Houseman, Ryan, Coull (2004, JASA)]
tsreg <- function (model, data, ord = 1, maxiter = 10, tol = sqrt(.Machine$double.eps)){
y <- model.extract(model.frame(model,data),"response") # Get response
X <- model.matrix(model,data) # Get covariate matrix
dd <- dim(X)
n <- dd[1]
p <- dd[2]
beta <- solve(t(X) %*% X, t(X) %*% y)
e <- y - as.vector(X %*% beta)
if(ord == 0){
sigma <- sqrt( sum(e*e)/(n-p) )
V.beta <- sigma * solve(t(X) %*% X)
V.Sigma <- matrix(2 * sigma^4 * n/(n-p)^2, 1, 1)
o <- list(beta = beta, rho = numeric(0), iters = 0,
ord = ord, V.beta = V.beta, V.Sigma = V.Sigma,
R.elems = NULL, dR.elems = NULL,
sigma = sigma, y = y, X = X, shifts = NULL, call = match.call())
class(o) <- "tsreg"
return(o)
}
shifts <- tsreg.makeShifts(n)
Z <- tsreg.makeZ(e,ord,n)
rho0 <- rho <- solve(t(Z) %*% Z, t(Z) %*% e[-(1:ord)])
for(i in 1:maxiter){
R.elems <- tsreg.makeCorrElems(rho, shifts)
R <- tsreg.makeCorr(R.elems)
R.qr <- qr(R)
Xt <- t(solve(R.qr,X))
XtX <- Xt %*% X
beta <- solve(XtX, Xt %*% y)
e <- y - as.vector(X %*% beta)
Z <- tsreg.makeZ(e,ord,n)
rho0 <- rho
ZtZ <- t(Z) %*% Z
rho <- solve(ZtZ, t(Z) %*% e[-(1:ord)])
if(max(abs(rho-rho0))<tol) break;
}
sigma <- as.vector(sqrt(e %*% solve(R.qr,e)/n))
sigma2 <- sigma*sigma
bnames <- rownames(beta)
beta <- as.vector(beta)
names(beta) <- bnames
V.beta <- sigma2*solve(XtX)
dR.elems <- tsreg.makeDCorrElems(rho, shifts)
RiR. <- list()
Info <- matrix(0, ord+1, ord+1)
for(j in 1:ord){
R. <- tsreg.makeCorr(dR.elems[[j]])
diag(R.) <- 0
RiR.[[j]] <- solve(R.qr, R.)
for(k in 1:j){
if(j == k) Info[j,j] <- sum(diag( RiR.[[j]] %*% RiR.[[j]] ))
else Info[j,k] <- Info[k,j] <- sum(diag(RiR.[[k]] %*% RiR.[[j]]))
}
Info[j,ord+1] <- Info[ord+1,j] <- sum(diag(RiR.[[j]])) / sigma2
}
Info[ord+1,ord+1] <- n / sigma2 / sigma2
V.Sigma <- 2 * solve(Info)
obj <- list(beta = beta, rho = as.vector(rho), iters = i,
ord = ord, V.beta = V.beta, V.Sigma = V.Sigma,
R.elems = R.elems, dR.elems = dR.elems,
sigma = sigma, y = y, X = X, shifts = shifts, call = match.call())
class(obj) <- "tsreg"
return(obj)
}
# Redefine chol as a method
chol <- function (x, ...) UseMethod("chol")
chol.matrix <- base::chol
# Summary for tsreg
summary.tsreg <- function (o, alpha = 0.05){
beta <- cbind(est = o$beta, SE = sqrt(diag(o$V.beta)))
Sigma <- cbind(est = c(o$rho,o$sigma^2), SE = sqrt(diag(o$V.Sigma)))
crit <- abs(qnorm(alpha/2))
beta <- cbind(beta, beta[,1] - crit*beta[,2])
beta <- cbind(beta, beta[,1] + crit*beta[,2])
Sigma <- cbind(Sigma, Sigma[,1] - crit*Sigma[,2])
Sigma <- cbind(Sigma, Sigma[,1] + crit*Sigma[,2])
sig.parm <- dim(Sigma)[1]
# Base CI's for sigma^2 on log transformation
Sigma[sig.parm,3] <- exp(log(Sigma[sig.parm,1])-crit*Sigma[sig.parm,2]/Sigma[sig.parm,1])
Sigma[sig.parm,4] <- exp(log(Sigma[sig.parm,1])+crit*Sigma[sig.parm,2]/Sigma[sig.parm,1])
colnames(beta)[3:4] <- paste((1-alpha)*100,"% CI ",c("lo","hi"),sep = "")
colnames(Sigma)[3:4] <- paste((1-alpha)*100,"% CI ",c("lo","hi"),sep = "")
if(o$ord == 0) rownames(Sigma) <- "sigma^2"
else rownames(Sigma) <- c(paste("rho[",1:length(o$rho),"]",sep = ""), "sigma^2")
so <- list(call = o$call, beta = beta, Sigma = Sigma)
class(so) <- "tsregsummary"
so
}
print.tsreg <- function (fit){
cat("Output from 'tsreg'\n")
print(fit$call)
cat("\nIterations:",fit$iter,"\n")
cat("\nRegression parameters:\n")
print(fit$beta)
cat("\nCorrelation parameters:\n")
print(fit$rho)
cat("\nCovariance scale (sigma^2)\n")
print(fit$sigma^2)
}
print.tsregsummary <- function (sfit){
cat("Summary of output from 'tsreg'\n")
print(sfit$call)
cat("\nRegression parameters:\n")
print(sfit$beta)
cat("\nCovariance parameters:\n")
print(sfit$Sigma)
cat("\n")
}
# Cholesky residuals for tsreg solution
chol.tsreg <- function (o){
if(o$ord == 0) R <- diag(length(o$y))
else R <- tsreg.makeCorr(o$R.elems)
R.qr <- qr(R)
R.inv <- solve(R.qr)
Lt <- chol(R.inv) / o$sigma
z <- as.vector(Lt %*% (o$y - o$X %*% o$beta))
ord <- o$ord
p <- length(o$beta)
n <- length(z)
ds2 <- matrix(0,n,ord+1)
ds2[,ord+1] <- -diag(R.inv)/o$sigma^2
if(ord>0) for(r in 1:ord){
R. <- tsreg.makeCorr(o$dR.elems[[r]])
diag(R.) <- 0
ds2[,r] <- - diag(Lt %*% R. %*% t(Lt))
}
ds2 <- ds2/o$sigma^2
dm <- Lt %*% o$X
aux <- new.env()
assign("dm", dm, env = aux)
assign("ds2", ds2, env = aux)
assign("V.beta", o$V.beta, env = aux)
assign("V.Sigma", o$V.Sigma, env = aux)
attr(z,"aux") <- aux
class(z) <- "tsregchol"
z
}
# Summary of Cholesky residuals for getting confidence bands
summary.tsregchol <- function (z, x = z, alpha = 0.05){
p <- pnorm(x)
n <- length(z)
dm <- apply(get("dm", env = attr(z,"aux")), 2, mean)
ds2 <- apply(get("ds2", env = attr(z,"aux")), 2, mean)
V.beta <- get("V.beta", env = attr(z,"aux"))
V.Sigma <- get("V.Sigma", env = attr(z,"aux"))
delta.m <- t(dm) %*% V.beta %*% dm
delta.s2 <- t(ds2) %*% V.Sigma %*% ds2
dnormx <- dnorm(x)
var.correction <- as.vector(dnormx*dnormx*delta.m +
x*x*dnormx*dnormx*delta.s2/4)
var.p <- p*(1-p)/n
se <- sqrt(var.p - var.correction)
z.alpha <- qnorm(alpha/2)
ci.lo <- p - z.alpha * se
ci.hi <- p + z.alpha * se
sort <- order(x)
cbind(p,se,ci.lo,ci.hi,zci.lo = qnorm(ci.lo),zci.hi = qnorm(ci.hi),index = sort)[sort,]
}
# QQ plot with confidence "bands"
qqnorm.tsregchol <- function (z,xlab = NULL,ylab = NULL,col = 2,...){
band <- suppressWarnings( summary(z) )
qqnorm(as.vector(z),xlab = xlab,ylab = ylab)
n <- length(z)
qz <- qnorm(band[,1])
lines(band[,5],qz,col = col,...)
lines(band[,6],qz,col = col,...)
}
qqnorm.tsreg <- function (fit,xlab = NULL,ylab = NULL,col = 2,...){
qqnorm(chol(fit),xlab = xlab,ylab = ylab,col = col,...)
}
#########################
#### SUPPORTING FUNCTIONS
#########################
# Make the covariate matrix for LS regression to get rho parameters
# [see Durbin (1960)]
tsreg.makeZ <- function (e, ord, n){
Z <- matrix(0,n-ord,ord)
for(j in 1:ord){
Z[,j] <- e[(ord-j+1):(n-j)]
}
Z
}
# Get shift matrices for constructing correlations
tsreg.makeShifts <- function (n){
E <- list()
E[[1]] <- diag(n)
for(i in 1:n) {
e <- rep(0,n)
if(i>1) e[i-1] <- 1
E[[i+1]] <- rbind(e,E[[i]][-n,])
}
E
}
# Make elements of a correlation matrix
tsreg.makeCorrElems <- function (rho, shifts){
shiftmat <- tsreg.makeShiftMat(rho, shifts)
solve(shiftmat$E0, shiftmat$Rho)
}
# Sum of shift matrices for constructing correlation matrix
tsreg.makeShiftMat <- function (rho, shifts){
E0 <- shifts[[1]]
n <- length(shifts)-1
Rho <- c(rho, rep(0,n-length(rho)))
for(i in 1:n){
E0 <- E0 - Rho[i]*shifts[[i+1]]
}
list(E0 = E0, Rho = Rho)
}
# Make a correlation matrix
tsreg.makeCorr <- function (elems){
n <- length(elems)
R <- matrix(0,n,n)
for(i in 1:(n-1)){
R[i,i+(1:(n-i))] <- elems[1:(n-i)]
}
R <- R + t(R)
diag(R) <- 1
R
}
# Make elements of correlation derivatives
tsreg.makeDCorrElems <- function (rho, shifts){
ord <- length(rho)
dR.elems <- list()
shiftmat <- tsreg.makeShiftMat(rho, shifts)
E0.qr <- qr(shiftmat$E0)
for(r in 1:ord){
E.dot <- shifts[[r+1]]
dR.elems[[r]] <- as.vector(solve(E0.qr)[,r] +
solve(E0.qr, E.dot) %*% solve(E0.qr, shiftmat$Rho))
}
dR.elems
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment