Final Update for V2

Modif.R

@@ -1,26 +1,45 @@
-Modif <- function(pred, data, co.r = NaN, pc.r = NaN, sc.1 = NaN, sc.2 = NaN){
-  # transform specified variable with parameters
-  # this function **return a dataset** with certain variable transformed
+modif <- function(pred, type, data, co.r=NULL, sc.1=NULL, sc.2=NULL, pc.r=NULL, object=NULL){
+  #---------------------------------
+  # transform specified variable with(out) parameters
+  # this function returns **a dataset** with selected variable transformed
   # make sure the raw data is already loaded into global environment
+  #---------------------------------
+  # pred: the name of variable to be transormed of class "character"
+  # type: transformation method:
+  #       1.1~1.2 media: carryover + s-curve / carryover + power curve
+  #       2.1.1~2.1.4 basic: +,-,*,/; corresponding to parameter "object"
+  #       2.2 basic: logarithm
+  #       2.3 basic: root
+  #       2.4 basic: exponent
+  #       2.5 basic: reciprocal
+  #       2.6 basic: time lag
+  #       the input for type should be strictly controlled to the list of option number above!!!
+  # data: data frame to store the variable transformed (original var. will be covered)
+  #---------------------------------
 
-  # source('BasicTrans.R')
-  df <- data
-
-  # check if the variable exists in the data frame
-  check <- pred %in% names(df)
-  if (check == FALSE){
-    stop(paste("The variable '", pred, "' does not exist in the database!", sep=""))
+  trans.collection <- c("co","sc","pc","bt")
+  if(!sum(sapply(trans.collection, existsFunction))==4){
+    # check if the functions co(), sc(), pc() & bt() exist or not, if not, source "odds.R"
+
+    if(file.exists("odds.R")){
+      source("odds.R")
+    }else{
+      source("https://raw.githubusercontent.com/elliott828/boulot-test/master/odds.R")
+    }
   }
-  pred <- as.character(pred) # ensure input like factor to be coerced to character
 
-  if(is.na(co.r)){
-    # if no parameter for co rate, then no transformation happens
-    message(paste("No transformation for variable '", pred, "'!", sep=""))
-  }else if(is.na(sc.1)){
-    # if sc parameters is NA, then do co+pc transformation
-    df[[pred]] <- pc(co(df[[pred]], co.r), pc.r)
+  df <- data
+  x <- data[[pred]]
+
+  if(type == 1.1){
+    df[[pred]] <- cs(data[[pred]], co.r, sc.1, sc.2)
+  }else if(type == 1.2){
+    df[[pred]] <- cp(data[[pred]], co.r, pc.r)
   }else{
-    df[[pred]] <- sc(co(df[[pred]], co.r), sc.1, sc.2)
+    opt <- substr(as.character(type),3,nchar(as.character(type)))
+    # turn 2.1.1~2.1.4 and 2.2~2.6 to 1.1~1.4 and 2~6
+    df[[pred]] <- bt(data[[pred]], opt, object)
   }
+
   return(df)
-}
+}

final.output.R

@@ -0,0 +1,53 @@
+final.output <- fucntion(resp, data, fit, prmt, contri) {    		
+    # prmt is a dataframe including 8 columns: 
+    #       variable, trans.meth, co.r, sc.1, sc.2, pc.r, oth, status
+
+    # Output: 
+
+
+    # Needed Packages: "MASS" (for stepAIC()), "car" (for vif())             
+
+    # 1.Transformation (trans.meth, co.r, sc.1, sc.2, pc.r, oth)
+    write.csv(prmt, paste(getwd(), "/prmt.history.csv", sep = ""))
+    message('Variable parameters history is exported to "prmt.csv" under default working directory')
+    cat(paste(rep("-+-",20),collapse=""))
+
+    # 2. Residuals
+    resid <- cbind(data[, resp], fit$fitted.values, summary(fit)$residuals)
+    rownames(resid) <- NULL
+    colnames(resid) <- c(resp, "Prediction", "Residuals")
+
+    write.csv(resid, paste(getwd(), "/residuals.csv", sep = ""))
+    message('Value of Response Variable, Prediction and Residuals
+            is exported to "residuals.csv" under default working directory')
+    cat(paste(rep("-+-",20),collapse=""))
+
+    # 3. Coefficients
+    coef <- coef(summary(fit))     # Estimate  Std. Error    t value   Pr(>|t|)
+
+    # 4. VIF
+    if(length(coef(fit)) > 1) {
+        vif <- vif(fit)
+    }
+
+    # Merge model information into one data frame
+    # Output csv
+    prmt.alive <- prmt[which(prmt$status == "alive"), ][, -8]
+
+    if(rownames(coef)[1] == "(Intercept)"){
+        prmt.alive <- rbind(NA, prmt.alive)
+        vif <- rbind(NA, as.matrix(vif))
+    }
+    model <- as.data.frame(cbind(prmt.alive, coef, contri, vif))
+    rownames(model) <- NULL
+    colnames(model)[, -(1:11)] <- c("contri.rate", "VIF")
+
+    write.csv(model, paste(getwd(), "/model.results.csv", sep = "")))
+    message('Information of the Model 
+            is exported to "model.results.csv" under default working directory')
+    cat(paste(rep("-+-",20),collapse=""))
+
+    # print(summary(stepAIC(fit, direction = "both", trace = 0)))        
+
+}
+

loop.output.R

@@ -0,0 +1,89 @@
+loop.output <- function(resp, data, fit) {    					
+    # Input: resp(name), data(modified), fit					
+    # Output: only on memory and screen, nothing to files
+
+    # Consists of 4 parts:  Part. I  Summary of Fit & MAPE
+    #                       Part. II  Plots
+    #                       Part. III  DW-test
+    #                       Part. IV  Contribution Rates
+
+    # Needed Packages: "lmtest" (for dwtest()), "zoo" (for library(lmtest))
+
+
+    #---------------------------------------------------------------
+    # Part. I  Summary of Fit & MAPE
+    #---------------------------------------------------------------
+
+    cat("\n")
+    readline("Part. I  Summary of Fit & MAPE")
+    print(summary(fit))	
+    mape <<- mean(abs(as.vector(fit$residuals))/data[, resp])
+    cat("MAPE of the model is ", mape, "\n")
+    cat("\n")
+
+
+    #---------------------------------------------------------------
+    # Part. II  Plots
+    #---------------------------------------------------------------
+
+    if(!is.null(dev.list())) invisible(dev.off()) # Clear the visible equipment
+    na <- readline("Part. II  Plots")
+    message("Please look at the Plots area!")
+    layout(matrix(c(1, 2,
+                    3, 3), nr = 2, byrow = T))
+
+    # 1. histogram for residuals
+    hist(summary(fit)$residuals, main = "Histogram of Residuals", xlab = "Residuals")
+
+    # 2. scatter points of residuals
+    plot(summary(fit)$residuals, type = "p", pch = 21, bg = "black", 
+        main = "Scatter of Residuals", ylab = "Residuals")
+    abline(h = 0)
+    par(new = FALSE)
+
+    # 3. Actual Data vs. Modeled Data
+    plot(data[, resp], type = "l", col = "blue", xlab = "", ylab = "",  axes = FALSE)
+
+    # legend(2, 5, c("Actual", "Predicted"), fill = c("blue", "red"))
+    # legend("topright",  c("Actual", "Predicted"), fill = c("blue", "red"), 
+    #        border = "white", box.col = "white", box.lty = NULL)
+    par(new = TRUE)
+
+    # prediction <- (as.matrix(data[, names(coef(fit))[-1]]) 
+    #    %*% as.vector(coef(fit)[-1]) + (coef(fit)[1]))
+    plot(as.vector(fit$fitted.values), type = "l", col = "red", 
+        xlab = "", ylab = "", main = "Actual Data vs. Predicted Data")
+    cat("\n")
+
+
+    #---------------------------------------------------------------
+    # Part. III  DW-test
+    #---------------------------------------------------------------
+
+    na <- readline("Part. III  DW-test")
+    print(dwtest(fit))
+    cat("\n")
+
+
+    #---------------------------------------------------------------
+    # Part. IV  Contribution Rates
+    #---------------------------------------------------------------    
+
+    na <- readline("Part. IV  Contribution Rates")
+
+    simulation <- cbind(coef(fit)[1], as.matrix(data[, names(coef(fit))[-1]]) * as.vector(coef(fit)[-1]))
+    colnames(simulation)[1] <- "(Intercept)"
+    contri <<- sapply(as.data.frame(simulation), sum)/sum(fit$fitted.values)
+    contri <- as.matrix(contri[order(contri, decreasing = T)])
+
+    contri.top10 <- as.matrix(head(contri[ - which(rownames(contri)=="(Intercept)")], 10))
+    rownames(contri.top10) <- rownames(contri)[ - which(rownames(contri)=="(Intercept)")]
+    contri.top10
+
+    # percent <- function(x, digits = 4, format = "f") {
+    #     paste0(formatC(100 * x, format = format, digits = digits), "%")
+    # }
+    # contri.top10 <- as.matrix(percent(head(contri.rate, 10)))
+    # rownames(contri.top10) <- rownames(contri.rate)[1:10]
+
+}