Sunday, December 14, 2014

Monthly Weather in Netherlands

When I downloaded the KNMI meteorological data, the intention was to do something which takes more than just the computers memory. While it is clearly not big data, at the very least 100 years of daily data is not small either. So I took along a load of extra variables to see what trouble I would run into. I did not run out of memory, but did make some figures.

Data

Data are acquired from KNMI. They have various sets of data, this page has a selection form which leads to the data used today. The data comes with a header explaining details, unfortunately in Dutch.

Plots

Just about everybody knows days are shorter in winter. What I never realized, even within that shorter day, we get less daylight. The short days are often so clouded, we don't get sun, meanwhile, in summer the sun does shine a bigger part of the daylight period.


In real hours sunshine this results in the following plot. December clearly has the darkest hours.
What we do get, not surprising since Dutch weather is fairly similar to English weather, is rain. Not continuous rain, most of the time it is dry, but still, autumn and winter do have days where it does not seem to stop. Autumn has the bad reputation for rain, but this plot makes winter look particular bad.
All this rain gives a load of humidity. This humidity in its turn, gives rise to a weather we name 'waterkoud'. It is above zero C but still quite cold outside. The humidity makes for air with a high heat capacity, hence one cools down quickly. Temperatures below zero can make for much nicer weather, but that can hamper traffic quite a lot. Most of the time it just doesn't freeze.
Code

library(plyr)
library(dplyr)
library(ggplot2)
r1 <- read.csv('KNMI_20141115.edited.txt')
Sys.setlocale(category = "LC_TIME", locale = "C")
r2 <- mutate(r1,
    date = as.Date(format(YYYYMMDD),'%Y%m%d'),
    month =factor(months(date,abbreviate=TRUE),
        levels=months(as.Date(
                paste('2014',
                    formatC(1:12,digits=2,width=2,flag='0'),
                    '01',sep='-')),
            abbreviate=TRUE)),
    yearf=factor(format(date,'%Y')),
    yearn=as.numeric(substr(YYYYMMDD,1,4)),
    day=format(date,'%e'))

g1 <- ggplot(r2,aes(x=month,y=SP))
g1 + geom_violin() +
        ylab('% of longest possible sunshine')

g1 <- ggplot(r2,aes(x=month,y=SQ/10))
g1 + geom_violin()  +
        ylab('Sunshine duration (h)')

g1 <- ggplot(r2,aes(x=month,y=DR/10))
g1 + geom_violin() +
        scale_y_continuous('Precipitation Duration (h)',
                breaks=c(0,6,12,18,24))

g1 <- ggplot(r2,aes(x=month,y=UG))
g1 + geom_violin() +
        ylab('Relative Humidity (%)')
 

g1 <- ggplot(r2,aes(x=month,y=TG/10))
g1 + geom_violin() +
        ylab('Temperature (C)')

Saturday, December 6, 2014

SAS PROC MCMC in R: Nonlinear Poisson Regression Models

In exercise 61.1 the problem is that the model has bad mixing. In the SAS manual the mixing is demonstrated after which a modified distribution is used to fix the model.
In this post the same problem is tackled in R; MCMCpack, RJags, RStan and LaplaceDemon. MCMCpack has quite some mixing problems, RStan seems to do best.

Data

To quote the SAS manual "This example shows how to analyze count data for calls to a technical support help line in the weeks immediately following a product release. (...) You can model the number of daily calls as a Poisson random variable, with the average number of calls modeled as a nonlinear function of the number of weeks that have elapsed since the product’s release. (...) During the first several weeks after a new product is released, the number of questions that technical support receives concerning the product increases in a sigmoidal fashion. The expression for the mean value in the classic Poisson regression involves the log link. There is some theoretical justification for this link, but with MCMC methodologies, you are not constrained to exploring only models that are computationally convenient. The number of calls to technical support tapers off after the initial release, so in this example you can use a logistic-type function to model the mean number of calls received weekly for the time period immediately following the initial release."
observed <- scan(text='
1 0 1 2 2 2 2 1 3 1 3 3
4 5 4 8 5 5 5 9 6 17 6 9
7 24 7 16 8 23 8 27',
    what=list(integer(),integer()),
    sep=' ',
)
names(observed) <- c('weeks','calls')
observed <- as.data.frame(observed)

Analysis

MCMCpack

