stanLogTransformed.R

temp <- c(11.9,14.2,15.2,16.4,17.2,18.1,18.5,19.4,22.1,22.6,23.4,25.1)
units <- c(185L,215L,332L,325L,408L,421L,406L,412L,522L,445L,544L,614L)
library(rstan)
stanmodel <- stan_model(model_code = stanLogTransformed)
fit <- sampling(stanmodel,
                data = list(N=length(units),
                            units=units,
                            temp=temp),
                iter = 1000, warmup=200)
stanoutput <- extract(fit)

## Extract generated posterior predictive quantities
Sims <- data.frame(stanoutput[["units_pred"]])

## Calculate summary statistics
SummarySims <- apply(Sims, 2, summary)
colnames(SummarySims) <- paste(icecream$temp,"ºC")

## Extract estimated parameters
(parms <- sapply(stanoutput[c("alpha", "beta", "sigma")], mean))
##      alpha       beta      sigma 
## 4.41543439 0.08186099 0.14909969

## Use parameters to predict median and mean
PredMedian <- exp(parms['alpha'] + parms['beta']*temp)
PredMean <- exp(parms['alpha'] + parms['beta']*temp + 0.5*parms['sigma']^2)

## Compare predictions based on parameters with simulation statistics
round(rbind(SummarySims, PredMedian, PredMean),1)
##            11.9 ºC 14.2 ºC 15.2 ºC 16.4 ºC 17.2 ºC 18.1 ºC
## Min.         116.1   101.6   123.6   140.6   164.0   209.5
## 1st Qu.      196.2   237.6   259.1   285.7   307.5   329.7
## Median       218.5   263.7   286.1   318.3   338.8   364.1
## Mean         221.9   267.2   290.2   321.0   342.6   369.3
## 3rd Qu.      243.5   292.2   317.5   350.4   372.7   403.4
## Max.         418.7   542.4   740.9   688.0   736.0   778.5
## PredMedian   218.5   264.1   286.8   316.5   338.1   364.1
## PredMean     220.9   267.0   289.9   320.0   341.8   368.1
##            18.5 ºC 19.4 ºC 22.1 ºC 22.6 ºC 23.4 ºC 25.1 ºC
## Min.         198.9   177.2   241.9   271.6   249.1   337.0
## 1st Qu.      339.6   368.7   456.5   476.9   508.7   577.9
## Median       375.2   405.2   506.2   527.4   562.9   645.2
## Mean         379.3   410.4   512.2   534.3   572.1   653.7
## 3rd Qu.      412.2   446.0   557.5   582.6   625.4   718.2
## Max.         678.1   984.9   944.6  1054.0  1282.0  1358.0
## PredMedian   376.3   405.3   506.2   527.4   563.4   648.0
## PredMean     380.4   409.6   511.6   533.1   569.5   655.0