Modeling Summary.Rmd

---
title: "Modeling spatial variation in energy use"
author: "Malcolm Morgan"
date: "11 June 2019"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Setup

bringing in the basic settings we need for R.

```{r p1, include=FALSE}
secure_path <- "E:/Users/earmmor/OneDrive - University of Leeds/CREDS Data/Tim Share"
#secure_path <- "D:/Users/earmmor/OneDrive - University of Leeds/CREDS Data/Tim Share"
```

## Input data

First load the data, some is in England and Wales only so we will start with just there. The datasets in use are:

1. LSOA boundaries
1. Annual Gas & Electricity domestic usage for LSOA
1. 2011 MOT data for miles driven by car (also have vans but not using)
1. Building age (based on counts of bands but estimating mean age)
1. Basic census demographics and building types (e.g. semi-detached)
1. Heating types from census
1. EPC average score (aggregated up from OA)
1. Population counts and population density
1. Number of rooms/bedrooms from census
1. Average emissions factors for driving

```{r p2, include=FALSE}
library(sf)
library(dplyr)
library(Hmisc)
bounds <- st_read("data-prepared/LSOA_generalised.gpkg")
elec <- read.csv("data-prepared/Electricty_2010-17.csv", stringsAsFactors = FALSE)
elec <- elec[!duplicated(elec$LSOA11),] 
gas <- read.csv("data-prepared/Gas_2010-17.csv", stringsAsFactors = FALSE)
mot <- read.csv(paste0(secure_path,"/From Tim/MOT Data RACv9.3/MOT Data RACv9.3 LSOAoutputs_2011.csv"), stringsAsFactors = FALSE)
age <- readRDS("data-prepared/age.Rds") # Not full UK
census <- readRDS("data-prepared/census_lsoa.Rds") # Not full UK
heating <- readRDS("data-prepared/central_heating.Rds") # Not full UK
EPC <- readRDS("data-prepared/EPC.Rds")
population <- readRDS("data-prepared/population.Rds") # Not full UK
density <- readRDS("data-prepared/populationDensity.Rds")
rooms <- readRDS("data-prepared/rooms.Rds")

dir.create("temp")
unzip("../Secure-Data/Excess/XSExpData1.zip", exdir = "temp")
tim <- read.csv("temp/XSExpData1.csv", stringsAsFactors = FALSE)
unlink("temp", recursive = T)

dir.create("temp")
unzip("../Secure-Data/Excess/Experian.zip", exdir = "temp")
experian <- read.csv("temp/UKDA-5738-csv/csv/2011-experian-data.csv", stringsAsFactors = FALSE)
unlink("temp", recursive = T)
names(experian) <- experian[4,]
experian <- experian[5:nrow(experian),]
experian <- experian[,c("GeographyValue","Median_(H) Household Income Value")]
names(experian) <- c("LSOA","median_household_income")
emissions <- readRDS("E:/Users/earmmor/OneDrive - University of Leeds/CREDS Data/github-secure-data/lsoa_emissions.Rds")
latitude <- st_centroid(bounds)
latitude <- cbind(st_drop_geometry(latitude), st_coordinates(latitude))
latitude <- latitude[,c("LSOA11","Y")]
names(latitude) <- c("LSOA11","northing")
```

We will subsets the data and join together into a single master table. When possible converting values to percentages to aid comparison and prevent model distortions.

