diff --git a/impetusModel/addBaselineWE.Rmd b/impetusModel/addBaselineWE.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..8c1009cab1d0dabb93c804fe6b8cf8e720a19589
--- /dev/null
+++ b/impetusModel/addBaselineWE.Rmd
@@ -0,0 +1,555 @@
+# 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}
+```
diff --git a/impetusModel/applyDroughtWEmodel.Rmd b/impetusModel/applyDroughtModel.Rmd
similarity index 94%
rename from impetusModel/applyDroughtWEmodel.Rmd
rename to impetusModel/applyDroughtModel.Rmd
index dc7d023dd7a1fd2149772f936e54ad86fa4c34aa..e2e687f2debf39838299f13a60329aa8cc2be7d0 100644
--- a/impetusModel/applyDroughtWEmodel.Rmd
+++ b/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