library(dplyr)
library(readr)
library(ggplot2)
library(tidyr)
library(tidyverse)
library(caret)
library(glmnet)

newspopData <- read_csv('OnlineNewsPopularity.csv')

# select specified data channel from params and drop other data channel columns
newspopData <- newspopData %>% filter(newspopData[params$data_channel] == 1) %>%
  select( -starts_with('data_channel_is_'), -url)


newspopData

summary(newspopData$shares)

pShares <- newspopData$shares
pShares <- sort(pShares)

quantileValues <- quantile(pShares, c(.25, .5, .75)) 

newspopData <- newspopData %>% mutate(article_popularity_index = if_else(shares<quantileValues[1], "4",
                                                                  if_else(shares<quantileValues[2], "3",
                                                                  if_else(shares<quantileValues[3],"2", "1"))))

newspopData$article_popularity_index <- as_factor(newspopData$article_popularity_index)

newspopData$scaledShare <- scale(newspopData$shares, center = T, scale = T)

#calculating correlation index
corrIndex <- round(cor(newspopData$scaledShare, newspopData$n_tokens_content),2)

# plotting the shares as a function of words in an article
ggplot(newspopData, aes(x= n_tokens_content, y = scaledShare)) +
  geom_point()+
  labs(subtitle = 'Word Count v/s Number of Shares',
      y = 'Number of Share (Scaled)',
      x = 'Number of Words in article', 
      title = toupper(str_replace_all(params$data_channel, "_", " ")),
      caption = 'Source: News Popularity Dataset') + 
      geom_text(x = max(newspopData$n_tokens_content)/2, 
                y = max(newspopData$scaledShare)/2, size = 4, 
                label = paste0("Correlation coefficient = ", corrIndex), color = 'red')

newspopData <- newspopData %>% select(-scaledShare)
newspopData <- newspopData %>% mutate(day = if_else(weekday_is_monday==   1, "Monday", 
                                            if_else(weekday_is_tuesday==  1, "Tuesday",
                                            if_else(weekday_is_wednesday==1, "Wednesday",
                                            if_else(weekday_is_thursday== 1, "Thursday",
                                            if_else(weekday_is_friday==   1, "Friday",
                                            if_else(weekday_is_saturday== 1,"Saturday",
                                            if_else(weekday_is_sunday==   1, "Sunday", 
                                                    "-"))))))))

newspopData$day <- as_factor(newspopData$day)

knitr::kable(newspopData %>% group_by(day) %>% summarize(Count_of_Articles = n()))

knitr::kable(newspopData %>% group_by(day, article_popularity_index) %>% summarize(Count_of_Articles = n()))

ggplot(newspopData, aes(x = day, y = shares/1000000)) +
  geom_bar(stat = 'identity', width = 0.5, fill = 'tomato3')+
  labs(subtitle = 'Number of Shares (Million) Vs Day of Week',
       caption = 'Source : News Popularity Dataset',
       y = 'Total Share count in Million',
       x = 'Day of the Week',
       title = toupper(str_replace_all(params$data_channel, "_", " ")),) +
       theme(axis.text.x = element_text(angle = 65, vjust = 0.6)) + 
       theme(plot.caption = element_text(size=9, color="red", face="italic", hjust = 1))

ggplot(newspopData, aes(x = global_rate_positive_words, y = global_sentiment_polarity)) +
       geom_point(aes(col = day, shape = article_popularity_index)) + 
       geom_smooth(aes(col = day), method = 'lm', se = F) + 
       scale_shape_discrete(name="Article Popularity Index", labels = c("4", "3", "2", "1"))+
       labs(subtitle = 'Positive Rate VS Sentiment Polarity Plot',
       x = 'Global Positive Word Rate',
       y = 'Global Sentiment Polarity',
       title = toupper(str_replace_all(params$data_channel, "_", " ")),
       caption = "Source : News Popularity Dataset",
       color = 'Days',
       size = 'Shares per Million') + 
       theme(plot.caption = element_text(size=9, color="red", face="italic", hjust = 1))

ggplot(newspopData, aes(x = global_rate_negative_words, y = global_sentiment_polarity)) + 
      geom_point(aes (col = day, shape = article_popularity_index)) + 
      geom_smooth(aes(col = day), method = 'lm', se = F) + 
      scale_shape_discrete(name="Article Popularity Index", labels = c("4", "3", "2", "1"))+
      labs(subtitle = 'Positive Rate VS Sentiment Polarity Plot',
       x = 'Global Negative Word Rate',
       y = 'Global Sentiment Polarity',
       title = toupper(str_replace_all(params$data_channel, "_", " ")),
       caption = "Source : News Popularity Dataset",
       color = 'Days',
       size = 'Shares per Million') + 
       theme(plot.caption = element_text(size=9, color="red", face="italic", hjust = 1))  

