Bayes models from SAS PROC MCMC in R, post 2

This is my second post in converting SAS's PROC MCMC examples in R. The task in his week is determining the transformation parameter in a Box-Cox transformation. SAS only determines Lambda, but I am not so sure about that. What I used to do was get an interval for Lambda, select an interpretable value (e.g. 1/2 is square root, -1 is inverse) and progress using this parameter. However, if that is my approach in a Bayesian context, maybe I should make my prior on lambda only those values I find acceptable. On the other hand, if I model lambda continuous, should not the uncertainty in Lambda also effect the uncertainty in the remainder of the model? In this post I will consider them all continuous. Hence the first step is to get SAS to give me those values, for which I added those parameters in the code and ran it through SAS Studio. Subsequently find the code for the various R approaches. For readability I stuck the code next to text, while output is at the bottom of the post. 

Note on previous post

Last week I failed to get some samplers. I now realize that some of these might have been solved with the ones trick or the zero trick. In addition, LaplacesDemon has dnormm() which is more elegant than my approach.

SAS

This is just the model from the SAS manual with extra parameters monitored.

proc mcmc data=boxcox nmc=50000 propcov=quanew seed=12567 monitor=(lda beta0 beta1 sd);
    ods select PostSumInt TADpanel;
    parms beta0 0 beta1 0 lda 1 s2 1;
    beginnodata;
    prior beta: ~ general(0);
    prior s2 ~ gamma(shape=3, scale=2);
    prior lda ~ unif(-2, 2);
    sd=sqrt(s2);
    endnodata;
    ys=(y**lda-1)/lda;
    mu=beta0+beta1*x;
    ll=(lda-1)*log(y)+lpdfnorm(ys, mu, sd);
    model general(ll);
run;

STAN

Rstan has the Jacobian needed for this exercise. It was fairly straightforward to fit it all together.
library(plyr)
library(dplyr)
# stan
library(rstan)
datain <- list(n=100,x=set2$x,y=set2$y)
parameters <- c('lda','beta0','beta1','std')
fit <- '
    data {
    int <lower=1> n;
    vector[n] y;
    vector[n] x;
    }
    parameters {
    real <lower=-2,upper=2> lda;
    real beta0;
    real beta1;
    real <lower=0> s2;
    }
    transformed parameters {
    vector[n] ys;
    vector[n] ll;
    vector[n] mu;
    real <lower=0> std;
    std <- sqrt(s2);
    for (i in 1:n) {
    ys[i] <- (pow(y[i],lda)-1)/lda;
    }
    mu <- beta0+beta1*x;
    ll <- (lda-1)*log(y);
    }
    model {
    s2 ~ gamma(3, 0.5);
    increment_log_prob(ll);
    increment_log_prob(normal_log(ys,mu,std));
    }
    ' %>%
stan(model_code = .,
    data = datain,
    pars=parameters)
fit

JAGS 

The JAGS code is not straightforward for two reasons.
The first reason is that Jags does not have the concept of Jacobian. However, upon some investigation I found the zeros trick. Which basically means that an arbitrary contribution to the deviance might be added by transforming that deviance and pretending to JAGS the values are from a Poisson distribution.
The second reason is that JAGS has directed acyclic graphs. At least, I think that is the reason. One cannot first calculate transformed values and then make them dependent on a distribution. My solution was to subtract transformed and estimated values and make the result normal around 0.
To get acceptable results in terms of n.eff and Rhat I increased the number of samples to 10000 per chain.
library(R2jags)
datain <- list(n=100,
    x=set2$x,
    y=set2$y,
    zero1 = rep(0,100),
    zero2 = rep(0,100))
parameters <- c('LDA','beta0','beta1','std')

jmodel <- function() {
    tau <- 1/s2
    std <- sqrt(s2)
    for (i in 1:n) {       
        mu[i] <- beta0+x[i]*beta1
        ys[i] <- (pow(y[i],LDA)-1)/LDA -mu[i]
        zero1[i] ~ dnorm(ys[i],tau)
        ll[i] <- -((LDA-1)*log(y[i]))+100
        zero2[i] ~ dpois(ll[i])
    }
    beta0 ~ dnorm(0,.0001)
    beta1 ~ dnorm(0,.0001)
    s2 ~ dgamma(3,.5)
    LDA ~ dunif(-2,2)
}
inits  <- function() {
    list(beta0=rnorm(1),
        beta1=rnorm(1),
        LDA=runif(1,-1,1),
        s2=runif(1,.1,3))
}
jj <- jags(model.file=jmodel,
    data=datain,
    inits=inits,
    n.iter=10000,
    parameters=parameters)
jj

MCMC

As such it is relatively easy to build the model for MCMC. Based on the demo vignette I was also able to set scale to obtain an acceptance rate about 0.2. After some additional calculations I came to batch length of 10000, and chose for 20 batches to keep computation time somewhat limited. All in all it took quite some attempts to get to these settings in terms of nbatch and batch length. MCMC only gives summary statistics based on each separate batch. These are subsequently summarized. The results seem a bit different from SAS and JAGS. In addition a small tweak was needed to get SD as output rather than internal variable S2.
mydens <- function(x,...) {
    if (x[4]<0 ) return(-Inf)
    std <- sqrt(x[4])
    mu <- x[2]+set2$x*x[3]
    ys <- (set2$y^x[1]-1)/x[1]
    ll <- (x[1]-1)*log(set2$y)
    sum(dnorm(ys,mu,std,log=TRUE)+ll) +dgamma(x[4],shape=3,scale=2)
}
sam <-  metrop(mydens,
    initial=c(.46,-2,1.8,1.3),
    nbatch=20,
    blen=10000,
    scale=.025,
    outfun=function(z, ...) {
        z[4] <- sqrt(z[4])
        c(z, z^2)}
)
sam$accept
foo <- apply(sam$batch, 2, mean)
mu <- foo[1:4]
sigmasq <- foo[5:8] - mu^2
data.frame(param=c('LDA','beta0','beta1','SD'),mu,sigmasq)

