# Opel Corsa Diesel Usage

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

I wanted to extend my car weight distribution calculation of June 16 from only 2000 to years 2000 to 2013. Unfortunately, come Sunday afternoon the code seemed too slow and not even the beginning of a post. So, I went on to another calculation I was working on; scraping Diesel usage data from the web and examining how it varied. The background of this is the repeated news items one reads, which state that actual fuel usage is much higher than reported. Since I may be in the market for an Opel Corsa Diesel, this was the car of choice.

### Sourcing the data

### Model

### Engines

### Months

### Drivers

### Factory results and driver results

### R code

library(R2jags)

data_list <- list(NDriver=nlevels(usage$Driver),

NBY=nlevels(usage$buildyear),

NEngine=nlevels(usage$engine),

NMonth=nlevels(usage$Month),

N=nrow(usage))

data_list$Driver <- (1:data_list$NDriver)[usage$Driver]

data_list$BY <- (1:data_list$NBY)[usage$buildyear]

data_list$Engine<- (1:data_list$NEngine)[usage$engine]

data_list$Month <- (1:data_list$NMonth)[usage$Month]

data_list$lp100 <- usage$lp100

uu <- unique(data.frame(engine=usage$engine,cDriver=usage$Driver,Driver=data_list$Driver,Engine=data_list$Engine))

data_list$uDriver <- uu$Driver

data_list$uEngine <- uu$Engine

data_list$Driver1 <- uu$Driver[uu$Engine==1]

data_list$Driver2 <- uu$Driver[uu$Engine==2]

data_list$Driver3 <- uu$Driver[uu$Engine==3]

data_list$NDriver1 <- length(data_list$Driver1)

data_list$NDriver2 <- length(data_list$Driver2)

data_list$NDriver3 <- length(data_list$Driver3)

mymod <- function() {

for (i in 1:N) {

fit[i] <- grandmean +

MMonth[Month[i]] +

MEngine[Engine[i]] +

MBY[BY[i]] +

MDriver[Driver[i]]

lp100[i] ~dnorm(fit[i],tauOfDriver[Driver[i]])

}

grandmean ~ dnorm(0,.001)

MMonth[1] <- 0

for (i in 2:NMonth) {

MMonth[i] ~ dnorm(0,tauMonth)

}

tauMonth ~ dgamma(0.001,0.001)

sdMonth <- sqrt(1/tauMonth)

MBY[1] <- 0

for (i in 2:NBY) {

MBY[i] ~ dnorm(0,tauBY)

}

tauBY ~ dgamma(0.001,0.001)

sdBY <- sqrt(1/tauBY)

MEngine[1] <- 0

for (i in 2:NEngine) {

MEngine[i] ~ dnorm(0,tauEngine)

}

tauEngine ~ dgamma(0.001,0.001)

sdEngine <- sqrt(1/tauEngine)

MDriver[1] <- 0

for (i in 2:NDriver) {

MDriver[i] ~ dnorm(0,tauDriver)

}

tauDriver ~ dgamma(0.001,0.001)

sdDriver <- sqrt(1/tauDriver)

for (i in 1:NDriver) {

tauOfDriver[i] <- pow(sdOfDriver[i],-2)

sdOfDriver[i] ~ dlnorm(mu.lsdD,tau.lsdD)

}

mu.lsdD ~ dnorm(0,.0001)

tau.lsdD <- pow(sigma.lsdD,-2)

sigma.lsdD ~ dunif(0,100)

for (i in 1:NDriver) {

meanDriver[i] <- grandmean + MDriver[uDriver[i]] + mean(MMonth[1:NMonth]) + mean(MBY[1:NBY]) + MEngine[uEngine[i]]

}

for (i in 1:NMonth) {

meanMonth[i] <- grandmean + mean(MDriver[1:NDriver]) + MMonth[i] + mean(MBY[1:NBY]) + mean(MEngine[1:NEngine])

}

for (i in 1:NBY) {

meanBY[i] <- grandmean + mean(MDriver[1:NDriver]) + mean(MMonth[1:NMonth]) +MBY[i]+ mean(MEngine[1:NEngine])

}

for ( i in 1:NDriver1) {

MDriver1[i] <- MDriver[Driver1[i]]

}

for ( i in 1:NDriver2) {

MDriver2[i] <- MDriver[Driver2[i]]

}

for ( i in 1:NDriver3) {

MDriver3[i] <- MDriver[Driver3[i]]

}

meanEngine[1] <- grandmean + mean(MDriver1) + mean(MMonth[1:NMonth])+ mean(MBY[1:NBY]) +MEngine[1]

meanEngine[2] <- grandmean + mean(MDriver2) + mean(MMonth[1:NMonth])+ mean(MBY[1:NBY]) +MEngine[2]

meanEngine[3] <- grandmean + mean(MDriver3) + mean(MMonth[1:NMonth])+ mean(MBY[1:NBY]) +MEngine[3]

}