ggplot(newspopData, aes(x=day, y =   shares/1000000)) + 
  geom_bar(aes(fill = as_factor(num_imgs)), stat="identity", position="dodge") + 
  scale_fill_discrete(name = 'Number of Images \n in Article') +
  labs(subtitle = 'Number of Images in article v/s number of shares',
       x = 'Number of Images',
       y = 'Number of shares (in millions)',
       title = toupper(str_replace_all(params$data_channel, "_", " ")),
       caption = "Source : News Popularity Dataset") + 
       theme(plot.caption = element_text(size=9, color="red", face="italic", hjust = 1))

ggplot(newspopData, aes(x =  num_keywords, y =  shares/1000000)) + 
  geom_bar(aes(fill = day), stat="identity") + 
  scale_x_discrete(limits=c(2:10)) +
  labs(subtitle = 'Shares based on number of keywords',
       x = 'Number of Keywords',
       y = 'Shares (in millions)',
       title = toupper(str_replace_all(params$data_channel, "_", " ")),
       caption = "Source : News Popularity Dataset",
       fill = 'Days') + 
       theme(plot.caption = element_text(size=9, color="red", face="italic", hjust = 1))

pols = c("Extremely Positive", "Positive", "Negative", "Extremely Negative")

newspopData <- newspopData %>% mutate(title_Sentiment_Class = 
                                        if_else(title_sentiment_polarity > 0.5, pols[1],
                                        if_else(title_sentiment_polarity > 0, pols[2],
                                        if_else(title_sentiment_polarity > -0.5, pols[3], pols[4]))))


newspopData$title_Sentiment_Class <- as_factor(newspopData$title_Sentiment_Class)

knitr::kable(newspopData %>% group_by(title_Sentiment_Class) %>% summarize(Count_of_articles = n()))

knitr::kable(newspopData %>% group_by(title_Sentiment_Class, article_popularity_index) %>% summarize(Count_of_articles = n()))

knitr::kable(newspopData %>% group_by(title_Sentiment_Class, day) %>% summarize(Count_of_articles = n()))

newspopData <- newspopData %>% mutate(content_Sentiment_Class = 
                                        if_else( global_sentiment_polarity > 0.5, pols[1],
                                        if_else( global_sentiment_polarity > 0, pols[2],
                                        if_else( global_sentiment_polarity > -0.5, pols[3], pols[4]))))


newspopData$content_Sentiment_Class <- as_factor(newspopData$content_Sentiment_Class)

# average article length 
knitr::kable(newspopData %>%
  group_by(content_Sentiment_Class, day)%>%
    select(day, content_Sentiment_Class, n_tokens_content) %>%
    summarise(Average_Content_Length = round(sum(n_tokens_content)/n())))

ggplot(newspopData, aes(x=title_Sentiment_Class, y=shares/1000000)) +
  geom_bar(aes(fill = day), stat="identity", position="dodge") +
  labs(subtitle = 'Shares based on title sentiment and day of the week',
       x = 'Title Sentiment Category',
       y = 'Shares (in millions)',
       title = toupper(str_replace_all(params$data_channel, "_", " ")),
       caption = "Source : News Popularity Dataset",
       fill = 'Days') +
      theme(plot.caption = element_text(size=9, color="red", face="italic", hjust = 1)) + 
       geom_text(x = 3, y = 0.6, size = 3, 
                label = paste("Title Sentiment Score range", "Extremely Positive = 0.5 <--> 1", 
                              'Positive = 0 <--> 0.5 ', 
                              'Negative = -0.5 <--> 0', 
                              'Extremely Negative = -1 <--> -0.5', 
                              sep = '\n' ), color = 'black')

ggplot(newspopData, aes(x= title_sentiment_polarity)) +
  geom_density(aes(fill=title_Sentiment_Class), alpha = 0.7) +
  facet_grid(. ~day)+
  scale_fill_discrete(name = "Title Sentiment") +
  labs(subtitle = 'Distribution of articles on title sentiment and day of the week',
       x = 'Title Sentiment Polarity',
       y = 'Density',
       title = toupper(str_replace_all(params$data_channel, "_", " ")),
       caption = "Source : News Popularity Dataset") + 
       theme(plot.caption = element_text(size=9, color="red", face="italic", hjust = 1))

ggplot(newspopData) + geom_density(aes(x = shares))+ 
  labs(title = "Shares Predictor Distribution in Dataset")

ggplot(newspopData) + geom_density(aes(x = log(shares)))+
  labs(title = "Logarithmic transformation of shares distribution")

set.seed(42)

# Removing mutated variables from the data set
lmData <- newspopData %>% select(-day, -title_Sentiment_Class, -content_Sentiment_Class, -article_popularity_index)

lmData$logShares <- log(lmData$shares)

