ampli/R/combmodels.R

# Combined models
# Continuation of clusterviz.R and weathmod.R

library(TSA)
library(caTools)
library(dplyr)
library(ggplot2)
library(reticulate)
library(tidyr)
library(MASS)
theme_set(theme_bw())
use_virtualenv("../venv/")

p <- import("pandas")
sns <- import("seaborn")
aggdf <- p$read_pickle("../data/9-clusters.agg.pkl")
aggdf$cluster <- factor(aggdf$cluster)
clusters <- levels(aggdf$cluster)
str(aggdf)
mtempdf <- read.csv("../data/weatherharm.csv", stringsAsFactors = FALSE) %>% 
    mutate(x = as.POSIXct(x, tz = "UTC")) %>%
    rename(read_time = x, rollingmin = y, fitmin = f, resmin = r)
str(mtempdf)

ntps <- length(unique(aggdf$read_time))

clus = "9"


harm.y <- ts(1:ntps, frequency = floor(48 * 365.25)) %>% harmonic(2)
harm.w <- ts(1:ntps, frequency = floor(48 * 7))      %>% harmonic(3)
harm.d <- ts(1:ntps, frequency = floor(48))          %>% harmonic(3)
colnames(harm.y) <- sprintf("%s.%s.%s", "year", rep(c("cos", "sin"), each = ncol(harm.y)/2), rep(1:(ncol(harm.y)/2), times = 2))
colnames(harm.w) <- sprintf("%s.%s.%s", "week", rep(c("cos", "sin"), each = ncol(harm.w)/2), rep(1:(ncol(harm.w)/2), times = 2))
colnames(harm.d) <- sprintf("%s.%s.%s", "day",  rep(c("cos", "sin"), each = ncol(harm.d)/2), rep(1:(ncol(harm.d)/2), times = 2))

clusdf <- filter(aggdf, cluster == clus) %>% 
    dplyr::select(read_time, kwh = kwh_tot_mean) %>% 
    left_join(mtempdf, by = "read_time") %>% cbind(harm.y, harm.w, harm.d)
str(clusdf)

ycols <- paste(colnames(harm.y), collapse = " + ") 
wcols <- paste(colnames(harm.w), collapse = " + ") 
dcols <- paste(colnames(harm.d), collapse = " + ") 

nform.full <- sprintf("kwh ~ %s + %s + %s + (%s):(%s) + (%s):(%s) + (%s):(%s) + resmin + resmin:(%s) + resmin:(%s) + resmin:(%s)", 
        ycols, wcols, dcols, ycols, wcols, ycols, dcols, wcols, dcols, ycols, wcols, dcols) %>% formula()
nform.comp <- sprintf("kwh ~ %s + %s + %s + (%s):(%s) + (%s):(%s) + resmin + resmin:(%s) + resmin:(%s) + resmin:(%s)", 
        ycols, wcols, dcols, ycols, dcols, wcols, dcols, ycols, wcols, dcols) %>% formula()
nform.now <- sprintf("kwh ~ %s + %s + %s + (%s):(%s) + (%s):(%s)", 
        ycols, wcols, dcols, ycols, dcols, wcols, dcols) %>% formula()
nform.min <- formula("kwh ~ 1")
nform.start <- sprintf("kwh ~ %s + %s + %s + resmin", 
        ycols, wcols, dcols) %>% formula()

# charmmod <- lm(kwh ~ resmin + harm.y * harm.w * harm.d + resmin:harm.y, data = clusdf)
charmmod <- lm(nform.comp, data = clusdf)
# charmmod <- lm(nform.full, data = clusdf)
# charmmod <- lm(kwh ~ ., data = clusdf)
summary(charmmod)

mean(abs(lm(nform.now,  data = clusdf)$residuals))
mean(abs(lm(nform.comp, data = clusdf)$residuals))
mean(abs(lm(nform.full, data = clusdf)$residuals))
sd(lm(nform.now,  data = clusdf)$residuals)
sd(lm(nform.comp, data = clusdf)$residuals)
sd(lm(nform.full, data = clusdf)$residuals)

cmdf <- data.frame(x = clusdf$read_time, y = clusdf$kwh, f = fitted(charmmod), r = resid(charmmod))
cmplot <-ggplot(cmdf, aes(x = x, y = y)) + geom_line(aes(y = f), color = "blue", size = 2) + geom_point() +
    geom_point(aes(y = r), color = "darkgreen")

cmplot

cmplot + coord_cartesian(xlim = c(as.POSIXct("2017-08-01", tz = "UTC"), as.POSIXct("2017-09-01", tz = "UTC")))

# sres <- stepAIC(charmmod, scope = list(upper = nform.full, lower = nform.min),
#                 direction = "both", steps = 300)