MCMCpack

Obviously, the deviance code is same as MCMC. Based on the ACF a thinning of 20 was chosen. I did not see an obvious reason to change from the default 20000 samples.
library(MCMCpack)
library(coda)
mydens <- function(x,...) {
    if (x[4]<0 ) return(-Inf)
    std <- sqrt(x[4])
    mu <- x[2]+set2$x*x[3]
    ys <- (set2$y^x[1]-1)/x[1]
    ll <- (x[1]-1)*log(set2$y)
    sum(dnorm(ys,mu,std,log=TRUE)+ll) +dgamma(x[4],shape=3,scale=2)
}
sam <- MCMCmetrop1R(mydens,
    theta.init=c(.46,-2,1.8,1.3),
    tune=1.5,thin=20)
summary(sam)

rcppbugs

Rcppbugs has many of the same limits as JAGS in terms of available distributions. The same tricks as used in JAGS are hence used in rcppbugs, with one exception, where JAGS has a Poison distribution, rcppbugs lacks it, hence the ones trick is used, which is the equivalent to zeros trick, but using Bernoulli distribution. The ACF gave quite a big lag, hence I upped both the number of samples and increased thinning.

library(rcppbugs)
x <- set2$x
y <- set2$y
one <- rep(1,100)
yzero <- rep(0,100)
beta <- mcmc.normal(c(.1,.1),mu=0,tau=0.0001)
LDA <- mcmc.uniform(.5,-2,2)
tau.y <- mcmc.gamma(sd(as.vector(y)),alpha=0.1,beta=0.1)
y.hat <- deterministic(function(x,beta,y,LDA) {
        ypred=beta[1]+beta[2]*x
        ytrans=(y^LDA-1)/LDA
        ypred-ytrans
    }, x, beta,y,LDA)
Jac <- deterministic(function(LDA,y)
        exp((LDA-1)*log(y))/100,LDA,y)
cJac <- mcmc.bernoulli(one,Jac,observed=TRUE)
y.lik <- mcmc.normal(yzero,mu=y.hat,tau=tau.y,observed=TRUE)

m <- create.model(beta,LDA, tau.y,
    y.hat,
    y.lik,
    Jac,
    cJac)
ans <- run.model(m, iterations=1e6L, burn=1e4L, adapt=1e3L, thin=1000L)

mysum <- function(x) c(mean=mean(x),sd=sd(x),quantile(x,c(.025,.5,.975)))
rbind(LDA=mysum(ans$LDA),
    beta0=mysum(ans$b[,1]),
    beta1=mysum(ans$b[,2]),
    sd=mysum(1/sqrt(ans$tau.y)))

LaplacesDemon

LaplacesDemon also needed a lot of samples. At a later attempt I tried the HARM sampler which worked much better. Getting the right sampler is important.
library('LaplacesDemon')
mon.names <- "LP"
parm.names <- c('beta0','beta1','LDA','std')
PGF <- function(Data) c(rnorm(2,0,1),runif(1,.4,.5),runif(1,.5,2))

MyData <- list(mon.names=mon.names,
    parm.names=parm.names,
    PGF=PGF,
    y=set2$y,
    x=set2$x)
Model <- function(parm, Data)
{
    ### Parameters
    std <- interval(parm[4], 1e-100, Inf)
    parm[4] <- std
    s2 <- parm[4]^2
    mu <- parm[1]+Data$x*parm[2]
    ys <- (Data$y^parm[3]-1)/parm[3]
    ll <- (parm[3]-1)*log(Data$y)
    LL <- sum(dnorm(ys,mu,std,log=TRUE)+ll)
    LP <- LL+dgamma(s2,shape=3,scale=2)
    Modelout <- list(LP=LP, Dev=-2*LL, Monitor=LP,
        yhat=1,
        parm=parm)
    return(Modelout)
}
# very dependent on initial values
Initial.Values <- GIV(Model, MyData, PGF=TRUE)
Fit <- LaplacesDemon(Model,
    Data=MyData,
    Algorithm='HARM',
    Specs=list(alpha.star=0.234, B=NULL),
    Initial.Values = Initial.Values,
    Iterations=40000,Status=2000,Thinning=30
)
Fit

Output

SAS

                              Standard 
Parameter     N     Mean      Deviation     95% HPD Interval
lda        50000    0.4701     0.0290     0.4142     0.5281
beta0      50000   -2.2064     0.4215    -3.0564    -1.4201
beta1      50000    1.9260     0.1806     1.5854     2.2826
sd         50000    1.3275     0.1482     1.0648     1.6328

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
lda      0.47    0.00 0.03    0.41    0.45    0.47    0.49    0.53   595 1.01
beta0   -2.20    0.02 0.40   -3.01   -2.45   -2.18   -1.94   -1.44   684 1.01
beta1    1.92    0.01 0.17    1.59    1.81    1.91    2.02    2.28   571 1.01
std      1.32    0.01 0.15    1.06    1.22    1.31    1.42    1.65   691 1.00
lp__  -309.43    0.05 1.45 -313.02 -310.13 -309.09 -308.36 -307.64   977 1.00

Samples were drawn using NUTS(diag_e) at Sat Oct  4 15:14:25 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).

JAGS 

Inference for Bugs model at "/tmp/RtmpSgjWNR/modelef330519173.txt", fit using jags,
 3 chains, each with 10000 iterations (first 5000 discarded), n.thin = 5
 n.sims = 3000 iterations saved
           mu.vect sd.vect      2.5%       25%       50%       75%     97.5%