# train test split
index <- createDataPartition(newspopData$shares, p = 0.7, list = FALSE)

# Removing the shares variable because we are using logShares as response
trainDf <- lmData[index, ] %>% select(-shares)
testDf  <- lmData[-index, ] %>% select(-shares)

modelLmVar <- lm(logShares~. ,data = trainDf)

# getting p-values 
pVal <- data.frame(summary(modelLmVar)$coefficients[,4])
pVal$row_names <- row.names(pVal)

# checking if slope of intercept is significant (i.e. <<< 0.05)
pVal <- pVal %>% filter(pVal[,1] < 0.05)
cols = c()

# remove intercept term and store all significant intercepts in "cols" variable 
if (pVal[,2][1] == "(Intercept)"){
  cols = pVal[,2][-1]
}else{
  cols = pVal[,2]
}
cols

trControl <- trainControl(method = 'repeatedcv', number = 10)

# model training
modelLm <- train(logShares ~ .,
              data = trainDf %>% 
              select(all_of(cols), logShares),
              method = "lm",
              preProcess = c('center', 'scale'),
              trControl = trControl)

title = toupper(str_replace_all(params$data_channel, "_", " "))

# variable importance for linear model
plot(varImp(modelLm), 
     main = paste("Variable importance Plot for ", "\n", title))

lmPredict <- predict(modelLm, newdata = testDf %>% select(all_of(cols), logShares))
lmRes <- postResample(lmPredict, obs = testDf$logShares)


# Adjusted R-square
n <- dim(testDf)[1]
p <- dim(testDf %>% select(all_of(cols)))[2]
adjR2Linear <-  1-(((1-lmRes[2]^2)*(n-1))/(n-p-1))  
 
lmRes

set.seed(43)

# train
fitControl <- trainControl(method = "repeatedcv",number = 5)
fitLASSO <- train(logShares ~ . ,
                  data=trainDf,
                  method = "lasso", 
                  trControl = fitControl,
                  preProcess = c("center", "scale")
                )


predLASSO <- predict(fitLASSO, newdata = testDf)

coeff <- predict(fitLASSO$finalModel, type = "coef", mode = "fraction",
                 s = fitLASSO$bestTune$fraction, testDf)


lassoRes <- postResample(predLASSO, obs = testDf$logShares)


# Adjusted R-square
n <- dim(testDf)[1]
p <- dim(testDf)[2]
adjR2Lasso <-  1-(((1-lassoRes[2]^2)*(n-1))/(n-p-1))

lassoRes

set.seed(200)

newspopData$logShares <- log(newspopData$shares)

trainIndex <- createDataPartition(newspopData$shares, p = 0.7, list = FALSE)

trainData <- newspopData[trainIndex, ]
testData <- newspopData[-trainIndex, ]

# Removing shares, logShares and mutated variables to form PCs
PCs <- prcomp(trainData %>% 
                select(-shares, -logShares, -day, 
                       -title_Sentiment_Class, 
                       -content_Sentiment_Class,
                       -article_popularity_index), 
                        center=TRUE, scale=TRUE)

# cumulative variance
var = cumsum(PCs$sdev**2/sum(PCs$sdev**2))

show(plot(var, xlab = "Principal Component", 
          ylab = "Cum. Prop of Variance Explained", ylim = c(0, 1), type = 'b'))

pt <- min(which(var > 0.80))

points(pt, var[pt], col="red", cex=2, pch=20)
text(8, 0.8, col="red", paste0(round(100*var[pt],2), "% variance \nis explained by ", pt, " \nnumber of PC(s)"))

# Creating training data PCs without the logShares variables
treeTrainData <- predict(PCs, newdata = trainData %>% select(-shares, -logShares)) %>% 
  as_tibble() %>% 
  select(PC1:paste0("PC", pt))

treeTrainData$logShares <- trainData$logShares

# This the final training data that will be used to train random forest model
treeTrainData

set.seed(400)
trctrl <- trainControl(method = "repeatedcv", number = 3)
randomForest <- train(logShares ~., data = treeTrainData, 
                      method = "rf",
                      trControl=trctrl,
                      preProcess = c("center", "scale"),
                      tuneGrid = expand.grid(mtry = 1:as.integer(pt/3)))

randomForest

# Creating test data PCs
treeTestData <- predict(PCs, newdata = testData %>% select(-shares, -logShares)) %>% 
                         as_tibble() %>% 
                         select(PC1:paste0("PC", pt))

# Completing the test data
treeTestData$logShares <- testData$logShares

rfPredict <- predict(randomForest, newdata = treeTestData)
randomFRes <- postResample(rfPredict, obs = treeTestData$logShares)


# Adjusted R-square
n <- dim(treeTestData)[1]
p <- dim(treeTestData)[2] - 1
adjR2Rf <-  1-(((1-randomFRes[2]^2)*(n-1))/(n-p-1))
 
