edits to plots for paper R2; re-run to html

Theme-1/changeOverTime/paper/changingPeakDemandMtus1974_2014_v2.0.Rmd

@@ -4,22 +4,19 @@ author: Ben Anderson (b.anderson@soton.ac.uk, `@dataknut`), Jacopo Torriti (j.to
   `@JTorriti`)
 date: 'Last run: `r Sys.time()`'
 output:
-  html_document:
-    df_print: paged
-    toc: yes
-    toc_depth: '4'
-  html_document2:
+  bookdown::html_document2:
+    code_folding: hide
     fig_caption: yes
-    keep_md: yes
     number_sections: yes
-    self_contained: no
+    self_contained: yes
     toc: yes
-    toc_depth: 4
-    toc_float: yes
-  pdf_document2:
+    toc_depth: 2
+    toc_float: TRUE
+  bookdown::pdf_document2:
+    fig_caption: yes
     number_sections: yes
     toc: yes
-    toc_depth: 4
+    toc_depth: 2
 ---
 
 ```{r knitr setup, include=FALSE}
@@ -77,6 +74,26 @@ makeHmsTime <- function(dt,timeChar ="halfHourChar"){
   return(dt)
 }
 
+addPeakPeriod <- function(plot,caption = myCaption){
+  # adds pae blue box to 16:00 - 20:00
+  newPlot <- plot + 
+  annotate("rect", xmin = hms::parse_hms("16:00:00"),
+             xmax = hms::parse_hms("20:00:00"), 
+             ymin = yMin, ymax = yMax*1.01, alpha = rectAlpha, fill = vLineCol, colour = vLineCol) +
+  scale_x_time(breaks = c(hms::parse_hms("03:00:00"),
+                          hms::parse_hms("06:00:00"),
+                          hms::parse_hms("09:00:00"),
+                           hms::parse_hms("12:00:00"),
+                           hms::parse_hms("15:00:00"),
+                          hms::parse_hms("18:00:00"),
+                          hms::parse_hms("21:00:00")
+                           ) #, date_labels = "%H:%M"
+                ) + 
+  labs(caption = paste0(caption, "\nPeak demand period shown shaded")) + # add to whatever comes in as a caption
+  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 0.5))
+  return(newPlot)
+}
+
 ba_createStackedPlot <- function(dt,fillVar,yVar,yLabel,facetVars){
   # intended for testing overall distributions
   ggplot(dt, aes(y = get(yVar), x = r_hmsTime, fill = get(fillVar))) + # forcats::fct_rev() to reverse categories
@@ -109,6 +126,22 @@ ba_createStepPlot <- function(dt,fillVar,yVar,yLabel,facetVars){
     ) 
 }
 
+ba_createLinePlot <- function(dt,fillVar,yVar,yLabel,facetVars){
+  # intended for testing overall distributions
+  ggplot(dt, aes(y = get(yVar), x = r_hmsTime, linetype = get(fillVar), colour = get(fillVar))) + # forcats::fct_rev() to reverse categories
+    geom_line() +
+    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 0.5)) + 
+    theme(strip.text.y = element_text(angle = 0, vjust = 0.5, hjust = 0.5)) + 
+    facet_grid(eval(facetVars)) +
+    #scale_fill_grey() +
+    theme(legend.title = element_blank()) +
+    theme(legend.position = "right") +
+    labs(caption = myCaption,
+         x = "Hour of the day",
+         y = eval(yLabel)
+    ) 
+}
+
 ba_createGreyStackedPlot <- function(dt,fillVar,yVar,yLabel,facetVars){
   # intended for testing overall distributions
   ggplot(dt, aes(y = get(yVar), x = r_hmsTime, fill = get(fillVar))) + # forcats::fct_rev() to reverse categories
@@ -405,26 +438,9 @@ halfHourlyNgPlotDT <- extractDT[, .( meanGW = mean(ENGLAND_WALES_DEMAND)/1000),
 The next two charts show demand levels in January 2006 & 2016 as mean MW per hour. Figure \ref(fig:compareTotal) shows mean MW across weekdays compared with each weekend day. Fig 1 in paper.
 
 ```{r compareTotal, fig.cap="NG demand totals over time", fig.height=6}