LDA          0.472   0.030     0.413     0.451     0.472     0.491     0.535
beta0       -2.229   0.443    -3.200    -2.512    -2.199    -1.920    -1.450
beta1        1.938   0.192     1.608     1.805     1.925     2.054     2.356
std          1.336   0.155     1.079     1.227     1.319     1.433     1.680
deviance 20620.595   3.365 20616.435 20618.142 20619.763 20622.197 20629.012
          Rhat n.eff
LDA      1.009  2900
beta0    1.003  1400
beta1    1.008  1500
std      1.002  3000
deviance 1.011   480

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 = 5.6 and DIC = 20626.2
DIC is an estimate of expected predictive error (lower deviance is better).

MCMC

  param         mu      sigmasq
1   LDA  0.4614172 0.0007537545
2 beta0 -2.0808097 0.1414263712
3 beta1  1.8693639 0.0267559710
4    SD  1.2832247 0.0180115043

MCMCpack

Note that MCMCpack has s2 rather than sd in output

Iterations = 501:20481
Thinning interval = 20
Number of chains = 1
Sample size per chain = 1000

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

        Mean      SD  Naive SE Time-series SE
[1,]  0.4657 0.02927 0.0009256      0.0009256
[2,] -2.1565 0.40274 0.0127358      0.0113576
[3,]  1.8988 0.17764 0.0056175      0.0052342
[4,]  1.7048 0.36930 0.0116782      0.0133882

2. Quantiles for each variable:

        2.5%     25%    50%     75%   97.5%
var1  0.4058  0.4479  0.466  0.4854  0.5215
var2 -2.9944 -2.4064 -2.146 -1.8799 -1.4346
var3  1.5606  1.7810  1.890  2.0099  2.2647
var4  1.1204  1.4473  1.652  1.9261  2.5331

rcppbugs

            mean         sd       2.5%        50%      97.5%
LDA    0.4597165 0.02875052  0.4050122  0.4592264  0.5207108
beta0 -2.0937414 0.39740763 -2.9511520 -2.0668905 -1.4015904
beta1  1.8647679 0.17411536  1.5640795  1.8566258  2.2509009
sd     1.2612126 0.14011464  1.0228554  1.2502009  1.5634482

LaplacesDemon

Call:
LaplacesDemon(Model = Model, Data = MyData, Initial.Values = Initial.Values,
    Iterations = 40000, Status = 2000, Thinning = 30, Algorithm = "HARM",
    Specs = list(alpha.star = 0.234, B = NULL))

Acceptance Rate: 0.22825
Algorithm: Hit-And-Run Metropolis
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
       beta0        beta1          LDA          std
0.1329756466 0.0283948376 0.0007961449 0.0174295808

Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
         All Stationary
Dbar 619.952    619.768
pD     9.531      3.530
DIC  629.483    623.298
Initial Values:
[1] -0.2820944  0.3190142  0.4778030  1.3380326

Iterations: 40000
Log(Marginal Likelihood): -313.4286
Minutes of run-time: 0.12
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: 532
Recommended Burn-In of Un-thinned Samples: 15960
Recommended Thinning: 31
Specs: (NOT SHOWN HERE)
Status is displayed every 2000 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 1333
Thinning: 30


Summary of All Samples
                 Mean        SD        MCSE       ESS           LB       Median
beta0      -2.0923727 0.3614098 0.091343057  17.08625   -2.8552455   -2.0861434
beta1       1.8711745 0.1631220 0.037952180  22.13164    1.5969508    1.8620763
LDA         0.4612769 0.0282230 0.004279672  34.00465    0.4086439    0.4616171
std         1.2653917 0.1320557 0.020805981  98.82246    1.0352842    1.2574597
Deviance  619.9523160 4.3660058 0.415674235 188.93047  616.3618665  619.2464350
LP       -309.9042345 2.1778219 0.206890192 190.00073 -313.2903260 -309.5529196
                   UB
beta0      -1.4732055
beta1       2.2047884
LDA         0.5171681
std         1.5495557
Deviance  626.7398118
LP       -308.1117717


Summary of Stationary Samples
                 Mean         SD        MCSE        ESS           LB
beta0      -2.0498690 0.36586611 0.118090993   8.836406   -2.9348999
beta1       1.8509120 0.16921671 0.050214255  15.798638    1.5936363
LDA         0.4577846 0.02886563 0.006775378  21.788913    0.4078682
std         1.2449204 0.12994746 0.027000671  55.838883    1.0332563
Deviance  619.7681778 2.65715568 0.326080727 121.032423  616.3311632
LP       -309.8148748 1.32309998 0.206890192 123.122008 -313.0799005
               Median           UB
beta0      -1.9902222   -1.5178152
beta1       1.8236876    2.2326418
LDA         0.4558953    0.5206561
std         1.2320277    1.5344502
Deviance  619.3065650  626.3929182
LP       -309.5806926 -308.1020531

Appendix: data used

