# Detrend Example

library(knitr)
library(nimble)
library(earlywarning)
library(ggplot2)
library(tidyr)
opts_chunk$set(dev='png', fig.width=5, fig.height=5, results='hide') some sample data from earlywarning: set.seed(123) data(ibms) plot(ibm_critical) raw <- as.data.frame(ibm_critical) names(raw) <- "x" Rather than explicitly modeling the trend element predicted by the linearization, let us simply remove it: N <- length(raw$x)
raw$t <- 1:N detrend <- loess(x ~ t, raw) data <- data.frame(x = detrend$residuals/sqrt(var(detrend$residuals))) qplot(raw$t, data$x, geom='line') ## LSN version Modify the LSN model to explicitly model the changing parameter as a hidden, stochastic variable lsn <- nimbleCode({ theta ~ dunif(-100.0, 100.0) sigma_x ~ dunif(1e-10, 100.0) sigma_y ~ dunif(1e-10, 100.0) m ~ dunif(-1e2, 1e2) x[1] ~ dunif(-100, 100) y[1] ~ dunif(-100, 100) for(i in 1:(N-1)){ mu_x[i] <- x[i] + y[i] * (theta - x[i]) x[i+1] ~ dnorm(mu_x[i], sd = sigma_x) mu_y[i] <- y[i] + m * t[i] y[i+1] ~ dnorm(mu_y[i], sd = sigma_y) } }) Constants in the model definition are the length of the dataset, $$N$$ and the time points of the sample. Note we’ve made time explicit, we’ll assume uniform spacing here. constants <- list(N = N, t = raw$t)

Initial values for the parameters

inits <- list(theta = 6, m = 0, sigma_x = 1, sigma_y = 1, y = rep(1,N))

and here we go:

Rmodel <- nimbleModel(code = lsn,
constants = constants,
data = data,
inits = inits)
Cmodel <- compileNimble(Rmodel)
mcmcspec <- configureMCMC(Rmodel, print=TRUE,thin=2e2)
Rmcmc <- buildMCMC(mcmcspec)
Cmcmc <- compileNimble(Rmcmc, project = Cmodel)
Cmcmc$run(1e6) NULL and examine results samples <- as.data.frame(as.matrix(Cmcmc$mvSamples))
dim(samples)
[1] 5000   84
samples <- samples[,1:4]
long <- gather(samples)
apply(samples, 2, mean)
            m       sigma_x       sigma_y         theta
0.0003790592  1.0792385676  0.1920288851 -0.0150533955 
ggplot(long) +
geom_line(aes(seq_along(value), value)) +
facet_wrap(~key, scale='free')
ggplot(long) +
geom_density(aes(value)) +
facet_wrap(~key, scale='free')
sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)

locale:
[1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8       LC_NAME=C
[9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] methods   stats     graphics  grDevices utils     datasets
[7] base

other attached packages:
[1] tidyr_0.2.0        ggplot2_1.0.0      earlywarning_0.0-1
[4] nimble_0.3         yaml_2.1.13        knitr_1.9

loaded via a namespace (and not attached):
[1] codetools_0.2-10 colorspace_1.2-4 deSolve_1.11
[4] digest_0.6.8     evaluate_0.5.5   formatR_1.0
[7] grid_3.1.2       gtable_0.1.2     igraph_0.7.1
[10] labeling_0.3     MASS_7.3-39      mnormt_1.5-1
[13] munsell_0.4.2    parallel_3.1.2   plyr_1.8.1
[16] proto_0.3-10     psych_1.5.1      Rcpp_0.11.4
[19] reshape2_1.4.1   scales_0.2.4     stringr_0.6.2
[22] tools_3.1.2