```{r  p3, include=FALSE}
bounds <- bounds[bounds$LSOA11 %in% census$CODE,]
elec <- elec[elec$LSOA11 %in% bounds$LSOA11, c("LSOA11","DomMet_11",
                                        "DomMet_17","MeanDomElec_11_kWh",
                                        "MeanDomElec_17_kWh",
                                        "TotDomElec_11_kWh","TotDomElec_17_kWh")]
gas <- gas[gas$LSOA11 %in% bounds$LSOA11, c("LSOA11","GasMet_11", "GasMet_17",
                                            "MeanDomGas_11_kWh","MeanDomGas_17_kWh",
                                            "TotDomGas_11_kWh", "TotDomGas_17_kWh")]
mot <- mot[mot$LSOA %in% bounds$LSOA11,]
names(mot) <- c("LSOA", "cars_total","cars_miles","pu5k","p5_12k","po12k","age_av","miles_av_u3",
                     "miles_av_o13","pcars_diesel","pmiles_diesel","vans_total","vans_miles",
                     "pmiles_car","pmiles_vans","cars_percap","miles_percap")
age <- age[age$lsoa %in% bounds$LSOA11, c("lsoa","BP_PRE_1900","BP_1900_1918","BP_1919_1929",
                                         "BP_1930_1939","BP_1945_1954","BP_1955_1964",
                                         "BP_1965_1972","BP_1973_1982","BP_1983_1992",
                                         "BP_1993_1999","BP_2000_2009","BP_2010_2015",
                                         "BP_UNKNOWN","ALL_PROPERTIES")]
census <- census[census$CODE %in% bounds$LSOA11,]
heating <- heating[heating$LSOA11 %in% bounds$LSOA11,]
EPC <- EPC[EPC$LSOA11 %in% bounds$LSOA11, c("LSOA11","Dwellings","Crr_EE","Ptn_EE")]
population <- population[population$LSOA11 %in% bounds$LSOA11,]
population$pop2011 <- as.numeric(population$pop2011)
population$pop2016 <- as.numeric(population$pop2016)
density$dense_2017 <- as.numeric(density$dense_2017)
density <- density[density$LSOA11 %in% bounds$LSOA11, ]



# Some values are more meaningful as percentages or averages
age$mean_house_age <- sapply(1:nrow(age), function(i){
  sub <- age[i,2:13]
  sub <- as.numeric(sub)
  res <- weighted.mean(x = c(1850,1909,1924,1934.5,1949.5,1959.5,1968.5,
                             1977.5,1987.5,1996,2004.5,2012.5) ,w = sub)
  return(res)
})

age$pP1900  <- round(age$BP_PRE_1900 / age$ALL_PROPERTIES * 100,2)
age$p1900_18 <- round(age$BP_1900_1918 / age$ALL_PROPERTIES * 100,2)
age$p1919_29 <- round(age$BP_1919_1929 / age$ALL_PROPERTIES * 100,2)
age$p1930_39 <- round(age$BP_1930_1939 / age$ALL_PROPERTIES * 100,2)
age$p1945_54 <- round(age$BP_1945_1954 / age$ALL_PROPERTIES * 100,2)
age$p1955_64 <- round(age$BP_1955_1964 / age$ALL_PROPERTIES * 100,2)
age$p1965_72 <- round(age$BP_1965_1972 / age$ALL_PROPERTIES * 100,2)
age$p1973_82 <- round(age$BP_1973_1982 / age$ALL_PROPERTIES * 100,2)
age$p1983_92 <- round(age$BP_1983_1992 / age$ALL_PROPERTIES * 100,2)
age$p1993_99 <- round(age$BP_1993_1999 / age$ALL_PROPERTIES * 100,2)
age$p2000_09 <- round(age$BP_2000_2009 / age$ALL_PROPERTIES * 100,2)
age$p2010_15 <- round(age$BP_2010_2015 / age$ALL_PROPERTIES * 100,2)
age$pUNKNOWN   <- round(age$BP_UNKNOWN / age$ALL_PROPERTIES * 100,2)

age <- age[,c("lsoa","mean_house_age","pP1900","p1900_18","p1919_29",
              "p1930_39","p1945_54","p1955_64","p1965_72","p1973_82",
              "p1983_92","p1993_99","p2000_09","p2010_15","pUNKNOWN")]

heating$pHeating_None <- round(heating$`No CH` / heating$All * 100,2)
heating$pHeating_Gas <- round(heating$Gas / heating$All * 100,2)
heating$pHeating_Electric <- round(heating$Electric / heating$All * 100,2)
heating$pHeating_Other <- round((heating$Oil + heating$`Solid fuel` + heating$Other)/ heating$All * 100,2)

heating <- heating[,c("LSOA11","pHeating_None","pHeating_Gas","pHeating_Electric","pHeating_Other")]

# Join Togther
all <- left_join(elec, gas, by = "LSOA11")
all <- left_join(all, mot, by = c("LSOA11" = "LSOA"))
all <- left_join(all, age, by = c("LSOA11" = "lsoa"))
all <- left_join(all, census, by = c("LSOA11" = "CODE"))
all <- left_join(all, heating, by = c("LSOA11" = "LSOA11"))
all <- left_join(all, EPC, by = c("LSOA11" = "LSOA11"))
all <- left_join(all, population, by = c("LSOA11" = "LSOA11"))
all <- left_join(all, density, by = c("LSOA11" = "LSOA11"))
all <- left_join(all, experian, by = c("LSOA11" = "LSOA"))
all <- left_join(all, rooms, by = c("LSOA11" = "geography.code"))
all <- left_join(all, emissions, by = c("LSOA11" = "LSOA"))
all <- left_join(all, latitude, by = c("LSOA11" = "LSOA11"))
rm(age, census, density, elec, EPC, gas, heating, mot, population, experian, rooms, emissions,latitude)

all$median_household_income <- as.numeric(all$median_household_income)
# Get average car emissions, by weighting miles by emissions

all$petrol_emissions <- all$cars_miles * 1.60934 * (1 - all$pmiles_diesel) * all$petrol_co2 /1000 
all$diesel_emissions <- all$cars_miles * 1.60934 * (all$pmiles_diesel) * all$diesel_co2 / 1000
# Convert from kgco2 to litres of fuel
all$petrol_litres <- all$petrol_emissions / 2.18943
all$diesel_litres <- all$diesel_emissions / 2.58935
# Convert to kWh
all$petrol_kwh <- all$petrol_litres * 9.42
all$diesel_kwh <- all$petrol_litres * 10.89
# Total
all$driving_kwh <- all$petrol_kwh + all$diesel_kwh
all$driving_kwh_percap <- all$driving_kwh / all$pop2011

```


