---
title: "Workshop 06: Regression and Time Series"
author: "oran.k@chula.ac.th"
date: "August 2021"
output: html_document
params:
  intab: TRUE
  
---

```{r,setup,echo=F,include=FALSE}
  #rm(list=ls())
  knitr::opts_chunk$set(echo=TRUE,eval=TRUE,result="show",warning = FALSE,message=FALSE)

  knitr::opts_template$set(hidden = list(
    fig.width = 10, fig.height = 6,
    fig.retina = 2, out.width = '100%',
    echo=!params$intab, eval=!params$intab,
    warning=FALSE,message=FALSE
  ))
  
  knitr::opts_template$set(badcode = list(
    echo=TRUE, eval=FALSE,
    warning=FALSE,message=FALSE
  ))
  
  #intab <- TRUE
```

```{r klippy, echo=FALSE, include=TRUE}
klippy::klippy(position = c('top', 'right'),color = 'darkred')
```

```{css, echo=F}
  .badcodeStyle { background-color: lightpink; font-family: courier }
  
```

```{css, echo=F}
  .hiddenStyle { background-color: lightblue; font-family: courier }
```

# Instruction  # `r if (params$intab) "{.tabset}"`
Read the [instruction](https://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/workshop06.pdf) carefully and think about how to develop R code to answer each questions. 




## Question 1 ## `r if (params$intab) "{.tabset}"`
Consider **`AirPassengers`** dataset and explore the following data mining procedure for a time series data

---

### 00: background `r if (params$intab) "{.tabset}"`

The time series in R are encoded as 'ts' class (an extension of 'vector') and performed through many libraries:

* __forecast:__ A unified library to manage and forecast with time series
* __ggplot2:__ A extension to help with customizing `autoplot()` and `autolayer()` 
* __fpp2:__ a library of data and method related to [Forecasting: Principles and Practice](https://otexts.com/fpp2/) by R.J. Hyndman and G.Athanasopoulos, one of the best online textbook about forecasting in R. 


#### Info 

A basic information about data in 'ts' class can be derived using

```{r}
  tsp(AirPassengers)
```

The command reveals `start()`, `end()`, and `frequency()` of time series data. Notice that 'ts' data are converted into fraction number, particular with interval of `deltat()`

```{r}
  start(AirPassengers)
  end(AirPassengers)
  frequency(AirPassengers)
  
    deltat(AirPassengers)
```

#### View

Full description of 'ts' can be done by:

```{r}
  time(AirPassengers)
  cycle(AirPassengers)
```

### 01: visual `r if (params$intab) "{.tabset}"`
##### plot the decomposition of the data and separate data into training set (1948-1958) and testing set (1959)

#### Slice
The 'ts' data slice can be done with `window()`

```{r}
  air.train <-  window(AirPassengers,start=1948.0,end=1958.917)
  air.test  <-  window(AirPassengers,start=1959.0,end=1959.917)
```

#### Plot
The 'ts' data plot can be done with `forecast::tsdisplay()` with option of select type of lower.right with arg `plot.type`. Some even ignores type of plot, using ``

```{r}
 forecast::tsdisplay(air.train)

 forecast::tsdisplay(air.train, plot.type = "histogram")  
 
 forecast::autoplot(air.train)
```

##### Decompose 

decompose can be done using `decompose()` with "additive" or "multiplicative". Each term can be extracted by `forecast::trendcycle()` `forecast::seasadj()` and `forecast::remainder()`, respectively.

```{r}
  forecast::autoplot(decompose(AirPassengers,type="additive"))

  air.decom <- decompose(air.train,type="multiplicative")
  
  # forecast::trendcycle(air.decom)
  # forecast::seasadj(air.decom)
  # forecast::remainder(air.decom)
```


##### __Note__ `plot()` is still work, but not popular. This is equivalent to `stat::ts.plot()`. `forecast::autoplot()` and `forecast::autolayer()` are more popular.  

```{r}
 plot(AirPassengers)

 ts.plot(air.train)
```


### 02: smoothing `r if (params$intab) "{.tabset}"`
##### develop the following smoothing forecasting models and compute accuracy 

#### Outlier

`forecast` package can be used to verify outliers using Friedman’s super smoother with command `forecast::tsoutliers()` 

```{r}
  ### SOLUTION TO QUESTION 1B ### 

  require(forecast)

  tsoutliers(air.train)  
```  

#### Smoothing Group

`forecast` package embedded many exponential smoothing functions:

* __forecast::ses()__ Simple exponential smoothing model 
* __forecast::holt()__ Holt’s linear trend model 
* __forecast::hw()__ Holt-Winters’ model 

They are parts of `ets()` a general exponential smoothing state space model [ETS](https://otexts.com/fpp2/estimation-and-model-selection.html). 

```{r}
  air.ses   <- ses(air.train,h=12)
  air.holt  <- holt(air.train,h=12)
  air.hw    <- hw(air.train,h=12)
  
  autoplot(air.ses)
  
  ##-- accuracy    
  accuracy(air.test,air.hw$mean)  
```

#### ARIMA Group

```{r}
  air.arima <- Arima(air.train,c(0,1,1),seasonal = list(order=c(0,1,1),period=12))
  

  ##-- auto optimize  
  air.arima <- auto.arima(air.train) 

  autoplot(forecast(air.arima,h=12))
  checkresiduals(air.arima)
  
```

##### __NOTE__ `forecast::Arima` overwrites `stat::arima()`. This is a reason, why we cannot use arg 'h' in `Arima()` and use `tsplot()` or `autoplot()` without `forecast()`.

```{r}  
  ##-- possible forecasing model
  arima(air.train,c(0,1,1),seasonal = list(order=c(0,1,1),period=12))

  plot(air.arima)
  
  ## ar()
  ## ma()
```

#### Other Group 

* __forecast::naive()__ naive model (last constant value)
* __forecast::meanf()__ mean forecast model 
* __forecast::rwf()__ random walk model 
* __forecast::thetaf()__ 	Theta forecast method


```{r}  
  plot(naive(air.train,h=12))
  plot(meanf(air.train,h=12))
  plot(rwf(air.train,h=12))
  plot(thetaf(air.train,h=12))
  plot(stlf(air.train,h=12))
```

#### Autolayer

Many forecasting can visually combine using `autolayer()` and `autoplot()`

```{r}
require(ggplot2)
air.naive <- naive(air.train,h=12)

hSize <- 12
autoplot(air.train) + 
  autolayer(forecast(air.naive,h=hSize),series="naive",PI=F) +  
  autolayer(forecast(air.holt ,h=hSize),series="Holt",PI=F)  +
  autolayer(forecast(air.arima,h=hSize),series="ARIMA",PI=F) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(limits = c(100,550),breaks=seq(100,500,50), expand = c(0, 0)) +
  xlab("Year") + ylab("Passengers") + theme_bw()  
```

### 03: API
#### develop DOS/Java application for exponential series forecasting that takes any csv file as an input (**hint:** this require seperate R script and basic on Batch file)

* export some data for API

```{r}
  ## -------------------------------------------------------------------------
  ## export airPassenger Data 
  ## -------------------------------------------------------------------------


  write.csv(AirPassengers,file = "air.csv")
```

* here is a code inside [fGUI.R](http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/algorithm/fGUI.R) that can execute with DOS and [batch file](http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/algorithm/runRScript.bat)

```{r,opts.label="badcode",class.source="badcodeStyle"}

  ## set path for r library and apply
  cmdArgs = commandArgs(trailingOnly = TRUE)
  libPath <- cmdArgs[1]
  .libPaths(c(libPath)) 
  
  
  ## load library (require .libPaths)
  
  require(tcltk)
  require(fgui) ## simple java script package
  require(forecast)
  
  
  ## create sub-function
  expoSeries  <-  function(period=12,opt="additive",alpha=0.5,
                           beta=0.05,csvFilename=NULL,predict=12)  {
    data <- read.csv( guiGetValue("csvFilename") ) [[2]]
    data.ts <- ts(data,start=1949,frequency=period)
    data.ht <- hw(data.ts,alpha=0.5,beta=0.01,seasonal=opt,h=predict)
    windows()
    plot(data.ht)
  }
  
  
  ## run main GUI
  gui(expoSeries,
       title="GUI: Exponental Smoothing Series",
       argSlider=list(alpha=c(0,1,0.1),
                      beta=c(0,.1,0.001)),
       argFilename = list(csvFilename=NULL),
       argOption = list(opt=c("additive","multiplicative"))
  )
```


Alternatively, you can use service of 'shinyapps.io' with [shinny](https://orankitti.shinyapps.io/shiny/) using the following source code [src](http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/algorithm/app.R) 

## Question 2 ## `r if (params$intab) "{.tabset}"`

Consider self-report weight/height and actual report of such measures of 200 participants

---

### 01: import
#### import/check for NA value separate into groups based on availability of data 

```{r}
  ### SOLUTION TO QUESTION 2A ### 

  fileURL <- "http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/data/selfReport.csv"

  self <- read.csv(url(fileURL),header = T)
  head(self) ## View(self)
```

* check for 'NA'
```{r}
  naSet<- as.data.frame(which(is.na(self),arr.ind = T))
  naSmy<- xtabs(~col+row,data=naSet)
  
  noRepWt<- as.numeric(names(naSmy[1,naSmy[1,]==0]))
  noRepHt<- as.numeric(names(naSmy[2,naSmy[2,]==0]))  
  noBoth <- as.numeric(names(naSmy[1,naSmy[1,]==1 & naSmy[2,]==1]))
  
  noNA   <- setdiff(1:nrow(self),c(noRepWt,noRepHt,noBoth))
```

* Alternatively, `complete.cases()` can be applied to eliminate all 'NA' 
```{r}
 self.nona <- self[complete.cases(self),]
```


### 02: separate 

* for participants who report both weight and height, separate data into testing data (15%) and training data (85%)

```{r }
  ### SOLUTION TO QUESTION 2B ### 

  set.seed(17)
  testingSet<- sample(noNA,size=round(length(noNA)*.15))
  self.test <- self[testingSet,]
  self.train<- self[setdiff(noNA,testingSet),] ### why not "self[-testingSet,]"? 

```

### 03: compare

* For training set, check difference between report and actual

```{r }
  ### SOLUTION TO QUESTION 2C ### 

  attach(self.train)
  self.train$gapHt <- height - repht
  self.train$gapWt <- weight - repwt  
  detach(self.train)
```

Does these commands work?


```{r}
  xtabs(~sex+gapHt,data=self.train) ## do you notice any error?  how to correct it?
  xtabs(~sex+gapWt,data=self.train) ## do you notice any error?  how to correct it?
```



```{r}
  head(self.train) ## fix(self.train)
  temp <- xtabs(~gapHt+gapWt,data=self.train)
  
  temp.mat <- as.matrix(temp)
  filled.contour(z=temp.mat,
                 x=as.numeric(rownames(temp.mat)),
                 y=as.numeric(colnames(temp.mat)),
                 nlevels=10)
  
```

### 04: predict 

* For training set, build a linear regression model to predict actual weight with actual height

```{r }
  ### SOLUTION TO QUESTION 2D ### 

  weight.lm <- lm(weight~height+sex,data=self.train)
  anova(weight.lm)
  confint(weight.lm)
  plot(weight.lm)  
  
  weight.predict <- predict(weight.lm,newdata=self.test)
  
  forecast::accuracy(weight.predict,self.test$weight)
  
  ##-- advance construction of lm
  
  ##-- combination of variable, i.e., lm(weight~height+sex+height:sex,data=self.train)
  # lm(weight~height*sex,data=self.train)

  ##-- interaction term   
  # lm(weight~height:sex,data=self.train)  
  
  ##-- no constant term
  # lm(weight~height:sex -1,data=self.train)  
  
  ##-- apply function 
  # lm(weight~I(sqrt(height)),data=self.train)
  
  ##-- all main term
  # lm(weight~.,data=self.train)
```

### 05: Assumption

* __L__: Linear relationship using `cor` or visual inspection  
* __I__: Independence or at least no AR(1) using `car::durbinWatsonText()` or `ols_vif_tol()` (VIF<5.0)
* __N__: Normality using `ols_test_normality()`
* __E__: weak Exogeneity ($x_j$ is not random variable). In Engineering, this assumption is usually ignored. 
* __S__: homoScedasticity (single variance) using `ols_plot_resid_fit()`

The main package to identify these assumption is [olsrr](https://olsrr.rsquaredacademy.com/articles/intro.html)

```{r}
  require(olsrr)

  # ols_plot_resid_hist(weight.lm)
  # ols_plot_resid_fit(weight.lm) + theme_bw()
  # ols_plot_cooksd_chart(weight.lm)  ## check for major influence points
  # ols_plot_added_variable(weight.lm)
  # ols_plot_response(weight.lm)
  
  ols_vif_tol(weight.lm)
  ols_test_normality(weight.lm)
  ols_test_correlation(weight.lm)
  ols_test_outlier(weight.lm)
```

```{r,opts.label="hidden",class.source="hiddenStyle"}
  ##--- solution 

  fileURL <- "http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/data/selfReport_clean.csv"

  require(data.table)
  self.DT <- fread(fileURL)
  
  self.DT[,sex:=factor(sex)]
  self.DT[,gapHt:= height-repht]
  self.DT[,gapWt:= weight-repwt]
  self.DT[,id:=1:.N]  
  
  self.DT <- self.DT[complete.cases(self.DT)]
  self.DT[,.N,by=sex]
  
  set.seed(17)
  strat.idx  <- self.DT[,sample(.I,10),by=sex]$V1
  self.train <- self.DT[id%in%strat.idx]
  self.test  <- self.DT[id%in%setdiff(self.DT$id,strat.idx)]
  
  set.seed(17)
  randm.idx  <- sample(self.DT$id,size=round(0.85*nrow(self.DT)))
  self.train <- self.DT[id%in%randm.idx]
  self.test  <- self.DT[id%in%setdiff(self.DT$id,randm.idx)]
  
  
  reg.model <- list()
  reg.model$lm <- lm(weight~height,data=self.train)
  reg.model$opt<- step(lm(weight~. -id-gapWt,data=self.train))
  reg.model$nav<- lm(weight~sex+height-1,data=self.train)
  
  require(olsrr)
  ols_step_all_possible(reg.model[[2]])
```


## Question 3 ## `r if (params$intab) "{.tabset}"`

Apply your knowledge of logistic regression into the following dataset

---

### 01: import

* observe pattern 

```{r }
  ### SOLUTION TO QUESTION 3A ### 

  fileURL <- "http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/data/height.csv"

  height <- read.csv(url(fileURL),header=T)
  height$Gender <- factor(height$Gender,levels = c(1,2)
                         ,labels = c("F","M"))  
  xtabs(~Gender+Height,data=height)
```

### 02: construct  

```{r}
  height.glm <- glm(Gender~Height,data=height,family = "binomial")
  
  
  prob <- data.frame(Gender=factor(rep("M",21)),Height=60:80)
  prob$prob<- predict(height.glm,newdata= prob,type = "response")
  plot(prob$Height,prob$prob,col="black")
```

* What is the cut-off threshold? 
* Can we apply other methods?



### 99: EXTRA

* Consider the previous dataset **`Davis`** or **`self.train`**, construct a model to predict **sex** of participants using other factors.

```{r,opts.label="badcode",class.source="badcodeStyle" }
  ### SOLUTION TO QUESTION 3B ### 

  fileURL <- "http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/data/selfReport.csv"

  self <- read.csv(file=fileURL,header = T)
  self$sex   <- factor(self$sex)
  self$gapHt <- self$height - self$repht
  self$gapWt <- self$weight - self$repwt 
  

  noNA <- setdiff(1:200,unique(which(is.na(self),arr.ind = T)[,1]))
  
  set.seed(17)
  testingSet<- sample(noNA,size=round(length(noNA)*.15))  
  self.test <- self[testingSet,]
  self.train<- self[setdiff(noNA,testingSet),]  ### why not "self[-testingSet,]"? 
  
  sex.log <- glm(sex~gapWt+gapHt,data=self.train,family="binomial")
  sex.best<- step(glm(sex~.,data=self.train,family="binomial"))
  
  anova(sex.log)
  anova(sex.best)
  
  self.test$probSex<- predict(sex.log,newdata=self.test,type = "response")
  self.test$bestSex<- predict(sex.best,newdata=self.test,type = "response")
  self$sex
```

## Quesiton 4 ## `r if (params$intab) "{.tabset}"`

Consider the [award](http://www.ats.ucla.edu/stat/data/poisson_sim.csv) dataset of difference student groups with different Mathematics score and follow the data mining steps


### 01: import

* import and analysis the data

```{r, message=FALSE}
  ### SOLUTION TO QUESTION 4A ### 
  fileURL <- "http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/data/poisson_sim.csv"

  require(data.table)
  award.DT <- as.data.table(read.csv(url(fileURL),header=T))
  award.DT[,prog:= factor(prog, levels=1:3, labels=c("General", "Academic", "Vocational"))]
  award.DT[,num_awards:=factor(num_awards)]
    
  require(ggplot2)
  require(dplyr)
  award.DT %>% ggplot( aes(x=num_awards,y=math,fill=prog) ) + 
    geom_violin() + facet_grid(~prog) 
```


```{r}

  award.DT[,.(min=min(math),avg=mean(math),max=max(math),sd=sd(math),.N),by=.(prog,num_awards)] -> temp.DT


  
  xtabs(avg~prog+num_awards,temp.DT)
  xtabs(N~prog+num_awards,temp.DT)
  
```

### 02: glm()
* use Poisson Regression to predict response


```{r}
  award.DT[,count_awards:=as.numeric(num_awards)]
  award.pois <- glm(count_awards ~ prog + math, family="poisson", data=award.DT)
  summary(award.pois)
  
  head(predict(award.pois,type="response"))
  
  award.grid <- expand.grid( math=seq(32,70,2)
                             ,prog=factor(1:3, labels=c("General", "Academic", "Vocational"))
                             ,num_awards=NA)
  
  award.grid$pois_awards <- predict(award.pois,newdata = award.grid,type="response")

  xtabs(pois_awards~prog+math,award.grid)
 
  ##-- 
  # plot(jitter(award.DT$math),jitter(award.DT$count_awards),col=award.DT$prog,pch=16,cex=0.5,xlab="Math",ylab="# awards")   
  
```

What is your conclusion?

### 03: compare

* compare the Poisson regression with simple linear model and calculate AIC of each model

```{r}
  ### SOLUTION TO QUESTION 4C ### 

  award.lm   <- lm(count_awards ~ prog + math, data=award.DT)
  summary(award.lm)

  award.grid$lm_awards <- predict(award.lm, newdata=award.grid,type="response")
  xtabs(lm_awards~prog+math,award.grid)  
  xtabs(pois_awards~prog+math,award.grid)  
  
  AIC(award.pois)
  AIC(award.lm)
```

## Question 5 ## `r if (params$intab) "{.tabset}"`

Consider the dataset [Professor Salaries](http://ie.eng.chula.ac.th/~oran/classes/CompMeth_resources/scripts/Salaries) and follow these instruction.

---

### 01: import
#### import data from MS Excel file with package **`xlsx`**


```{r,eval=FALSE}
  ### SOLUTION TO QUESTION 5A ### 
 
  require(openxlsx)
  fileURL <- "http://ie.eng.chula.ac.th/~oran/classes/compMeth_resources/workshop/scripts/data/Salaries.xlsx" 

  ##-- import sheet by name
  prof <- read.xlsx(xlsxFile=fileURL,sheet = "Salaries") 
  
  ##-- import sheet by index
  prof <- read.xlsx(xlsxFile=fileURL, sheet=1)   
  glimpse(prof) 
    
  ##-- import specific row/col of sheet index
  wb       <- loadWorkbook("Salaries.xlsx") ##-- where is your file?  
  sheets   <- getSheets(wb)
  sheet    <- sheets[[1]] ##-- access to first sheet
  nLastRow <- sheet$getLastRowNum()
  prof     <- as.data.frame(readRows(sheet,startRow=2,endRow=nLastRow,startColumn = 1)) 
  colnames(prof) <- readRows(sheet,startRow=1,endRow=1,startColumn = 1)                          
  glimpse(prof)  ##-- this leads to ERROR; can you explain WHY? 
  
  
  require(dplyr)
  require(lattice)
  ##-- convert data into appropiate type
  prof <- mutate(prof,rank=factor(rank,labels=c("AsstProf","AssocProf","Prof")),
                 discipline=factor(discipline),
                 yrs.since.phd=as.numeric(yrs.since.phd) ,
                 yrs.service=as.numeric(yrs.service) ,
                 sex=factor(sex,labels=c("M","F")),
                 salary=as.numeric(salary) ) 
  
  ##-- play with data
  alys01 <- filter(prof,discipline=="B") %>% group_by(rank,sex) %>% summarise(avg=mean(salary),n=length(salary))
  xtabs(avg~sex+rank,data=alys01)  
```

### 02:  observe  
#### observe the pattern of `**salary**` and relationship of each column 

```{r, eval=FALSE} 
  ### SOLUTION TO QUESTION 5B ### 

  glimpse(prof) ## data contain quantitative (number) // qualitative (description)
  
  ##-- analysis of qualitative data (description)
  bwplot(salary~rank|discipline:sex,data=prof) ##-- is there any diff in terms of rank, discipline, and sex 
  anova(lm(salary~rank*sex*discipline,data=prof))

  ##-- analysis of quantitative (number) data
  cor(prof[,c(3,4,6)])
  anova(lm(salary~yrs.since.phd+yrs.service,data=prof))
  pairs(prof[,c(3,4,6)],col=prof$discipline)
  
  ##-- analysis of both  
  xyplot(salary~yrs.since.phd|(discipline:rank),data=prof)
  xyplot(salary~yrs.service|(discipline:rank),data=prof)  
  prof.lm  <- step( lm(salary~.,data=prof)) ## find the best model
  anova(prof.lm)
```  

### 02: `hist()`
#### visualize histogram of salary base on main factors

```{r, eval=FALSE} 
  ### SOLUTION TO QUESTION 5C ### 

  prof.hist <- hist(prof$salary,levels=20,plot = F)
  prof.B    <- filter(prof,discipline=="B")
  prof.Asst <- filter(prof,rank=="AsstProf")
  prof.nFull<- filter(prof,rank!="Prof")
  
  hist(prof$salary,breaks=prof.hist$breaks,col="blue")
  hist(prof.B$salary,breaks=prof.hist$breaks,col="red",add=T)
  rug(prof$salary,col="blue")
  rug(prof.B$salary,col="red")  
  
  hist(prof$salary,breaks=prof.hist$breaks,col="blue")
  hist(prof.nFull$salary,breaks=prof.hist$breaks,col="red",add=T)
  hist(prof.Asst$salary,breaks=prof.hist$breaks,col="orange",add=T)  
  rug(prof$salary,col="blue")
  rug(prof.nFull$salary,col="red")
  rug(prof.Asst$salary,col="orange")      

```

### 03: OLR
#### predict salary with ordinary linear regression model

```{r, eval=FALSE}       
  ### SOLUTION TO QUESTION 5D ### 

  prof.ctree <- ctree(salary~.,data=prof)
  plot(prof.ctree)
  summary(prof.ctree)
  
```

### 04: DT
#### use decision tree (**`ctree`** in **`party`** package or **`rpart`** in **`rpart`** package) to analysis data 

```{r, eval=FALSE}       
  ### SOLUTION TO QUESTION 5E ### 

  prof.ctree <- ctree(salary~.,data=prof)
  plot(prof.ctree)
  summary(prof.ctree)
  
```



<!----------- BEGIN: STANDARD OKC FOOTER  -------------------------------->
#

---

Copyright 2019 &ensp; Oran Kittithreerapronchai.  &emsp; All Rights Reserved.
&emsp; __Last modified:__ `r format(Sys.time(), '%Y-%M-%d, ')` 
<!----------- END: STANDARD OKC FOOTER  ---------------------------------->



<!-------- BEGIN: STANDARD OKC PRINT OUT with SYNTEX  -------------------->
```{r, echo=F,eval=F,include=F} 

##--- CODE PRINING

curPath <- dirname(rstudioapi::getSourceEditorContext()$path)
fileName<- basename(rstudioapi::getSourceEditorContext()$path)
outName <- paste( c(fileName,".html"),collapse="")

setwd(curPath)
rmarkdown::render(fileName,params=list(intab=FALSE),output_file = outName)
# library(knitr)
# stitch(
#   script = fileName,
#   system.file("misc", "knitr-template.Rhtml", package = "knitr"),
#   output = outName
# )

#---
```
<!--------- ENDED: STANDARD OKC PRINT OUT with SYNTEX  ------------------->