R Notebook

#load needed packages
library(tidyverse)
library(MASS)
library(forecast)
library(fpp)
library(corrplot)
library(magrittr)
library(zoo)
library(RColorBrewer)
library(gridExtra)
library(Amelia)

#bring in data

train_features = read.csv('dengue_features_train.csv')
train_labels   = read.csv('dengue_labels_train.csv')
densubformat = read.csv('submission_format.csv')

#separate data by city
sj_train_features = train_features %>% filter(city == 'sj')
sj_train_labels   = train_labels   %>% filter(city == 'sj')

iq_train_features = train_features %>% filter(city == 'iq')
iq_train_labels   = train_labels   %>% filter(city == 'iq')

#data shape for each city
cat('\nSan Juan\n',
    '\t features: ', sj_train_features %>% ncol, 
    '\t entries: ' , sj_train_features %>% nrow,
    '\t labels: '  , sj_train_labels %>% nrow)

San Juan
     features:  24   entries:  936   labels:  936

cat('\nIquitos\n',
    '\t features: ', iq_train_features %>% ncol, 
    '\t entries: ' , iq_train_features %>% nrow,
    '\t labels: '  , iq_train_labels %>% nrow)

Iquitos
     features:  24   entries:  520   labels:  520

#Look at Data 
head(sj_train_features[1:7])

tail(sj_train_features)

head(iq_train_features[1:7])

tail(iq_train_features)

#Remove "week_start_date"
sj_train_features %<>% dplyr::select(-week_start_date)
iq_train_features %<>% dplyr::select(-week_start_date)

apply(sj_train_features, 2, function(x) 
  round(100 * (length(which(is.na(x))))/length(x) , digits = 1)) %>%
  as.data.frame() %>%
  `names<-`('Percent of Missing Values')

#missing values

missmap(sj_train_features)

#example plot of variable with out the NA fix

sj_train_features %>%
  mutate(index = as.numeric(row.names(.))) %>%
  ggplot(aes(index, ndvi_ne)) + 
  geom_line(colour = 'dodgerblue') +
  ggtitle("Vegetation Index over Time")

#fill in the NA 
sj_train_features$ndvi_ne %<>% na.locf(fromLast = TRUE)
iq_train_features$ndvi_ne %<>% na.locf(fromLast = TRUE)

#distrubution of labels SJ
cat('\nSan Juan\n',
    '\t total cases mean: ',      sj_train_labels$total_cases %>% mean(), 
    '\t total cases variance: ' , sj_train_labels$total_cases %>% var() )

San Juan
     total cases mean:  34.18056     total cases variance:  2640.045

#distrubution of labels IQ
cat('\nIquitos\n',
    '\t total cases mean: ',      iq_train_labels$total_cases %>% mean(), 
    '\t total cases variance: ' , iq_train_labels$total_cases %>% var() )

Iquitos
     total cases mean:  7.565385     total cases variance:  115.8955

#Total cases of disease
rbind(iq_train_labels, sj_train_labels) %>% 
  ggplot(aes(x = total_cases,fill = ..count..)) + 
  geom_histogram(bins = 12, colour = 'black') + ggtitle('Total Cases of Dengue') +
  scale_y_continuous(breaks = seq(0,700,100)) + facet_wrap(~city)

#correlations ?
sj_train_features %<>% mutate('total_cases' = sj_train_labels$total_cases)
iq_train_features %<>% mutate('total_cases' = iq_train_labels$total_cases)

#san juan corr
sj_train_features %>% 
  dplyr::select(-city, -year, -weekofyear) %>%
  cor(use = 'pairwise.complete.obs') -> M1

corrplot(M1, type="lower", method="color",
         col=brewer.pal(n=8, name="RdBu"),diag=FALSE)

#iquitos corr
iq_train_features %>% 
  dplyr::select(-city, -year, -weekofyear) %>%
  cor(use = 'pairwise.complete.obs') -> M2

corrplot(M2, type="lower", method="color",
         col=brewer.pal(n=8, name="RdBu"),diag=FALSE)

#correlation bar
sort(M1[21,-21]) %>%  
  as.data.frame %>% 
  `names<-`('correlation') %>%
  ggplot(aes(x = reorder(row.names(.), -correlation), y = correlation, fill = correlation)) + 
  geom_bar(stat='identity', colour = 'black') + scale_fill_continuous(guide = FALSE) + scale_y_continuous(limits =  c(-.15,.25)) +
  labs(title = 'San Juan\n Correlations', x = NULL, y = NULL) + coord_flip() -> cor1