## Examining Correlations

Before getting into the detail lets Let see which variables are correlated with energy use. This gives an overall idea of what matters.

### Gas

```{r p4, include=FALSE}
all_matrix <- all[,2:ncol(all)]
all_matrix <- data.matrix(all_matrix)

correlations <- rcorr(x = all_matrix, type="pearson")
correlations_gas <- data.frame(R = correlations$r[,"MeanDomGas_11_kWh"], P = correlations$P[,"MeanDomGas_11_kWh"])
correlations_gas <- correlations_gas[order(correlations_gas$R),]
correlations_gas <- correlations_gas[correlations_gas$R > 0.3 | correlations_gas$R < -0.3,]
```

A table of the top correlations, showing persons R and P-values. Correlating with 2011 gas consumption as most data (e.g. census if from that year)

```{r  p5, echo=FALSE}
correlations_gas
```

Some expected results

```{r  p5a, echo=FALSE}
all$northing_km <- all$northing / 1000

ggplot(all, aes(northing_km, MeanDomGas_11_kWh)) +
  geom_point() +
  geom_smooth(method = "lm") +
  xlab("Kilometres North") +
  ylab("Mean annual gas consumption kWh (2011)") +
  geom_text(x = 400, y = 40000, label = "y = 241 + 0.002496x,  r2 = 0.004", parse = FALSE) +
  ggsave("plots/northing_gas_consumption.png")
```



The strongest correlation is with the number of rooms. This is unsurprising as bigger houses require more heating. The next strongest correlation is with the number of bedrooms, but as the number of bedrooms is closely related to the number of rooms this is not very informative.

Other strong correlations are % of AB social grade people being strongly positive and % DE Social grade people being strongly negative. Also note that working from home and being self-employed is predictive of higher gas usage, perhaps showing households that are heated all day?

We could try a model of the top 5 variaibles
```{r p6}
#gas_lm0 = lm(MeanDomGas_11_kWh ~ mean_rooms + SocGrade_DE. + mean_bedrooms + SocGrade_AB. + median_household_income, data=all, na.action = na.exclude)
# summary(gas_lm0)

gas_lmoa = lm(all$TotDomGas_11_kWh ~ all$mean_rooms)
gas_lmob = lm(all$TotDomGas_11_kWh ~ poly(all$mean_rooms, 2))
summary(gas_lmoa)
summary(gas_lmob)
```

It would be more interesting to remove the effect of the number of rooms and see what is strongly correlated with the residuals.

```{r  p7, echo=FALSE}
plot(all$mean_rooms, all$MeanDomGas_11_kWh, 
     ylab = "Mean Domestic Gas consumption 2011 (kWh)",
     xlab = "Mean number of rooms per house")

ggplot(all, aes(mean_rooms, MeanDomGas_11_kWh)) +
  geom_point() +
  geom_smooth(method = "lm") +
  xlab("Mean number of rooms per house") +
  ylab("Mean Domestic Gas consumption 2011 (kWh)") +
  ggsave("plots/rooms_gas_consumption.png")
```

Let's look at what best predicts the residuals after accounting for the number of rooms.

