Updated Townsend scores and regression
Ellen Talbot
20/03/2018
Packages
Data prep
Load data for Townsend calculation
Load data for postcode sector reweighting
- Load data for accomodation analysis
Accomodation variable as count per PCS to calculate largest group
census <- read.csv("~/Google_Drive/PhD/Data/DEP/Output_data/Census_clean.csv")[,1:6] #Just accomodation variables before proportions were calculated
census$total <- rowSums(census[,2:6])
results <- as.data.frame(t(colSums(census[,2:7], na.rm=TRUE)))
#as.data.frame(t(results)) -> results
results %<>%
mutate(detpc = ACCOM_unshared_house_detached/total *100) %>%
mutate(semipc = ACCOM_unshared_house_semi/total *100) %>%
mutate(terrpc = ACCOM_unshared_house_terrace/total *100) %>%
mutate(flatpc = ACCOM_flat/total *100) %>%
mutate(sharepc = ACCOM_shared/total *100)
results <- as.list(results[,7:11])- Load energy data
Calculations
Calculate postcode sector proportions for townsend score by merging PCS weights and townsend scores
#Merge weights onto census data
Census_and_props <- as.data.frame(merge(PCS_weights,td, by.x="OA_CODE", by.y="GEO_CODE", all.x=TRUE))
#Create new dataset containing only proportions
Weighted_variables <- Census_and_props
for(colnum in c(6:9)){
Weighted_variables[,colnum] <- round((Census_and_props[,colnum]*Census_and_props$prop_OA),2)
}
PCS_TD <- Weighted_variables %>%
group_by(PCS) %>%
summarise_each(funs(sum), -OA_CODE, -prop_OA, -perc_OA)## `summarise_each()` is deprecated.
## Use `summarise_all()`, `summarise_at()` or `summarise_if()` instead.
## To map `funs` over a selection of variables, use `summarise_at()`Calculate Z score
PCS_TD$zCar <- scale(PCS_TD$carpc, center = TRUE, scale = TRUE)
PCS_TD$zOvc <- scale(PCS_TD$ovcpc, center = TRUE, scale = TRUE)
PCS_TD$zTen <- scale(PCS_TD$tenpc, center = TRUE, scale = TRUE)
PCS_TD$zEau <- scale(PCS_TD$eaupc, center = TRUE, scale = TRUE)
# Combine z-scores in to one score
PCS_TD$z <- rowSums(PCS_TD[c("zCar", "zOvc", "zTen", "zEau")])
# Drop unnecessary items
PCS_TD <- select(PCS_TD, PCS, z)Merge energy data and weighted TD score
ENERGY_TD <- merge(energy, PCS_TD, by="PCS", all.x=TRUE)
ENERGY_TD <- ENERGY_TD %>%
select(-readdate,-smart_meter,-day,-month, -row_mean) %>%
group_by(PCS) %>%
summarise_all(mean) # this also means that the z-score is now a mean of the original.
ENERGY_TD$rowtot <- rowSums(ENERGY_TD[,2:49])
ENERGY_TD$rowmean <- rowMeans(ENERGY_TD[,2:49])Townsend score regression model
ENERGY_TD1 <- ENERGY_TD
ENERGY_TD1$PCS <- NULL
hist(ENERGY_TD1$rowtot, breaks = 100)
hist(ENERGY_TD1$z, breaks = 100)
TD.mod1 <- lm(rowtot ~ z , ENERGY_TD1, na.action = na.exclude) #Daily total ~ Townsend score
summary(TD.mod1)##
## Call:
## lm(formula = rowtot ~ z, data = ENERGY_TD1, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6095.0 -501.5 -137.3 287.0 11320.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5934.405 11.989 495.00 <2e-16 ***
## z -56.223 3.752 -14.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1042 on 7559 degrees of freedom
## (1159 observations deleted due to missingness)
## Multiple R-squared: 0.02885, Adjusted R-squared: 0.02872
## F-statistic: 224.6 on 1 and 7559 DF, p-value: < 2.2e-16TD.mod2 <- lm(rowmean ~ z -1, ENERGY_TD1, na.action = na.exclude) #Daily total ~ Townsend score with no intercept
summary(TD.mod2)##
## Call:
## lm(formula = rowmean ~ z - 1, data = ENERGY_TD1, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4 112.9 120.5 129.3 361.9
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## z -0.5486 0.4517 -1.214 0.225
##
## Residual standard error: 125.5 on 7560 degrees of freedom
## (1159 observations deleted due to missingness)
## Multiple R-squared: 0.000195, Adjusted R-squared: 6.276e-05
## F-statistic: 1.475 on 1 and 7560 DF, p-value: 0.2247Townsend score and accomodation variable regression model
TD_ACCOM <- merge(ENERGY_TD, subset, by.x="PCS", by.y="RMSect", all.x=TRUE, na.rm=TRUE)
TD_AC.mod1 <- lm(rowtot ~ z + detached_pc + semi_pc + terrace_pc + flat_pc + shared_pc, TD_ACCOM, na.action = na.exclude) #Daily total ~ Townsend score and all accomodation variables without intercept
summary(TD_AC.mod1)##
## Call:
## lm(formula = rowtot ~ z + detached_pc + semi_pc + terrace_pc +
## flat_pc + shared_pc, data = TD_ACCOM, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6062.8 -501.7 -91.7 317.2 11108.0
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5499.711 56.744 96.921 < 2e-16 ***
## z -23.678 3.648 -6.490 9.1e-11 ***
## detached_pc 1843.845 75.059 24.565 < 2e-16 ***
## semi_pc -87.500 90.404 -0.968 0.3331
## terrace_pc -239.762 100.208 -2.393 0.0168 *
## flat_pc NA NA NA NA
## shared_pc 3093.735 1320.952 2.342 0.0192 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 972 on 7555 degrees of freedom
## (1159 observations deleted due to missingness)
## Multiple R-squared: 0.1559, Adjusted R-squared: 0.1554
## F-statistic: 279.2 on 5 and 7555 DF, p-value: < 2.2e-16As discussed with Dani - Semi-detached removed to create it as reference category
TD_AC.mod2 <- lm(rowtot ~ z + detached_pc + terrace_pc + flat_pc + shared_pc, TD_ACCOM, na.action = na.exclude) #Daily total ~ Townsend score and all accomodation variables without intercept
summary(TD_AC.mod2)##
## Call:
## lm(formula = rowtot ~ z + detached_pc + terrace_pc + flat_pc +
## shared_pc, data = TD_ACCOM, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6062.8 -501.7 -91.7 317.2 11108.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5412.211 64.950 83.329 < 2e-16 ***
## z -23.678 3.648 -6.490 9.1e-11 ***
## detached_pc 1931.346 102.497 18.843 < 2e-16 ***
## terrace_pc -152.262 113.280 -1.344 0.1790
## flat_pc 87.500 90.404 0.968 0.3331
## shared_pc 3093.735 1320.952 2.342 0.0192 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 972 on 7555 degrees of freedom
## (1159 observations deleted due to missingness)
## Multiple R-squared: 0.1559, Adjusted R-squared: 0.1554
## F-statistic: 279.2 on 5 and 7555 DF, p-value: < 2.2e-16Also as discussed with Dani - combine flat and shared into one category and model ran again - removing semi detached as the biggest group to control for housing type
PCS_level_DT[, -c(2:50,56:80)] -> subset2
subset2 %<>%
mutate(ACCOM_TOTAL = ACCOM_unshared_house_detached+ACCOM_unshared_house_semi+ACCOM_unshared_house_terrace+ACCOM_flat) %>%
mutate(detached_pc = ACCOM_unshared_house_detached/ACCOM_TOTAL) %>%
mutate(semi_pc = ACCOM_unshared_house_semi/ACCOM_TOTAL) %>%
mutate(terrace_pc = ACCOM_unshared_house_terrace/ACCOM_TOTAL) %>%
mutate(flatandshare_pc = (ACCOM_flat+ACCOM_shared)/ACCOM_TOTAL) # combine flat and shared into one category due to small group sizes
subset2[,c(1,7,9:12)] -> subset2
TD_ACCOM2 <- merge(ENERGY_TD, subset2, by.x="PCS", by.y="RMSect", all.x=TRUE, na.rm=TRUE)
TD_AC2.mod1 <- lm(rowtot ~ z + detached_pc + terrace_pc + flatandshare_pc, TD_ACCOM2, na.action = na.exclude) #Daily total ~ Townsend score and all accomodation variables without intercept
summary(TD_AC2.mod1)##
## Call:
## lm(formula = rowtot ~ z + detached_pc + terrace_pc + flatandshare_pc,
## data = TD_ACCOM2, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6060.5 -501.7 -91.7 313.6 11097.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5393.185 64.394 83.753 < 2e-16 ***
## z -23.527 3.648 -6.448 1.2e-10 ***
## detached_pc 1953.272 102.042 19.142 < 2e-16 ***
## terrace_pc -135.597 113.058 -1.199 0.2304
## flatandshare_pc 179.383 80.301 2.234 0.0255 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 972.3 on 7556 degrees of freedom
## (1159 observations deleted due to missingness)
## Multiple R-squared: 0.1554, Adjusted R-squared: 0.1549
## F-statistic: 347.5 on 4 and 7556 DF, p-value: < 2.2e-16Spatial models
# PCT stands for postcode town e.g. L or LA.
TD_ACCOM2$pct. <-substr(TD_ACCOM2$PCS, 1, 2)
TD_ACCOM2$pct. <- gsub('[0-9]+', "", TD_ACCOM2$pct.)
# Drop 3rd character from small postcode sectors in inner London
# Count the number of occurances in each postcode area
sample_size <- as.data.frame(TD_ACCOM2$pct.) %>%
group_by(TD_ACCOM2$pct.) %>%
summarise(count=n())
print("Highest - 352, lowest - 10")## [1] "Highest - 352, lowest - 10"TD_AC2.mod2 <- lm(rowtot ~ pct. + z + detached_pc + terrace_pc + flatandshare_pc -1, TD_ACCOM2, na.action = na.exclude) #Daily total ~ Townsend score and all accomodation variables without intercept
summary(TD_AC2.mod2)##
## Call:
## lm(formula = rowtot ~ pct. + z + detached_pc + terrace_pc + flatandshare_pc -
## 1, data = TD_ACCOM2, na.action = na.exclude)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6218.7 -448.4 -84.2 274.3 11069.4
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## pct.AL 5698.228 171.335 33.258 < 2e-16 ***
## pct.B 5822.464 88.544 65.758 < 2e-16 ***
## pct.BA 5735.909 133.807 42.867 < 2e-16 ***
## pct.BB 5638.697 137.952 40.874 < 2e-16 ***
## pct.BD 5733.605 118.196 48.509 < 2e-16 ***
## pct.BH 5255.641 129.676 40.529 < 2e-16 ***
## pct.BL 5499.152 155.004 35.477 < 2e-16 ***
## pct.BN 5425.537 116.601 46.531 < 2e-16 ***
## pct.BR 5899.638 183.030 32.233 < 2e-16 ***
## pct.BS 5550.410 105.178 52.772 < 2e-16 ***
## pct.CA 5475.920 134.473 40.721 < 2e-16 ***
## pct.CB 5513.511 142.811 38.607 < 2e-16 ***
## pct.CF 5238.597 108.788 48.154 < 2e-16 ***
## pct.CH 5431.123 94.661 57.375 < 2e-16 ***
## pct.CM 5596.910 125.482 44.603 < 2e-16 ***
## pct.CO 5140.006 141.486 36.329 < 2e-16 ***
## pct.CR 6059.143 175.599 34.506 < 2e-16 ***
## pct.CT 5274.154 131.140 40.218 < 2e-16 ***
## pct.CV 5434.131 121.608 44.686 < 2e-16 ***
## pct.CW 5359.866 164.226 32.637 < 2e-16 ***
## pct.DA 6010.425 149.182 40.289 < 2e-16 ***
## pct.DE 4949.381 123.225 40.165 < 2e-16 ***
## pct.DH 5031.296 158.860 31.671 < 2e-16 ***
## pct.DL 5234.711 138.996 37.661 < 2e-16 ***
## pct.DN 4957.389 112.079 44.231 < 2e-16 ***
## pct.DT 5181.921 179.149 28.925 < 2e-16 ***
## pct.DY 5442.688 142.523 38.188 < 2e-16 ***
## pct.E 6215.465 136.839 45.422 < 2e-16 ***
## pct.EC 4455.134 307.627 14.482 < 2e-16 ***
## pct.EN 5954.866 178.684 33.326 < 2e-16 ***
## pct.EX 5314.832 116.998 45.427 < 2e-16 ***
## pct.FY 5290.776 177.181 29.861 < 2e-16 ***
## pct.GL 5593.018 114.836 48.704 < 2e-16 ***
## pct.GU 5584.145 113.301 49.286 < 2e-16 ***
## pct.HA 6275.844 148.568 42.242 < 2e-16 ***
## pct.HD 5465.755 184.030 29.700 < 2e-16 ***
## pct.HG 5356.460 209.587 25.557 < 2e-16 ***
## pct.HP 5629.152 138.759 40.568 < 2e-16 ***
## pct.HR 5292.329 200.097 26.449 < 2e-16 ***
## pct.HU 5188.342 145.032 35.774 < 2e-16 ***
## pct.HX 5484.177 186.942 29.336 < 2e-16 ***
## pct.IG 6581.686 185.572 35.467 < 2e-16 ***
## pct.IP 5552.097 118.933 46.682 < 2e-16 ***
## pct.KT 5861.809 130.755 44.830 < 2e-16 ***
## pct.L 5524.549 81.939 67.423 < 2e-16 ***
## pct.LA 5201.225 141.330 36.802 < 2e-16 ***
## pct.LD 5553.970 278.511 19.942 < 2e-16 ***
## pct.LE 5132.009 111.771 45.915 < 2e-16 ***
## pct.LL 5531.925 110.818 49.919 < 2e-16 ***
## pct.LN 4764.413 164.262 29.005 < 2e-16 ***
## pct.LS 5538.963 111.033 49.886 < 2e-16 ***
## pct.LU 5663.283 189.383 29.904 < 2e-16 ***
## pct.M 5643.027 97.633 57.798 < 2e-16 ***
## pct.ME 5560.847 127.901 43.478 < 2e-16 ***
## pct.MK 5500.551 126.613 43.444 < 2e-16 ***
## pct.N 6155.712 132.266 46.540 < 2e-16 ***
## pct.NE 5543.130 94.148 58.877 < 2e-16 ***
## pct.NG 4844.627 102.985 47.042 < 2e-16 ***
## pct.NN 5355.583 125.835 42.560 < 2e-16 ***
## pct.NP 5155.152 132.450 38.922 < 2e-16 ***
## pct.NR 5277.248 116.932 45.131 < 2e-16 ***
## pct.NW 6218.729 137.511 45.224 < 2e-16 ***
## pct.OL 5578.503 138.909 40.160 < 2e-16 ***
## pct.OX 5873.946 116.322 50.497 < 2e-16 ***
## pct.PE 5258.030 112.330 46.809 < 2e-16 ***
## pct.PL 5567.602 125.497 44.364 < 2e-16 ***
## pct.PO 5265.282 112.989 46.600 < 2e-16 ***
## pct.PR 5472.327 133.046 41.131 < 2e-16 ***
## pct.RG 5856.277 109.517 53.474 < 2e-16 ***
## pct.RH 5532.393 132.601 41.722 < 2e-16 ***
## pct.RM 5885.324 142.120 41.411 < 2e-16 ***
## pct.S 4943.358 89.672 55.127 < 2e-16 ***
## pct.SA 5431.932 107.954 50.317 < 2e-16 ***
## pct.SE 5743.867 122.741 46.797 < 2e-16 ***
## pct.SG 5786.448 144.061 40.167 < 2e-16 ***
## pct.SK 5512.899 118.405 46.560 < 2e-16 ***
## pct.SL 5897.219 150.833 39.098 < 2e-16 ***
## pct.SM 5997.183 205.209 29.225 < 2e-16 ***
## pct.SN 5786.982 133.868 43.229 < 2e-16 ***
## pct.SO 5396.686 115.813 46.598 < 2e-16 ***
## pct.SP 5981.240 167.800 35.645 < 2e-16 ***
## pct.SR 5222.390 168.354 31.020 < 2e-16 ***
## pct.SS 5553.610 139.497 39.812 < 2e-16 ***
## pct.ST 5061.671 123.401 41.018 < 2e-16 ***
## pct.SW 6097.899 124.111 49.132 < 2e-16 ***
## pct.SY 5528.997 132.442 41.746 < 2e-16 ***
## pct.TA 5446.481 142.975 38.094 < 2e-16 ***
## pct.TD 5340.816 540.506 9.881 < 2e-16 ***
## pct.TF 4957.372 178.714 27.739 < 2e-16 ***
## pct.TN 5575.589 113.198 49.255 < 2e-16 ***
## pct.TQ 5213.079 161.646 32.250 < 2e-16 ***
## pct.TR 5882.025 149.326 39.390 < 2e-16 ***
## pct.TS 5166.699 115.394 44.774 < 2e-16 ***
## pct.TW 6007.944 134.791 44.572 < 2e-16 ***
## pct.UB 6256.302 175.227 35.704 < 2e-16 ***
## pct.W 6611.679 142.557 46.379 < 2e-16 ***
## pct.WA 5388.648 121.343 44.408 < 2e-16 ***
## pct.WC 6448.086 362.404 17.793 < 2e-16 ***
## pct.WD 5979.122 164.849 36.270 < 2e-16 ***
## pct.WF 5240.496 129.670 40.414 < 2e-16 ***
## pct.WN 5247.505 177.281 29.600 < 2e-16 ***
## pct.WR 5374.230 163.633 32.843 < 2e-16 ***
## pct.WS 5416.984 137.695 39.340 < 2e-16 ***
## pct.WV 5443.592 143.884 37.833 < 2e-16 ***
## pct.YO 5383.120 113.918 47.254 < 2e-16 ***
## z -20.862 3.928 -5.312 1.12e-07 ***
## detached_pc 1983.692 109.079 18.186 < 2e-16 ***
## terrace_pc -119.111 114.522 -1.040 0.298
## flatandshare_pc -440.659 91.547 -4.813 1.51e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 929.1 on 7452 degrees of freedom
## (1159 observations deleted due to missingness)
## Multiple R-squared: 0.9766, Adjusted R-squared: 0.9762
## F-statistic: 2849 on 109 and 7452 DF, p-value: < 2.2e-16# Effect of z score on energy usage by Postcode area
TD_ACCOM2$one <- 1
TD_AC2.mod3 <- lm('log(rowtot) ~ 0 +(one + z):(pct.)', TD_ACCOM2)
summary(TD_AC2.mod3)##
## Call:
## lm(formula = "log(rowtot) ~ 0 +(one + z):(pct.)", data = TD_ACCOM2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3462 -0.0667 -0.0047 0.0579 1.0529
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## one:pct.AL 8.720e+00 3.543e-02 246.146 < 2e-16 ***
## one:pct.B 8.705e+00 1.202e-02 724.121 < 2e-16 ***
## one:pct.BA 8.726e+00 2.107e-02 414.077 < 2e-16 ***
## one:pct.BB 8.682e+00 2.375e-02 365.586 < 2e-16 ***
## one:pct.BD 8.682e+00 1.876e-02 462.754 < 2e-16 ***
## one:pct.BH 8.682e+00 1.981e-02 438.191 < 2e-16 ***
## one:pct.BL 8.650e+00 2.564e-02 337.323 < 2e-16 ***
## one:pct.BN 8.648e+00 1.707e-02 506.566 < 2e-16 ***
## one:pct.BR 8.717e+00 3.155e-02 276.275 < 2e-16 ***
## one:pct.BS 8.679e+00 1.501e-02 578.172 < 2e-16 ***
## one:pct.CA 8.675e+00 2.556e-02 339.401 < 2e-16 ***
## one:pct.CB 8.715e+00 2.338e-02 372.770 < 2e-16 ***
## one:pct.CF 8.617e+00 1.571e-02 548.409 < 2e-16 ***
## one:pct.CH 8.641e+00 1.897e-02 455.501 < 2e-16 ***
## one:pct.CM 8.723e+00 1.947e-02 448.112 < 2e-16 ***
## one:pct.CO 8.674e+00 2.344e-02 370.067 < 2e-16 ***
## one:pct.CR 8.765e+00 3.279e-02 267.348 < 2e-16 ***
## one:pct.CT 8.641e+00 2.067e-02 418.061 < 2e-16 ***
## one:pct.CV 8.670e+00 1.969e-02 440.287 < 2e-16 ***
## one:pct.CW 8.705e+00 2.809e-02 309.909 < 2e-16 ***
## one:pct.DA 8.716e+00 2.582e-02 337.594 < 2e-16 ***
## one:pct.DE 8.632e+00 2.113e-02 408.479 < 2e-16 ***
## one:pct.DH 8.585e+00 2.943e-02 291.771 < 2e-16 ***
## one:pct.DL 8.634e+00 2.390e-02 361.313 < 2e-16 ***
## one:pct.DN 8.612e+00 1.707e-02 504.557 < 2e-16 ***
## one:pct.DT 8.671e+00 3.277e-02 264.575 < 2e-16 ***
## one:pct.DY 8.685e+00 2.377e-02 365.349 < 2e-16 ***
## one:pct.E 8.625e+00 3.209e-02 268.786 < 2e-16 ***
## one:pct.EC 7.709e+00 3.140e-01 24.552 < 2e-16 ***
## one:pct.EN 8.737e+00 3.124e-02 279.686 < 2e-16 ***
## one:pct.EX 8.682e+00 1.810e-02 479.722 < 2e-16 ***
## one:pct.FY 8.543e+00 3.534e-02 241.755 < 2e-16 ***
## one:pct.GL 8.712e+00 1.730e-02 503.714 < 2e-16 ***
## one:pct.GU 8.733e+00 1.773e-02 492.610 < 2e-16 ***
## one:pct.HA 8.757e+00 2.504e-02 349.788 < 2e-16 ***
## one:pct.HD 8.683e+00 3.463e-02 250.753 < 2e-16 ***
## one:pct.HG 8.678e+00 3.658e-02 237.225 < 2e-16 ***
## one:pct.HP 8.725e+00 2.209e-02 394.890 < 2e-16 ***
## one:pct.HR 8.736e+00 4.006e-02 218.083 < 2e-16 ***
## one:pct.HU 8.600e+00 2.387e-02 360.239 < 2e-16 ***
## one:pct.HX 8.621e+00 3.172e-02 271.778 < 2e-16 ***
## one:pct.IG 8.796e+00 3.291e-02 267.247 < 2e-16 ***
## one:pct.IP 8.706e+00 1.789e-02 486.607 < 2e-16 ***
## one:pct.KT 8.751e+00 2.287e-02 382.608 < 2e-16 ***
## one:pct.L 8.580e+00 1.935e-02 443.426 < 2e-16 ***
## one:pct.LA 8.629e+00 2.476e-02 348.559 < 2e-16 ***
## one:pct.LD 8.691e+00 5.444e-02 159.629 < 2e-16 ***
## one:pct.LE 8.633e+00 1.830e-02 471.613 < 2e-16 ***
## one:pct.LL 8.689e+00 1.639e-02 530.226 < 2e-16 ***
## one:pct.LN 8.648e+00 2.919e-02 296.278 < 2e-16 ***
## one:pct.LS 8.655e+00 1.685e-02 513.571 < 2e-16 ***
## one:pct.LU 8.665e+00 4.155e-02 208.571 < 2e-16 ***
## one:pct.M 8.629e+00 1.421e-02 607.315 < 2e-16 ***
## one:pct.ME 8.685e+00 2.011e-02 431.816 < 2e-16 ***
## one:pct.MK 8.701e+00 2.011e-02 432.668 < 2e-16 ***
## one:pct.N 8.708e+00 2.245e-02 387.931 < 2e-16 ***
## one:pct.NE 8.628e+00 1.330e-02 648.712 < 2e-16 ***
## one:pct.NG 8.584e+00 1.562e-02 549.608 < 2e-16 ***
## one:pct.NN 8.674e+00 1.969e-02 440.535 < 2e-16 ***
## one:pct.NP 8.635e+00 2.097e-02 411.733 < 2e-16 ***
## one:pct.NR 8.692e+00 1.866e-02 465.752 < 2e-16 ***
## one:pct.NW 8.688e+00 2.099e-02 413.883 < 2e-16 ***
## one:pct.OL 8.653e+00 2.342e-02 369.546 < 2e-16 ***
## one:pct.OX 8.725e+00 1.932e-02 451.587 < 2e-16 ***
## one:pct.PE 8.684e+00 1.752e-02 495.698 < 2e-16 ***
## one:pct.PL 8.707e+00 1.898e-02 458.755 < 2e-16 ***
## one:pct.PO 8.640e+00 1.655e-02 522.137 < 2e-16 ***
## one:pct.PR 8.685e+00 2.179e-02 398.486 < 2e-16 ***
## one:pct.RG 8.747e+00 1.687e-02 518.635 < 2e-16 ***
## one:pct.RH 8.709e+00 2.201e-02 395.759 < 2e-16 ***
## one:pct.RM 8.687e+00 2.358e-02 368.385 < 2e-16 ***
## one:pct.S 8.579e+00 1.239e-02 692.193 < 2e-16 ***
## one:pct.SA 8.698e+00 1.541e-02 564.377 < 2e-16 ***
## one:pct.SE 8.612e+00 2.333e-02 369.065 < 2e-16 ***
## one:pct.SG 8.718e+00 2.348e-02 371.324 < 2e-16 ***
## one:pct.SK 8.681e+00 1.832e-02 473.846 < 2e-16 ***
## one:pct.SL 8.756e+00 2.455e-02 356.686 < 2e-16 ***
## one:pct.SM 8.721e+00 3.539e-02 246.428 < 2e-16 ***
## one:pct.SN 8.744e+00 2.105e-02 415.307 < 2e-16 ***
## one:pct.SO 8.698e+00 2.059e-02 422.505 < 2e-16 ***
## one:pct.SP 8.766e+00 2.865e-02 305.922 < 2e-16 ***
## one:pct.SR 8.444e+00 2.915e-02 289.721 < 2e-16 ***
## one:pct.SS 8.691e+00 2.289e-02 379.686 < 2e-16 ***
## one:pct.ST 8.614e+00 1.949e-02 441.928 < 2e-16 ***
## one:pct.SW 8.635e+00 1.749e-02 493.781 < 2e-16 ***
## one:pct.SY 8.725e+00 2.192e-02 398.028 < 2e-16 ***
## one:pct.TA 8.696e+00 2.325e-02 373.989 < 2e-16 ***
## one:pct.TD 8.660e+00 9.923e-02 87.278 < 2e-16 ***
## one:pct.TF 8.647e+00 3.048e-02 283.684 < 2e-16 ***
## one:pct.TN 8.716e+00 1.640e-02 531.519 < 2e-16 ***
## one:pct.TQ 8.650e+00 2.718e-02 318.254 < 2e-16 ***
## one:pct.TR 8.767e+00 2.381e-02 368.226 < 2e-16 ***
## one:pct.TS 8.599e+00 1.903e-02 451.936 < 2e-16 ***
## one:pct.TW 8.711e+00 2.274e-02 382.975 < 2e-16 ***
## one:pct.UB 8.741e+00 3.234e-02 270.282 < 2e-16 ***
## one:pct.W 8.728e+00 2.148e-02 406.377 < 2e-16 ***
## one:pct.WA 8.666e+00 2.068e-02 419.129 < 2e-16 ***
## one:pct.WC 8.963e+00 1.213e-01 73.921 < 2e-16 ***
## one:pct.WD 8.715e+00 3.678e-02 236.959 < 2e-16 ***
## one:pct.WF 8.608e+00 2.217e-02 388.284 < 2e-16 ***
## one:pct.WN 8.622e+00 3.177e-02 271.420 < 2e-16 ***
## one:pct.WR 8.699e+00 2.850e-02 305.274 < 2e-16 ***
## one:pct.WS 8.672e+00 2.306e-02 376.007 < 2e-16 ***
## one:pct.WV 8.675e+00 2.449e-02 354.224 < 2e-16 ***
## one:pct.YO 8.672e+00 1.746e-02 496.645 < 2e-16 ***
## z:pct.AL -3.672e-04 1.373e-02 -0.027 0.978665
## z:pct.B -4.016e-04 3.350e-03 -0.120 0.904594
## z:pct.BA -2.877e-02 8.461e-03 -3.400 0.000678 ***
## z:pct.BB -5.993e-03 8.451e-03 -0.709 0.478222
## z:pct.BD 4.482e-03 6.789e-03 0.660 0.509178
## z:pct.BH -6.049e-03 8.057e-03 -0.751 0.452826
## z:pct.BL -3.527e-04 9.274e-03 -0.038 0.969659
## z:pct.BN -1.017e-02 6.190e-03 -1.643 0.100525
## z:pct.BR -7.428e-03 9.594e-03 -0.774 0.438811
## z:pct.BS -1.167e-02 4.935e-03 -2.364 0.018105 *
## z:pct.CA -2.133e-02 9.979e-03 -2.138 0.032575 *
## z:pct.CB -1.261e-02 5.884e-03 -2.143 0.032115 *
## z:pct.CF -1.185e-02 6.795e-03 -1.743 0.081308 .
## z:pct.CH -1.245e-02 5.743e-03 -2.167 0.030260 *
## z:pct.CM -8.878e-03 8.385e-03 -1.059 0.289720
## z:pct.CO -1.870e-02 8.595e-03 -2.176 0.029592 *
## z:pct.CR -8.899e-03 5.082e-03 -1.751 0.079978 .
## z:pct.CT -1.019e-02 6.926e-03 -1.471 0.141374
## z:pct.CV -7.396e-03 6.318e-03 -1.171 0.241834
## z:pct.CW -1.442e-02 8.431e-03 -1.710 0.087232 .
## z:pct.DA -1.317e-02 1.045e-02 -1.261 0.207526
## z:pct.DE -5.945e-03 4.633e-03 -1.283 0.199478
## z:pct.DH -7.469e-03 1.978e-02 -0.378 0.705788
## z:pct.DL -2.120e-02 1.064e-02 -1.992 0.046387 *
## z:pct.DN -4.499e-03 6.580e-03 -0.684 0.494216
## z:pct.DT -2.645e-02 1.741e-02 -1.519 0.128775
## z:pct.DY -1.471e-02 8.499e-03 -1.731 0.083570 .
## z:pct.E 3.138e-03 2.992e-03 1.049 0.294316
## z:pct.EC -1.462e-01 8.364e-02 -1.747 0.080600 .
## z:pct.EN 5.775e-03 1.186e-02 0.487 0.626367
## z:pct.EX -6.994e-03 8.257e-03 -0.847 0.396987
## z:pct.FY 2.885e-02 1.141e-02 2.527 0.011522 *
## z:pct.GL -1.452e-03 6.756e-03 -0.215 0.829855
## z:pct.GU -1.272e-02 7.271e-03 -1.750 0.080205 .
## z:pct.HA 2.994e-05 1.036e-02 0.003 0.997694
## z:pct.HD -1.213e-02 9.013e-03 -1.346 0.178317
## z:pct.HG -1.813e-02 1.684e-02 -1.076 0.281800
## z:pct.HP -5.338e-03 8.638e-03 -0.618 0.536661
## z:pct.HR -2.522e-02 1.692e-02 -1.491 0.136070
## z:pct.HU 2.687e-03 9.807e-03 0.274 0.784116
## z:pct.HX 2.347e-02 1.099e-02 2.136 0.032684 *
## z:pct.IG -3.804e-03 1.014e-02 -0.375 0.707486
## z:pct.IP -1.300e-02 7.601e-03 -1.711 0.087182 .
## z:pct.KT -2.453e-03 1.059e-02 -0.232 0.816754
## z:pct.L -2.137e-02 6.028e-03 -3.545 0.000395 ***
## z:pct.LA -8.135e-03 1.097e-02 -0.742 0.458208
## z:pct.LD -5.059e-02 2.766e-02 -1.829 0.067445 .
## z:pct.LE -1.529e-03 4.290e-03 -0.356 0.721519
## z:pct.LL -3.907e-02 5.797e-03 -6.740 1.71e-11 ***
## z:pct.LN -1.845e-02 9.450e-03 -1.952 0.050953 .
## z:pct.LS -3.938e-03 4.826e-03 -0.816 0.414566
## z:pct.LU 2.605e-03 6.725e-03 0.387 0.698493
## z:pct.M 9.219e-04 4.370e-03 0.211 0.832909
## z:pct.ME -8.168e-03 8.271e-03 -0.988 0.323423
## z:pct.MK -2.655e-03 6.928e-03 -0.383 0.701615
## z:pct.N -1.322e-02 4.642e-03 -2.848 0.004417 **
## z:pct.NE -2.253e-02 4.111e-03 -5.480 4.39e-08 ***
## z:pct.NG -4.669e-03 5.052e-03 -0.924 0.355424
## z:pct.NN -1.333e-02 6.048e-03 -2.205 0.027518 *
## z:pct.NP -1.167e-02 7.247e-03 -1.611 0.107271
## z:pct.NR -2.072e-02 7.233e-03 -2.864 0.004189 **
## z:pct.NW -1.008e-02 5.458e-03 -1.846 0.064884 .
## z:pct.OL 9.338e-04 8.782e-03 0.106 0.915327
## z:pct.OX -2.784e-02 7.321e-03 -3.803 0.000144 ***
## z:pct.PE 2.272e-03 4.575e-03 0.497 0.619463
## z:pct.PL -1.908e-02 7.003e-03 -2.725 0.006453 **
## z:pct.PO -1.316e-02 6.324e-03 -2.081 0.037488 *
## z:pct.PR -8.216e-03 7.783e-03 -1.056 0.291202
## z:pct.RG -1.749e-02 6.506e-03 -2.688 0.007212 **
## z:pct.RH -1.308e-03 1.153e-02 -0.113 0.909685
## z:pct.RM -4.426e-03 6.528e-03 -0.678 0.497774
## z:pct.S -2.355e-03 3.884e-03 -0.606 0.544276
## z:pct.SA -2.465e-02 6.433e-03 -3.832 0.000128 ***
## z:pct.SE -5.226e-03 4.177e-03 -1.251 0.210932
## z:pct.SG -2.051e-02 1.103e-02 -1.859 0.063074 .
## z:pct.SK -6.484e-03 7.515e-03 -0.863 0.388221
## z:pct.SL 1.149e-03 7.878e-03 0.146 0.884088
## z:pct.SM -6.845e-03 1.213e-02 -0.564 0.572610
## z:pct.SN -1.194e-02 7.682e-03 -1.554 0.120315
## z:pct.SO 9.358e-03 8.171e-03 1.145 0.252122
## z:pct.SP -3.713e-02 1.294e-02 -2.870 0.004116 **
## z:pct.SR 1.313e-01 1.500e-02 8.751 < 2e-16 ***
## z:pct.SS -1.454e-03 8.360e-03 -0.174 0.861909
## z:pct.ST 1.155e-02 9.098e-03 1.269 0.204377
## z:pct.SW 5.466e-03 5.939e-03 0.920 0.357418
## z:pct.SY -2.864e-02 9.799e-03 -2.923 0.003476 **
## z:pct.TA -1.815e-02 1.125e-02 -1.613 0.106690
## z:pct.TD -3.524e-03 1.813e-02 -0.194 0.845915
## z:pct.TF -9.340e-03 1.185e-02 -0.788 0.430475
## z:pct.TN -1.435e-02 5.781e-03 -2.482 0.013085 *
## z:pct.TQ -1.284e-02 9.949e-03 -1.290 0.197025
## z:pct.TR -2.354e-02 1.386e-02 -1.698 0.089546 .
## z:pct.TS -5.562e-03 6.888e-03 -0.807 0.419441
## z:pct.TW -2.022e-03 9.591e-03 -0.211 0.833066
## z:pct.UB -1.587e-03 6.421e-03 -0.247 0.804783
## z:pct.W -2.100e-02 5.506e-03 -3.813 0.000138 ***
## z:pct.WA -8.819e-03 7.755e-03 -1.137 0.255479
## z:pct.WC 1.284e-01 4.722e-02 2.720 0.006543 **
## z:pct.WD -1.577e-02 1.459e-02 -1.080 0.280053
## z:pct.WF 1.093e-02 8.677e-03 1.260 0.207758
## z:pct.WN -4.103e-04 1.104e-02 -0.037 0.970354
## z:pct.WR -1.177e-02 9.805e-03 -1.201 0.229967
## z:pct.WS 2.307e-04 6.478e-03 0.036 0.971593
## z:pct.WV -5.339e-03 8.095e-03 -0.659 0.509597
## z:pct.YO -2.468e-02 7.245e-03 -3.407 0.000661 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1697 on 7351 degrees of freedom
## (1159 observations deleted due to missingness)
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9996
## F-statistic: 9.406e+04 on 210 and 7351 DF, p-value: < 2.2e-16# This allows us to get a seperate constant term and estimate of the impact of the townsend score on energy usage for every postcode town.
#Keep coefficients for the z scores and wrangle into DT
require(broom)## Loading required package: broomrequire(stringr)
mod3.data <- tidy(TD_AC2.mod3)
mod3.data %<>%
filter(str_detect(term, "^z")) #Keep just z score coefs
mod3.data$id <- sub('.*\\.', '', mod3.data$term) # Extract postcode areas for joining to shapefiles
### Visualisation
#Load in shapefile of Postcode town boundaries
require(raster)## Loading required package: raster##
## Attaching package: 'raster'## The following object is masked from 'package:magrittr':
##
## extract## The following object is masked from 'package:data.table':
##
## shift## The following object is masked from 'package:dplyr':
##
## select## The following object is masked from 'package:tidyr':
##
## extractpostalareas <- shapefile("~/Google_Drive/PhD/Data/DEP/Raw_data/PostalSector/PostalArea.shp")
# Join the data to the shapefile
postalareas@data <- merge(postalareas@data, mod3.data, by.x="PostArea", by.y="id", all.x=TRUE, na.rm=TRUE)
# Write the data for plotting in QGIS
#shapefile(postalareas, '~/Google_Drive/PhD/Data/DEP/postal.energy.shp')GOR_PA <- read.csv("~/Google_Drive/PhD/Data/DEP/Output_data/GOR_postalarea.csv")[,c(3,18,19)]
TD_ACCOM3 <- merge(TD_ACCOM2, GOR_PA, by.x="pct.", by.y="PostArea", all.x=TRUE)
TD_AC3.mod1 <- lm('log(rowtot) ~ 0 +(one + z):(label)', TD_ACCOM3)
summary(TD_AC3.mod1)##
## Call:
## lm(formula = "log(rowtot) ~ 0 +(one + z):(label)", data = TD_ACCOM3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.9445 -0.0697 -0.0048 0.0609 1.1264
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## one:labelA 8.597178 0.009877 870.385 < 2e-16 ***
## one:labelB 8.648424 0.005799 1491.481 < 2e-16 ***
## one:labelD 8.633723 0.006078 1420.579 < 2e-16 ***
## one:labelE 8.637250 0.008024 1076.389 < 2e-16 ***
## one:labelF 8.678326 0.007177 1209.160 < 2e-16 ***
## one:labelG 8.702669 0.006565 1325.595 < 2e-16 ***
## one:labelH 8.685583 0.006604 1315.143 < 2e-16 ***
## one:labelJ 8.703050 0.005092 1709.157 < 2e-16 ***
## one:labelK 8.702972 0.006082 1430.839 < 2e-16 ***
## z:labelA -0.011521 0.003483 -3.308 0.000946 ***
## z:labelB -0.003850 0.001950 -1.974 0.048389 *
## z:labelD -0.004037 0.002113 -1.910 0.056160 .
## z:labelE -0.007384 0.002266 -3.258 0.001126 **
## z:labelF -0.002539 0.002292 -1.108 0.267967
## z:labelG -0.008068 0.002135 -3.779 0.000159 ***
## z:labelH -0.004456 0.001297 -3.437 0.000593 ***
## z:labelJ -0.011282 0.001950 -5.787 7.49e-09 ***
## z:labelK -0.013358 0.002420 -5.520 3.52e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1752 on 6862 degrees of freedom
## (1840 observations deleted due to missingness)
## Multiple R-squared: 0.9996, Adjusted R-squared: 0.9996
## F-statistic: 9.363e+05 on 18 and 6862 DF, p-value: < 2.2e-16require(broom)
require(stringr)
mod1.data <- tidy(TD_AC3.mod1)
mod1.data %<>%
filter(str_detect(term, "^z")) #Keep just z score coefs
mod1.data$id <- sub('.*\\.', '', mod1.data$term) # Extract postcode areas for joining to shapefiles
mod1.data$id <-substr(mod1.data$term, 8, 8)
### Visualisation
#Load in shapefile of Postcode town boundaries
require(raster)
GORareas <- shapefile("~/Google_Drive/PhD/Data/DEP/Raw_data/England_gor_2001_clipped/england_gor_2001_clipped.shp")
# Join the data to the shapefile
GORareas@data <- merge(GORareas@data, mod1.data, by.x="label", by.y="id", all.x=TRUE, na.rm=TRUE)
# Write the data for plotting in QGIS
#shapefile(GORareas, '~/Google_Drive/PhD/Data/DEP/GOR.energy.shp')