The MCMCpack solution is derived from LaplacesDemon solution below. It is placed as first because it shows some of the problems with the mixing. As a change from LaplacesDemon, gamma could get negative, for which I made a -Inf likelihood.
library(MCMCpack)
MCMCfun <- function(parm) {
    names(parm) <- c('alpha','beta','gamma')
    if (parm['gamma']<0) return(-Inf)
    mu <-parm['gamma']*
        LaplacesDemon::invlogit(parm['alpha']+parm['beta']*observed$weeks)
    LL <- sum(dpois(observed$calls,mu,log=TRUE))
    LP <- LL+ dgamma(parm['gamma'],shape=3.4,scale=12,log=TRUE) +
        dnorm(parm['alpha'],-5,sd=.25,log=TRUE) +
        dnorm(parm['beta'],0.75,.5,log=TRUE)
    return(LP)
}
mcmcout <- MCMCmetrop1R(MCMCfun,
    c(alpha=-4,beta=1,gamma=2))

The figures show bad mixing. Parameters, especially Beta and Gamma, get stuck. There is quite some autocorrelation.
plot(mcmcout)
acf(mcmcout)
The cause is a nasty correlation between Beta and Gamma
pairs(as.data.frame(mcmcout))

LaplacesDemon

I chose an adaptive algorithm for LaplacesDemon. For the specs, the numbers are derived from the standard deviation of a previous run. Stepsize keeps reasonably constant through the latter part of run. The samples look much better than MCMCpack, although the mixing is not ideal either.
At a later stage I tried this same analysis with reflective Slice Sampler, however, that did was quite a bit more difficult and the end result was not better than this.
library('LaplacesDemon')
mon.names <- "LP"
parm.names <- c('alpha','beta','gamma')
PGF <- function(Data) c(rnorm(3,0,1))
N <-1
MyData <- list(mon.names=mon.names,
    parm.names=parm.names,
    PGF=PGF,
    calls=observed$calls,
    weeks=observed$weeks)
Model <- function(parm, Data)
{
    mu <-parm['gamma']*
        invlogit(parm['alpha']+parm['beta']*Data$weeks)
    LL <- sum(dpois(Data$calls,mu,log=TRUE))
    LP <- LL+ dgamma(parm['gamma'],shape=3.4,scale=12,log=TRUE) +
        dnorm(parm['alpha'],-5,sd=.25,log=TRUE) +
        dnorm(parm['beta'],0.75,.5,log=TRUE)
    Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LP,
        yhat=mu,
        parm=parm)
    return(Modelout)
}

Initial.Values <- c(alpha=-4,beta=1,gamma=2) #GIV(Model, MyData, PGF=TRUE)
Fit1 <- LaplacesDemon(Model,
    Data=MyData,
    Initial.Values = Initial.Values,
    Algorithm = "AHMC",
    Specs = list(epsilon = c(.23,.2,13.5)/4,
        L = 2, Periodicity = 10),
    Iterations=40000,Status=2000   
)
BurnIn <- Fit1$Rec.BurnIn.Thinned
plot(Fit1, BurnIn, MyData, PDF=FALSE)

Jags

I do not think using one chain is advisable, especially since Jags makes more chains so easy. But in the spirit of this analysis I am using one. Coda plots are used since they are a bit more compact.
library(R2jags)
datain <- list(
    calls=observed$calls,
    weeks=observed$weeks,
    n=nrow(observed))
parameters <- c('alpha','beta','gamma')

jmodel <- function() {
    for (i in 1:n) {       
        mu[i] <- gamma*ilogit(alpha+beta*weeks[i])
        calls[i] ~ dpois(mu[i])
    }
    alpha ~ dnorm(-5,1/(.25*.25))
    gamma ~ dgamma(3.4,1/12)
    beta ~ dnorm(.75,1/(.5*.5))
}

jj <- jags(model.file=jmodel,
    data=datain,
    parameters=parameters,
    n.chains=1,
    inits=list(list(alpha=-4,beta=1,gamma=2))
    )

cc <- as.mcmc(jj$BUGSoutput$sims.array[,1,])

plot(cc)
acfplot(cc)

Stan