```{r  p8, echo=FALSE}
gas_lm1 = lm(MeanDomGas_11_kWh ~ mean_rooms, data=all, na.action = na.exclude)
#gas_lm1 = lm(TotDomGas_11_kWh ~ poly(mean_rooms, 2), data=all, na.action = na.exclude)
summary(gas_lm1)
plot(predict(gas_lm1),all$MeanDomGas_11_kWh,
     xlab="predicted",ylab="actual")
abline(a=0,b=1, col = "red")

all$gas_lm1_res = residuals(gas_lm1, na.action = na.exclude)

all_matrix <- all[,2:ncol(all)]
all_matrix <- data.matrix(all_matrix)

correlations <- rcorr(x = all_matrix, type="pearson")
correlations_gas_res <- data.frame(R = correlations$r[,"gas_lm1_res"], P = correlations$P[,"gas_lm1_res"])
correlations_gas_res <- correlations_gas_res[order(correlations_gas_res$R),]
correlations_gas_res <- correlations_gas_res[correlations_gas_res$R > 0.3 | correlations_gas_res$R < -0.3,]
correlations_gas_res
```

Some of these are more related to car use (e.g. T2W_Car is the % of people travelling to work by car). Of the ones about houses, SocialGrade_C2 is the strongest correlation. This is interesting as it seems that by accounting for the number of rooms we have removed some but not all the of the social grade effect. Could C2 class people be slightly different from the other social classes?

Plots of the relationship between the proportion of people with class C2 and gas consumption. And the correlation between residuals for the model and social grade C2.

```{r  p9, echo=FALSE}
plot(all$SocGrade_C2., all$MeanDomGas_11_kWh, 
     ylab = "Mean Domestic Gas consumption 2011 (kWh)",
     xlab = "Social Grade C2")

plot(all$gas_lm1_res, all$SocGrade_C2., 
     ylab = "Number of rooms residuals",
     xlab = "Social Grade C2")
```

So let try a new model with the two variables. This improved the model from an  R squared of 0.54 to 0.63

```{r  p10, echo=FALSE}
gas_lm2 = lm(MeanDomGas_11_kWh ~ mean_rooms + SocGrade_C2., data=all, na.action = na.exclude)
summary(gas_lm2)

plot(predict(gas_lm2),all$MeanDomGas_11_kWh,
     xlab="predicted",ylab="actual")
abline(a=0,b=1, col = "red")
```


We can repeat the process, of looking at the residuals.

```{r p11, echo=FALSE}
all$gas_lm2_res = residuals(gas_lm2, na.action = na.exclude)

all_matrix <- all[,2:ncol(all)]
all_matrix <- data.matrix(all_matrix)

correlations <- rcorr(x = all_matrix, type="pearson")
correlations_gas_res <- data.frame(R = correlations$r[,"gas_lm2_res"], P = correlations$P[,"gas_lm2_res"])
correlations_gas_res <- correlations_gas_res[order(correlations_gas_res$R),]
correlations_gas_res <- correlations_gas_res[correlations_gas_res$R > 0.2 | correlations_gas_res$R < -0.2,]
correlations_gas_res
```

Hear some other appropriate variables such as the EPC rating (Crr_EE) and building age (mean_house_age) stand out as we the proportion of self-employed people.

Let try a final model with all these added variables.

```{r p12, echo=FALSE}
gas_lm3 = lm(MeanDomGas_11_kWh ~ mean_rooms + SocGrade_C2. + Crr_EE + Self.Emp. + mean_house_age, data=all, na.action = na.exclude)
summary(gas_lm3)

plot(predict(gas_lm3),all$MeanDomGas_11_kWh,
     xlab="predicted",ylab="actual")
abline(a=0,b=1, col = "red")
```

So we have a model that can explain 68.9% of the variation in gas usage at LSOA level based on just 5 variables. Gas is tricky due to the off-gas grid areas which we have not properly captured, although the type of heating has not been very predictive.

### Kitchen Sink Approach

Let's throw the kitchen sink at the data that is even slightly correlated and see what is the best model we can get.