randomFRes

# Training the model
boostedTree <- train(logShares ~., data = trainData %>% select(-shares), method = "gbm",
                     trControl = trainControl(method = "repeatedcv", number = 3),
                     preProcess = c("center", "scale"),
                     tuneGrid = expand.grid(n.trees = c(25, 50, 100, 150, 200),
                                            interaction.depth = 1:4,
                                            shrinkage = 0.1,
                                            n.minobsinnode = 10),
                     verbose = FALSE)

# Best combination of tuning hyperparameters
print(boostedTree$bestTune)

boostedTPredict <- predict(boostedTree, newdata = testData %>% select(-shares))

boostedTRes <- postResample(boostedTPredict, obs = (testData %>% select(-shares))$logShares)

# plot
plot(boostedTree, xlab = "Max Tree Depth ", 
     ylab = "Root Mean Square Error (Found through CV)", 
     main = "All Predictors - Minimizing RMSE")

# Adjusted R-square
n <- dim(testData)[1]
p <- dim(testData)[2] - 1
adjR2Boosted <-  1-(((1-boostedTRes[2]^2)*(n-1))/(n-p-1))
 

boostedTRes

model_types = c("Linear", "Ensemble")

modelComparisonDf = tibble(model = c("Linear Model"), model_type = c(model_types[1]), RMSE = c(round(lmRes[[1]],3)), R_Squared = c(round(lmRes[[2]], 4)), Adj_R2 = c(round(adjR2Linear,4)) )

modelComparisonDf <- rbind(modelComparisonDf, c("LASSO", model_types[1], round(lassoRes[[1]],3), round(lassoRes[[2]],4), round(adjR2Lasso,4)))
modelComparisonDf <- rbind(modelComparisonDf, c("Random Forest", model_types[2], round(randomFRes[[1]],3), round(randomFRes[[2]],4), round(adjR2Rf,4)))
modelComparisonDf <- rbind(modelComparisonDf, c("Boosted Trees", model_types[2], round(boostedTRes[[1]],3), round(boostedTRes[[2]],4), round(adjR2Boosted,4)))

knitr::kable(modelComparisonDf)

bestLinearModel = ""
bestLinearRMSE = ""
bestLinearR2 = ""

bestEnsembleModel = ""
bestEnsembleRMSE = ""
bestEnsembleR2 = ""

bestOverallModel = ""
bestOverallRMSE = ""
bestOverallR2 = ""

lineardf <- modelComparisonDf %>% filter(model_type == "Linear")

ensembleDf <- modelComparisonDf %>% filter(model_type == "Ensemble")

if (length(unique(lineardf$R_Squared)) > 1){
  bestLinearModel <- lineardf[which.min(lineardf$RMSE),][[1]]
  bestLinearRMSE <- lineardf[which.min(lineardf$RMSE),][[3]]
  bestLinearR2 <- lineardf[which.min(lineardf$RMSE),][[4]]
} else{
  bestLinearModel <- lineardf[which.max(lineardf$Adj_R2),][[1]]
  bestLinearRMSE <- lineardf[which.max(lineardf$Adj_R2),][[3]]
  bestLinearR2 <- lineardf[which.max(lineardf$Adj_R2),][[4]]
}


if (length(unique(ensembleDf$R_Squared)) > 1){
  bestEnsembleModel <- ensembleDf[which.min(ensembleDf$RMSE),][[1]]
  bestEnsembleRMSE <- ensembleDf[which.min(ensembleDf$RMSE),][[3]]
  bestEnsembleR2 <- ensembleDf[which.min(ensembleDf$RMSE),][[4]]
} else{
  bestEnsembleModel <- ensembleDf[which.max(ensembleDf$Adj_R2),][[1]]
  bestEnsembleRMSE <- ensembleDf[which.max(ensembleDf$Adj_R2),][[3]]
  bestEnsembleR2 <- ensembleDf[which.max(ensembleDf$Adj_R2),][[4]]
}


if (length(unique(modelComparisonDf$R_Squared)) > 1){
  bestOverallModel <- modelComparisonDf[which.min(modelComparisonDf$RMSE),][[1]]
  bestOverallRMSE <- modelComparisonDf[which.min(modelComparisonDf$RMSE),][[3]]
  bestOverallR2 <- modelComparisonDf[which.min(modelComparisonDf$RMSE),][[4]]
} else{
  bestOverallModel <- modelComparisonDf[which.max(modelComparisonDf$Adj_R2),][[1]]
  bestOverallRMSE <- modelComparisonDf[which.max(modelComparisonDf$Adj_R2),][[3]]
  bestOverallR2 <- modelComparisonDf[which.max(modelComparisonDf$Adj_R2),][[4]]
}