mages

library(lattice)
key <- list(
  rep=FALSE, 
  lines=list(col=c("#00526D", "blue"), type=c("p","l"), pch=1),
  text=list(lab=c("Observation","Estimate")),
  rectangles = list(col=adjustcolor("yellow", alpha.f=0.5), border="grey"),
  text=list(lab="95% Prediction credible interval"))
xyplot(l.95..CI + u.95..CI + Estimate + Units_sold ~ Temperature | Model, 
       data=modelData, as.table=TRUE, main="Ice cream model comparision",
       xlab="Temperatures (C)", ylab="Units sold", 

mages

modelData <- data.frame(
  Model=factor(c(rep("Linear model", n), 
                 rep("Log-transformed LM", n),
                 rep("Poisson (log)",n),
                 rep("Binomial (logit)",n)),  
               levels=c("Linear model", 
                        "Log-transformed LM",
                        "Poisson (log)",
                        "Binomial (logit)"), 
               ordered = TRUE),

mages

library(brms)
# Linear Gaussian model
lin.mod <- brm(units ~ temp, family="gaussian")
# Log-transformed Linear Gaussian model
log.lin.mod <- brm(log_units ~ temp, family="gaussian")
# Poisson model
pois.mod <- brm(units ~ temp, family="poisson")
# Binomial model
bin.mod <- brm(units | trials(market.size) ~ temp, family="binomial")

mages

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)
log_units <- log(units)
n <- length(units)
market.size <- rep(800, n)

mages

glmModelPlot <- function(x, y, xlim,ylim, meanPred,  LwPred, UpPred, 
                         plotData, main=NULL){
                           ## Based on code by Arthur Charpentier:
                           ## http://freakonometrics.hypotheses.org/9593
                           par(mfrow=c(1,1))
                           n <- 2
                           N <- length(meanPred)
                           zMax <- max(unlist(sapply(plotData, "[[", "z")))*1.5
                           mat <- persp(xlim, ylim, matrix(0, n, n), main=main,
                                        zlim=c(0, zMax), theta=-30, 

mages

meanPred <- apply(Sims, 2, mean)
LwPred <- apply(Sims, 2, quantile, 0.05)
UpPred <- apply(Sims, 2, quantile, 0.95)
xlim <- c(min(icecream$temp)*0.95, max(temp)*1.05)
ylim <- c(floor(min(units)*0.95),
          ceiling(max(units)*1.05))
plotStan <- lapply(
  seq(along=temp),
  function(i){
    stp = 251

mages

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)

mages

stanLogTransformed <-"
data {
  int N;
  vector[N] units;
  vector[N] temp;
}
transformed data {  
  vector[N] log_units;        
  log_units = log(units);
}

mages

# Markus Gesmann
library(arm) # for 'display' function only
icecream <- data.frame(
  # http://www.statcrunch.com/5.0/viewreport.php?reportid=34965&groupid=1848
  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)
)
basicPlot <- function(...){

mages

library(latex2exp)
x <- seq(-4, 4, len = 101)
y <- cbind(sin(x), cos(x))

op=par(mfrow=c(2,1))
## plotmath
matplot(x, y, type = "l", xaxt = "n",
        main = expression(paste("plotmath: ", plain(sin) * phi, "  and  ",
                                plain(cos) * phi)),
        ylab = expression("sin" * phi, "cos" * phi), # only 1st is taken