R: Basic Codes

SAS has been the leader in statistical softwares in corporates and industries for a long time, while SPSS and STATA have been on top in academics. However, with the huge influx of start-ups, R and Python are coming up quickly to the top, and have challenged SAS for the leading position.

In this blog, I will share a basic set of R commands and codes that will be helpful to start working with R. These are easily available online, and this blog is just a small step to consolidate the important codes in one place. If someone really wants to learn R, it is recommended to browse through coursera and edx.

In this blog, we will go through the following commands:

1. Reading Data in R
2. List of columns in a dataset
3. Different joins in R
4. If conditions and multiple if conditions
5. Appending two dataset
6. Group by command
7. Finding distinct values in a vector

# 1. Reading Data in R

SampleData <- read.csv('C:/Users/Desktop/Sample_R_Data.csv',header = TRUE,dec=".")

SampleData1 <- read.table('C:/Users/Desktop/Sample_R_Data.csv',sep=",",header = TRUE,dec=".")

# 2. List of columns in the dataset

colnames(SampleData)

colist <- colnames(SampleData)

data1 <- subset(SampleData, select=c("cust_i", "y_var"))

data2 <- SampleData[,c(2,3,4)]

data2_1 <- SampleData[c("cust_i", "y_var")]

data3 <- SampleData[1:10,c(2,3,4)]

# 3. Joins in R

 data_LeftJoin <- merge(x = data1, y = data2, by = "cust_i", all.x = TRUE)

data_RightJOin <- merge(x = data1, y = data2, by = "cust_i", all.y = TRUE)

data_OuterJoin <- merge(x = data1, y = data2, by = "cust_i", all = TRUE)

data_CrossJoin <- merge(x = data1, y = data2, by = NULL)

data_InnerJoin <- merge(x = data1, y = data3, by = "cust_i")

# 4a. If conditions

data4 <- subset(SampleData, y_var==1)

data5 <- SampleData[SampleData$y_var==0,]

# install.packages("sqldf")

library(sqldf)

data6 <- sqldf("select * from SampleData where y_var==1")

# 4b. Multiple If condition (SUBSET is inefficient)

data7 <- subset(SampleData, y_var==1 & SampleData <= 300)

data8 <- SampleData[SampleData$y_var==0 & SampleData$str_trnx_gap <= 300,]

data9 <- SampleData[SampleData$y_var==0 | SampleData$str_trnx_gap <= 300,1:4]

# 5. Append

data10 <- rbind(data4,data5)

# 6. Group By Summary Statistics in R

# Using TABLE or AGGREGATE or TAPPLY

# Note the use of USER_DEFINED function

table(SampleData$y_var, responseName= "sum")

data.frame(table(SampleData$y_var))

aggregate(SampleData$cust_i, by=list(Category=SampleData$y_var), FUN=function(x){NROW(x)})

tapply(SampleData$y_var, SampleData$y_var, FUN=function(x){NROW(x)})

tapply(SampleData$y_var, SampleData$y_var, FUN=function(x){NROW(x)/NROW(SampleData)})

# 7. Get distinct values of a vector

unique(SampleData$y_var)

Outlier Treatment in SAS

A very important step for any type of analysis is the outlier treatment. An outlier is an observation which lies at a distant from the majority of the observations, may be because of some exceptional cases, or due to issues in data storage and management. For a regression model to be robust, it is essential to have a sanity check of the data, to remove the presence of any such anomaly.

A very basic way to remove the ouliers are to delete such extreme observations completely from the data. Another way to handle it is called capping and flooring, where the upper extreme values and the lower extreme values are replaced by a pre-determined threshold. Another common way to do this is to implement both of these together.

When the number of variables are huge, it might not be possible to check the distribution of every variable one by one, instead a rule can be created to treat the outliers. Let us create the following rule:

For all variables, we will delete the values which are above the 99 percentile or below the 1 percentile. After that, we would cap all the values that lie above 95 percentile to the 95th percentile value, and floor all the values that are less than 5 percentile to the 5th percentile value.

The following SAS code has been automated to implement the above rule and create a outlier-treated dataset. The rule can be modified as per scenario, or as per the analysts' discretion.

libname dataloc "C:/Desktop/Data";  /* MODEL DATASET LOCATION  */

%let inset=Base_Data;               /* MODEL DATASET NAME      */

%let libout=C:/Desktop/Output;

libname outlib "&libout.";          /* OUTPUT LOCATION         */

%let outset=Clean_Data;             /* OUTPUT DATASET          */

%let upper_del_threshold = 99;      /* UPPER LIMIT TO DELETE   */

%let lower_del_threshold = 1;       /* LOWER LIMIT TO DELETE   */

%let upper_cap_threshold = 95;      /* UPPER LIMIT TO CAP      */

%let lower_floor_threshold = 5;     /* LOWER LIMIT TO FLOOR    */

%let varlist = X1 X2 X3;            /* LIST OF VARIABLES       */

data inset1;

set dataloc.&inset.;

run;

