#Set up

•rm(list=ls(all=TRUE))

•library(xlsReadWrite)

•library(ggplot2)

•theme_set(theme_bw())

•library(quantreg)

•library(gplots)

•setwd("I:/Project/2017 LD GHG Rule/Footprint_analysis/Outputs")

•fleetorig <- xlsReadWrite::read.xls("I:/Project/2017 LD GHG 

Rule/Footprint_analysis/fleet/2010_20120509_fleet.xls")

#Drop Electrics and Diesels, keep hybrids like DOT did in NPRM.

fleet<-fleetorig[fleetorig$Engine.Cycle %in% c("H","G","R"),]

##Calculate Weights - Don't remove vehicles after this, as these are based on total 

volumes

•a<-subset(fleet,NHTSA.Defined.New.NHTSA.Car.Truck == "C")

•car.vol<-sum(a$Volume.2021)

•b<-subset(fleet,NHTSA.Defined.New.NHTSA.Car.Truck == "T")

•truck.vol<-sum(b$Volume.2021)

•fleet$Weight<-ifelse(fleet$NHTSA.Defined.New.NHTSA.Car.Truck=="C",fleet$Volume.2021/car.v

ol,fleet$Volume.2021/truck.vol)

•

annot<-function(a)

•{paste(signif(coef(a)[[1]],digits=3),"+",signif(coef(a)[[2]],digits=3),"x : 

R2=",signif(summary(a)[[8]],digits=3),sep = "")}

•

annotrq<-function(a)

•{paste(signif(coef(a)[[1]],digits=3),"+",signif(coef(a)[[2]],digits=3),"x",sep = "")}

•

analysis<-function(x) {

•#Create Fleet

•

•fleet1<-subset(fleet,NHTSA.Defined.New.NHTSA.Car.Truck == x)

••

•#Diagnostics

•#HP/WT ~ Footprint

•lm.1<-lm(fleet1$Horsepower/fleet1$Curb.Weight ~ fleet1$Footprint, weights=fleet1$Weight)

•#HP/WT ~ Footprint Squared

•lm.2<-lm(fleet1$Horsepower/fleet1$Curb.Weight ~ poly(fleet1$Footprint,2,raw=TRUE), 

weights=fleet1$Weight)

•#WT/FP ~ FP

•lm.12<-lm(I(fleet1$Curb.Weight/fleet1$Footprint) ~ fleet1$Footprint, weights = 

fleet1$Weight)

•#HP ~ WT

•lm.3<-lm(fleet1$Horsepower~poly(fleet1$Curb.Weight,1, raw=TRUE), weights=fleet1$Weight)

•#HP ~WT Squared

•lm.4<-lm(fleet1$Horsepower~poly(fleet1$Curb.Weight,2, raw=TRUE), weights=fleet1$Weight)

•#CO2 ~ Footprint

•lm.5<-lm(fleet1$CO2.Package.1 ~ fleet1$Footprint, weights=fleet1$Weight)

•#CO2 ~ HP/WT + WT

•lm.8<-lm(fleet1$CO2.Package.1 ~ 

I(fleet1$Horsepower/fleet1$Curb.Weight)+I(fleet1$Curb.Weight),weights = fleet1$Weight)

•#Create CO2 w/ Adjustments

•#HP/WT Adjusted CO2

•fleet1$hp.wt.avg<-sum(fleet1$Horsepower/fleet1$Curb.Weight*fleet1$Weight)

•fleet1$hp.wt.CO2<-(fleet1$CO2.Package.1-(fleet1$Horsepower/fleet1$Curb.Weight-fleet1$hp.w

t.avg)*coef(lm.8)[2])

•#WT/FP adjusted CO2

•fleet1$crb.wt.fp.avg<-sum(fleet1$Curb.Weight/fleet1$Footprint*fleet1$Weight)

•fleet1$crbwt.norm.CO2<-(fleet1$CO2.Package.1-coef(lm.8)[3]*(fleet1$Curb.Weight-fleet1$crb

.wt.fp.avg*fleet1$Footprint))

•#HP/WT and WT/FP Adjusted CO2
•fleet1$hpwt.crbwt.norm.CO2<-(fleet1$hp.wt.CO2-coef(lm.8)[3]*(fleet1$Curb.Weight-fleet1$cr

b.wt.fp.avg*fleet1$Footprint))

•#Return Answers

•ans<-list(fleet1,lm.1,lm.2,lm.12,lm.3,lm.4,lm.5,lm.8)

•invisible(return(ans))

}

#Function for the regressions needed on each of the CO2 datasets.

regressions<-function(dataset){

•datasetb<-data.frame(dataset[1])

•#List the datasets here

•CO2set<-data.frame(datasetb$CO2, datasetb$No.Tech.CO2, datasetb$CO2.Package.1, 

datasetb$hp.wt.CO2, datasetb$crbwt.norm.CO2, datasetb$hpwt.crbwt.norm.CO2)

••runr<-function(dsb){

•••weighted.ols•<-lm(dsb~datasetb$Footprint,weights = datasetb$Weight)

•••unweighted.ols•<-lm(dsb~datasetb$Footprint,weights = NULL)

•••weighted.mad•<-rq(dsb~datasetb$Footprint, tau=.5, weights = datasetb$Weight)

•••unweighted.mad•<-rq(dsb~datasetb$Footprint, tau=.5,weights = NULL)

••invisible(return(regr<-list(weighted.ols,unweighted.ols,weighted.mad,unweighted.mad)))

••}

•invisible(return(j<-lapply(CO2set,runr)))

•}

#Create the Datasets, and excute the regressions, then merge the car and truck fleets

cars.egr<-analysis("C")

car.regressions<-regressions(cars.egr)

trucks.egr<-analysis("T")

truck.regressions<-regressions(trucks.egr)

fleet4<-data.frame(rbind(cars.egr[[1]], trucks.egr[[1]]))

•

car_base_equation<-annot(cars.egr[[8]])

truck_base_equation<-annot(trucks.egr[[8]])

••

summarize_regressions<-function(name,number){••

•name<-data.frame(rbind(

••cbind(name,"weighted 

OLS",annot(car.regressions[[number]][[1]]),annot(truck.regressions[[number]][[1]])),

••cbind(name,"unweighted 

OLS",annot(car.regressions[[number]][[2]]),annot(truck.regressions[[number]][[2]])),

••cbind(name,"weighted 

MAD",annotrq(car.regressions[[number]][[3]]),annotrq(truck.regressions[[number]][[3]])),

••cbind(name,"unweighted 

MAD",annotrq(car.regressions[[number]][[4]]),annotrq(truck.regressions[[number]][[4]]))

••))

•

•return(name)

•}

•

all_regressions<-rbind(•

summarize_regressions("2010",1),

summarize_regressions("No.Tech",2),

summarize_regressions("Tech.Normalized",3),

summarize_regressions("HP/WT-Tech-Normalized",4),

summarize_regressions("WT/FP-Tech-Normalized",5),

summarize_regressions("HP/WT&WT/FP-Tech-Normalized",6)

)

colnames(all_regressions)<-c("Dataset","Regression","Car Equation", "Truck Equation")

all_regressions

dchart<-function(tit,jj,mult){

•fleet4$RCO2<-fleet4[,jj]

•base.chart<-ggplot(data=fleet4, aes(y=RCO2, x=Footprint, size=Volume.2021, 
weight=Weight, group=NHTSA.Defined.New.NHTSA.Car.Truck))+

••geom_point(aes(shape=1))+

••labs(x="Footprint", y="CO2", size="2021 Projected Sales")+

••opts(title=paste(tit,"CO2 v. Footprint"))+

••facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+

••scale_x_continuous(limits=c(35,80))+

••#weighted lines

••geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE) +

••geom_quantile(quantiles=.5, colour="red")+

••#unweighted lines

••geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE, 

weight=I(1),linetype=5) +

••geom_quantile(quantiles=.5, colour="red",weight=1,linetype=5)+

••#no legend

••opts(legend.position = "none")

•return(base.chart)

}

•

•physics<-function(){

#Plot of WT/ v. FP - Weighted OLS

p2<-ggplot(data=fleet, aes(y=I(fleet$Curb.Weight), x=fleet$Footprint, size=Volume.2021, 

weight=fleet$Weight, shape=1))+

•geom_point(aes(shape=1))+

•labs(x="Footprint", y="WT", colour="C/T", size="2021 Sales")+

•geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE)+

•geom_smooth(span=1,colour="red")+

•facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+

•opts(title="WT v. FP - Weighted OLS and Loess Fit")+

•opts(legend.position = "none")

#Plot of WT/FP v. FP - Weighted OLS

p1<-ggplot(data=fleet, aes(y=I(fleet$Curb.Weight/fleet$Footprint), x=fleet$Footprint, 

size=Volume.2021, weight=fleet$Weight))+

•geom_point(aes(shape=1))+

•labs(x="Footprint", y="WT/FP", colour="C/T", size="2021 Sales")+

•geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE)+

•geom_smooth(span=1,colour="red")+

•facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+

•opts(title="WT/FP v. FP - Weighted OLS and Loess Fit")+

•opts(legend.position = "none")

•

# Plot of HP/WT v. FP - Weighted OLS

p4<-ggplot(data=fleet, aes(y=I(fleet$Horsepower/fleet$Curb.Weight), x=fleet$Footprint, 

size=Volume.2021, weight=fleet$Weight))+

•geom_point(aes(shape=1))+

•geom_smooth(method=lm, formula= y~ poly(x,1, raw=TRUE), se=FALSE)+

•facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck)+

•labs(x="Footprint", y="HP/WT", colour="C/T", size="2021 Sales")+

•opts(title="HP/WT v. FP - Weighted OLS")+

•opts(legend.position = "none")

return(list(p1,p2,p4))

}

#Footprint Histogram•

fphist<-ggplot(data=fleet,aes(x=Footprint, weight=Volume.2021))+

••geom_bar(color="black", fill="white", binwidth=1)+

••facet_wrap(~NHTSA.Defined.New.NHTSA.Car.Truck,ncol=1)+

••labs(x="Footprint",y="MY 2021 Projected Sales")+

••opts(title=paste("Footprint Distribution by Car and Truck"))+

••geom_vline(xintercept=41, linetype=2)+

••geom_vline(xintercept=56, linetype=2)+

••geom_vline(xintercept=74, linetype=6)+

••geom_vline(xintercept=41)
••

#Plot by vehicle class

vcplot<-function(){

man.sales<-aggregate((fleet$Volume.2021), by=list(fleet$Our.Class), sum)

colnames(man.sales)<- c("Our.Class", "Man.2021.Sales")

man.footprint<-stack(by(fleet,list(fleet$Our.Class),function(subset)weighted.mean(subset$F

ootprint, subset$Volume.2021), simplify=TRUE))

colnames(man.footprint)<- c("Avg.Footprint","Our.Class")

•tyu2<-merge(man.sales,fleet)•

•tyu2$man.proportion<-tyu2$Volume.2021/tyu2$Man.2021.Sales

•tyu2$wt.fp.proportion<-tyu2$Curb.Weight/tyu2$Footprint*tyu2$man.proportion

•manwtfp<-aggregate((tyu2$wt.fp.proportion), by=list(tyu2$Our.Class), sum)

•colnames(manwtfp) <- c("Our.Class", "Avg.Density")

•manwtfp<- transform(manwtfp, Our.Class=reorder(Our.Class,Avg.Density))

•testmerge<-merge(manwtfp,man.sales)

•manwtfp<-merge(testmerge, man.footprint)

•#wt/fp averages

•

ct<-signif(sum(cars.egr[[1]]$Curb.Weight/cars.egr[[1]]$Footprint*cars.egr[[1]]$Weight),dig

its=3)

dt<-signif(sum(trucks.egr[[1]]$Curb.Weight/trucks.egr[[1]]$Footprint*trucks.egr[[1]]$Weigh

t), digits=3)

•ggplot(data=manwtfp,aes(x=Avg.Footprint, y=Avg.Density,label=Our.Class))+

••geom_point(size=1)+

••geom_text(hjust=-.05, angle=30, size=3)+

••labs(x="Footprint",y="WT/FP")+

••opts(title=paste("WT/FP by Vehicle Class"))+

••opts(axis.text.x=theme_text(hjust=1))+

••geom_hline(yintercept=ct, linetype=2, colour="red")+

••annotate(geom="text", x = 65, y= 72, colour="red",label=paste("Car fleet avg.:",ct), 

size=3)+

••geom_hline(yintercept=dt, linetype=2, colour="blue")+

••annotate(geom="text", x = 65, y= 84, colour="blue",label=paste("Truck fleet avg.:",dt), 

size=3)+

••scale_x_continuous(limits=c(40,70))+

••scale_y_continuous(limits=c(65,100))

••}

tester3<-function(df){lm(Curb.Weight ~ Footprint, weights=Volume.2021, data=df)}

phys1a<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester3)[[1]]

phys1b<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester3)[[2]]

tester2<-function(df){lm(Curb.Weight/Footprint ~ Footprint, weights=Volume.2021, data=df)}

phys2a<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester2)[[1]]

phys2b<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester2)[[2]]

tester<-function(df){lm(Horsepower/Curb.Weight ~ Footprint, weights=Volume.2021, data=df)}

phys3a<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester)[[1]]

phys3b<-dlply(fleet4,.(NHTSA.Defined.New.NHTSA.Car.Truck),tester)[[2]]

lapply(list(phys1a,phys1b,phys2a,phys2b,phys3a,phys3b),summary)

pdf(paste("2010_curve_analysis-",Sys.Date( ),".pdf"))

textplot(all_regressions,cex=0.5,halign=c("left"))

fphist

dchart("2010","CO2",0)•

dchart("No.Tech","No.Tech.CO2",1)

dchart("Tech-Normalized","CO2.Package.1",2)

dchart("HP/WT-Tech-Normalized","hp.wt.CO2",3)

dchart("WT/FP-Tech-Normalized","crbwt.norm.CO2",4)

dchart("HP/WT&WT/FP-Tech-Normalized","hpwt.crbwt.norm.CO2",5)

physics()

vcplot()

dev.off()