Stan probably does best handling this typical distribution. Again, one chain is only in the context of this posting.
library(rstan)
smodel <- '
    data {
    int <lower=1> n;
    int calls[n];
    real  weeks[n];
    }
    parameters {
    real Alpha;
    real Beta;
    real Gamma;
    }
    transformed parameters {
    vector[n] mu;
    for (i in 1:n) {
       mu[i] <- Gamma*inv_logit(Alpha+Beta*weeks[i]);
    }
    }
    model {
    calls ~ poisson(mu);
    Gamma ~ gamma(3.4,1./12.);
    Beta ~ normal(.75,.5);
    Alpha ~ normal(-5,.25);
    }
    '
fstan <- stan(model_code = smodel,
    data = datain,
    pars=c('Alpha','Beta','Gamma'),
    chains=1,
    init=list(list(Alpha=-4,Beta=1,Gamma=2)))

traceplot(fstan,inc_warmup=FALSE)

smc <- as.mcmc(as.matrix(fstan))
acfplot(smc)






Sunday, November 30, 2014

Change in temperature in Netherlands over the last century

I read a post 'race for the warmest year' at sargasso.nl. They used a plot, originating from Ed Hawkins to see how 2014 progressed to be warmest year. Obviously I wanted to make the same plot using R. In addition, I wondered which parts of the year had most increased in temperature.

Data

Similar to last week, data are acquired from KNMI. They have various sets of data, this page has a selection form which leads to the data used today. The data comes with a header explaining details, unfortunately in Dutch. Of relevance for this post is TG, average and minimum temperature in 0.1 C. Station used is de Bilt, where they got most data. Prior to reading the data into R, the explanatory text header was removed.
The data input is completed by converting YYYYMMDD to date, year, month and dayno variables. Prior to analysis for simplicity a leap day was removed. I chose merge() rather than a dplyr merge function since the latter releveled my moth factor. The data.frame mylegend is used later on to label the first of the month.
library(plyr)
library(dplyr)
library(ggplot2)
library(locfit)
library(plotrix)
r1 <- read.csv('KNMI_20141115.edited.txt')
Sys.setlocale(category = "LC_TIME", locale = "C")
r2 <- mutate(r1,
    date = as.Date(format(YYYYMMDD),'%Y%m%d'),
    month =factor(months(date,abbreviate=TRUE),
        levels=months(as.Date(
                paste('2014',
                    formatC(1:12,digits=2,width=2,flag='0'),
                    '01',sep='-')),
            abbreviate=TRUE)),
    yearf=factor(format(date,'%Y')),
    yearn=as.numeric(substr(YYYYMMDD,1,4)),
    day=format(date,'%e'))

days <- filter(r2,yearn==1901) %>%
    mutate(.,dayno=format(date,'%j') ) %>%
    select(.,month,day,dayno)

r3 <- merge(r2,days,all=TRUE) %>%
    filter(.,!grepl('0229',YYYYMMDD)) %>%
    mutate(.,daynon=as.numeric(dayno))

mylegend <- filter(days,day==' 1') %>%
    mutate(.,daynon=as.numeric(dayno))

Plots

Cumulative average by year

Each line is a separate year. For this plot is is needed to have year a factor. Unfortunately, I was not able to get a colourbar legend for colors, that required a continuous year variable. Green is beginning last century, pinkish is recent, the fat black line is 2014.

r4 <- group_by(r3,yearf) %>%
    mutate(.,cmtemp = cummean(TG/10))

g1 <- ggplot(r4,aes(x=daynon,y=cmtemp,
        col=yearf))
g1 + geom_line(alpha=.4,show_guide=FALSE)  +
    scale_x_continuous('Day',
        breaks=mylegend$daynon,
        labels=mylegend$month,
        expand=c(0,0)) +
    scale_y_continuous('Temperature (C)')  +
    geom_line(data=r4[r4$yearf=='2014',],
        aes(x=daynon,y=cmtemp),
        col='black',
        size=2)

2014 with average of 30 years

To get a better idea how 2014 compares to previous years, the average of 30 years has been added. We had warm year, except for August, which suggested an early spring. In hindsight, second half of August had colder days than beginning April or end October.


r3$Period <- cut(r3$yearn,c(seq(1900,2013,30),2013,2014),
    labels=c('1901-1930','1931-1960',
        '1961-1990','1991-2013','2014'))
