<file_sep>###############################################
## PS 5625 - Applied Statistical Programming
## Problem Set 5
## Author: <NAME>
## Tasks: 1 - Develop the package FitStatisticsPack

rm(list = ls())

# Load libraries
library(devtools)
library(roxygen2)

# Define WD
setwd("~/Google Drive/2. Academicos/6. Doutorado/Class/2017/Spring/PS 5625 - Applied Statistical Programming/Homeworks/Problem Sets/PS5/")

# Generate the Package Directory:
if(!dir.exists("FitStatisticsPack")){
   create("FitStatisticsPack", rstudio = FALSE)
}

# Add the function into R file.
file.copy(from = "PS5_CreateFunction.R" , to = "FitStatisticsPack/R/") 
# Rename the function file.
file.rename(from = "FitStatisticsPack/R/PS5_CreateFunction.R", to = "FitStatisticsPack/R/FitStats.R")

# Package the code with the new DESCRIPTION file
current.code <- as.package("FitStatisticsPack")

# Generate the Documentation 
document(current.code)

# Edit the DESCRIPTION FILE
descriptionfile <- read.dcf("FitStatisticsPack/DESCRIPTION")
descriptionfile[1, 2] <- "Calculate fit statistics" 
descriptionfile[1, 3] <- "0.1"
descriptionfile[1, 4] <- "person(\"Patrick\", \"<NAME>\", email = \"<EMAIL>\", role = c(\"aut\", \"cre\"))"
descriptionfile[1, 5] <- "Calculates fit statistics for statistical models."
descriptionfile[1, 6] <- c("\nR (>= 3.3.2), \n stats")
descriptionfile[1, 7] <- "GPL (>= 2)"
# Add the file "Imports" to DESCRIPTION FILE.
Imports <- "testthat"
descriptionfile <- cbind(descriptionfile, Imports)

# Save the new version of the  DESCRIPTION FILE
write.dcf(descriptionfile, "FitStatisticsPack/DESCRIPTION")

# Generate tests folder
if(!dir.exists("FitStatisticsPack/tests/")){
   dir.create("FitStatisticsPack/tests/")
}

# Generate testthat folder
if(!dir.exists("FitStatisticsPack/tests/testthat")){
   dir.create("FitStatisticsPack/tests/testthat")
}

# Add testall to test folder
file.copy(from = "PS5_test_all.R" , to = "FitStatisticsPack/tests/") 

# Rename testall file.
file.rename(from = "FitStatisticsPack/tests/PS5_test_all.R", 
            to = "FitStatisticsPack/tests/test-all.R")

# Add tests into testthat folder
file.copy(from = "PS5_test_class.R", to = "FitStatisticsPack/tests/testthat")
file.copy(from = "PS5_test_error_Pclass.R", to = "FitStatisticsPack/tests/testthat")

# Rename files 
file.rename(from = "FitStatisticsPack/tests/testthat/PS5_test_class.R", 
            to = "FitStatisticsPack/tests/testthat/test_class.R")
file.rename(from = "FitStatisticsPack/tests/testthat/PS5_test_error_Pclass.R", 
            to = "FitStatisticsPack/tests/testthat/test_Pclass.R")

# Package (Again) the code with the new DESCRIPTION file
current.code <- as.package("FitStatisticsPack")

# Check if the package content
check(current.code)

# Build the package
build("FitStatisticsPack", path = getwd())

# Install the package
install("FitStatisticsPack", reload = TRUE)


<file_sep>###############################################
## PS 5625 - Applied Statistical Programming
## Problem Set 5
## Author: <NAME>
## Tasks: 1 - Import dataset
##        2 - Clean and recode data
##        3 - Estimate the models and generate the predicted values.
##        4 - Generate a vector with naive predictions
##        5 - Calculate the Fit Statistics

rm(list = ls())

# Define WD
setwd("~/Google Drive/2. Academicos/6. Doutorado/Class/2017/Spring/PS 5625 - Applied Statistical Programming/Homeworks/Problem Sets/PS5/")

# Load the libraries
library(randomForest)
library(foreign)
library(AER)
library(FitStatisticsPack)

# Import dataset
anesTS2012 <- read.dta("anes_timeseries_2012_stata12.dta")

# Select the variables of interest
keep_var <- c("ft_dpc", "gender_respondent_x", "ecblame_pres", "econ_unpast", 
              "dem_raceeth_x", "libcpre_self", "dem_edu")

# Keep only the variables of interest
anesTS2012sub <- anesTS2012[,is.element(colnames(anesTS2012), keep_var)]

# Recode variables
# Female is recoded as 1
anesTS2012sub$female <- ifelse(anesTS2012sub$gender_respondent_x=="2. Female", 1, 0)

# Is president blamed by the economy?
# "A great deal", "a lot", and "A moderate amount" are recoded as 1.
anesTS2012sub$pres_blamed <- ifelse(anesTS2012$ecblame_pres=="-9. Refused" | 
                                        anesTS2012$ecblame_pres=="-8. Don't know", NA, 
                                             ifelse(anesTS2012$ecblame_pres=="4. A little" |
                                                       anesTS2012$ecblame_pres=="5. Not at all", 0, 1))

