Should I use premium Diesel? Result: No

[This article was first published on 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.

A while ago I had a post: ‘Should I use premium Diesel? Setup. Since that time the data has been acquired. This post describes the results.

Data

Data is registered by me in 2014 and 2015. 2014 has standard Diesel, while 2015 has premium. Both are from the same months of the year, so should have approximately the same weather and traffic.

Results

The figure shows the main result: While the 2015 has some better results, it is not by much. It is not in the ballpark of 5% which was originally defined.

Results from distribution fitting and t.test confirm the visual observation.

fitdistr

Standard
      mean          sd    
  3.59517971   0.19314598 
 (0.04828649) (0.03414371)
Premium 
      mean          sd    
  3.54015573   0.23905047 
 (0.05976262) (0.04225855)

t.test

Welch Two Sample t-test

data:  usage by kind
t = -0.6934, df = 28.732, p-value = 0.4936
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.2173822  0.1073343
sample estimates:
 mean in group Premium mean in group Standard 
              3.540156               3.595180 

Posterior density for premium Diesel

Previously I defined a prior for the analysis of the data. This was a mixed distribution from three parts. Like Standard fuel, 5% better and something else in the neighborhood of 3.6. Given the simplicity of the model and the non-standard prior I used MCMCpack’s MCMCmetrop1R function. This allowed me to write it all as R code.
The result is a density for mean usage quite close to 3.6.

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,] 3.5409 0.06413 0.0004535       0.001722
[2,] 0.2699 0.05353 0.0003785       0.001448

2. Quantiles for each variable:

       2.5%   25%    50%    75%  97.5%
var1 3.4147 3.494 3.5465 3.5884 3.6523
var2 0.1873 0.232 0.2629 0.2998 0.3953

Conclusion

While premium Diesel had a few occasions with remarkable fuel usage, it did not result in a relevant decrease in fuel usage overall. In addition, while the advised prices for standard and premium are 5% apart, the achievable price difference is 10%. The real price fighters do not do premium, hence the pressure to discount is less. This 10% price increase is more than I initially expected. Hence the data do not show benefit for using premium Diesel.

Data and code

library(MASS)
library(ggplot2)
library(MCMCpack)

r1 <- read.table(text='l km
28.25 710.6
22.93 690.4
28.51 760.5
23.22 697.9
31.52 871.2
24.68 689.6
30.85 826.9
23.04 699
29.96 845.3
30.16 894.7
25.71 696
23.6 669.8
28.57 739
27.23 727.4
18.31 499.9
24.28 689.5′,header=TRUE)
r1$usage=100*r1$l/r1$km
r1$Fuel=’Standard’

r2 <- read.table(text='l km
23.97 696
25.33 693.5
24.19 698.5
23.27 689.6
25.78 707.8
15.84 510.1
32.23 860.2
23.29 674.6
25.33 678.5
26.44 686
36.43 913.2
26.46 850.1
24.38 678.8
34.85 1001.3
28.34 858.5
25.85 698.8′
,header=TRUE)
r2$usage=100*r2$l/r2$km
r2$Fuel=’Premium’

rr <- rbind(r1,r2)
rr$Fuel <- factor(rr$Fuel)
ggplot(rr,aes(Fuel,usage))+geom_violin()+scale_y_continuous(‘usage (l/100 km)’)

fitdistr(r1$usage,’normal’)
fitdistr(r2$usage,’normal’)

t.test(usage ~ kind, data=rr)

priormeandens <- function(usage) {
  dens <- (dnorm(usage,3.6,.05)+
      dnorm(usage,3.6/1.05,.05)+
      dnorm(usage,(3.6+3.6/1.05)/2,.15))/3
   log(dens)
 }
 mydens <- function(theta,x) {
   musage=theta[1]
   sdusage=theta[2]
   datdens <- sum(dnorm(x,musage,sdusage,log=TRUE))
   priormean <- priormeandens(musage)
   priorsd <- dunif(sdusage,0,1,log=TRUE)
   datdens+priormean+priorsd
 } 
mcmc1 <- MCMCmetrop1R(mydens,
    theta.init=c(3.8,.2),
    seed=as.integer((Sys.time())),
    x=r2$usage)
plot(mcmc1)
summary(mcmc1)

To 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.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)