```{r p13, echo=FALSE, eval=FALSE}
all_gas <- all[,rownames(correlations_gas)]
all_gas <- all_gas[, !names(all_gas) %in% c("MeanDomElec_17_kWh","MeanDomElec_11_kWh","MeanDomGas_17_kWh",
                                            "TotDomGas_17_kWh","TotDomGas_11_kWh","electric diesel_n")]

all_gas_sub <- dplyr::sample_n(all_gas, 10000)



all_gas_sub <- as.data.frame(lapply(all_gas_sub, scale))
gas_lm4 = lm(MeanDomGas_11_kWh ~ (.)^2, data=all_gas_sub, na.action = na.exclude)
coeff <- summary(gas_lm4)$coefficients
coeff <- as.data.frame(coeff)
coeff <- coeff[coeff$`Pr(>|t|)` > 0.75,]

#plot(all$TotDomElec_11_kWh, all$diesel_co2)

```


## Electricity

Let's do the same for electricity.

```{r p14}
all_matrix <- all[,2:ncol(all)]
all_matrix <- data.matrix(all_matrix)

correlations <- rcorr(x = all_matrix, type="pearson")
correlations_elec <- data.frame(R = correlations$r[,"MeanDomElec_11_kWh"], P = correlations$P[,"MeanDomElec_11_kWh"])
correlations_elec <- correlations_elec[order(correlations_elec$R),]
correlations_elec <- correlations_elec[correlations_elec$R > 0.5 | correlations_elec$R < -0.5,]
correlations_elec
```

Here we see the top factors are (excluding car related ones) Number of rooms, Self Employed, and % gas heating, detached houses, and % other heating. Again income matters but is not the top variable. That gas heating reduces electricity demand is clear, but why does electric heating not increase demand? A simple model based on a few top variables gets an R squared of 0.74.

```{r p15}
elec_lm0 = lm(MeanDomElec_11_kWh ~ mean_rooms + Self.Emp. + pHeating_Gas + Whole_House_Detached. + pHeating_Other, 
              data=all, na.action = na.exclude)
summary(elec_lm0)
```

What about the residuals after accounting for the number of rooms?

```{r  p16}
elec_lm1 = lm(MeanDomElec_11_kWh ~ mean_rooms, data=all, na.action = na.exclude)
summary(elec_lm1)

```

```{r p17, echo=FALSE}
plot(predict(elec_lm1),all$MeanDomElec_11_kWh,
     xlab="predicted",ylab="actual")
abline(a=0,b=1, col = "red")

all$elec_lm1_res = residuals(elec_lm1, na.action = na.exclude)

all_matrix <- all[,2:ncol(all)]
all_matrix <- data.matrix(all_matrix)

correlations <- rcorr(x = all_matrix, type="pearson")
correlations_elec_res <- data.frame(R = correlations$r[,"elec_lm1_res"], P = correlations$P[,"elec_lm1_res"])
correlations_elec_res <- correlations_elec_res[order(correlations_elec_res$R),]
correlations_elec_res <- correlations_elec_res[correlations_elec_res$R > 0.3 | correlations_elec_res$R < -0.3,]
correlations_elec_res
```

Here the number of rooms is much less predictive than for gas. This suggests that other behaviours matter more than house size. Also, note there is a lot less variation between LSOAs in terms of their electricity usage.

The top residuals are % of different heating types, self-employed and working from home.

```{r p18}
elec_lm2 = lm(MeanDomElec_11_kWh ~ mean_rooms + Self.Emp. + pHeating_Gas + pHeating_Electric + pHeating_Other + T2W_Home., 
              data=all, na.action = na.exclude)
summary(elec_lm2)
```

We can get an R squared of 0.75 which is probably pushing the upper limit of what is possible, considering the quality of the data and that there will be some inherent randomness.

## Driving

Finally, let's look at driving. The top correlations are.

```{r p19, echo=FALSE}
all_matrix <- all[,2:ncol(all)]
all_matrix <- data.matrix(all_matrix)

correlations <- rcorr(x = all_matrix, type="pearson")
correlations_cars <- data.frame(R = correlations$r[,"miles_percap"], P = correlations$P[,"miles_percap"])
correlations_cars <- correlations_cars[order(correlations_cars$R),]
correlations_cars <- correlations_cars[correlations_cars$R > 0.5 | correlations_cars$R < -0.5,]
correlations_cars
```

The top correlations are quite broad but include household without cars, unemployment, and population density, car ownership per capita, number of rooms and bedrooms, and % travel to work by car.

```{r p20}
cars_lm0 = lm(miles_percap ~ cars_percap + T2W_Car. + NoCarsHH + Unemployed. + dense_2011, 
              data=all, na.action = na.exclude)
summary(cars_lm0)

plot(all$cars_percap, all$miles_percap,
     xlab="Cars per person",ylab="Miles per person",
     xlim = c(0,2), ylim = c(0,20000))


plot(predict(cars_lm0),all$miles_percap,
     xlab="predicted",ylab="actual", xlim = c(0,20000), ylim = c(0,20000))
abline(a=0,b=1, col = "red")
```