# can use ncol(M1) instead of 21 to generalize the code
sort(M2[21,-21]) %>%  
  as.data.frame %>% 
  `names<-`('correlation') %>%
  ggplot(aes(x = reorder(row.names(.), -correlation), y = correlation, fill = correlation)) + 
  geom_bar(stat='identity', colour = 'black') + scale_fill_continuous(guide = FALSE) + scale_y_continuous(limits =  c(-.15,.25)) +
  labs(title = 'Iquitos\n Correlations', x = NULL, y = NULL) + coord_flip() -> cor2

grid.arrange(cor1, cor2, nrow = 1)

preprocessData <- function(data_path, labels_path = NULL)
{
  # load data 
  df <- read.csv(data_path)
  
  # features we want
  features = c("reanalysis_specific_humidity_g_per_kg", "reanalysis_dew_point_temp_k",
               "station_avg_temp_c", "station_min_temp_c")
  
  # fill missing values
  df[features] %<>% na.locf(fromLast = TRUE) 
  
  # add city if labels data aren't provided
  if (is.null(labels_path)) features %<>% c("city", "year", "weekofyear")
  
  # select features we want
  df <- df[features]
  
  # add labels to dataframe
  if (!is.null(labels_path)) df %<>% cbind(read.csv(labels_path))
  
  # filter by city
  df_sj <- filter(df, city == 'sj')
  df_iq <- filter(df, city == 'iq')
  
  # return a list with the 2 data frames 
  return(list(df_sj, df_iq))
}


#preprocess files
preprocessData(data_path = 'dengue_features_train.csv', labels_path = 'dengue_labels_train.csv') -> trains
sj_train <- trains[[1]]; iq_train <- as.data.frame(trains[2])

summary(sj_train)

 reanalysis_specific_humidity_g_per_kg reanalysis_dew_point_temp_k
 Min.   :11.72                         Min.   :289.6              
 1st Qu.:15.23                         1st Qu.:293.8              
 Median :16.83                         Median :295.4              
 Mean   :16.54                         Mean   :295.1              
 3rd Qu.:17.85                         3rd Qu.:296.4              
 Max.   :19.44                         Max.   :297.8              
 station_avg_temp_c station_min_temp_c city          year        weekofyear   
 Min.   :22.84      Min.   :17.80      iq:  0   Min.   :1990   Min.   : 1.00  
 1st Qu.:25.81      1st Qu.:21.70      sj:936   1st Qu.:1994   1st Qu.:13.75  
 Median :27.21      Median :22.80               Median :1999   Median :26.50  
 Mean   :27.00      Mean   :22.59               Mean   :1999   Mean   :26.50  
 3rd Qu.:28.18      3rd Qu.:23.90               3rd Qu.:2003   3rd Qu.:39.25  
 Max.   :30.07      Max.   :25.60               Max.   :2008   Max.   :53.00  
  total_cases    
 Min.   :  0.00  
 1st Qu.:  9.00  
 Median : 19.00  
 Mean   : 34.18  
 3rd Qu.: 37.00  
 Max.   :461.00  

summary(iq_train)

 reanalysis_specific_humidity_g_per_kg reanalysis_dew_point_temp_k
 Min.   :12.11                         Min.   :290.1              
 1st Qu.:16.12                         1st Qu.:294.6              
 Median :17.43                         Median :295.9              
 Mean   :17.10                         Mean   :295.5              
 3rd Qu.:18.19                         3rd Qu.:296.6              
 Max.   :20.46                         Max.   :298.4              
 station_avg_temp_c station_min_temp_c city          year        weekofyear   
 Min.   :21.40      Min.   :14.70      iq:520   Min.   :2000   Min.   : 1.00  
 1st Qu.:27.00      1st Qu.:20.60      sj:  0   1st Qu.:2003   1st Qu.:13.75  
 Median :27.57      Median :21.35               Median :2005   Median :26.50  
 Mean   :27.52      Mean   :21.20               Mean   :2005   Mean   :26.50  
 3rd Qu.:28.09      3rd Qu.:22.00               3rd Qu.:2007   3rd Qu.:39.25  
 Max.   :30.80      Max.   :24.20               Max.   :2010   Max.   :53.00  
  total_cases     
 Min.   :  0.000  
 1st Qu.:  1.000  
 Median :  5.000  
 Mean   :  7.565  
 3rd Qu.:  9.000  
 Max.   :116.000  

#split into train and test
sj_train_subtrain <- head(sj_train, 800)
sj_train_subtest  <- tail(sj_train, nrow(sj_train) - 800)

