split drought model so baseline WE done seperately (first) then drought WE & TUB

impetusModel/addBaselineWE.Rmd

+# Introduction
+
+Now we add the basline water efficiency uptake to each of the two models
+
+## Add baseline WE uptake
+
+This has to be done after the expansion to months as it is a monthly uptake model. how this works:
+
+ * Add baseline WE uptake rates of dual flush WC & low flow shower head uptake and adjust l/day/hh for WCs and showers accordingly
+
+### Backwards estimation of water efficiency uptake
+The first step is to allocate households to dual-flush/no dual-flush status and low-flow/no low-flow shower heads for all years. The logic applied is as follows:
+
+  * dual-flush/no dual-flush WC:
+    * linearly interpolate rate for a given year and month using the 2% yearly adoption rate
+    * apply these rates randomly to households. We acknowledge that those with meters may be more likely to also have low-flush WCs but we do not model this association at this time.
+  * low-flow/normal flow shower head:
+    * as for WC but stratified by top & bottom 50% shower litres/day within household size
+
+> somehow this process needs to take account of the 'existing' uptake levels in 2011 in model v2 (based on SPRG survey data)
+
+```{r Backcast WE uptake rates}
+# start with a collapsed year/month spine
+
+hhFinalDataComboExpandedDT <- hhFinalDataComboExpandedDT[, obsDate := ymd(paste0(currYear, "-", currMon, "-15"))]  # make it the middle
+
+yearMonthSpineDT <- hhFinalDataComboExpandedDT[, .(nObs = .N), by = .(currYear, currMon, obsDate)]
+
+
+minYear <- min(yearMonthSpineDT$currYear)
+maxYear <- max(yearMonthSpineDT$currYear)
+
+# reduce occupancy levels to 1->5, 6+ for this purpose
+hhFinalDataComboExpandedDT <- hhFinalDataComboExpandedDT[, occRed := ifelse(occupancy >= 6,
+                                                                                  6,
+                                                                                  occupancy)
+                                                               ]
+
+hhFinalDataComboExpandedDT$loWC <- 0
+hhFinalDataComboExpandedDT$loSh <- 0
+
+minOcc <- min(hhFinalDataComboExpandedDT$occRed)
+maxOcc <- max(hhFinalDataComboExpandedDT$occRed)
+
+wcCollectorDT <- NULL
+shCollectorDT <- NULL
+
+
+models <- c("v1_3","v2_0") # by model (as usage distributions vary)
+for(mod in models){
+  countYear = 0 # reset for the start of each model
+  for(y in maxYear:minYear){
+    print(paste0("Updating water efficiency rates for model: ",mod," in ",y))
+    for(m in 12:1){ # run backwards so we can work backwards from the 2011 rates
+      yearMonthSpineDT <- yearMonthSpineDT[currYear == y & currMon == m, 
+                                           dfWCRate := estdfRate - (countYear * (dfWCAdopt/12))]
+      yearMonthSpineDT <- yearMonthSpineDT[currYear == y & currMon == m, 
+                                           lfShowerRate := estlfRate - (countYear * (lfShowerAdopt/12))]
+      countYear <- countYear + 1
+      
+      monthlyWCRate <- yearMonthSpineDT[currYear == y & currMon == m, dfWCRate]
+      monthlyShRate <- yearMonthSpineDT[currYear == y & currMon == m, lfShowerRate]
+      
+      tempDT <- hhFinalDataComboExpandedDT[currYear == y & currMon == m & model == mod,]
+      if(nrow(tempDT) > 0){ # check we have household data for this combination of month, year & model
+        #print(paste0("Model: ", mod , " Year: ", y , " Month: ", m , " N rows: ", nrow(tempDT)))
+        # for WCs
+        #print(paste0("monthlyWCRate: ", monthlyWCRate))
+        # sample households within model/year/month according to the monthly rate
+        # use sample_n as it randomly selects n, use floor to force rounding so n = an integer
+        dualFlushDT <- as.data.table(sample_n(tempDT[, 
+                                                     .(currYear,currMon,obsDate, hhid, model)], 
+                                              floor(nrow(tempDT)*monthlyWCRate)
+        )
+        )
+        #print(paste0("WC: Selected ", nrow(dualFlushDT), " using monthlyWCRate = ", round(monthlyWCRate,3)))
+        dualFlushDT <- dualFlushDT[, dualFlushWC := "Dual flush"]
+        
+        wcCollectorDT <- rbind(wcCollectorDT, dualFlushDT) # add to the data collector
+        
+        # same for showers
+        showerCut <- tempDT[, median(Shower.baseline.madj, na.rm = TRUE)]
+        # sample households within model/year/month according to the monthly rate
+        # use sample_n as it randomly selects n, use floor to force rounding so n = an integer
+        lowFlowShDT <- as.data.table(sample_n(tempDT[, 
+                                                     .(currYear,currMon,obsDate, hhid, model)], 
+                                              floor(nrow(tempDT)*monthlyShRate)
+        )
+        )
+        #print(paste0("Shower: Selected ", nrow(lowFlowShDT), " using monthlyShRate = ", round(monthlyShRate,3)))
+        lowFlowShDT <- lowFlowShDT[, loFlowShower := "Low flow"]
+        shCollectorDT <- rbind(shCollectorDT, lowFlowShDT) # add to the data collector
+      }
+    }
+  }
+}
+
+wcCollectorDT <- as.data.table(wcCollectorDT)
+setkey(wcCollectorDT, currYear,currMon,hhid,model,obsDate)
+setkey(hhFinalDataComboExpandedDT, currYear,currMon,hhid,model,obsDate)
+# merge the dualFlushWC indicator back on to main data file
+hhFinalDataComboExpandedDT <- merge(hhFinalDataComboExpandedDT, wcCollectorDT, all.x = TRUE) # keep all records
+# set NA to 0 (not low flush/flow)
+hhFinalDataComboExpandedDT <- hhFinalDataComboExpandedDT[,dualFlushWC := ifelse(is.na(dualFlushWC), "Single flush", dualFlushWC)]
+          
+# repeat for showers
+shCollectorDT <- as.data.table(shCollectorDT)
+setkey(shCollectorDT, currYear,currMon,hhid,model,obsDate)
+setkey(hhFinalDataComboExpandedDT, currYear,currMon,hhid,model,obsDate)
+# merge the low flow shower indicator back on to main data file
+hhFinalDataComboExpandedDT <- merge(hhFinalDataComboExpandedDT, shCollectorDT, all.x = TRUE) # keep all records
+hhFinalDataComboExpandedDT <- hhFinalDataComboExpandedDT[, 
+                                                                 loFlowShower := ifelse(is.na(loFlowShower), "Normal flow", loFlowShower)]
+
+#summary(yearMonthSpineDT)
+```
+
+```{r check uptake rates}
+# check imputed update rates
+ggplot(yearMonthSpineDT, aes(x=obsDate)) +
+  geom_line(aes(y = 100*dfWCRate, colour = "Dual flush WC")) +
+  geom_line(aes(y = 100*lfShowerRate, colour = "Low flow shower")) +
+  theme(legend.title = element_blank()) +
+  labs(title = "Modelled uptake rates",
+       y = "%",
+       x = "Date") 
+
+# check imputed uptake
+dt <- hhFinalDataComboExpandedDT[, .(nHHs = .N,
+                                         pcHH = 100*(.N/1800)
+                                         ),
+                                     by=.(obsDate, model, dualFlushWC)]
+ggplot(dt, aes(x=obsDate, colour = factor(dualFlushWC))) +
+  geom_line(aes(y = pcHH)) +
+  theme(legend.title = element_blank()) +
+  facet_grid(model ~ .) +
+  labs(title = "Dual Flush WCs: Modelled household uptake",
+       y = "%",
+       x = "Date") 
+
+# check imputed uptake
+dt <- hhFinalDataComboExpandedDT[, .(nHHs = .N,
+                                         pcHH = 100*(.N/1800)
+                                         ), 
+                                     by=.(obsDate, model, loFlowShower)]
+
+ggplot(dt, aes(x=obsDate, colour = factor(loFlowShower))) +
+  geom_line(aes(y = pcHH)) +
+  theme(legend.title = element_blank()) +
+  facet_grid(model ~ .) +
+  labs(title = "Low flow showers: Modelled household uptake",
+       y = "%",
+       x = "Date") 
+```
+
+Next we use the water use reduction values given at the start of the section to update the l/day/hh for those who have the dual flush WCs or low flow showers.
+
+```{r update baseline water consumption}
+# wc ----
+hhFinalDataComboExpandedDT <- hhFinalDataComboExpandedDT[, WC.baseline.madj.we := ifelse(dualFlushWC == "Dual flush", WC.baseline.madj * dfWCReduction, WC.baseline.madj)
+                                                               ]
+
+# shower ----
+hhFinalDataComboExpandedDT <- hhFinalDataComboExpandedDT[, Shower.baseline.madj.we := ifelse(loFlowShower == "Low flow", Shower.baseline.madj * lfShowerReduction, Shower.baseline.madj)
+                                                               ]
+
+# update total ----
+hhFinalDataComboExpandedDT <- hhFinalDataComboExpandedDT[, sumDaily.baseline.madj.we := Basin.baseline.madj + 
+                                                           Bath.baseline.madj + 
+                                                           Dishwasher.baseline.madj + 
+                                                           External.baseline.madj + 
+                                                           KitchenSink.baseline.madj + 
+                                                           Shower.baseline.madj.we + 
+                                                           WC.baseline.madj.we + 
+                                                           WashingMachine.baseline.madj
+                                                         ]
+```
+
+
+### Compare effects
+Check the effects on each model
+
+```{r plot we adjusted shower and WC use}
+plotDT <- hhFinalDataComboExpandedDT[, .(meanSh = mean(Shower.baseline.madj, na.rm = TRUE),
+                                            meanShWE = mean(Shower.baseline.madj.we, na.rm = TRUE),
+                                            meanWC = mean(WC.baseline.madj, na.rm = TRUE),
+                                            meanWCWE = mean(WC.baseline.madj.we, na.rm = TRUE)
+                                            ),
+                                        by = .(obsDate, model, metered)]
+
+ggplot(plotDT, aes(x = obsDate,)) +
+  geom_line(aes(y = meanSh, colour = "Shower")) + 
+  geom_line(aes(y = meanShWE, colour = "Shower with modeled low flow shower uptake")) +
+  facet_grid(model ~ metered) +
+  theme(legend.title = element_blank()) +
+  theme(legend.position = "bottom") +
+  labs(title = "Low flow showers: impact of modeled uptake",
+       y = "Mean l/hh/day",
+       x = "Date")
+
+ggplot(plotDT, aes(x = obsDate,)) +
+  geom_line(aes(y = meanWC, colour = "WC")) + 
+  geom_line(aes(y = meanWCWE, colour = "WC with modeled dual flush uptake")) +
+  facet_grid(model ~ metered) +
+  theme(legend.title = element_blank()) +
+  theme(legend.position = "bottom") +
+  labs(title = "Dual flush WCs: impact of modeled uptake",
+       y = "Mean l/hh/day",
+       x = "Date")
+```
+
+Well that seems to have an effect!
+
+## Model 1 (synthetic) results
+
+Now re-draw final baseline charts for papers
+
+```{r final model v1_3 2012 by month}
+myTitle <- "All uses (2012 only)"
+myCaption <- "IMPETUS model: synthetic households (n = 1800 per month)\nModel v1\nBaseline water efficiency uptake"
+
+dt <- hhFinalDataComboExpandedDT[currYear == 2012 & model == "v1_3", .(Basin = mean(Basin.baseline.madj),
+                                 Bath = mean(Bath.baseline.madj),
+                                 KitchenSink = mean(KitchenSink.baseline.madj),
+                                 Dishwasher = mean(Dishwasher.baseline.madj),
+                                 External = mean(External.baseline.madj),
+                                 Shower = mean(Shower.baseline.madj.we),
+                                 WC = mean(WC.baseline.madj.we),
+                                 WashingMachine = mean(WashingMachine.baseline.madj)
+                                 ), by = .(metered, currMon)]
+# recast dt to make plotting easier
+plotDT <- melt(dt, id.vars = c("currMon", "metered"))
+plotDT <- plotDT[, Usage := variable]
+
+myPlot <- ggplot(plotDT, aes(x = factor(currMon), y = value, group = Usage)
+                 ) +
+  geom_line(aes(colour = Usage)) +
+  facet_grid(. ~ metered, scales = "free") +
+  theme(legend.title = element_blank()) +
+  theme(legend.position = "bottom") +
+  labs(y = "Mean litres/day",
+       x = "Month",
+       title = myTitle,
+       caption = myCaption)
+
+myPlot
+
+# version with linetype for use in bw fig
+myPlot <- ggplot(plotDT, aes(x = factor(currMon), y = value, group = Usage)
+                 ) +
+  geom_line(aes(linetype = Usage)) +
+  facet_grid(. ~ metered, scales = "free") +
+  theme(legend.title = element_blank()) +
+  theme(legend.position = "bottom") +
+  labs(y = "Mean litres/day",
+       x = "Month",
+       caption = myCaption)
+
+myPlot
+
+# Grey scale version if required
+myPlot <- myPlot + theme_bw()
+
+# Figure for IWA Bath final paper (http://ws.iwaponline.com/content/early/2018/02/13/ws.2018.035)
+ggsave(paste0("plots_v1/Fig3_Final_model_v1_3_2012_by_month.pdf"), plot = myPlot, dpi = 400)
+
+```
+
+```{r test model v1_3 by component}
+myTitle <- "External: mean daily total (no WE adjustment)"
+myCaption <- "Model: v1_3 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v1_3"], "External.baseline.madj")
+myPlot
+
+myTitle <- "Shower: mean daily total (WE adjustment)"
+myCaption <- "Model: v1_3 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v1_3"], "Shower.baseline.madj.we")
+myPlot
+
+myTitle <- "WC: mean daily total (WE adjustment)"
+myCaption <- "Model: v1_3 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v1_3"], "WC.baseline.madj.we")
+myPlot
+
+myTitle <- "Dishwasher: mean daily total (no WE adjustment)"
+myCaption <- "Model: v1_3 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v1_3"], "Dishwasher.baseline.madj")
+myPlot
+
+myTitle <- "Bath: mean daily total (no WE adjustment)"
+myCaption <- "Model: v1_3 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v1_3"], "Bath.baseline.madj")
+myPlot
+
+myTitle <- "Washing machine: mean daily total (no WE adjustment)"
+myCaption <- "Model: v1_3 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v1_3"], "WashingMachine.baseline.madj")
+myPlot
+
+myTitle <- "Kitchen sink: mean daily total (no WE adjustment)"
+myCaption <- "Model: v1_3 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v1_3"], "KitchenSink.baseline.madj")
+myPlot
+
+```
+
+```{r final model v1_3 all dates}
+myTitle <- "All uses (all years)"
+myCaption <- "IMPETUS model: synthetic households (n = 1800 per month)\nModel v1\nBaseline water efficiency uptake"
+
+dt <- hhFinalDataComboExpandedDT[model == "v1_3", .(Basin = mean(Basin.baseline.madj),
+                                 Bath = mean(Bath.baseline.madj),
+                                 KitchenSink = mean(KitchenSink.baseline.madj),
+                                 Dishwasher = mean(Dishwasher.baseline.madj),
+                                 External = mean(External.baseline.madj),
+                                 Shower = mean(Shower.baseline.madj.we),
+                                 WC = mean(WC.baseline.madj.we),
+                                 WashingMachine = mean(WashingMachine.baseline.madj)
+                                 ), by = .(metered, obsDate)]
+
+myPlot <- ggplot(dt, aes(x = obsDate, group = metered)
+                 ) +
+  #geom_boxplot()  +
+  geom_line(aes(y = Basin, colour = "Basin")) +
+  geom_line(aes(y = Bath, colour = "Bath")) +
+  geom_line(aes(y = KitchenSink, colour = "Kitchen Sink")) +
+  geom_line(aes(y = Dishwasher, colour = "Dishwasher")) +
+  geom_line(aes(y = External, colour = "External")) +
+  geom_line(aes(y = Shower, colour = "Shower")) +
+  geom_line(aes(y = WC, colour = "WC")) +
+  geom_line(aes(y = WashingMachine, colour = "Washing Machine")) +
+  facet_grid(. ~ metered, scales = "free") +
+  theme(legend.title = element_blank()) +
+  theme(legend.position = "bottom") +
+  labs(y = "Mean litres/day",
+       x = "Month",
+       title = myTitle,
+       caption = myCaption)
+
+myPlot
+```
+
+```{r test model v1_3 total}
+myTitle <- "Total: mean daily total"
+myCaption <- "Baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v1_3"], "sumDaily.baseline.madj.we")
+myPlot
+```
+
+```{r model v1_3 with ci}
+# all uses mean daily
+dt <- hhFinalDataComboExpandedDT[model == "v1_3"]
+
+# create a monthly mean & 95% CI plot using the dt & var
+  plotDT <- dt[, .(meanDaily = mean(sumDaily.baseline.madj, na.rm = TRUE),
+                   sd = sd(sumDaily.baseline.madj, na.rm = TRUE),
+                   nObs = .N), by = .(metered, obsDate, currYear)]
+  plotDT <- plotDT[, ciLower := meanDaily - qnorm(0.975) * (sd/sqrt(nObs))]
+  plotDT <- plotDT[, ciUpper := meanDaily + qnorm(0.975) * (sd/sqrt(nObs))]
+  library(lubridate)
+  plotDT[, month := lubridate::month(obsDate, label = TRUE, abbr = TRUE)]
+  
+  ggplot(plotDT, aes(x = month, y = meanDaily, 
+                     linetype = metered,
+                     colour = as.factor(currYear),
+                     group = as.factor(currYear))) +
+    geom_line() +
+    geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), width = 0.2) +
+    theme(legend.title = element_blank()) +
+    facet_grid(. ~ metered) +
+    theme(legend.position = "bottom") +
+    labs(y = "Mean litres/day",
+         x = "Month")
+  
+  ggsave(paste0("plots_v1/baselineMadjWeMonthlyByYearMetering.pdf"), dpi = 480)
+```
+
+
+## Model 2 (SPRG practices) results
+
+Now redraw final baseline charts for papers
+
+```{r final model v2_0 2012 by month}
+myTitle <- "All uses (2012 only)"
+myCaption <- "IMPETUS model: SPRG households (n = 1800)\nModel v2\nBaseline water efficiency uptake"
+
+# have to include na.rm
+dt <- hhFinalDataComboExpandedDT[currYear == 2012 & model == "v2_0", .(Basin = mean(Basin.baseline.madj),
+                                 Bath = mean(Bath.baseline.madj, na.rm = TRUE),
+                                 KitchenSink = mean(KitchenSink.baseline.madj, na.rm = TRUE),
+                                 Dishwasher = mean(Dishwasher.baseline.madj, na.rm = TRUE),
+                                 External = mean(External.baseline.madj, na.rm = TRUE),
+                                 Shower = mean(Shower.baseline.madj.we, na.rm = TRUE),
+                                 WC = mean(WC.baseline.madj.we, na.rm = TRUE),
+                                 WashingMachine = mean(WashingMachine.baseline.madj, na.rm = TRUE)
+                                 ), by = .(metered, currMon)]
+# recast dt to make plotting easier
+plotDT <- melt(dt, id.vars = c("currMon", "metered"))
+plotDT <- plotDT[, Usage := variable]
+
+myPlot <- ggplot(plotDT, aes(x = factor(currMon), y = value, group = Usage)
+                 ) +
+  geom_line(aes(colour = Usage)) +
+  facet_grid(. ~ metered, scales = "free") +
+  theme(legend.title = element_blank()) +
+  theme(legend.position = "bottom") +
+  labs(y = "Mean litres/day",
+       x = "Month",
+       title = myTitle,
+       caption = myCaption)
+
+myPlot
+
+# version with linetype for use in bw fig
+myPlot <- ggplot(plotDT, aes(x = factor(currMon), y = value, group = Usage)
+                 ) +
+  geom_line(aes(linetype = Usage)) +
+  facet_grid(. ~ metered, scales = "free") +
+  theme(legend.title = element_blank()) +
+  theme(legend.position = "bottom") +
+  labs(y = "Mean litres/day",
+       x = "Month",
+       caption = myCaption)
+
+myPlot
+
+# Grey scale version if required
+myPlot <- myPlot + theme_bw()
+
+# Figure for IWA Bath final paper (http://ws.iwaponline.com/content/early/2018/02/13/ws.2018.035)
+ggsave(paste0("plots_v2/Final_model_v2_0_2012_by_month.pdf"), plot = myPlot, dpi = 400)
+
+```
+
+```{r test model v2 by component}
+myTitle <- "External: mean daily total (no WE adjustment)"
+myCaption <- "Model: v2_0 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v2_0"], "External.baseline.madj")
+myPlot
+
+myTitle <- "Shower: mean daily total (WE adjustment)"
+myCaption <- "Model: v2_0 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v2_0"], "Shower.baseline.madj.we")
+myPlot
+
+myTitle <- "WC: mean daily total (WE adjustment)"
+myCaption <- "Model: v2_0 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v2_0"], "WC.baseline.madj.we")
+myPlot
+
+myTitle <- "Dishwasher: mean daily total (no WE adjustment)"
+myCaption <- "Model: v2_0 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v2_0"], "Dishwasher.baseline.madj")
+myPlot
+
+myTitle <- "Bath: mean daily total (no WE adjustment)"
+myCaption <- "Model: v2_0 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v2_0"], "Bath.baseline.madj")
+myPlot
+
+myTitle <- "Washing machine: mean daily total (no WE adjustment)"
+myCaption <- "Model: v2_0 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v2_0"], "WashingMachine.baseline.madj")
+myPlot
+
+myTitle <- "Kitchen sink: mean daily total (no WE adjustment)"
+myCaption <- "Model: v2_0 baseline with water efficiency"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v2_0"], "KitchenSink.baseline.madj")
+myPlot
+```
+
+```{r final model v2_0 all dates}
+myTitle <- "All uses (all years)"
+myCaption <- "IMPETUS model: SPRG households (n = 1800)\nModel v2\nBaseline water efficiency uptake"
+
+dt <- hhFinalDataComboExpandedDT[model == "v2_0", .(Basin = mean(Basin.baseline.madj, na.rm = TRUE),
+                                 Bath = mean(Bath.baseline.madj, na.rm = TRUE),
+                                 KitchenSink = mean(KitchenSink.baseline.madj, na.rm = TRUE),
+                                 Dishwasher = mean(Dishwasher.baseline.madj, na.rm = TRUE),
+                                 External = mean(External.baseline.madj, na.rm = TRUE),
+                                 Shower = mean(Shower.baseline.madj.we, na.rm = TRUE),
+                                 WC = mean(WC.baseline.madj.we, na.rm = TRUE),
+                                 WashingMachine = mean(WashingMachine.baseline.madj, na.rm = TRUE)
+                                 ), by = .(metered, obsDate)]
+
+myPlot <- ggplot(dt, aes(x = obsDate, group = metered)
+                 ) +
+  #geom_boxplot()  +
+  geom_line(aes(y = Basin, colour = "Basin")) +
+  geom_line(aes(y = Bath, colour = "Bath")) +
+  geom_line(aes(y = KitchenSink, colour = "Kitchen Sink")) +
+  geom_line(aes(y = Dishwasher, colour = "Dishwasher")) +
+  geom_line(aes(y = External, colour = "External")) +
+  geom_line(aes(y = Shower, colour = "Shower")) +
+  geom_line(aes(y = WC, colour = "WC")) +
+  geom_line(aes(y = WashingMachine, colour = "Washing Machine")) +
+  facet_grid(. ~ metered, scales = "free") +
+  theme(legend.title = element_blank()) +
+  theme(legend.position = "bottom") +
+  labs(y = "Mean litres/day",
+       x = "Month")
+
+myPlot
+ggsave(paste0("plots_v2/baselineMadjWeMonthlyComponentByMetering.pdf"), dpi = 480)
+```
+
+```{r test model v2 total}
+myTitle <- "Total: mean daily total"
+myCaption <- "Model: v2_0 baseline with water efficiency for shower & WC"
+myPlot <- ba_IMPETUSmakeYearMonthPlot(hhFinalDataComboExpandedDT[model == "v2_0"], "sumDaily.baseline.madj.we")
+myPlot
+
+# recreate without the ribbons etc
+plotDT <- hhFinalDataComboExpandedDT[model == "v2_0", .(mean = mean(sumDaily.baseline.madj.we, na.rm = TRUE),
+                   nObs = .N), by = .(metered, obsDate, model)]
+  ggplot(plotDT, aes(x = obsDate, y = mean, colour = metered, group = metered)) +
+    geom_line() +
+    theme(legend.title = element_blank()) +
+    theme(legend.position = "bottom") +
+    labs(y = "Mean litres/day",
+         x = "Year")
+  ggsave(paste0("plots_v2/baselineMadjWeMonthlyTotalByMetering.pdf"), dpi = 480)
+```
+
+```{r model v2 with ci}
+# all uses mean daily
+dt <- hhFinalDataComboExpandedDT[model == "v2_0"]
+
+# create a monthly mean & 95% CI plot using the dt & var
+  plotDT <- dt[, .(meanDaily = mean(sumDaily.baseline.madj, na.rm = TRUE),
+                   sd = sd(sumDaily.baseline.madj, na.rm = TRUE),
+                   nObs = .N), by = .(metered, obsDate, currYear)]
+  plotDT <- plotDT[, ciLower := meanDaily - qnorm(0.975) * (sd/sqrt(nObs))]
+  plotDT <- plotDT[, ciUpper := meanDaily + qnorm(0.975) * (sd/sqrt(nObs))]
+  library(lubridate)
+  plotDT[, month := lubridate::month(obsDate, label = TRUE, abbr = TRUE)]
+  
+  ggplot(plotDT, aes(x = month, y = meanDaily, 
+                     linetype = metered,
+                     colour = as.factor(currYear),
+                     group = as.factor(currYear))) +
+    geom_line() +
+    geom_errorbar(aes(ymin = ciLower, ymax = ciUpper), width = 0.2) +
+    theme(legend.title = element_blank()) +
+    facet_grid(. ~ metered) +
+    theme(legend.position = "bottom") +
+    labs(y = "Mean litres/day",
+         x = "Month")
+  
+  ggsave(paste0("plots_v2/baselineMadjWeMonthlyByYearMetering.pdf"), dpi = 480)
+```
+
+
+```{r run baseline WE update}
+```