g1 <- ggplot(r3[r3$yearn<2014,],aes(x=daynon,y=TG/10,col=Period))
g1 + geom_smooth(span=.15,method='loess',size=1.5)  +
    scale_x_continuous('Day',
        breaks=mylegend$daynon,
        labels=mylegend$month,
        expand=c(0,0)) +
    geom_line(#aes(x=daynon,y=TG/10),
        data=r3[r3$yearn==2014,]) +
    scale_y_continuous('Temperature (C)')

Change by year

Finally, a plot showing how temperature changed within the years. To obtain this plot, I needed a day corrected base temperature. The baseline temperature is smoothed over days for years 1901 to 1924. The baseline was used to get a corrected baseline, which was subsequently smoothed over years and days.
Smoothers have edge effects, to remove these from the visual part, January and December have been added as extra to the data. Hence within the year there are only minimal edge effects.
The plot shows that middle last century, some parts of the year actually had a drop in temperature. In contrast, November has gradually been getting warmer since middle last century. The new century has seen quite an increase. 
myyears <- r3[r3$yearn<1925,]
m13 <- filter(myyears,daynon<30) %>%
    mutate(.,daynon=daynon+365)
m0 <- filter(myyears,daynon>335) %>%
    mutate(.,daynon=daynon-365)
myyears <- rbind_list(m0,myyears,m13)

nn <- .2
mymod <- locfit(TG ~ lp(daynon,nn=nn),
    data=myyears)
topred <- data.frame(daynon=1:365)
topred$pp <- predict(mymod,topred)
#plot(pp~ daynon,data=topred)

r5 <- merge(r3,topred) %>%
    mutate(.,tdiff=(TG-pp)/10) %>%
    select(.,tdiff,daynon,yearn)
m13 <- filter(r5,daynon<30) %>%
    mutate(.,daynon=daynon+365,
        yearn=yearn-1)
m0 <- filter(r5,daynon>335) %>%
    mutate(.,daynon=daynon-365,
        yearn=yearn+1)
r6 <- rbind_list(m0,r5,m13)

topred <- expand.grid(
    daynon=seq(1:365),
    yearn=1901:2014)
topred$pp2 <- locfit(
        tdiff ~ lp(yearn,daynon,nn=nn),
        data=r6) %>%
    predict(.,topred)
#topred <- arrange(topred,daynon,yearn)

myz <- matrix(topred$pp2,ncol=365)
zmin <- floor(min(topred$pp2)*10)/10
zmax <- ceiling(max(topred$pp2)*10)/10
myseq <- seq(zmin,zmax,.1)
par(mar=c(5,4,4,6))
image(myz,useRaster=TRUE,
    axes=FALSE,frame.plot=TRUE,
    col=colorRampPalette(c('blue','red'))(length(myseq)-1),
    breaks=myseq)
axis((seq(10,114,by=10)-1)/113,labels=seq(1910,2010,by=10),side=1)
axis((mylegend$daynon-1)/365,labels=mylegend$month,side=2)
color.legend(1.1,0,1.2,1,legend=c(zmin,zmax),gradient='y',
    rect.col=colorRampPalette(c('blue','red'))(length(myseq)-1))

Sunday, November 23, 2014

When should I change to snow tires in Netherlands

The Royal Netherlands Meteorological Institute has weather information by day for a number of Dutch stations. In this post I want to use those data for a practical problem: when should I switch to winter tires? (or is that snow tires? In any case nails or metal studs are forbidden in Netherlands). Netherlands does not have laws prescribing winter tires, it does not get bad enough, though there are some regulations in neighboring countries, which is of relevance for people driving to winter sports. For others, it is up to themselves.

Data

Data are acquired from KNMI. They have various sets of data, this page has a selection form which leads to the data used today. The data comes with a header explaining details, unfortunately in Dutch. Of relevance for this post are TG and TN, average and minimum temperature in 0.1 C. Station used is de Bilt, where they got most data.

Analysis 

It seems under 7 C  there is advantage to using winter tires. Temperature varies within the day, but in general the coldest part should be just before dawn which is morning rush hour. Warmest part around noon, with evening rush hour already cooling down. Based on this I decided that days with an average temperature of 7 C or less, or a minimum temperature of 0 C or less benefit from winter tires. In addition I chose the period 1980 to 2013.

Result

It seems October is way too early to switch. Second half of November is an appropriate time.


Code