model.file <- file.path(tempdir(),'mijnmodel.txt')

write.model(mymod,con=model.file)

inits <- function() list(

grandmean = rnorm(1,5,1),

MDriver = c(NA,rnorm(data_list$NDriver-1)) ,

MBY = c(NA,rnorm(data_list$NBY-1)) ,

MEngine = c(NA,rnorm(data_list$NEngine-1)) ,

MMonth = c(NA,rnorm(data_list$NMonth-1)) ,

tauDriver = runif(1,1,3),

tauBY = runif(1,1,3),

tauEngine=runif(1,1,3),

tauMonth=runif(1,1,3),

sdOfDriver = runif(data_list$NDriver,1,3)

)

parameters <- c('sdOfDriver','meanMonth','meanDriver','meanBY',

‘meanEngine’,’sdDriver’,’sdMonth’,’sdEngine’)

jagsfit <- jags(data=data_list,inits=inits,model.file=model.file,

parameters.to.save=parameters,n.chains=4,DIC=FALSE,n.iter=20000,

progress.bar=’none’)

png(‘jagsfit.png’)

plot(jagsfit)

dev.off()

jagsfit.mc <- as.mcmc(jagsfit)

fitsummary <- summary(jagsfit.mc)

stats <- fitsummary$statistics

quants <- fitsummary$quantiles

library(ggplot2)

Engines <- as.data.frame(cbind(stats[grepl('meanEngine',rownames(stats)),],

quants[grepl(‘meanEngine’,rownames(quants)),]))

Engines$type <- levels(usage$engine)

Engines

limits <- aes(ymax = `25%`, ymin=`75%`)

p <- ggplot(Engines, aes(y=Mean, x=type))

dodge <- position_dodge(width=0.9)

png(‘engines.png’)

p + geom_point(position=dodge) +

geom_errorbar(limits, position=dodge, width=0.25) +

coord_flip()

dev.off()

Months <- as.data.frame(cbind(stats[grepl('meanMonth',rownames(stats)),],

quants[grepl(‘meanMonth’,rownames(quants)),]))

Months$Month <- levels(usage$Month)

limits <- aes(ymax = `25%`, ymin=`75%`)

p <- ggplot(Months, aes(y=Mean, x=Month))

dodge <- position_dodge(width=0.9)

png(‘Months.png’)

p + geom_point(position=dodge) +

geom_errorbar(limits, position=dodge, width=0.25) +

coord_flip()

dev.off()

Driver <- as.data.frame(cbind(stats[grepl('meanDriver',rownames(stats)),],

quants[grepl(‘meanDriver’,rownames(quants)),]))

Driver <- cbind(uu,Driver)

limits <- aes(ymax = `25%`, ymin=`75%`)

p <- ggplot(Driver, aes(y=Mean, x=Driver,colour=engine))

dodge <- position_dodge(width=0.9)

png(‘driver.png’)

p + geom_point(position=dodge) +

geom_errorbar(limits, position=dodge, width=0.25) +

coord_flip()

dev.off()

ag <- aggregate(usage$lp100,by=list(Driver=usage$Driver),mean)

min(ag$x)

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