impetusModel/applyDroughtWEmodel.Rmd → impetusModel/applyDroughtModel.Rmd

@@ -305,7 +305,7 @@ Temporary use bans:
 
 We only run the TUB processes for each input model separately as they use an external consumption threshold which will be different across models.
 
-The others are run on the pooled model households as all processes are the same and we want as large a pool of 'uncovnerted; households as possible for the random selection processes.
+The others are run on the pooled model households as all processes are the same and we want as large a pool of 'unconverted; households as possible for the random selection processes.
 
 > Need to find ways to speed up these loops
 
@@ -319,7 +319,7 @@ Add the additional dual flush WCs.
 
 # sample households with single flush at appropriate rate
 # save the resulting data.table and add more rows for each drought phase (setting dualFlushWCdr as we go) then merge back onto original file to add dualFlushWCdr
-# how do I get it to do this _within_ the table?
+# how do I get it to do this _within_ the table?
 
 
 minYear <- min(hhFinalDataComboExpandedCEHDT$currYear)
@@ -337,14 +337,14 @@ for(y in minYear:maxYear){
   for(m in 1:12){
     #print(paste0("Month: ", m))
     notNormal <- 0 # assume normal
-    weMulti <- 1 # set to 1 here so that in normal years it has no effect - it is a multiplier used to increase the % who adopt in a non-nomral drough phase (See below)
+    weMulti <- 1 # set to 1 here so that in normal years it has no effect - it is a multiplier used to increase the % who adopt in a non-normal drought phase (See below)
     droughtPhase <- "1. Normal"
     monthYearDT <- hhFinalDataComboExpandedCEHDT[currYear == y & currMon == m]
     if(nrow(monthYearDT[Colne == "2. Developing",]) > 0){
       # we're in Developing
       droughtPhase <- "2. Developing"
       notNormal <- 1
-      weMulti <- 2/12 # used to multiply the annual dfWCAdopt rate (set in aparameters above).
+      weMulti <- 2/12 # used to multiply the annual dfWCAdopt rate (set in parameters above).
       # We need a monthly uptake in Developing drought to reflect increased WE efforts by water companies (2/12 could be a bit low?)
     }
     if(nrow(monthYearDT[Colne == "3. Drought",]) > 0){
@@ -686,7 +686,7 @@ hhFinalDataComboExpandedCEHDT <- hhFinalDataComboExpandedCEHDT[, sumDaily.baseli
 ```
 
 ## Test intervention effects by model
-i.e. compare to baseline without interventions. 
+i.e. compare baseline madj with we to baseline madj with we and dr
 
 ### WC
 
@@ -807,7 +807,22 @@ ba_IMPETUSaddDroughtPhases(myPlot, minY, maxY)
 
 ### Total
 
-```{r Compare sum}
+hhFinalDataComboExpandedCEHDT contains all the data
+
+we have two models
+
+`r table(hhFinalDataComboExpandedCEHDT$model)`
+
+ * v1_3 = original we uptake
+ * v2_0 = enhanced we uptake
+ 
+for both we have madj (monthly adjusted) we (water efficiency uptake) and dr (drought) scenarios:
+
+```{r testMeans}
+hhFinalDataComboExpandedCEHDT[, .(mean = mean(sumDaily.baseline.madj.we.dr, na.rm = TRUE)), keyby = .(model, currMonS)]
+```
+
+```{r selectData}
 plotTotalDT <- hhFinalDataComboExpandedCEHDT[, .(
   meanBaselineWeSum = mean(sumDaily.baseline.madj.we, na.rm = TRUE),
   sdBaselineWeSum = sd(sumDaily.baseline.madj.we, na.rm = TRUE),
@@ -817,23 +832,32 @@ plotTotalDT <- hhFinalDataComboExpandedCEHDT[, .(
                                               )
                                               , by = .(obsDate, model, metered)
                                           ]
+# calculate CIs if we want them
+
 plotTotalDT <- plotTotalDT[, yminWE := meanBaselineWeSum - qnorm(0.975)*(sdBaselineWeSum/sqrt(nObs))]
 plotTotalDT <- plotTotalDT[, ymaxWE := meanBaselineWeSum + qnorm(0.975)*(sdBaselineWeSum/sqrt(nObs))]
-
 plotTotalDT <- plotTotalDT[, yminDr := meanBaselineDrSum - qnorm(0.975)*(sdBaselineDrSum/sqrt(nObs))]
 plotTotalDT <- plotTotalDT[, ymaxDr := meanBaselineDrSum + qnorm(0.975)*(sdBaselineDrSum/sqrt(nObs))]
+```
+
+Let's compare the baseline we with baseline we + dr results for each version of the model
+
+### Per-model reporting
 
-# 1.3 ----
+#### Model v1_3 (synthetic)
+
+Model 1 - synthetic
+
+```{r compare_v_1.3, fig.height=4}
 myPlot <- ggplot(plotTotalDT[model == "v1_3"], aes(x = obsDate)) +
-  geom_ribbon(aes(ymin = yminWE, ymax = ymaxWE, fill = "Baseline 95% CI", group = metered),alpha = 0.5) +
+  # don't use CI
+#  geom_ribbon(aes(ymin = yminWE, ymax = ymaxWE, fill = "Baseline 95% CI", group = metered),alpha = 0.5) +
   geom_line(aes(y = meanBaselineWeSum, colour = "Baseline", linetype = metered, group = metered)) +
-  geom_ribbon(aes(ymin = yminDr, ymax = ymaxDr, fill = "Drought model 95% CI", group = metered),alpha = 0.5) +
+#  geom_ribbon(aes(ymin = yminDr, ymax = ymaxDr, fill = "Drought model 95% CI", group = metered),alpha = 0.5) +
   geom_line(aes(y = meanBaselineDrSum, colour = "Drought model", linetype = metered)) +
   theme(legend.title = element_blank()) +
-  facet_grid(model ~ .) +
   theme(legend.position = "bottom") +
-  labs(title = "Mean total usage",
-       y = "Mean l/hh/day",
+  labs(y = "Mean l/hh/day",
        x = "Date") 
 
 minDate <- plotTotalDT[, min(obsDate)]
@@ -842,8 +866,14 @@ maxDate <- plotTotalDT[, max(obsDate)]
 maxY <- max(plotTotalDT[model == "v1_3"]$ymaxWE)
 minY <- min(plotTotalDT[model == "v1_3"]$yminDr)
 myDrPlot <- ba_IMPETUSaddDroughtPhases(myPlot, minY, maxY)
-myDrPlot
+myDrPlot + 
+  guides(colour=guide_legend(ncol=1)) +
+ # guides(fill=guide_legend(ncol=1)) +
+  guides(linetype=guide_legend(ncol=1)) +
+ggsave(paste0("plots_v1/compareBaselineWEDroughtWE_v1.3.pdf"), height = 3, dpi = 480)
+```
 
+```{r iwaBath_paperFig}
 # Figure for IWA Bath final paper (http://ws.iwaponline.com/content/early/2018/02/13/ws.2018.035)
 # take the mean of the baseline & WE models (i.e. ignore metering as not mentioned in the text & caption etc)
 # also do not display drought phases (too complex to create key!)
@@ -864,23 +894,27 @@ minY <- min(paperPlotDT$meanBaselineDrSum)
 
 paperPlot <- ba_IMPETUSaddDroughtPhases(paperPlot, minY, maxY)
 paperPlot
-ggsave("Fig4_Compare_sum_model_v1_3.pdf", dpi = 400)
+ggsave(paste0("plots_v1/Fig4_Compare_sum_model_v1_3.pdf"), dpi = 400)
 
 # Grey scale version if required
 #myDrPlot <- myDrPlot + theme_bw()
 #ggsave("Fig4_Compare_sum_model_v1_3_gs.pdf", plot = myDrPlot, dpi = 400)
+```
 
-# 2.0 ----
+Now compare model v2 - enhanced we uptake rates
+
+#### Model v2 (SPRG based)
+
+```{r compare_v_2_0, fig.height=4}
 myPlot <- ggplot(plotTotalDT[model == "v2_0"], aes(x = obsDate)) +
-  geom_ribbon(aes(ymin = yminWE, ymax = ymaxWE, fill = "Baseline 95% CI", group = metered),alpha = 0.5) +
+  # don't use CI
+#  geom_ribbon(aes(ymin = yminWE, ymax = ymaxWE, fill = "Baseline 95% CI", group = metered),alpha = 0.5) +
   geom_line(aes(y = meanBaselineWeSum, colour = "Baseline", linetype = metered, group = metered)) +
-  geom_ribbon(aes(ymin = yminDr, ymax = ymaxDr, fill = "Drought model 95% CI", group = metered),alpha = 0.5) +
+#  geom_ribbon(aes(ymin = yminDr, ymax = ymaxDr, fill = "Drought model 95% CI", group = metered),alpha = 0.5) +
   geom_line(aes(y = meanBaselineDrSum, colour = "Drought model", linetype = metered)) +
   theme(legend.title = element_blank()) +
-  facet_grid(model ~ .) +
   theme(legend.position = "bottom") +
-  labs(title = "Mean total usage",
-       y = "Mean l/hh/day",
+  labs(y = "Mean l/hh/day",
        x = "Date") 
 
 minDate <- plotTotalDT[, min(obsDate)]
@@ -889,35 +923,22 @@ maxDate <- plotTotalDT[, max(obsDate)]
 maxY <- max(plotTotalDT[model == "v2_0"]$ymaxWE)
 minY <- min(plotTotalDT[model == "v2_0"]$yminDr)
 myDrPlot <- ba_IMPETUSaddDroughtPhases(myPlot, minY, maxY)
-myDrPlot
+myDrPlot + 
+  guides(colour=guide_legend(ncol=1)) +
+ # guides(fill=guide_legend(ncol=1)) +
+  guides(linetype=guide_legend(ncol=1)) +
+ggsave(paste0("plots_v2/compareBaselineWEDroughtWE_v2_0.pdf"), height = 4, dpi = 480)
 
-myPlot <- ggplot(plotTotalDT, aes(x = obsDate)) +
-  geom_ribbon(aes(ymin = yminWE, ymax = ymaxWE, fill = "Baseline 95% CI", group = metered),alpha = 0.5) +
-  geom_line(aes(y = meanBaselineWeSum, colour = "Baseline", linetype = metered, group = metered)) +
-  geom_ribbon(aes(ymin = yminDr, ymax = ymaxDr, fill = "Drought model 95% CI", group = metered),alpha = 0.5) +
-  geom_line(aes(y = meanBaselineDrSum, colour = "Drought model", linetype = metered)) +
-  theme(legend.title = element_blank()) +
-  facet_grid(model ~ .) +
-  theme(legend.position = "bottom") +
-  labs(title = "Mean total usage",
-       y = "Mean l/hh/day",
-       x = "Date") 
-
-minDate <- plotTotalDT[, min(obsDate)]
-maxDate <- plotTotalDT[, max(obsDate)]
-
-maxY <- max(plotTotalDT$ymaxWE)
-minY <- min(plotTotalDT$yminDr)
-myDrPlot <- ba_IMPETUSaddDroughtPhases(myPlot, minY, maxY)
-myDrPlot
+```
 
-plotTotalDT <- plotTotalDT[, pcSaved := 100*((meanBaselineWeSum - meanBaselineDrSum)/meanBaselineWeSum)] 
 
-```
+### Compare models
 
 The table below shows the mean and maximum savings over the period for each model.
 
 ```{r reduction tables}
+plotTotalDT[, pcSaved := 100*(meanBaselineWeSum - meanBaselineDrSum)/meanBaselineWeSum]
+
 t <- plotTotalDT[, .("Min reduction" = round(min(pcSaved),2),
                      "Mean reduction" = round(mean(pcSaved),2),
                      "Max reduction" = round(max(pcSaved),2)),
@@ -948,7 +969,7 @@ This is clearly shown in the chart below where:
  * sevDroughtCol <- "red" - Severe Drought
  * recoveringCol <- "yellowgreen" - Recovering
 
-```{r Chart savings}
+```{r model v1_3 savings}
 myPlot <- ggplot(plotTotalDT[model == "v1_3"], aes(x = obsDate, y = pcSaved, colour = metered)) + 
   geom_point() +
   theme(legend.title = element_blank()) +
@@ -982,38 +1003,25 @@ minY <- 0
 
 myPlot <- ba_IMPETUSaddDroughtPhases(myPlot, minY, maxY)
 myPlot
-ggsave("Fig5_Chart_savings_v1_3_by_month.pdf", plot = myPlot, dpi = 400)
+ggsave(paste0("plots_v1/Fig5_Chart_savings_v1_3_by_month.pdf"), plot = myPlot, dpi = 400)
 
 # Grey scale version if required
 #myPlot <- myPlot + theme_bw()
 #ggsave("Fig5_Chart_savings_v1_3_by_month_gs.pdf", plot = myPlot, dpi = 400)
+```
 
+```{r model v2_0 savings}
 myPlot <- ggplot(plotTotalDT[model == "v2_0"], aes(x = obsDate, y = pcSaved, colour = metered)) + 
-  geom_point() +
-  theme(legend.title = element_blank()) +
-  theme(legend.position = "bottom") +
-  facet_grid(model ~ .) +
-  labs(title = "% saving in total l/hh/day",
-       y = "%",
-       x = "Date") 
-
-maxY <- max(plotTotalDT$pcSaved)
-minY <- 0
-ba_IMPETUSaddDroughtPhases(myPlot, minY, maxY)
-
-myPlot <- ggplot(plotTotalDT, aes(x = obsDate, y = pcSaved, colour = metered)) + 
-  geom_point() +
+  geom_point() +s
   theme(legend.title = element_blank()) +
   theme(legend.position = "bottom") +
-  facet_grid(model ~ .) +
-  labs(title = "% saving in total l/hh/day",
-       y = "%",
+  labs(y = "%",
        x = "Date") 
 
 maxY <- max(plotTotalDT$pcSaved)
 minY <- 0
-
 ba_IMPETUSaddDroughtPhases(myPlot, minY, maxY)
+ggsave(paste0("plots_v2/baselineWeDroughtModelPcSavings.pdf"), dpi = 400)
 ```
 
 Interesting - savings under model v2 are larger. Why?
@@ -1043,7 +1051,7 @@ ggplot(droughtPhaseDT, aes(x = phase, y = V1, fill = phase)) +
   scale_fill_manual(values = c(develCol, droughtCol, sevDroughtCol, recoveringCol)
                     )
 
-ggsave("Fig4_5_DroughtPlotKey.pdf", dpi = 400)
+ggsave(paste0("plots_v1/Fig4_5_DroughtPlotKey.pdf"), dpi = 400)
 ```
 
 ## Extract drought & WE-adjusted hot water volumes (for BECC 2017 paper)
@@ -1054,5 +1062,5 @@ Use: hhFinalDataComboExpandedCEHDT
 summary(hhFinalDataComboExpandedCEHDT)
 ```
 
-```{r Run WE Model v1_0}
+```{r Run drought model}
 ```
\ No newline at end of file