iq_train_subtrain <- head(iq_train, 400)
iq_train_subtest  <- tail(iq_train, nrow(sj_train) - 400)

#negative binomial model replication


mae <- function(error) return(mean(abs(error)) )

get_bst_model <- function(train, test)
{
  
  # Step 1: specify the form of the model
  form <- "total_cases ~ 1 +
            reanalysis_specific_humidity_g_per_kg +
            reanalysis_dew_point_temp_k + 
            station_avg_temp_c +
            station_min_temp_c"
  
  grid = 10 ^(seq(-8, -3,1))
  
  best_alpha = c()
  best_score = 1000
  
  # Step 2: Find the best hyper parameter, alpha
  for (i in grid)
  {
    model = glm.nb(formula = form,
                   data = train,
                   init.theta = i)
    
    results <-  predict(model, test)
    score   <-  mae(test$total_cases - results)
    
    if (score < best_score) {
      best_alpha <- i
      best_score <- score
      cat('\nbest score = ', best_score, '\twith alpha = ', best_alpha)
    }
  }
  
  # Step 3: refit on entire dataset
  combined <- rbind(train, test)
  combined_model = glm.nb(formula=form,
                          data = combined,
                          init.theta = best_alpha)
  
  return (combined_model)
}
sj_model <- get_bst_model(sj_train_subtrain, sj_train_subtest)

best score =  21.01167  with alpha =  1e-08

iq_model <- get_bst_model(iq_train_subtrain, iq_train_subtest)

best score =  6.421811  with alpha =  1e-08
best score =  6.421811  with alpha =  1e-06
best score =  6.421811  with alpha =  1e-05

#plot san juan

sj_train$fitted = predict(sj_model, sj_train, type = 'response')
sj_train %>% 
  mutate(index = as.numeric(row.names(.))) %>%
  ggplot(aes(x = index)) + ggtitle("San Juan") +
  geom_line(aes(y = total_cases, colour = "total_cases")) + 
  geom_line(aes(y = fitted, colour = "fitted"))

#plot Iquitos
iq_train$fitted = predict(iq_model, iq_train, type = 'response')
iq_train %>% 
  mutate(index = as.numeric(row.names(.))) %>%
  ggplot(aes(x = index)) + ggtitle("Iquitos") + 
  geom_line(aes(y = total_cases, colour = "total_cases")) + 
  geom_line(aes(y = fitted, colour = "fitted"))

#submit predictions
tests <- preprocessData('dengue_features_test.csv')
sj_test <- tests[[1]]; iq_test <- tests[[2]]

sj_test$predicted = predict(sj_model , sj_test, type = 'response')
iq_test$predicted = predict(iq_model , iq_test, type = 'response')

submissions = read.csv('submission_format.csv')
inner_join(submissions, rbind(sj_test,iq_test)) %>%
  dplyr::select(city, year, weekofyear, total_cases = predicted) ->
  predictions

Joining, by = c("city", "year", "weekofyear")

predictions$total_cases %<>% round()
write.csv(predictions, 'negbinomialpredictions.csv', row.names = FALSE)

#time series plot 

sj.ts<-ts(sj_train_features$total_cases, frequency=52, start=c(1990, 18))
plot(sj.ts)

iq.ts<-ts(iq_train_features$total_cases, frequency=52, start=c(2000, 26))
plot(iq.ts)

#train and test sets
sj.train.ts<-window(sj.ts, end=c(2004,37))
sj.test.ts<-window(sj.ts,start=c(2004,38))

iq.train.ts<-window(iq.ts, end=c(2008,27))
iq.test.ts<-window(iq.ts,start=c(2008,28))

#neural network models

sj.fit.nnet<-nnetar(sj.train.ts)
sj.fit.nnet

Series: sj.train.ts 
Model:  NNAR(14,1,8)[52] 
Call:   nnetar(y = sj.train.ts)

Average of 20 networks, each of which is
a 15-8-1 network with 137 weights
options were - linear output units 

sigma^2 estimated as 44.45

iq.fit.nnet<-nnetar(iq.train.ts)
iq.fit.nnet

Series: iq.train.ts 
Model:  NNAR(3,1,2)[52] 
Call:   nnetar(y = iq.train.ts)

Average of 20 networks, each of which is
a 4-2-1 network with 13 weights
options were - linear output units 

sigma^2 estimated as 34.12

#neural net forecast

sj.fcast.nnet<-forecast(sj.fit.nnet,h=260)
plot(sj.fcast.nnet) 
lines(sj.test.ts, col="red")  

legend("topright",lty=1,col=c("red","blue"),c("actual values","forecast"))