# It Unemployment better or worse than 1 year ago? 
# Worse and Same are recoded as 0
anesTS2012sub$unemployment <- ifelse(anesTS2012$econ_unpast=="-9. Refused" | 
                                        anesTS2012$econ_unpast=="-8. Don't know", NA, 
                                          ifelse(anesTS2012$econ_unpast=="1. Better", 1, 0))

# Race and ethnicity group
# "White and non-Hispanic" is recoded as 1
anesTS2012sub$white <- ifelse(anesTS2012sub$dem_raceeth_x == "-9. Refused" |
                                 anesTS2012$dem_raceeth_x == "-8. Don't know", NA, 
                                     ifelse(anesTS2012sub$dem_raceeth_x == "1. White non-Hispanic", 1, 0)) 


# 7pt scale Liberal/conservative self-placement
anesTS2012sub$lcscale <- ifelse(anesTS2012sub$libcpre_self == "-9. Refused" |
                                anesTS2012$libcpre_self == "-8. Don't know" |
                                anesTS2012$libcpre_self == "-2. Haven't thought much about this", NA, 
                                ifelse(anesTS2012sub$libcpre_self=="1. Extremely liberal", 1,
                                    ifelse(anesTS2012sub$libcpre_self=="2. Liberal", 2, 
                                       ifelse(anesTS2012sub$libcpre_self=="3. Slightly liberal", 3,
                                          ifelse(anesTS2012sub$libcpre_self=="4. Moderate; middle of the road", 4,
                                             ifelse(anesTS2012sub$libcpre_self=="5. Slightly conservative", 5,
                                                ifelse(anesTS2012sub$libcpre_self=="6. Conservative", 6, 7)))))))

# Education
# 95. Other {SPECIFY} is recoded as NA
anesTS2012sub$education <- ifelse(anesTS2012sub$dem_edu == "-9. Refused" |
                                     anesTS2012$dem_edu == "-8. Don't know" |
                                     anesTS2012$dem_edu == "95. Other {SPECIFY}", 
                                          NA, as.numeric(anesTS2012$dem_edu))
j <- 0
for(i in sort(unique(anesTS2012sub$education))){
   anesTS2012sub$education[anesTS2012sub$education==as.numeric(i)] <- j
   j <- j+1
}

# Obama Thermometer
anesTS2012sub$obamather <- ifelse(anesTS2012sub$ft_dpc>=0, anesTS2012sub$ft_dpc, NA)

# Remove the vectors j, i, and keep_var.
rm(j, i, keep_var)

# Split the sample
training_set <- row.names(anesTS2012sub) %in% sample(1:nrow(anesTS2012sub),
                                                     nrow(anesTS2012sub)/2, replace = FALSE)
# Estimate the models
# OLS
model_OLS <- lm(obamather ~ female + pres_blamed + unemployment + white + education
             + lcscale, data = anesTS2012sub[training_set,])
# Tobit
model_Tobit <- tobit(obamather ~ female + pres_blamed + unemployment + white + education
                     + lcscale, data=anesTS2012sub[training_set,], left = 0, right = 100)
# Random Forest
model_RM <- randomForest(obamather ~ female + pres_blamed + unemployment + white + education
                         + lcscale, data = na.omit(anesTS2012sub[training_set,]), trees = 1000)

# Generate the predicted values
# OLS
predOLS <- predict(model_OLS, newdata = anesTS2012sub[!training_set,], type="response")
# Tobit
predTobit <- predict(model_Tobit, newdata = anesTS2012sub[!training_set,], type="response")
# Random Trees
predRM <- predict(model_RM, newdata = anesTS2012sub[!training_set,], type="response")
# Add them into a matrix
pred_matrix <- matrix(c(predOLS, predTobit, predRM), ncol = 3)
# Add names to the matrix colums
colnames(pred_matrix) <- c("OLS", "Tobit", "RM")

# Generate the vector R with the naive predictions
R <- sample(1:100, nrow(anesTS2012sub[!training_set,]), replace = TRUE)

# Estimate the Fit Statistics 
# It's going to return Warning messages because r is not specified, and Tobit's
# predictions vector contains values less than zero.
FitStats(P = pred_matrix, y = anesTS2012sub[!training_set,"obamather"],
              statistics = c("rmse", "mad", "mape", "meape", "rmsle", "mrae"))
# It's going to return Warning messages because Tobit's 
# predictions vector contains values less than zero.
FitStats(P = pred_matrix, y = anesTS2012sub[!training_set,"obamather"],
              statistics = c("rmse", "mad", "mape", "meape", "rmsle"))
# It calculates only MRAE
FitStats(P = pred_matrix, y = anesTS2012sub[!training_set,"obamather"],
              r = R, statistics = c("mrae"))
# We can observe that the models have similar fit statistics. However,
# the random forest model is the one that has the lowest statistics, except by 
# MAPE (Tobit is the one with the lowest value for this statistic). Moreover,
# we also have that Tobit produces some negative predicted values and, therefore,
# the RMSLE is not possible to be calculated to this model.

