# A JAGS calculation on pattern of rain January 1906-1915 against 2003-2012

**Wiekvoet**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Two weeks ago I showed rain data from six stations in Netherlands years 1906 till now. Last week I showed that frequency of days with and without rain differed between December 1906-1915 and December 2003-2012. This week I am considering the same data as t-distributed with a left truncation at 0 mm rain. This can be modeled very easy in JAGS. This model gives the posterior distribution of the amount of rain moved 1 mm of rain upwards and half a mm narrower, but both these intervals include 0.

### The data

### The model

The model is almost off the shelf. Truncation, two t-distributions with different means and standard deviations. Priors there taken from *Data Analysis Using Regression and Multilevel/Hierarchical Models*, Gelman and Hill.

sel1 <- all[(all$year <1916 | all$year>2003)

& all$Month==’Jan’

& all$location==’DE-BILT’

& all$day %in% seq(1,31,by=5),]

sel1$decade <- factor(c('1906-1915','2004-2013')[1+(sel1$year>1950)])

datain <- list(

isCensored = 1-as.numeric(sel1$RD==0),

N = nrow(sel1),

RD = ifelse(sel1$RD==0,NA,sel1$RD/10),

period = 1+(sel1$year>1950)

)

model1 <- function() {

for ( i in 1:N ) {

isCensored[i] ~ dinterval( RD[i] , 0 )

RD[i] ~ dt( mu[period[i]] , tau[period[i]],nu )

}

for (i in 1:2) {

tau[i] <- 1/pow(sigma[i],2)

sigma[i] ~ dlnorm(sigmamu,sigmatau)

mu[i] ~ dnorm(mumu,mutau)

}

mumu ~ dnorm(0,1e-6)

mutau <- 1/pow(musigma,2)

musigma ~ dunif(0,100)

sigmamu ~ dnorm(0,1e-6)

sigmatau <- 1/pow(sigmasigma,2)

sigmasigma ~ dunif(1,100)

mudif <- mu[2]-mu[1]

sigmadif <- sigma[2]-sigma[1]

nu <- 1/nuinv

nuinv ~ dunif(0,.5)

}

params <- c('mu','sigma','mudif','sigmadif','nu')

inits <- function() {

list(mumu = rnorm(1,0),

musigma = runif(1,1,2),

sigmamu = rnorm(1,0),

sigmasigma = runif(1,1,2),

nuinv = runif(1,0,.5))

}

jagsfit <- jags(datain, model=model1, inits=inits,

parameters=params,progress.bar=”gui”,

n.chains=1,n.iter=10000,n.burnin=3000,n.thin=3)

jagsfit

### Interpretation

**leave a comment**for the author, please follow the link and comment on their blog:

**Wiekvoet**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.