proc means data=inset1 StackODSOutput P&upper_del_threshold. P&lower_del_threshold. P&upper_cap_threshold. P&lower_floor_threshold.;

var &varlist;

ods output summary=LongPctls;

run;

data LongPctls;

set LongPctls;

run;

data _null_;

set LongPctls;

call symput(compress("var"||_n_),compress(variable));

call symput(compress("up_del"||_n_),compress(P&upper_del_threshold.));

call symput(compress("low_del"||_n_),compress(P&lower_del_threshold.));

call symput(compress("up_ext"||_n_),compress(P&upper_cap_threshold.));

call symput(compress("low_ext"||_n_),compress(P&lower_floor_threshold.));

call symput(compress("n"),compress(_n_));

run;

%macro dataset;

data outset1;

set inset1;

%do i = 1 %to &n.;

if &&var&i..  gt &&up_del&i.. then delete;

else if &&var&i..  lt &&low_del&i.. then delete;

else if &&var&i..  gt &&up_ext&i.. then &&var&i.. = &&up_ext&i..;

else if &&var&i..  lt &&low_ext&i.. then &&var&i.. = &&low_ext&i..;

%end;

run;

%mend;

%dataset;

data outlib.&outset.;

set outset1;

run;

/*A quick check in the univariate analysis to validate the outlier treatment*/

%macro chck;

%do i = 1 %to &n.;

proc univariate data=inset1;

var &&var&i..;

run;

proc univariate data=outset1;

var &&var&i..;

run;

%end;

%mend;

%chck;

Logistic Regression in SAS

Logistic Regression is one of the most used technique in the analytics world, and for every propensity modelling, risk modelling etc., this is one of the most important as well as well-accepted steps. The main difference between the logistic regression and the linear regression is that the Dependent variable (or the Y variable) is a continuous variable in linear regression, but is a dichotomous or categorical variable in a logistic regression. You can read more about logistic regression here or the wiki page.

While validating a logistic model, we try to see some of the statistics like Concordance and Discordance, Sensitivity and Specificity, Precision and Recall, Area under the ROC curve. One of the most important aspect is the Precision and Recall. The problem faced by the analysts is how to balance between the two. If we try to increase one of them, the other reduces. A measure that tries to balance both of them is called as F-ratio or F-measure, which is calculated as the Harmonic Mean of precision and recall. However, based on the project requirement, one can calculate an adjusted F-ratio, which calculates the harmonic mean of precision and recall after giving a higher weight on one of them.

The following SAS code is an attempt to simplify the SAS code, and it has been automated for future use. A detailed documentation about the Logistic regression output is given here. The various outputs like parameter estimate, concordance-discordance, classification table etc. will be stored as tables. The html output contains the regular stuffs, along with the ROC curve for the training data as well as the ROC curve of the validation data. The score statement gives the classification table for the validation data, and also scores the validation data which can be used for calculating other validation statistics like Kolmogorov-Smirnov etc.

%let train = TRAINING_DATA;           /* TRAINING DATA */

%let validate = VALIDATION_DATA       /* VALIDATION DATA    */

%let targetvar = Y;                   /* DEPENDENT BINARY VARIABLE */

%let varlist = X1 X2 X3 X4;           /* INDEPENDENT VARIABLES   */

%let binvarlist = X2 X3;              /* LIST OF BINARY INDEPENDENT VARIABLES */

ods graphics on;

ods html;

ods output           

parameterestimates = TBL_ParamEst     /* PARAMETER ESTIMATES     */

OddsRatios =TBL_OddRatio              /* ODD RATIOS              */

LackFitPartition = TBL_HLpartition    /* HOSMER-LEMESHOW PARTITIONS    */

LackFitChiSq = TBL_HLstatistic        /* HOSMER-LEMESHOW STATISTIC     */

Association = TBL_Association         /* CONCORDANCE DISCORDANCE ETC   */

FitStatistics = TBL_FitStatistic      /* AIC, -2LOG, SC, ETC           */

GlobalTests = TBL_GlobalTests         /* WALD, LOGLIKELIHOOD, ETC      */

Classification = TBL_Classification   /* CLASSIFICATION TABLE    */

;

proc logistic data= &train. descending outest=LogisticTest outmodel=LogisticModel plots(only)=(roc) PLOTS(MAXPOINTS=NONE);

class      &binvarlist.    /param=reference ref=first;

model &targetvar.(event='1')=

&varlist.

/lackfit ctable selection=stepwise slentry= .05 sls=.05 ridging=none;

score data=&train. out= Scored_Training_Data outroc= ROC_TABLE_TRAINING_Data;

score data=&validate. out= Scored_Validation_Data outroc= ROC_TABLE_VALIDATION_DATA;

run;

ods html close;

ods graphics off;

ods output close;

proc sort data= TBL_ParamEst;

by descending waldchisq;

run;

/* The above code sorts the significant estimates based on the Wald Chi Square */

data TBL_Classification (drop = B);

set TBL_classification;

precision = correctevents/(correctevents + Incorrectevents);

recall = correctevents/(correctevents + Incorrectnonevents);

F_stat1 = harmean(precision,recall);