library(plyr)
library(dplyr)
library(ggplot2)
# data prepared without text heading
r1 <- read.csv('KNMI_20141115.edited.txt')
Sys.setlocale(category = "LC_TIME", locale = "C") # Not NL months
r2 <- mutate(r1,
    date = as.Date(format(YYYYMMDD),'%Y%m%d'),
    month =factor(months(date,abbreviate=TRUE),
        levels=months(as.Date(
                paste('2014',
                    formatC(1:12,digits=2,width=2,flag='0'),
                    '01',sep='-')),
            abbreviate=TRUE)),
    yearf=factor(format(date,'%Y')),
    yearn=as.numeric(substr(YYYYMMDD,1,4)),
    day=format(date,'%e'))
# days number 1 to 365, using numbers from 1901

days <- filter(r2,yearn==1901) %>%
    mutate(.,dayno=format(date,'%j') ) %>%
    select(.,month,day,dayno)

# select correct years and Months, remove leap day to sync day numbers
# and flag selected days
r3 <- merge(r2,days,all=TRUE) %>%
    filter(.,!grepl('0229',YYYYMMDD)) %>%
    mutate(.,daynon=as.numeric(dayno)) %>%
    filter(.,yearn>1980 & yearn<2014 & month %in% c('Oct','Nov','Dec')) %>%
    mutate(.,wt = TN<0 | TG<70 ) %>%
    mutate(.,month=factor(month,levels=c('Oct','Nov','Dec')))

# count the days
r4 <- as.data.frame(xtabs(wt ~ dayno,data=r3)/length(unique(r3$yearn)))
r4$daynon <- as.numeric(as.character(r4$dayno))

# plot

mylegend <- select(r3,day,month,daynon,dayno) %>%
    unique(.) %>%
    filter(.,day ==' 1')
g1 <- ggplot(r4,aes(x=daynon,y=Freq))
g1 + geom_smooth()  +
    geom_point() +
    scale_x_continuous('Date',breaks=mylegend$daynon,labels=mylegend$month) +
    ylab('Proportion wintery')  

Sunday, November 16, 2014

SAS PROC MCMC example in R; Poisson Regression

In this post I will try to copy the calculations of SAS's PROC MCMC example 61.5 (Poisson Regression) into the various R solutions. In this post Jags, RStan, MCMCpack, LaplacesDemon solutions are shown. Compared to the first post in this series, rcppbugs and mcmc are not used. Rcppbugs has no poisson distribution and while I know how to go around that (ones trick, second post in series) I just don't think that should be the approach for a standard distribution such as poisson. Mcmc has a weird summary which I dislike.

Data and model

Data are insurance claims. To quote SAS manual 'The variable n represents the number of insurance policy holders and the variable c represents the number of insurance claims. The variable car is the type of car involved (classified into three groups), and it is coded into two levels. The variable age is the age group of a policy holder (classified into two groups)'.
There is a nice trick in the model, again to quote; 'The logarithm of the variable n is used as an offset—that is, a regression variable with a constant coefficient of 1 for each observation. Having the offset constant in the model is equivalent to fitting an expanded data set with 3000 observations, each with response variable y observed on an individual level. The log link relates the mean and the factors car and age'.
There was an nice SAS trick in the data. The raw data has car as one variable and a design matrix is needed, for which PROC TRANSREG is used. In R using model.matrix() is at least a bit more transparant.
insure <- read.table(text='
n c car  age
500  42 small 0
1200 37 medium 0
100  1 large 0
400 101 small 1
500  73 medium 1
300  14 large 1',
skip=1,header=TRUE)

insure$ln=log(insure$n)
insure$car  <- relevel(insure$car,'small')
input_insure <-  model.matrix(~ car + age,data=insure)

Model

The models are all straightforward implementations of the SAS code.

MCMCpack

Probably the shortest code is for MCMCpack.
library(MCMCpack)
MCMCfun <- function(parm) {
    mu <- input_insure %*% parm
    LL <- sum(dnorm(parm,0,1e3,log=TRUE))
    yhat <- exp(mu+insure$ln)
    LP <- LL+sum(dpois(insure$c,yhat,log=TRUE))
    return(LP)
}
mcmcout <- MCMCmetrop1R(MCMCfun,rep(0,4))
summary(mcmcout)

LaplacesDemon