iq.fcast.nnet<-forecast(iq.fit.nnet,h=156)
plot(iq.fcast.nnet) 
lines(iq.test.ts, col="red")  

legend("topright",lty=1,col=c("red","blue"),c("actual values","forecast"))

#accuracy
acc.sj.nnet <- accuracy(sj.fcast.nnet, sj.test.ts)
acc.sj.nnet

                      ME      RMSE       MAE  MPE MAPE      MASE        ACF1
Training set 0.006646099  6.667078  4.885579 -Inf  Inf 0.1225886 -0.01376727
Test set     7.338365616 32.377485 18.987895 -Inf  Inf 0.4764430  0.93389023
             Theil's U
Training set        NA
Test set             0

tsdisplay(residuals(sj.fcast.nnet))

acc.iq.nnet <- accuracy(iq.fcast.nnet,iq.test.ts)
acc.iq.nnet

                      ME      RMSE      MAE  MPE MAPE      MASE         ACF1
Training set -0.01527385  5.840966 3.543559 -Inf  Inf 0.3871471 -0.001184902
Test set      5.28822909 12.685205 7.300727 -Inf  Inf 0.7976317  0.844531021
             Theil's U
Training set        NA
Test set             0

tsdisplay(residuals(iq.fcast.nnet))

#solutions for neural net model
NNET_sj_sol=data.frame(densubformat[1:260, -4],total_cases=round(sj.fcast.nnet$mean))
NNET_iq_sol=data.frame(densubformat[261:416,-4],total_cases=round(iq.fcast.nnet$mean))

Dengue_NNET_solution=bind_rows(NNET_sj_sol,NNET_iq_sol)

Vectorizing 'ts' elements may not preserve their attributesVectorizing 'ts' elements may not preserve their attributes

write.csv(Dengue_NNET_solution,file="Dengue_NNET_predictions.csv",row.names=F)

#arima model SJ

sj.arima=Arima(sj.train.ts,order=c(1,1,2),seasonal=c(0,0,0))
sj.arima

Series: sj.train.ts 
ARIMA(1,1,2) 

Coefficients:
          ar1     ma1     ma2
      -0.7820  0.9503  0.1979
s.e.   0.2153  0.2147  0.0347

sigma^2 estimated as 192.5:  log likelihood=-3023.21
AIC=6054.43   AICc=6054.48   BIC=6072.89

sj.arima.fc=forecast(sj.arima,h=260)
acc.sj.arima = accuracy(sj.arima.fc,sj.test.ts)
acc.sj.arima

                       ME     RMSE       MAE  MPE MAPE      MASE        ACF1
Training set  0.006577048 13.83894  8.379156  NaN  Inf 0.2102492 0.007353076
Test set     14.740977080 34.48529 18.234742 -Inf  Inf 0.4575449 0.932173942
             Theil's U
Training set        NA
Test set             0

plot(sj.arima.fc)

tsdisplay(residuals(sj.arima.fc))

#arima modle IQ
iq.arima=Arima(iq.train.ts,order=c(0,1,2),seasonal=c(0,0,1))
iq.arima

Series: iq.train.ts 
ARIMA(0,1,2)(0,0,1)[52] 

Coefficients:
          ma1      ma2    sma1
      -0.2731  -0.3054  0.0385
s.e.   0.0469   0.0479  0.0437

sigma^2 estimated as 53.12:  log likelihood=-1418.7
AIC=2845.39   AICc=2845.49   BIC=2861.52

iq.arima.fc=forecast(iq.arima,h=156)
acc.iq.arima = accuracy(iq.arima.fc,iq.test.ts)
acc.iq.arima

                      ME      RMSE      MAE  MPE MAPE      MASE        ACF1
Training set 0.008711802  7.253488 3.808308  NaN  Inf 0.4160719 -0.00938541
Test set     7.897047615 13.832998 8.248390 -Inf  Inf 0.9011674  0.84182636
             Theil's U
Training set        NA
Test set             0

plot(iq.arima.fc)

tsdisplay(residuals(iq.arima.fc))

#solution for arima model
Arima_sj_sol=data.frame(densubformat[1:260,-4],total_cases=round(sj.arima.fc$mean))
Arima_iq_sol=data.frame(densubformat[261:416,-4],total_cases=round(iq.arima.fc$mean))

Dengue_Arima_solution=bind_rows(Arima_sj_sol,Arima_iq_sol)

Vectorizing 'ts' elements may not preserve their attributesVectorizing 'ts' elements may not preserve their attributes

write.csv(Dengue_Arima_solution,file="Dengue_Arima_predictions.csv",row.names=F)