There is a very strong correlation between cars per person and miles driven per person, suggesting that a car owner is a car driven. It is perhaps unsurprising that people do not own a lot of cars they don't need, conversely it is impossible to drive a non-existent car. So let's remove the car ownership to get a clearer idea of what else matters.

```{r p21}
cars_lm1 = lm(miles_percap ~  T2W_Car. + NoCarsHH + Unemployed. + dense_2011, 
              data=all, na.action = na.exclude)
summary(cars_lm1)

plot(predict(cars_lm1),all$miles_percap,
     xlab="predicted",ylab="actual", xlim = c(0,20000), ylim = c(0,20000))
abline(a=0,b=1, col = "red")
```

Still a reasonably good fit (R^2 of 0.79), but there is some clear non-linerality to this relationship.

```{r  p22}
pairs(all[all$miles_percap < 30000,
          c("miles_percap", "T2W_Car.", "NoCarsHH" , "Unemployed.", "dense_2011")])
```

It seems that density has a non-linear effect on driving.

```{r  p23}
plot(all$dense_2011, all$miles_percap,
     xlab="People per km ^2",ylab="Miles per person",
     ylim = c(0, 20000))
# fit non-linear model
cars_lm2a = lm(miles_percap ~  dense_2011, data=all, na.action = na.exclude)
cars_lm2b = nls(miles_percap ~ a * exp(b * dense_2011), data = all, start = list(a = 8000, b = -7e-5),
                na.action = na.exclude)
cars_lm2c = lm(miles_percap ~ exp(-7.25e-5 * dense_2011), data=all, na.action = na.exclude)
summary(cars_lm2a)
summary(cars_lm2c)
```

If we put that back into the original model.

```{r  p24}
cars_lm3 = lm(miles_percap ~  T2W_Car. + NoCarsHH + Unemployed. +                  exp(-7.25e-5 * dense_2011), 
              data=all, na.action = na.exclude)
summary(cars_lm3)

plot(predict(cars_lm3),all$miles_percap,
     xlab="predicted",ylab="actual", xlim = c(0,20000), ylim = c(0,20000))
abline(a=0,b=1, col = "red")
```

It didn't help very much :(

So try a polynomial on density and No car households

```{r  p25}
cars_lm4 = lm(miles_percap ~  T2W_Car. + poly(NoCarsHH, 2) + Unemployed. + poly(dense_2011, 2), 
              data=all, na.action = na.exclude)
summary(cars_lm4)

plot(predict(cars_lm4),all$miles_percap,
     xlab="predicted",ylab="actual", xlim = c(0,20000), ylim = c(0,20000))
abline(a=0,b=1, col = "red")
```

Much better but still seeing an S-bend to the data?

## Driving Energy consumption

How about looking at energy consumption rather than miles, does that make a difference.

```{r  p26, echo=FALSE}
all_matrix <- all[,2:ncol(all)]
all_matrix <- data.matrix(all_matrix)

correlations <- rcorr(x = all_matrix, type="pearson")
correlations_carskwh <- data.frame(R = correlations$r[,"driving_kwh_percap"], P = correlations$P[,"driving_kwh_percap"])
correlations_carskwh <- correlations_carskwh[order(correlations_carskwh$R),]
correlations_carskwh <- correlations_carskwh[correlations_carskwh$R > 0.5 | correlations_carskwh$R < -0.5,]
correlations_carskwh
```

Similar pattern to miles per capita, probably because variation in miles is much greater than variation in fuel efficnency. Travel to Work by Bus comes out as 


```{r  p27}
cars_lm5 = lm(driving_kwh_percap ~ miles_percap +  cars_percap + T2W_Car. + NoCarsHH + Unemployed. + dense_2011, 
              data=all, na.action = na.exclude)
summary(cars_lm5)

plot(all$miles_percap, all$driving_kwh_percap,
     xlab="Miles per person",ylab="kWh per person",
     xlim = c(0,20000), ylim = c(0, 12000)
     )


plot(predict(cars_lm5),all$driving_kwh_percap,
     xlab="predicted",ylab="actual", xlim = c(0,12000), ylim = c(0,12000))
abline(a=0,b=1, col = "red")
```