LaplacesDemon has a bit more code around essentially the same likelihood code as MCMCpack. In addition one needs to chose the sampling regime. I chose RWM with a first run providing values for parameter variances/covariances as a input for the final run. Of the four R based methods,  LaplacesDemon() is the function which creates most output.
library('LaplacesDemon')
mon.names <- "LP"
parm.names <- c('alpha','beta_car1','beta_car2','beta_age')
PGF <- function(Data) c(rnorm(4,0,1))

MyData <- list(mon.names=mon.names,
    parm.names=parm.names,
    PGF=PGF,
    y=insure$c,
    x=input_insure,
    ln=insure$ln)
Model <- function(parm, Data) {
    mu <- Data$x %*% parm
    LL <- sum(dnorm(parm,0,1e3,log=TRUE))
    yhat <- exp(mu+Data$ln)
    LP <- LL+sum(dpois(Data$y,yhat,log=TRUE))
    Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LP,
        yhat=yhat,
        parm=parm)
    return(Modelout)
}
Initial.Values <- GIV(Model, MyData, PGF=TRUE)
Fit1 <- LaplacesDemon(Model,
    Data=MyData,
    Algorithm='RWM',
    Covar=rep(.1,4),
    Initial.Values = Initial.Values
)
Fit2 <- LaplacesDemon(Model,
    Data=MyData,
    Algorithm='RWM',
    Covar=var(Fit1$Posterior2),
    Initial.Values = apply(Fit1$Posterior2,2,median)
)
Fit2

JAGS

The code is so simple, no comments are needed.

library(R2jags)
datain <- list(
    car1=input_insure[,2],
    car2=input_insure[,3],
    agev=insure$age,
    ln=insure$ln,
    c=insure$c,
    n=nrow(insure))
parameters <- c('a0','c1','c2','age')

jmodel <- function() {
    for (i in 1:n) {       
        mu[i] <- a0+car1[i]*c1+car2[i]*c2+agev[i]*age
        y[i] <- exp(mu[i]+ln[i])
        c[i] ~ dpois(y[i])
    }
    a0 ~ dnorm(0,.0001)
    c1 ~ dnorm(0,.0001)
    c2 ~ dnorm(0,.0001)
    age ~ dnorm(0,.0001)
}

jj <- jags(model.file=jmodel,
    data=datain,
    parameters=parameters)
jj

Stan

The code is again straightforward. As in previous posts, datain and parameters from Jags are re-used. Since I am still struggling a bit with Rstan I am using dplyr to combine model and simulation in one function. Stan is the slowest of all R based methods, the overhead of compiling is relative large for such a simple model 
library(dplyr)
library(rstan)
fit <- '
    data {
    int <lower=1> n;
    int c[n];
    vector[n] agev;
    vector[n] ln;
    vector[n] car1;
    vector[n] car2;
    }
    parameters {
    real a0;
    real c1;
    real c2;
    real age;
    }
    transformed parameters {
    vector[n] y;
    vector[n] mu;
    mu <- a0+car1*c1+car2*c2+agev*age;
    y <- exp(mu+ln);
    }
    model {
    c ~ poisson(y);
    a0 ~ normal(0,100);
    c1 ~ normal(0,100);
    c2 ~ normal(0,100);
    age ~ normal(0,100);
    }
    ' %>%
    stan(model_code = .,
        data = datain,
        pars=parameters)
fit

Results

MCMCpack

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
The Metropolis acceptance rate was 0.37195
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
> summary(mcmcout)

Iterations = 501:20500
Thinning interval = 1
Number of chains = 1
Sample size per chain = 20000

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

       Mean     SD  Naive SE Time-series SE
[1,] -2.643 0.1304 0.0009223       0.003448
[2,] -1.787 0.2736 0.0019350       0.007368
[3,] -0.696 0.1273 0.0008999       0.003315
[4,]  1.327 0.1335 0.0009440       0.003497

2. Quantiles for each variable:

        2.5%     25%     50%     75%   97.5%
var1 -2.9082 -2.7309 -2.6388 -2.5552 -2.3911
var2 -2.3606 -1.9639 -1.7772 -1.5944 -1.2809
var3 -0.9456 -0.7789 -0.6962 -0.6091 -0.4504
var4  1.0712  1.2363  1.3268  1.4188  1.5914

LaplacesDemon

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = apply(Fit1$Posterior2, 2, median), Covar = var(Fit1$Posterior2), Algorithm = "RWM")