<file_sep>library(testthat)
test_check("FitStatisticsPack")<file_sep>#' Calculate Fit Statistics
#'
#' @param y A vector with the observable values in the dataset.
#' @param P A matrix with the the predicted values.
#' @param r A vector with naive predictions (optional).
#' @param statistics A list of statistics to be calculated. 
#' @return A matrix with the selected statistics.
#' @author <NAME>
#' @details The user can specify the following statistics: RMSE, MAD, MAPE, MEAPE, RMSLE, and MRAE.
#'  If "mrae" is specified by the user, a vector r must be provided. Otherwise, MRAE. will be not calculated and
#'  a warning messenger will be shown (The remaining statistics will be calculated in any case).
#'  If P or y contain negative numbers, RMSLE will be not calculated and a warning messenger will be 
#'  shown.
#' @examples
#' x <- rpois(100, lambda = 100)
#' Y <- 4 + 3 * x + rnorm(100)
#' P <- matrix(predict(lm(Y ~ x)), ncol = 1)
#' r <- sample(min(Y):max(Y), 100, replace = TRUE)
#' FitStats(y = Y, P = P, r = r)
#' @rdname FitStats
#' @import stats
#' @export
FitStats <- function(y, P, r = NULL, 
                          statistics = c("rmse", "mad", "mape", "meape", "rmsle", "mrae")){
   # If y=0, it adds a small number to not generate Inf in the next 
   #  in the next steps.
   y[y==0 & !is.na(y)] <- y[y==0 & !is.na(y)]+0.01
   
   # Calculate the error term
   error <- apply(P, MARGIN = 2, function(x) abs(x-y))
   
   # Calculate the percentage of error
   perc_error <- (error/abs(y)) * 100
   
   # Select the observations in y that do not have NA in error
   yn <- y[complete.cases(error)]
   
   # Select the observation in P that do not have NA in error
   Pn <- matrix(P[complete.cases(error),], ncol = ncol(P))
   
   # Keep only the complete cases of error and perc_error (Non-NAs)
   e <- matrix(error[complete.cases(error),], ncol = ncol(error))
   a <- matrix(perc_error[complete.cases(perc_error),], ncol = ncol(perc_error))
   
   # Generate the number of non-NAs rows in e 
   n <- nrow(e)
   
   # Generate a matrix to store the results
   results <- matrix(1:ncol(Pn), nrow = ncol(Pn), ncol = 1)
   # Rename first matrix rows
   rownames(results) <- colnames(P)
   
   # Calculate the fit statistics
   # Calculate RMSE
   if("rmse" %in% statistics){
      rmse <- apply(e, MARGIN = 2, function(x) sqrt(sum(x^2)/n))
      results <- cbind(results, RMSE = round(rmse, 4))
   }
   # Calculate MAD
   if("mad" %in% statistics){
      mad <- apply(e, MARGIN = 2, median)
      results <- cbind(results, MAD = round(mad, 4))
   }
   # Calculate MAPE
   if("mape" %in% statistics){
      mape <- apply(a, MARGIN = 2, function(x) sum(x)/n)
      results <- cbind(results, MAPE = round(mape, 4))
   }
   # Calculate MEAPE
   if("meape" %in% statistics){
      meape <- apply(a, MARGIN = 2, median)
      results <- cbind(results, MEAPE = round(meape, 4))
   }
   # Calculate RMSLE
   if("rmsle" %in% statistics){
      rmsle <- apply(Pn, MARGIN = 2, function(x) sqrt(sum((log(x+1)-log(yn+1))^2)/n))
      results <- cbind(results, RMSLE = round(rmsle, 4))
      # Test if P or y contain negative numbers and generate a warning messeger
      neg_numbers <- apply(Pn, MARGIN = 2, function(x) (min(x)<0 | min(yn)<0))
      if(TRUE %in% neg_numbers){
         warning("NaNs introduced because P or x contain negative numbers")
      }
   }
   if("mrae" %in% statistics){
      if(is.null(r)){
         results[,-1]
         warning("r must be specified when mrae is selected.")
      }else if (!is.null(r)){
         # If r=0, it adds a small number to not generate Inf in the next 
         #  in the next steps.
         r[r==0 & !is.na(r)] <- r[r==0 & !is.na(r)]+0.01  
         # Calculate b
         b <- abs(r-y)
         # Keep only the observations that are not NA's in error:
         bn <- b[complete.cases(error)]
         # Calculate MRAE
         mrae <- apply(e, MARGIN = 2, function(x) median(x/bn))
         results <- cbind(results, MRAE = round(mrae, 4))
      }
   }
   results[,-1]
}
<file_sep>context("Checking if errors an error when P is not a matrix")
test_that("Error when P is not matrix", 
          expect_error(FitStats(y = seq(1, 100), 
                             P = seq(1, 100),
                             r = seq(3, 200, length.out = 100)))
)<file_sep>context("Checking object class")
test_that("Correct object is produced", 
          expect_is(FitStats(y = seq(1, 100), 
                             P = matrix(c(seq(1, 100), seq(2, 101), seq(2, 1010, length.out = 100)), ncol = 3),
                             r = seq(3, 200, length.out = 100)),  "matrix")
          )