B = 2;

F_stat2 = (((1 + B*B) * (precision * recall))/((B*B*precision) + recall));

run;

proc sort data= LOG_Classification;

by descending F_Stat;

run;

/* The above code uses the Classification Table and calculates the precision, recall and sorts the table based on the F ratio*/

/* It also calculates the adjusted F ratio where higher weightage is given on Recall */

/* The probability value where the F ratio (or adjusted F ratio) is maximum should be treated as the threshold probability during validation */

data ROC_TABLE_TRAINING_Data(drop = B);

set ROC_TABLE_TRAINING_Data;

precision = _POS_/(_POS_ + _FALPOS_);

recall = _POS_/(_POS_ + _FALNEG_);

F_stat = harmean(precision,recall);

B = 2;

F_stat2 = (((1 + B*B) * (precision * recall))/((B*B*precision) + recall));

run;

proc sort data= ROC_TABLE_TRAINING_Data;

by descending F_Stat;

run;

/* This code should give the same precision and recall value as the above table. Instead of using the classification table, here the ROC table made on the training data is used to calculate precision and recall     */

data prob_threshold;

set ROC_TABLE_TRAINING_Data;

if _n_ = 1 then output;

run;

data _null_;

set prob_threshold;

call symput(“prob_thresh”,compress(_prob_);

run;

/* Saving the threshold value of probability for demarcating between 1 and 0 */

data ROC_TABLE_VALIDATION_DATA (drop = B);

set ROC_TABLE_VALIDATION_DATA;

precision = _POS_/(_POS_ + _FALPOS_);

recall = _POS_/(_POS_ + _FALNEG_);

F_stat = harmean(precision,recall);

B = 2;

F_stat2 = (((1 + B*B) * (precision * recall))/((B*B*precision) + recall));

run;

data Validation_Check;

set ROC_TABLE_VALIDATION_DATA;

if _prob_ ge &prob_thresh. then output;     /* Use the Probability_Threshold from the step above */

run;

proc sort data= Validation_Check;

by _prob_;

run;

/* This part of the code calculates the precision and recall in the validation data at the pre-fixed level of probability threshold value. */

Automated SAS code for variable reduction: Variance Inflation Factor (VIF)

In statistics (or econometrics), the variance inflation factor (VIF) calculates incidence and severity of multicollinearity among the independent variables in an ordinary least squares (OLS) regression analysis. One can read more about problems of multicollinearity here and about VIF here. In a linear regression analysis, it is important to run the VIF test to remove the multicollinearity among the independent variables. As a rule of thumb, the VIF value should not be more than 2 for better modelling. However, the final decision depends on the analyst’s discretion.

While running the linear regression analysis, one should not remove all the variables which have VIF more than the pre-decided threshold value (in this case, say 5). Instead, the analyst should remove the variable having the highest VIF value, and then re-calculate the VIF values. The process of calculation and removal of variable should continue till the highest VIF comes lower than the threshold level, only variable being removed at a time.

The process is not a difficult one, but might turn to be cumbersome process if the number of independent variables is very high. The following SAS code is an automated code to solve the problem multiple iterations, and the final datasets gives the list of retained variables as well as removed variables. The SAS code uses proc reg as the only statistical procedure to calculate the VIF automatically. The iterations are used to remove one variable at a time.

libname dataloc “/Desktop/Model";   /* MODEL DATASET LOCATION      */

%let inset=MODEL_DATA;        /* MODEL DATASET NAME    */

%let target= Y_VAR;           /* DEPENDENT VARIABLE    */

libname outlib “/Desktop/Output"; /* OUTPUT LOCATION       */

%let VIF_limit = 2;           /* VIF LIMIT             */
%let VIF_val = 100;

/* VARIABLE LIST */

%let varlist =

X1

X2

X3

…

X100

;

data inset;

set dataloc.&inset.;

run;

/* TO CREATE A BLANK TABLE FOR REMOVED VARIABLE */

ods output "Parameter Estimates"=vif;
proc reg data=inset ;

model &target. =

&varlist.

/VIF;

run;

data outlib.removed_variable_list;

set vif;

if _n_ = 0 then output outlib.removed_variable_list; 
run;

/* LOOP FOR ITERATIONS */

%macro vif_automated;

%do %while (%sysevalf(&vif_val. > &VIF_limit.));

ods output "Parameter Estimates"=vif;

proc reg data=inset ;

model &target. =

&varlist.

/VIF;

run;

proc sort data=     vif;

by descending VarianceInflation;

run;

data vif_top vif_others;

set vif;

if _n_ = 1 then output vif_top;

if _n_ gt 1 then output vif_others;
run;

data vif_top;

set vif_top;

call symput( compress("vif_val"),compress(VarianceInflation));

run;

data outlib.removed_variable_list;

set vif_top outlib.removed_variable_list;

proc sql;

select distinct variable into: varlist separated by " "

from vif_others

where variable ^= "Intercept"

;

quit;

%end;

data outlib.final_variable_vif;

set vif;

run;

%mend;

%vif_automated;