Acceptance Rate: 0.3294
Algorithm: Random-Walk Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
     alpha  beta_car1  beta_car2   beta_age
0.02421653 0.09052011 0.01687834 0.02480825

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
        All Stationary
Dbar 62.614     62.614
pD    0.000      0.000
DIC  62.614     62.614
Initial Values:
     alpha  beta_car1  beta_car2   beta_age
-2.6482937 -1.7698730 -0.6848875  1.3175961

Iterations: 10000
Log(Marginal Likelihood): -33.58429
Minutes of run-time: 0.03
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 4
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 30
Specs: (NOT SHOWN HERE)
Status is displayed every 100 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 1000
Thinning: 10


Summary of All Samples
                 Mean           SD         MCSE      ESS          LB     Median
alpha      -2.6488552 1.342421e-01 6.565916e-03 582.9301  -2.9160940  -2.650964
beta_car1  -1.8137509 2.822791e-01 1.474688e-02 508.9384  -2.3852812  -1.815035
beta_car2  -0.6856638 1.283016e-01 6.234377e-03 701.7594  -0.9401333  -0.684194
beta_age    1.3241569 1.353041e-01 6.180607e-03 752.9705   1.0582629   1.324364
Deviance   62.6135632 1.379071e-06 6.585127e-08 643.3371  62.6135607  62.613563
LP        -49.8707670 1.445016e+00 6.295406e-02 648.5900 -53.5791639 -49.554267
                   UB
alpha      -2.4029950
beta_car1  -1.2906070
beta_car2  -0.4358381
beta_age    1.5855214
Deviance   62.6135661
LP        -48.0070303


Summary of Stationary Samples
                 Mean           SD         MCSE      ESS          LB     Median
alpha      -2.6488552 1.342421e-01 6.565916e-03 582.9301  -2.9160940  -2.650964
beta_car1  -1.8137509 2.822791e-01 1.474688e-02 508.9384  -2.3852812  -1.815035
beta_car2  -0.6856638 1.283016e-01 6.234377e-03 701.7594  -0.9401333  -0.684194
beta_age    1.3241569 1.353041e-01 6.180607e-03 752.9705   1.0582629   1.324364
Deviance   62.6135632 1.379071e-06 6.585127e-08 643.3371  62.6135607  62.613563
LP        -49.8707670 1.445016e+00 6.295406e-02 648.5900 -53.5791639 -49.554267
                   UB
alpha      -2.4029950
beta_car1  -1.2906070
beta_car2  -0.4358381
beta_age    1.5855214
Deviance   62.6135661
LP        -48.0070303

Jags

Inference for Bugs model at "/tmp/RtmpgHy1qo/modelc8f275c94f1.txt", fit using jags,
 3 chains, each with 2000 iterations (first 1000 discarded)
 n.sims = 3000 iterations saved
         mu.vect sd.vect   2.5%    25%    50%    75%  97.5%  Rhat n.eff
a0        -2.649   0.131 -2.908 -2.737 -2.647 -2.559 -2.394 1.001  2700
age        1.327   0.136  1.063  1.235  1.323  1.419  1.588 1.001  3000
c1        -1.795   0.280 -2.393 -1.976 -1.787 -1.595 -1.295 1.003   990
c2        -0.696   0.128 -0.948 -0.783 -0.695 -0.610 -0.436 1.002  1400
deviance  36.929   2.817 33.407 34.867 36.332 38.308 43.995 1.001  3000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = var(deviance)/2)
pD = 4.0 and DIC = 40.9
DIC is an estimate of expected predictive error (lower deviance is better).

Stan

Inference for Stan model: __LHS.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.

       mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
a0    -2.65    0.00 0.13  -2.93  -2.74  -2.64  -2.56  -2.40  1248    1
c1    -1.79    0.01 0.28  -2.34  -1.98  -1.78  -1.60  -1.27  1706    1
c2    -0.69    0.00 0.13  -0.94  -0.78  -0.69  -0.61  -0.44  1859    1
age    1.33    0.00 0.14   1.06   1.23   1.32   1.42   1.60  1411    1
lp__ 835.44    0.04 1.43 831.88 834.75 835.76 836.50 837.20  1294    1

Samples were drawn using NUTS(diag_e) at Sat Nov 15 14:15:39 2014.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).