-addPeakPeriod <- function(plot){
-  newPlot <- plot + 
-  annotate("rect", xmin = hms::parse_hms("16:00:00"),
-             xmax = hms::parse_hms("20:00:00"), 
-             ymin = yMin, ymax = yMax*1.01, alpha = rectAlpha, fill = vLineCol, colour = vLineCol) +
-  scale_x_time(breaks = c(hms::parse_hms("03:00:00"),
-                          hms::parse_hms("06:00:00"),
-                          hms::parse_hms("09:00:00"),
-                           hms::parse_hms("12:00:00"),
-                           hms::parse_hms("15:00:00"),
-                          hms::parse_hms("18:00:00"),
-                          hms::parse_hms("21:00:00")
-                           ) #, date_labels = "%H:%M"
-                ) + 
-  labs(caption = paste0(dataNote, "\nPeak demand period shown shaded")) + 
-  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 0.5))
-  return(newPlot)
-}
 
-dataNote <- "Source: National Grid half-hourly demand data (England & Wales) 2006-2016"
+
+myCaption <- "Source: National Grid half-hourly demand data (England & Wales) 2006-2016"
 
 rectAlpha <- 0.1
 vLineAlpha <- 0.4
@@ -435,7 +451,7 @@ rawPlot <- ggplot(halfHourlyNgPlotDT, aes(x = r_hmsTime, y = meanGW, colour = as
   geom_line() +
   facet_grid(r_month ~ wDay) +
   #ylim(0,NA) +
-    labs(caption = dataNote,
+    labs(caption = myCaption,
        x = "Time of Day",
        y = "Mean GW") +
   theme(legend.title=element_blank()) #http://www.cookbook-r.com/Graphs/Legends_%28ggplot2%29/#hiding-the-legend-title
@@ -443,12 +459,11 @@ rawPlot <- ggplot(halfHourlyNgPlotDT, aes(x = r_hmsTime, y = meanGW, colour = as
 rawPlot
 yMax <- max(halfHourlyNgPlotDT$meanGW)
 yMin <- min(halfHourlyNgPlotDT$meanGW)
-myPlot <- addPeakPeriod(rawPlot) + scale_colour_grey(start = 0.8, end = 0.2)
+myPlot <- addPeakPeriod(rawPlot, caption = myCaption) + scale_colour_grey(start = 0.8, end = 0.2)
 myPlot + 
-  labs(caption = paste0(dataNote, "\nPeak demand period shown shaded")) + 
   theme(legend.position = "bottom")
 
-ggsave("ngMeanGW.png")
+ggsave("ngMeanGW_R2Fig1.png", height = 6)
 ```
 
 Figure \ref(fig:compareNormalised) shows the same data but normalised to the overall mean for the year, month and day of the week. This shows the relative distribution of demand rather than absolute and illustrates the increased 'peakiness' as well as shift to later demand. Fig 2 in paper.
@@ -479,12 +494,11 @@ vLineAlpha <- 0.4
 vLineCol <- "#0072B2" # http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
 myTextSize <- 3
 
-myPlot <- addPeakPeriod(normPlot) + scale_colour_grey(start = 0.8, end = 0.2)
+myPlot <- addPeakPeriod(normPlot, caption = myCaption) + scale_colour_grey(start = 0.8, end = 0.2)
 myPlot +
-  labs(caption = paste0(dataNote, "\nPeak demand period shown shaded")) + 
   theme(legend.position = "bottom")
 
-ggsave("ngNormalisedMeanGW.png")
+ggsave("ngNormalisedMeanGW_R2Fig2.png", height = 6)
 ```
 
 
@@ -1059,10 +1073,10 @@ myTextSize <- 3
 # ba_createStackedPlot(plotDT, yVar = "pcAnyHHMainWt", fillVar = "actLabel", yLabel = "% halfhours (main act)", facetVars = "ba_survey ~ inHome") +
 #   scale_fill_hue()
 
-ba_createStepPlot(plotDT, yVar = "pcAnyHHMainWt", fillVar = "actLabel", yLabel = "% halfhours (main act)", facetVars = "ba_survey ~ inHome") +
+ba_createLinePlot(plotDT, yVar = "pcAnyHHMainWt", fillVar = "actLabel", yLabel = "% halfhours (main act)", facetVars = "ba_survey ~ inHome") +
   scale_fill_hue()
 
-ggsave("allYearsLinePlot.png", height = 8)
+ggsave("allYearsLinePlot_R2Fig3.png", height = 8)
 ```
 
 Separate out weekdays from weekend days and select specific variables to plot...
@@ -1123,6 +1137,9 @@ changePlotDT <- dtChange1974_2014WtDT[actLabel %like% "Personal" |
                                       actLabel %like% "Media", .(r_hmsTime, inHome, r_wday, actLabel, actLabelChar, pcAnyHHMainWt.1974, pcAnyHHMainWt.2014, deltaMain)
                                 ]
 # with(changePlotDT, table(inHome, actLabelChar))
+
+# recode actLabelChar to get plot order more sensibe (there are other ways...)
+changePlotDT <- changePlotDT[actLabelChar %like% "Personal", actLabelChar := "Home/personal care"]
 # fiddle around so we can avoid facet by inHome
 changePlotDT <- changePlotDT[, actLabelChar := ifelse(actLabel %like% "Work" & 
                                                         inHome %like% "Not", 
@@ -1144,8 +1161,6 @@ ba_plotChangeDow <- function(dt,actVar,yLabel){ # change plot function
                        linetype = get(actVar))) +
     geom_line() + 
     theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 0.5)) + 
-    #scale_colour_grey() +
-    facet_grid(r_wday ~ .) +
     theme(legend.title = element_blank()) +
     theme(legend.position = "bottom") +
     labs(caption = myCaption,
@@ -1157,15 +1172,23 @@ ba_plotChangeDow <- function(dt,actVar,yLabel){ # change plot function
 ```
 
 
-Fig \ref(fig:selectedActsDeltaPlot) drawa % point change plot 1974 - 2014 for selected variables. (paper v2.0 Fig 4)
+Fig \ref(fig:selectedActsDeltaPlot) shows % point change plot 1974 - 2014 for selected variables. (paper v2.0 Fig 4)
 
 > NB: ignores hours < 06:00 to avoid noisy data
 
-All activities at except Travel or where labelled.
+All activities at home except Travel or where labelled.
 
 ```{r selectedActsDeltaPlot, fig.height=6}
-myPlot <- ba_plotChangeDow(changePlotDT[plotMe == 1 & lubridate::hour(r_hmsTime) >= 6], # <- select after 05:59
-                           actVar = "actLabelChar", 
+# set label to create two paels and keep colours constant for this and subsequent plots
+changePlotDT <- changePlotDT[, plotPanel := "Panel B (Household)"]
+changePlotDT <- changePlotDT[ actLabelChar == "Travel" |
+                              actLabelChar %like% "Work"
+                              , plotPanel := "Panel A (Work and travel)"]
+
+
+myCaption <- "Synthetic MTUS 1974 - 2014, weighted"
+myPlot <- ba_plotChangeDow(changePlotDT[plotMe == 1 & plotPanel %like% "A" & lubridate::hour(r_hmsTime) >= 6], # <- select after 05:59
+                           actVar = "actLabelChar", # need to use this as it has been adjusted for work (not) at home
                         yLabel = "% point change 1974 - 2014") # plot as normal ggplot (easier for paper writing)
 yMax <- max(changePlotDT[plotMe == 1 & lubridate::hour(r_hmsTime) >= 6]$deltaMain)
 yMin <- min(changePlotDT[plotMe == 1 & lubridate::hour(r_hmsTime) >= 6]$deltaMain)
@@ -1174,9 +1197,25 @@ vLineAlpha <- 0.4
 vLineCol <- "#0072B2" # http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
 myTextSize <- 3
 
-addPeakPeriod(myPlot) + geom_hline(yintercept = 0)
+panelAColours <- c("orange","grey","black")
 
-ggsave("selectedActsDeltaPlot.png", height = 8)
+addPeakPeriod(myPlot, caption = myCaption) + geom_hline(yintercept = 0) + 
+  facet_grid(r_wday ~ plotPanel) +
+  scale_colour_manual(values = panelAColours)
+
+ggsave("selectedActsDeltaPlot_R2Fig4a.png", width = 4)
+
+myPlot <- ba_plotChangeDow(changePlotDT[plotMe == 1 & plotPanel %like% "B" & lubridate::hour(r_hmsTime) >= 6], # <- select after 05:59
+                           actVar = "actLabel", 
+                        yLabel = "% point change 1974 - 2014") # plot as normal ggplot (easier for paper writing)
+
+finalColours <- c("red","blue","darkgreen")
+
+addPeakPeriod(myPlot, caption = myCaption) + geom_hline(yintercept = 0) + 
+  facet_grid(r_wday ~ plotPanel) +
+  scale_colour_manual(values = finalColours)
+
+ggsave("selectedActsDeltaPlot_R2Fig4b.png", width = 4)
 ```
 
 This charts shows:
@@ -1323,7 +1362,6 @@ changePlotDT <- changePlotDT[, plotMe := ifelse(inHome %like% "At" |
 16-64 - All activities at home except Travel or where labelled (Figure \ref(fig:deltaPlotDow16_64) - Fig 5 in paper).
 
 ```{r deltaPlotDow16_64, fig.height=8, fig.cap="deltaPlotDow16_64"}
-
 plotDT <- changePlotDT[plotMe == 1 & ba_workingAge == "16-64" &
                                           lubridate::hour(r_hmsTime) >= 6] # <- select after 05:59
 myPlot <- ba_plotChangeDow(plotDT, 
@@ -1337,10 +1375,9 @@ vLineCol <- "#0072B2" # http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-col
 myTextSize <- 3
 addPeakPeriod(myPlot) + geom_hline(yintercept = 0) + facet_grid(r_wday ~ ba_sex)
 
-ggsave("deltaPlotDow16_64.png", height = 8)
 ```
 
-65+ - All activities at home except Travel or where labelled (Figure \ref(fig:deltaPlotDow64m) - Fig 6 in paper).
+65+ - All activities at home except Travel or where labelled (Figure \ref(fig:deltaPlotDow64m)).
 
 ```{r deltaPlotDow64m, fig.height=8, fig.cap="deltaPlotDow64m"}
 
@@ -1357,7 +1394,6 @@ vLineCol <- "#0072B2" # http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-col
 myTextSize <- 3
 addPeakPeriod(myPlot) + geom_hline(yintercept = 0) + facet_grid(r_wday ~ ba_sex)
 
-ggsave("deltaPlotDow64m.png", height = 8)
 ```
 
 ## Detailed Activities: Descriptive analysis
@@ -1372,9 +1408,9 @@ Test this using two poisson models - one for 'early' eating and one for 'late' e
 
 Also use two logit to test p(any reported) - this tests predictors of reporting any activity in the time period and may be easier to interpret.
 
-Create a change plot - Fig 7 in paper.
+Create a change plot - Figure 5 in paper.
 
-```{r plotChange}
+```{r deltaPlotWeekDaysAgeByGender}
 deltaPlotDT <- dtChange1974_2014_workingAgeWtDT[r_hmsTime >= hms::as.hms("07:00:00") & 
                                                   (actLabel == "Food" | actLabel == "Media" | 
                                       actLabel %like% "Personal") &
@@ -1395,9 +1431,10 @@ myPlot <- ggplot2::ggplot(deltaPlotDT, aes(x = r_hmsTime)) +
   facet_grid(ba_sex ~ ba_workingAge)
 
 dataNote <- "MTUS 1974 - 2014"
-addPeakPeriod(myPlot)
+addPeakPeriod(myPlot) +
+  scale_colour_manual(values = finalColours)
 
-ggsave("deltaPlotWeekDaysByGender.png", height = 8)
+ggsave("deltaPlotWeekDaysAgeByGender_R2Fig5.png", height = 6)
 ```
 
 NB: we use interaction terms for survey year & working age to understand how older people's reported activities have changed over time.
@@ -2048,9 +2085,10 @@ myPlot <- ggplot2::ggplot(compareJunDT, aes(x = meanNormDiffGW, y = deltaMain))
   labs(x = "Mean change in normalised GW",
        y = "% point change in recorded half hours")
 addPlotStuff(myPlot) + geom_smooth(colour = "black" , se=TRUE)
+ggsave("plotRelativeGWDeltaTU_R2Fig6.png", height = 6)
 ```
 
-Repeat for relative GW (Figure \@ref(fig:plotRelativeGWDeltaTUNoFit)) without best fit lines (for paper).
+Repeat for relative GW (Figure \@ref(fig:plotRelativeGWDeltaTUNoFit)) without best fit lines.
 
 ```{r plotRelativeGWDeltaTUNoFit, fig.cap="Comparing relative GW change"}
 myPlot <- ggplot2::ggplot(compareJunDT, aes(x = meanNormDiffGW, y = deltaMain)) +