s1 <- scan(what=list(double(),double()),text='
10.0  3.0  72.6  8.3  59.7  8.1  20.1  4.8  90.1  9.8   1.1  0.9
78.2  8.5  87.4  9.0   9.5  3.4   0.1  1.4   0.1  1.1  42.5  5.1
57.0  7.5   9.9  1.9   0.5  1.0 121.1  9.9  37.5  5.9  49.5  6.7
8.3  1.8   0.6  1.8  53.0  6.7 112.8 10.0  40.7  6.4   5.1  2.4
73.3  9.5 122.4  9.9  87.2  9.4 121.2  9.9  23.1  4.3   7.1  3.5
12.4  3.3   5.6  2.7 113.0  9.6 110.5 10.0   3.1  1.5  52.4  7.9
80.4  8.1   0.6  1.6 115.1  9.1  15.9  3.1  56.5  7.3  85.4  9.8
32.5  5.8  43.0  6.2   0.1  0.8  21.8  5.2  15.2  3.5   5.2  3.0
0.2  0.8  73.5  8.2   4.9  3.2   0.2  0.3  69.0  9.2   3.6  3.5
0.2  0.9 101.3  9.9  10.0  3.7  16.9  3.0  11.2  5.0   0.2  0.4
80.8  9.4  24.9  5.7 113.5  9.7   6.2  2.1  12.5  3.2   4.8  1.8
80.1  8.3  26.4  4.8  13.4  3.8  99.8  9.7  44.1  6.2  15.3  3.8
2.2  1.5  10.3  2.7  13.8  4.7  38.6  4.5  79.1  9.8  33.6  5.8
9.1  4.5  89.3  9.1   5.5  2.6  20.0  4.8   2.9  2.9  82.9  8.4
7.0  3.5  14.5  2.9  16.0  3.7  29.3  6.1  48.9  6.3   1.6  1.9
34.7  6.2  33.5  6.5  26.0  5.6  12.7  3.1   0.1  0.3  15.4  4.2
2.6  1.8  58.6  7.9  81.2  8.1  37.2  6.9
')
set2 <- data.frame(x=s1[[2]],y=s1[[1]])

Bayesian models in R

There are many ways to run general Bayesian calculations in or from R. The best known are JAGS, OpenBUGS and STAN. Then some time ago Rasmus Bååth had a post Three ways to run Bayesian models in R in which he mentioned LaplacesDemon (not on CRAN) on top of those. A check of the Bayes task view gives 'MCMCpack (...) contains a generic Metropolis sampler that can be used to fit arbitrary models', 'The mcmc package consists of an R function for a random-walk Metropolis algorithm for a continuous random vector' and 'Note that rcppbugs is a package that attempts to provide a pure R alternative to using OpenBUGS/WinBUGS/JAGS for MCMC'. 
Clearly there is a wealth of approaches and each will have its strengths and weaknesses. To get an idea on their comparison I decided to run a number of calculations through all of them. The source of these calculations will be the fine SAS manual, where PROC MCMC has 21 examples. The first example, which I will try this week, is drawing samples from a number of distributions, from which I selected the mixture of three normal distributions.

Task

Draw samples from a mixture of three normal distributions. Let me start by stating that I would not use anything as complex as an MCMC sampler for this task, my personal SAS solution would be:
data p1;
    do i=1 to 1000;
        s=rand('table', .3, .4);
        select (s);
            when (1) r=rand('normal', -3, 2);
            when (2) r=rand('normal', 2, 1);
            when (3) r=rand('normal', 10, 4);
        end;
        output;
    end;
    keep r;
run;

proc sgplot data=p1;
    density r / type=kernel (c=.5);
run;

The equivalent R code:
library(plyr)
library(dplyr)
nu <- c(-3,2,10)
sd <- c(2,1,4)
tau <- sd^-2
wgt <- c(.3,.4,.3)

sampls <- sample(1:3,1000,replace=TRUE,wgt)
sampls <- rnorm(1000,nu[sampls],sd[sampls])
density(sampls) %>%    plot

MCMCpack

Just having a density was enough for MCMCpack. The samples have quite a long autocorrelation but do follow the distribution(not shown).
library(MCMCpack)
mydens <- function(x,...) {
    step1 <- sum(wgt*dnorm(x,nu,sd))
    log(step1)
}
mcmcpack <- MCMCmetrop1R(mydens,theta.init=1)
acf(mcmcpack)

mcmc

mcmc has the property that its result is not the samples, but rather summary statistics of the samples. Hence it does not give samples from the desired distribution.

JAGS

The JAGS samples had no autocorrelation and followed the distribution.
library(R2jags)
jmodel <- function() {
    i ~ dcat(wgt)
    nunu <- (i==1)*nu[1]+
            (i==2)*nu[2]+
            (i==3)*nu[3]
        tautau <- (i==1)*tau[1]+
            (i==2)*tau[2]+
            (i==3)*tau[3]
       
    p ~ dnorm(nunu,tautau)
}

datain <- list(nu=nu,tau=tau,wgt=wgt)
parameters <- c('i','p')

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

STAN

I did not manage to get STAN to give me the samples.

rccpbugs

I did not manage to get rccpbugs to give me the samples.

LaplacesDemon

LaplacesDemon feels like the most complex of direct R solutions.  
library('LaplacesDemon')
mon.names <- "LP"
parm.names <- 'x'
MyData <- list(mon.names=mon.names, parm.names=parm.names)
N <-1
Model <- function(parm, Data)
{
    x <- parm[1]
    LL <- sum(wgt*dnorm(x,nu,sd)) %>% log
    Modelout <- list(LP=LL, Dev=-2*LL, Monitor=LL,
        yhat=x, parm=parm)
    return(Modelout)
}
Initial.Values <- 1
Fit <- LaplacesDemon(Model, Data=MyData, Initial.Values)

Trying dplyr on triathon data

There was a triathlon in Almere last week, like every year since 1983. I pulled the data of all years to get some idea how things have changed in that sport. To get a visual I decided to plot the best 10% athletes. Then later I decided this was an ideal moment to look at plyr and dplyr again, so rewrote everything using those tools. I must say I like them very much and intend to use them again.

Data

Data from triathlon-uitslagen.nl, this years website www.challenge-almere.com and wikipedia which told me there was one lady in 1983. The data is in a bunch of spreadsheets (.xls and .xlsx), with different columns used, detail level and even those merged cells which make things look nice. Some years had besides finishers also those who did not, broken down by time of certain milestones. Some years had teams next to individuals. It was a bit too much variation to be coded in R, I did my first cleaning, removing spurious columns and rows, in libreoffice and made relatively uniform .xls files.
During importing R coded most times as date times some as character. The date times had a date (1899-12-31), so I made time to difftimes and retained the original as character. The actual year of the match is pulled from the spreadsheet name.
library(XLConnect)
library(plyr)
library(dplyr)
library(lattice)
setwd('C:\\Users\\Kees\\Documents\\r hobby\\almeretriatlon')
dd <- dir(,'.xls')

readyear <- function(filename) {
  print(filename)
  wb <- loadWorkbook(filename)
  rw <- readWorksheet(wb,sheet=1)
  names(rw) <- tolower(names(rw))
  names(rw)[names(rw)=='total'] <- 'time'
  mutate(rw,
      timec=if (class(rw$time)[1]=='character') {
            time
          } else format(time,'%H:%M:%S') ,
      time=as.difftime(timec),
      year=as.numeric(gsub('[A-Za-z\\. ]','',filename))
  )
}
years0 <- lapply(dd,readyear)
The resulting data may have a gender variable, and/or, a code variable. The code variable may be character and contain gender plus age bracket as digits (or team, or PRO or ELI), but it can also only contain the age bracket and hence be numerical. Gender can be D, H, M, V, W, interpreted as dame, heer, man, vrouw, woman respectively (female, male, male, female, female).
For plotting purposes a two digit Year factor was added, running from 83 to 14.
triatlon <- ldply(years0,function(x) {
# first unify cat info & gender 
          x <- if ('cat' %in% names(x)) {
                if (class(x$cat)[1]=='character') {
                  mutate(x,
                      gender=substr(cat,1,1),
                      cat=substring(cat,2))
                } else  mutate(x,
                      cat=as.character(cat))
              } else mutate(x,cat=NA)
        }) %>% 
    filter(.,!(gender %in% c('T','P'))) %>%
    mutate(.,
        timen=as.numeric(time) , 
        gender=factor(
            ifelse(gender %in% c('H','M'),
                'Male','Female')) ,
        Year=factor(substr(format(year),3,4),
            levels=substr(format(
                    seq(min(year),max(year))),3,4))) %>%
    select(.,
        gender,
        cat,
        year,
        Year,
        time,
        timec,
        timen)  

Plot 

This is actually where plyr/dplyr worked best to get nice clean code, it was as if it was made for this purpose, define groups, define selection and plot. The choice for lattice plots was made before I went to recode and I did not feel like changing it.

The plot shows times steadily decreasing through the previous century and stabilizing around turn of the century. 2006 was an European championship and 2008 world championship, and these years do dot stick out as being faster. Hence these times should be fairly representative of what can be achieved in the Netherlands. There is quite some year to year variation, but that can be attributed to weather and tracks used.  

group_by(triatlon,year,gender) %>% 

    filter(.,percent_rank(time)<0.1) %>%

    bwplot(timen ~ Year | gender,

        data=.,

        ylab='Time (h)',

        xlab='Year',

        horizon=FALSE,

        layout=c(1,2,1),

        drop.unused.levels=FALSE,

        scales=list(cex=.7))

Trying a prefmap

Preference mapping is a key technique in sensory and consumer research. It links the sensory perception on products to the liking of products and hence provides clues to the development of new, well tasting, products. Even though it is a key technique, it is also a long standing problem how to perform such an analysis. In R the SensoMineR package provides a prefmap procedure. This post attempts to create such an analysis with Stan.

Data

Data are the coctail data from the SensoMineR package. After conversion to a scale 0 to 10 with 10 most liked, the product means are:
   means
1   5.03
2   3.02
3   5.42
4   3.55
5   5.67
6   5.74
7   3.84
8   3.75
9   4.17
10  4.26
11  3.20
12  3.88
13  5.98
14  3.95
15  6.47
16  4.90

Model

The model is build upon my post of last week: Mapping products in a space . What is added is a consumer section. Each consumer's preference is modeled as a ideal point, where liking is maximum, with points further away liked less and less. In mathematical terms the distance dependent function is max_liking * e^{-distance*scale}. The ideal point is captured by three numbers; its liking and its coordinates. The scale function is, for now, common for all consumers. In addition there is a lot of code for administration of all parameters.
model1 <- "
data {
        int<lower=0> ndescriptor;
        int<lower=0> nproduct;
        matrix[nproduct,ndescriptor] profile;
    int<lower=0> nconsumer;
    matrix[nproduct,nconsumer] liking;
}
parameters {
    row_vector[nproduct] prodsp1;
    row_vector[nproduct] prodsp2;
    real<lower=0,upper=1> sigma1;
    real<lower=0,upper=1> sigma2;
    matrix [nconsumer,2] optim;
    vector <lower=0> [nconsumer] maxima;
    real <lower=0> scale;
    real <lower=0> sliking;
}
transformed parameters {
   vector [ndescriptor] descrsp1;
   vector [ndescriptor] descrsp2;
     matrix[nproduct,ndescriptor] expected1;  
     matrix[nproduct,ndescriptor] expected2;  
     matrix[nproduct,ndescriptor] residual1;  
     vector[nproduct] distances;
     matrix[nproduct,nconsumer] likepred;


   descrsp1 <- profile'*prodsp1';
   expected1 <- (descrsp1*prodsp1)';
   residual1 <- profile-expected1;
   descrsp2 <- profile'*prodsp2';
   expected2 <- (descrsp2*prodsp2)';
   for (i in 1:nconsumer) {
      for (r in 1:nproduct) {
      distances[r] <- sqrt(square(prodsp1[r]-optim[i,1])+
                           square(prodsp2[r]-optim[i,2]));
      likepred[r,i] <- maxima[i]*exp(-distances[r]*scale);
      }
   }
}
model {  
     for (r in 1:nproduct) {
        prodsp1[r] ~ normal(0,1);
        prodsp2[r] ~ normal(0,1);
        for (c in 1:ndescriptor) {
           profile[r,c] ~ normal(expected1[r,c],sigma1);
           residual1[r,c] ~ normal(expected2[r,c],sigma2);
        }
        for (i in 1:nconsumer) {
           liking[r,i] ~ normal(likepred[r,i],sliking);
           optim[i,1] ~ normal(0,2);
           optim[i,2] ~ normal(0,2);
        }
    scale ~ normal(1,.1);
    maxima ~ normal(5,2);
    sliking ~ normal(2,1);
    }
}
generated quantities {
   vector [ndescriptor] descrspace1;
   vector [ndescriptor] descrspace2;
   row_vector [nproduct] prodspace1;
   row_vector [nproduct] prodspace2;
   matrix [nconsumer,2] optima;

   prodspace1 <-(
                     ((prodsp1[12]>0)*prodsp1)-
                     ((prodsp1[12]<0)*prodsp1)
                  );
   prodspace2 <-(
                     ((prodsp2[12]>0)*prodsp2)-
                     ((prodsp2[12]<0)*prodsp2)
                  ); 
   descrspace1 <-(
                     ((prodsp1[12]>0)*descrsp1)-
                     ((prodsp1[12]<0)*descrsp1)
                  );
   descrspace2 <-(
                     ((prodsp2[12]>0)*descrsp2)-
                     ((prodsp2[12]<0)*descrsp2)
                  ); 
   for (i in 1:nconsumer) {
      optima[i,1] <- (
                        ((prodsp1[12]>0)*optim[i,1])-
                        ((prodsp1[12]<0)*optim[i,1])
                     );
      optima[i,2] <- (
                        ((prodsp2[12]>0)*optim[i,2])-
                        ((prodsp2[12]<0)*optim[i,2])
                     );
   }
}
"

Analysis results

Sensominer's result

For comparative reasons the plot resulting from SensoMineR's carto() function. I have followed the parameter settings from the SensoMineR package to get this plot. Color is liking, numbered dots are products. The blue zone is best liked, as can be seen from the products with highest means residing there.

New method

In the plot the blue dots are samples of ideal points, the bigger black numbers are locations of products and the smaller red numbers are consumer's ideal points.
This is different from the SensoMineR map , the consumers have pulled well liked products such as 13 and 15 to the center. In a way, I suspect that in this analysis the consumer's preference has overruled most information from the sensory space. Given that, I will be splitting consumers.

Three groups of consumers

Three groups of consumers were created via k-means clustering. From sensory and consumer insight point of view the clusters may describe three different ways to experience the particular products. Obviously a clustering upon demographics or marketing segments may be equally valid, but I don't have that information. The cluster sizes are 15, 52 and 33 respectively.

Cluster 1

This cluster is characterized by liking for products 8 to 11. Compared to the original space, this cluster does not like products 13 and 15 so much, does not dislike product 4 and 12 so much.

 

Cluster 2

These are the bulk of the consumers and the result of all consumers is more pronounced. However, product 1 has shifter quite a distance to liked.

Cluster 3

This plot is again fairly similar to the all consumer plot. What is noticeable here is that there is a void in the center. The center of the most liked region is not occupied.

Next Steps

There are still some things to improve in this approach. Better tuning of the various priors in the model. Modeling the range of consumer's liking rather than solely their maximum. It may be possible to have the scale parameter subject dependent. Perhaps a better way to extract the dimensions from sensory space, thereby avoiding the Jacobian warning and using estimated standard deviations of the sensory profiling data. Finally, improved graphics.

Code

# Reading and first map

# senso.cocktail
# hedo.cocktail
library(SensoMineR)
data(cocktail)
res.pca <- PCA(senso.cocktail,graph=FALSE)
# SensoMineR does a dev.new for each graph, hence captured like this.
dev.new <- function() png('carto.png')
res.carto <- carto(res.pca$ind$coord[,1:2],
    graph.tree=FALSE,
    graph.corr=FALSE,
    hedo.cocktail)
dev.off()
# reset default graph settings
rm(dev.new)
dev.new()

# model

library(rstan)
nprod <- 16
ndescr <- 13
nconsumer <- 100
sprofile <- as.matrix(scale(senso.cocktail))
datain <- list(
    nproduct=nprod,
    ndescriptor=ndescr,
    profile=sprofile,
    nconsumer=nconsumer,
    liking = as.matrix(10-hedo.cocktail[,1:nconsumer])
)
data.frame(means=rowMeans(10-hedo.cocktail)  )

model1 <- "
data {
        int<lower=0> ndescriptor;
        int<lower=0> nproduct;
        matrix[nproduct,ndescriptor] profile;
    int<lower=0> nconsumer;
    matrix[nproduct,nconsumer] liking;
}
parameters {
    row_vector[nproduct] prodsp1;
    row_vector[nproduct] prodsp2;
    real<lower=0,upper=1> sigma1;
    real<lower=0,upper=1> sigma2;
    matrix [nconsumer,2] optim;
    vector <lower=0> [nconsumer] maxima;
    real <lower=0> scale;
    real <lower=0> sliking;
}
transformed parameters {
   vector [ndescriptor] descrsp1;
   vector [ndescriptor] descrsp2;
     matrix[nproduct,ndescriptor] expected1;  
     matrix[nproduct,ndescriptor] expected2;  
     matrix[nproduct,ndescriptor] residual1;  
     vector[nproduct] distances;
     matrix[nproduct,nconsumer] likepred;


   descrsp1 <- profile'*prodsp1';
   expected1 <- (descrsp1*prodsp1)';
   residual1 <- profile-expected1;
   descrsp2 <- profile'*prodsp2';
   expected2 <- (descrsp2*prodsp2)';
   for (i in 1:nconsumer) {
      for (r in 1:nproduct) {
      distances[r] <- sqrt(square(prodsp1[r]-optim[i,1])+
                           square(prodsp2[r]-optim[i,2]));
      likepred[r,i] <- maxima[i]*exp(-distances[r]*scale);
      }
   }
}
model {  
     for (r in 1:nproduct) {
        prodsp1[r] ~ normal(0,1);
        prodsp2[r] ~ normal(0,1);
        for (c in 1:ndescriptor) {
           profile[r,c] ~ normal(expected1[r,c],sigma1);
           residual1[r,c] ~ normal(expected2[r,c],sigma2);
        }
        for (i in 1:nconsumer) {
           liking[r,i] ~ normal(likepred[r,i],sliking);
           optim[i,1] ~ normal(0,2);
           optim[i,2] ~ normal(0,2);
        }
    scale ~ normal(1,.1);
    maxima ~ normal(5,2);
    sliking ~ normal(2,1);
    }
}
generated quantities {
   vector [ndescriptor] descrspace1;
   vector [ndescriptor] descrspace2;
   row_vector [nproduct] prodspace1;
   row_vector [nproduct] prodspace2;
   matrix [nconsumer,2] optima;

   prodspace1 <-(
                     ((prodsp1[12]>0)*prodsp1)-
                     ((prodsp1[12]<0)*prodsp1)
                  );
   prodspace2 <-(
                     ((prodsp2[12]>0)*prodsp2)-
                     ((prodsp2[12]<0)*prodsp2)
                  ); 
   descrspace1 <-(
                     ((prodsp1[12]>0)*descrsp1)-
                     ((prodsp1[12]<0)*descrsp1)
                  );
   descrspace2 <-(
                     ((prodsp2[12]>0)*descrsp2)-
                     ((prodsp2[12]<0)*descrsp2)
                  ); 
   for (i in 1:nconsumer) {
      optima[i,1] <- (
                        ((prodsp1[12]>0)*optim[i,1])-
                        ((prodsp1[12]<0)*optim[i,1])
                     );
      optima[i,2] <- (
                        ((prodsp2[12]>0)*optim[i,2])-
                        ((prodsp2[12]<0)*optim[i,2])
                     );
   }
}
"

pars <- c('prodspace1','prodspace2','optima','scale','maxima')

fit <- stan(model_code = model1,
    data = datain,
    pars=pars)

# plotting

fitsamps <- as.data.frame(fit)

combiplot <- function(fitsamps,datain,labs) {
    prod <- reshape(fitsamps,
        drop=names(fitsamps)[33:ncol(fitsamps)],
        direction='long',
        varying=list(names(fitsamps)[1:16],
            names(fitsamps)[17:32]),
        timevar='sample',
        times=1:16,
        v.names=c('PDim1','PDim2')
    )
        sa <- sapply(1:16,function(x)
            c(sample=x,
                Dim1=mean(prod$PDim1[prod$sample==x]),
                Dim2=mean(prod$PDim2[prod$sample==x])))
    sa <- as.data.frame(t(sa))
   
    optimindex <- grep('optima',names(fitsamps))
    noptim <- datain$nconsumer
    loc <- reshape(fitsamps,
        drop=names(fitsamps)[(1:ncol(fitsamps))[-optimindex]],
        direction='long',
        varying=list(names(fitsamps)[optimindex[1:noptim]],
            names(fitsamps)[optimindex[(1:noptim)+noptim]]),
        timevar='subject',
        times=1:noptim,
        v.names=c('Dim1','Dim2')
    )
    locx <- loc[sample(nrow(loc),60000),]
    plot(locx$Dim1,locx$Dim2,
        col='blue',
        pch=46,
        cex=2,
        xlim=c(-1,1)*.7,
        ylim=c(-1,1)*.7)
    sa2 <- sapply(1:noptim,function(x)
            c(sample=x,
                Dim1=mean(loc$Dim1[loc$subject==x]),
                Dim2=mean(loc$Dim2[loc$subject==x])))
    sa2 <- as.data.frame(t(sa2))
    text(x=sa2$Dim1,y=sa2$Dim2,labels=labs,cex=.8,col='red')
    text(x=sa$Dim1,y=sa$Dim2,labels=sa$sample,cex=1.5)
    invisible(fitsamps)
}

combiplot(fitsamps,datain,1:100)

# three clusters

tlik <- t(scale(hedo.cocktail))
km <- kmeans(tlik,centers=3)
table(km$cluster)

datain1 <- list(
    nproduct=nprod,
    ndescriptor=ndescr,
    profile=sprofile,
    nconsumer=sum(km$cluster==1),
    liking = as.matrix(10-hedo.cocktail[,km$cluster==1])
)
fit1 <- stan(model_code = model1,
    data = datain1,
    fit=fit,
    pars=pars)

fitsamps1 <- as.data.frame(fit1)
#

datain2 <- list(
    nproduct=nprod,
    ndescriptor=ndescr,
    profile=sprofile,
    nconsumer=sum(km$cluster==2),
    liking = as.matrix(10-hedo.cocktail[,km$cluster==2])
)
fit2 <- stan(model_code = model1,
    data = datain2,
    fit=fit,
    pars=pars)

fitsamps2 <- as.data.frame(fit2)
##
datain3 <- list(
    nproduct=nprod,
    ndescriptor=ndescr,
    profile=sprofile,
    nconsumer=sum(km$cluster==3),
    liking = as.matrix(10-hedo.cocktail[,km$cluster==3])
)
fit3 <- stan(model_code = model1,
    data = datain3,
    fit=fit,
    pars=pars)

fitsamps3 <- as.data.frame(fit3)
combiplot(fitsamps1,datain1,which(km$cluster==1))
combiplot(fitsamps2,datain2,which(km$cluster==2))
combiplot(fitsamps3,datain3,which(km$cluster==3))

Mapping products in a space

I have read about people doing a Bayesian PCA at some points and always wondered how that would work. Then, at some point I thought of a way to do so. As ideas evolved my interest became not PCA as such, but rather in a prefmap. As a first step in that this post contains the mapping from a sensory space to a two dimensional space. For prefmap this step is commonly done via a PCA.

Data

Data are the coctail data from sensominer package.

Algorithm

The calculation is mostly inspired by the many PLS algorithms to which I was exposed when I was doing chemometrics. Scores and loadings may be obtained from each other by multiplying with the data matrix. In this case it means I just take a set of product scores and obtain the associated descriptors via a simple matrix multiplication. The resulting product and descriptor vectors can be used to reconstruct the original matrix; the best solution minimizes difference between the constructed and original data. For dimension two subtract reconstructed data from original data and repeat on residuals.

Scaling

PCA has the possibility to have unit length scores or loadings, or R and Q mode if that is your favorite jargon. If one has a more singular value decomposition look, it is just where the eigenvalues go. At this point I made the choice to do that in the variable space. 

Unique solution

PCA is known not to have one unique solution; each solution is equivalent to its mirror image. It seemed most elegant to do this completely at the end, after inspection of the data it seemed the location of product 12 was suitable for making the solution unique, since it was extreme on both dimensions. The final step (generated quantities) forces the location to be top right quadrant for data reported.

Code

library(rstan)
nprod <- 16
ndescr <- 13
sprofile <- as.matrix(scale(senso.cocktail,scale=FALSE))
datain <- list(
    nproduct=nprod,
    ndescriptor=ndescr,
    profile=sprofile
    )
   
model1 <- "
data {
        int<lower=0> ndescriptor;
        int<lower=0> nproduct;
        matrix[nproduct,ndescriptor] profile;
}
parameters {
    row_vector[nproduct] prodsp1;
    row_vector[nproduct] prodsp2;
    real<lower=0,upper=1> sigma1;
    real<lower=0,upper=1> sigma2;
}
transformed parameters {
   vector [ndescriptor] descrsp1;
   vector [ndescriptor] descrsp2;
   matrix[nproduct,ndescriptor] expected1;  
   matrix[nproduct,ndescriptor] expected2;  
   matrix[nproduct,ndescriptor] residual1;  

   descrsp1 <- profile'*prodsp1';
   expected1 <- (descrsp1*prodsp1)';
   residual1 <- profile-expected1;
   descrsp2 <- profile'*prodsp2';
   expected2 <- (descrsp2*prodsp2)';
}
model {  
     for (r in 1:nproduct) {
        prodsp1[r] ~ normal(0,1);
        prodsp2[r] ~ normal(0,1);
        for (c in 1:ndescriptor) {
           profile[r,c] ~ normal(expected1[r,c],sigma1);
           residual1[r,c] ~ normal(expected2[r,c],sigma2);
        }
     }
}
generated quantities {
   vector [ndescriptor] descrspace1;
   vector [ndescriptor] descrspace2;
   row_vector [nproduct] prodspace1;
   row_vector [nproduct] prodspace2;
   prodspace1 <-(
                     ((prodsp1[12]>0)*prodsp1)-
                     ((prodsp1[12]<0)*prodsp1)
                  );
   prodspace2 <-(
                     ((prodsp2[12]>0)*prodsp2)-
                     ((prodsp2[12]<0)*prodsp2)
                  ); 
   descrspace1 <-(
                     ((prodsp1[12]>0)*descrsp1)-
                     ((prodsp1[12]<0)*descrsp1)
                  );
   descrspace2 <-(
                     ((prodsp2[12]>0)*descrsp2)-
                     ((prodsp2[12]<0)*descrsp2)
                  ); 
}
"
pars <- c('prodspace1','prodspace2','descrspace1','descrspace2')
fit1 <- stan(model_code = model1,
    data = datain,
    pars=pars)

Results

For comparison, first a standard biplot.

Product space

It is not difficult to extract the samples and plot them. See end of post. One notable property of the plot is that the products are in ellipses with the minor axis towards the center. Apparently part of variation between MCMC samples is rotational freedom between dimensions. Other than that the solution is actually pretty close to the PCA

Descriptor space

The rotational freedom is even more clear here.

Additional code

data

library(SensoMineR)
data(cocktail)

biplot

pr <- prcomp(senso.cocktail) 
plot(pr)
biplot(pr)

product plot

fit1samps <- as.data.frame(fit1)

prod <- reshape(fit1samps,
    drop=names(fit1samps)[33:59],
    direction='long',
    varying=list(names(fit1samps)[1:16],
        names(fit1samps)[17:32]),
    timevar='sample',
    times=1:16,
    v.names=c('PDim1','PDim2')
)
   
prod <- prod[order(prod$PDim1),]
plot(prod$PDim1,prod$PDim2,
    col=c(2,17,3,4,6,5,7:10,13,12,11,14:16)[prod$sample],
    pch=46,
    cex=2,
    xlim=c(-1,1)*.75,
    ylim=c(-1,1)*.75)
sa <- sapply(1:16,function(x)
        c(sample=x,
            Dim1=mean(prod$PDim1[prod$sample==x]),
            Dim2=mean(prod$PDim2[prod$sample==x])))
sa <- as.data.frame(t(sa))
text(x=sa$Dim1,y=sa$Dim2,labels=sa$sample,cex=1.5)

descriptor plot

descr <- reshape(fit1samps,
    drop=names(fit1samps)[c(1:32,59)],
    direction='long',
    varying=list(names(fit1samps)[33:45],
        names(fit1samps)[46:58]),
    timevar='sample',
    times=1:13,
    v.names=c('DDim1','DDim2')
)

descr <- descr[order(descr$DDim1),]
plot(descr$DDim1,descr$DDim2,
    col=c(2,1,3:13)[descr$sample],
    pch=46,
    cex=2,
    xlim=c(-1,1)*9,
    ylim=c(-1,1)*9)
sa <- sapply(1:13,function(x)
        c(sample=x,
            Dim1=mean(descr$DDim1[descr$sample==x]),
            Dim2=mean(descr$DDim2[descr$sample==x])))
sa <- as.data.frame(t(sa))
text(x=sa$Dim1,y=sa$Dim2,labels=names(senso.cocktail))