# K Nearest Neighbor in R

K-nearest neighbor is one of many nonlinear algorithms that can be used in machine learning. By non-linear I mean that a linear combination of the features or variables is not needed in order to develop decision boundaries. This allows for the analysis of data that naturally does not meet the assumptions of linearity.

KNN is also known as a “lazy learner”. This means that there are known coefficients or parameter estimates. When doing regression we always had coefficient outputs regardless of the type of regression (ridge, lasso, elastic net, etc.). What KNN does instead is used K nearest neighbors to give a label to an unlabeled example. Our job when using KNN is to determine the number of K neighbors to use that is most accurate based on the different criteria for assessing the models.

In this post, we will develop a KNN model using the “Mroz” dataset from the “Ecdat” package. Our goal is to predict if someone lives in the city based on the other predictor variables. Below is some initial code.

library(class);library(kknn);library(caret);library(corrplot)
library(reshape2);library(ggplot2);library(pROC);library(Ecdat)
data(Mroz)
str(Mroz)
## 'data.frame':    753 obs. of  18 variables:
##  $work : Factor w/ 2 levels "yes","no": 2 2 2 2 2 2 2 2 2 2 ... ##$ hoursw    : int  1610 1656 1980 456 1568 2032 1440 1020 1458 1600 ...
##  $child6 : int 1 0 1 0 1 0 0 0 0 0 ... ##$ child618  : int  0 2 3 3 2 0 2 0 2 2 ...
##  $agew : int 32 30 35 34 31 54 37 54 48 39 ... ##$ educw     : int  12 12 12 12 14 12 16 12 12 12 ...
##  $hearnw : num 3.35 1.39 4.55 1.1 4.59 ... ##$ wagew     : num  2.65 2.65 4.04 3.25 3.6 4.7 5.95 9.98 0 4.15 ...
##  $hoursh : int 2708 2310 3072 1920 2000 1040 2670 4120 1995 2100 ... ##$ ageh      : int  34 30 40 53 32 57 37 53 52 43 ...
##  $educh : int 12 9 12 10 12 11 12 8 4 12 ... ##$ wageh     : num  4.03 8.44 3.58 3.54 10 ...
##  $income : int 16310 21800 21040 7300 27300 19495 21152 18900 20405 20425 ... ##$ educwm    : int  12 7 12 7 12 14 14 3 7 7 ...
##  $educwf : int 7 7 7 7 14 7 7 3 7 7 ... ##$ unemprate : num  5 11 5 5 9.5 7.5 5 5 3 5 ...
##  $city : Factor w/ 2 levels "no","yes": 1 2 1 1 2 2 1 1 1 1 ... ##$ experience: int  14 5 15 6 7 33 11 35 24 21 ...

We need to remove the factor variable “work” as KNN cannot use factor variables. After this, we will use the “melt” function from the “reshape2” package to look at the variables when divided by whether the example was from the city or not.

Mroz$work<-NULL mroz.melt<-melt(Mroz,id.var='city') Mroz_plots<-ggplot(mroz.melt,aes(x=city,y=value))+geom_boxplot()+facet_wrap(~variable, ncol = 4) Mroz_plots From the plots, it appears there are no differences in how the variable act whether someone is from the city or not. This may be a flag that classification may not work. We now need to scale our data otherwise the results will be inaccurate. Scaling might also help our box-plots because everything will be on the same scale rather than spread all over the place. To do this we will have to temporarily remove our outcome variable from the data set because it’s a factor and then reinsert it into the data set. Below is the code. mroz.scale<-as.data.frame(scale(Mroz[,-16])) mroz.scale$city<-Mroz$city We will now look at our box-plots a second time but this time with scaled data. mroz.scale.melt<-melt(mroz.scale,id.var="city") mroz_plot2<-ggplot(mroz.scale.melt,aes(city,value))+geom_boxplot()+facet_wrap(~variable, ncol = 4) mroz_plot2 This second plot is easier to read but there is still little indication of difference. We can now move to checking the correlations among the variables. Below is the code mroz.cor<-cor(mroz.scale[,-17]) corrplot(mroz.cor,method = 'number') There is a high correlation between husband’s age (ageh) and wife’s age (agew). Since this algorithm is non-linear this should not be a major problem. We will now divide our dataset into the training and testing sets set.seed(502) ind=sample(2,nrow(mroz.scale),replace=T,prob=c(.7,.3)) train<-mroz.scale[ind==1,] test<-mroz.scale[ind==2,] Before creating a model we need to create a grid. We do not know the value of k yet so we have to run multiple models with different values of k in order to determine this for our model. As such we need to create a ‘grid’ using the ‘expand.grid’ function. We will also use cross-validation to get a better estimate of k as well using the “trainControl” function. The code is below. grid<-expand.grid(.k=seq(2,20,by=1)) control<-trainControl(method="cv") Now we make our model, knn.train<-train(city~.,train,method="knn",trControl=control,tuneGrid=grid) knn.train ## k-Nearest Neighbors ## ## 540 samples ## 16 predictors ## 2 classes: 'no', 'yes' ## ## No pre-processing ## Resampling: Cross-Validated (10 fold) ## Summary of sample sizes: 487, 486, 486, 486, 486, 486, ... ## Resampling results across tuning parameters: ## ## k Accuracy Kappa ## 2 0.6000095 0.1213920 ## 3 0.6368757 0.1542968 ## 4 0.6424325 0.1546494 ## 5 0.6386252 0.1275248 ## 6 0.6329998 0.1164253 ## 7 0.6589619 0.1616377 ## 8 0.6663344 0.1774391 ## 9 0.6663681 0.1733197 ## 10 0.6609510 0.1566064 ## 11 0.6664018 0.1575868 ## 12 0.6682199 0.1669053 ## 13 0.6572111 0.1397222 ## 14 0.6719586 0.1694953 ## 15 0.6571425 0.1263937 ## 16 0.6664367 0.1551023 ## 17 0.6719573 0.1588789 ## 18 0.6608811 0.1260452 ## 19 0.6590979 0.1165734 ## 20 0.6609510 0.1219624 ## ## Accuracy was used to select the optimal model using the largest value. ## The final value used for the model was k = 14. R recommends that k = 16. This is based on a combination of accuracy and the kappa statistic. The kappa statistic is a measurement of the accuracy of a model while taking into account chance. We don’t have a model in the sense that we do not use the ~ sign like we do with regression. Instead, we have a train and a test set a factor variable and a number for k. This will make more sense when you see the code. Finally, we will use this information on our test dataset. We will then look at the table and the accuracy of the model. knn.test<-knn(train[,-17],test[,-17],train[,17],k=16) #-17 removes the dependent variable 'city table(knn.test,test$city)
##
## knn.test  no yes
##      no   19   8
##      yes  61 125
prob.agree<-(15+129)/213
prob.agree
## [1] 0.6760563

Accuracy is 67% which is consistent with what we found when determining the k. We can also calculate the kappa. This done by calculating the probability and then do some subtraction and division. We already know the accuracy as we stored it in the variable “prob.agree” we now need the probability that this is by chance. Lastly, we calculate the kappa.

prob.chance<-((15+4)/213)*((15+65)/213)
kap<-(prob.agree-prob.chance)/(1-prob.chance)
kap
## [1] 0.664827

A kappa of .66 is actual good.

The example we just did was with unweighted k neighbors. There are times when weighted neighbors can improve accuracy. We will look at three different weighing methods. “Rectangular” is unweighted and is the one that we used. The other two are “triangular” and “epanechnikov”. How these calculate the weights is beyond the scope of this post. In the code below the argument “distance” can be set to 1 for euclidean and 2 for absolute distance.

kknn.train<-train.kknn(city~.,train,kmax = 25,distance = 2,kernel = c("rectangular","triangular",
"epanechnikov"))
plot(kknn.train)

kknn.train
##
## Call:
## train.kknn(formula = city ~ ., data = train, kmax = 25, distance = 2,     kernel = c("rectangular", "triangular", "epanechnikov"))
##
## Type of response variable: nominal
## Minimal misclassification: 0.3277778
## Best kernel: rectangular
## Best k: 14

If you look at the plot you can see which value of k is the best by looking at the point that is the lowest on the graph which is right before 15. Looking at the legend it indicates that the point is the “rectangular” estimate which is the same as unweighted. This means that the best classification is unweighted with a k of 14. Although it recommends a different value for k our misclassification was about the same.

Conclusion

In this post, we explored both weighted and unweighted KNN. This algorithm allows you to deal with data that does not meet the assumptions of regression by ignoring the need for parameters. However, because there are no numbers really attached to the results beyond accuracy it can be difficult to explain what is happening in the model to people. As such, perhaps the biggest drawback is communicating results when using KNN.

# Elastic Net Regression in R

Elastic net is a combination of ridge and lasso regression. What is most unusual about elastic net is that it has two tuning parameters (alpha and lambda) while lasso and ridge regression only has 1.

In this post, we will go through an example of the use of elastic net using the “VietnamI” dataset from the “Ecdat” package. Our goal is to predict how many days a person is ill based on the other variables in the dataset. Below is some initial code for our analysis

library(Ecdat);library(corrplot);library(caret);library(glmnet)
data("VietNamI")
str(VietNamI)
## 'data.frame':    27765 obs. of  12 variables:
##  $pharvis : num 0 0 0 1 1 0 0 0 2 3 ... ##$ lnhhexp  : num  2.73 2.74 2.27 2.39 3.11 ...
##  $age : num 3.76 2.94 2.56 3.64 3.3 ... ##$ sex      : Factor w/ 2 levels "female","male": 2 1 2 1 2 2 1 2 1 2 ...
##  $married : num 1 0 0 1 1 1 1 0 1 1 ... ##$ educ     : num  2 0 4 3 3 9 2 5 2 0 ...
##  $illness : num 1 1 0 1 1 0 0 0 2 1 ... ##$ injury   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $illdays : num 7 4 0 3 10 0 0 0 4 7 ... ##$ actdays  : num  0 0 0 0 0 0 0 0 0 0 ...
##  $insurance: num 0 0 1 1 0 1 1 1 0 0 ... ##$ commune  : num  192 167 76 123 148 20 40 57 49 170 ...
##  - attr(*, "na.action")=Class 'omit'  Named int 27734
##   .. ..- attr(*, "names")= chr "27734"

We need to check the correlations among the variables. We need to exclude the “sex” variable as it is categorical. Code is below.

p.cor<-cor(VietNamI[,-4])
corrplot.mixed(p.cor)

No major problems with correlations. Next, we set up our training and testing datasets. We need to remove the variable “commune” because it adds no value to our results. In addition, to reduce the computational time we will only use the first 1000 rows from the data set.

VietNamI$commune<-NULL VietNamI_reduced<-VietNamI[1:1000,] ind<-sample(2,nrow(VietNamI_reduced),replace=T,prob = c(0.7,0.3)) train<-VietNamI_reduced[ind==1,] test<-VietNamI_reduced[ind==2,] We need to create a grid that will allow us to investigate different models with different combinations of alpha ana lambda. This is done using the “expand.grid” function. In combination with the “seq” function below is the code grid<-expand.grid(.alpha=seq(0,1,by=.5),.lambda=seq(0,0.2,by=.1)) We also need to set the resampling method, which allows us to assess the validity of our model. This is done using the “trainControl” function” from the “caret” package. In the code below “LOOCV” stands for “leave one out cross-validation”. control<-trainControl(method = "LOOCV") We are no ready to develop our model. The code is mostly self-explanatory. This initial model will help us to determine the appropriate values for the alpha and lambda parameters enet.train<-train(illdays~.,train,method="glmnet",trControl=control,tuneGrid=grid) enet.train ## glmnet ## ## 694 samples ## 10 predictors ## ## No pre-processing ## Resampling: Leave-One-Out Cross-Validation ## Summary of sample sizes: 693, 693, 693, 693, 693, 693, ... ## Resampling results across tuning parameters: ## ## alpha lambda RMSE Rsquared ## 0.0 0.0 5.229759 0.2968354 ## 0.0 0.1 5.229759 0.2968354 ## 0.0 0.2 5.229759 0.2968354 ## 0.5 0.0 5.243919 0.2954226 ## 0.5 0.1 5.225067 0.2985989 ## 0.5 0.2 5.200415 0.3038821 ## 1.0 0.0 5.244020 0.2954519 ## 1.0 0.1 5.203973 0.3033173 ## 1.0 0.2 5.182120 0.3083819 ## ## RMSE was used to select the optimal model using the smallest value. ## The final values used for the model were alpha = 1 and lambda = 0.2. The output list all the possible alpha and lambda values that we set in the “grid” variable. It even tells us which combination was the best. For our purposes, the alpha will be .5 and the lambda .2. The r-square is also included. We will set our model and run it on the test set. We have to convert the “sex” variable to a dummy variable for the “glmnet” function. We next have to make matrices for the predictor variables and a for our outcome variable “illdays” train$sex<-model.matrix( ~ sex - 1, data=train ) #convert to dummy variable
test$sex<-model.matrix( ~ sex - 1, data=test ) predictor_variables<-as.matrix(train[,-9]) days_ill<-as.matrix(train$illdays)
enet<-glmnet(predictor_variables,days_ill,family = "gaussian",alpha = 0.5,lambda = .2)

We can now look at specific coefficient by using the “coef” function.

enet.coef<-coef(enet,lambda=.2,alpha=.5,exact=T)
enet.coef
## 12 x 1 sparse Matrix of class "dgCMatrix"
##                         s0
## (Intercept)   -1.304263895
## pharvis        0.532353361
## lnhhexp       -0.064754000
## age            0.760864404
## sex.sexfemale  0.029612290
## sex.sexmale   -0.002617404
## married        0.318639271
## educ           .
## illness        3.103047473
## injury         .
## actdays        0.314851347
## insurance      .

You can see for yourself that several variables were removed from the model. Medical expenses (lnhhexp), sex, education, injury, and insurance do not play a role in the number of days ill for an individual in Vietnam.

With our model developed. We now can test it using the predict function. However, we first need to convert our test dataframe into a matrix and remove the outcome variable from it

test.matrix<-as.matrix(test[,-9])
enet.y<-predict(enet, newx = test.matrix, type = "response", lambda=.2,alpha=.5)

Let’s plot our results

plot(enet.y)

This does not look good. Let’s check the mean squared error

enet.resid<-enet.y-test$illdays mean(enet.resid^2) ## [1] 20.18134 We will now do a cross-validation of our model. We need to set the seed and then use the “cv.glmnet” to develop the cross-validated model. We can see the model by plotting it. set.seed(317) enet.cv<-cv.glmnet(predictor_variables,days_ill,alpha=.5) plot(enet.cv) You can see that as the number of features are reduce (see the numbers on the top of the plot) the MSE increases (y-axis). In addition, as the lambda increases, there is also an increase in the error but only when the number of variables is reduced as well. The dotted vertical lines in the plot represent the minimum MSE for a set lambda (on the left) and the one standard error from the minimum (on the right). You can extract these two lambda values using the code below. enet.cv$lambda.min
## [1] 0.3082347
enet.cv$lambda.1se ## [1] 2.874607 We can see the coefficients for a lambda that is one standard error away by using the code below. This will give us an alternative idea for what to set the model parameters to when we want to predict. coef(enet.cv,s="lambda.1se") ## 12 x 1 sparse Matrix of class "dgCMatrix" ## 1 ## (Intercept) 2.34116947 ## pharvis 0.003710399 ## lnhhexp . ## age . ## sex.sexfemale . ## sex.sexmale . ## married . ## educ . ## illness 1.817479480 ## injury . ## actdays . ## insurance . Using the one standard error lambda we lose most of our features. We can now see if the model improves by rerunning it with this information. enet.y.cv<-predict(enet.cv,newx = test.matrix,type='response',lambda="lambda.1se", alpha = .5) enet.cv.resid<-enet.y.cv-test$illdays
mean(enet.cv.resid^2)
## [1] 25.47966

A small improvement.  Our model is a mess but this post served as an example of how to conduct an analysis using elastic net regression.

# Lasso Regression in R

In this post, we will conduct an analysis using the lasso regression. Remember lasso regression will actually eliminate variables by reducing them to zero through how the shrinkage penalty can be applied.

We will use the dataset “nlschools” from the “MASS” packages to conduct our analysis. We want to see if we can predict language test scores “lang” with the other available variables. Below is some initial code to begin the analysis

library(MASS);library(corrplot);library(glmnet)
data("nlschools")
str(nlschools)
## 'data.frame':    2287 obs. of  6 variables:
##  $lang : int 46 45 33 46 20 30 30 57 36 36 ... ##$ IQ   : num  15 14.5 9.5 11 8 9.5 9.5 13 9.5 11 ...
##  $class: Factor w/ 133 levels "180","280","1082",..: 1 1 1 1 1 1 1 1 1 1 ... ##$ GS   : int  29 29 29 29 29 29 29 29 29 29 ...
##  $SES : int 23 10 15 23 10 10 23 10 13 15 ... ##$ COMB : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...

We need to remove the “class” variable as it is used as an identifier and provides no useful data. After this, we can check the correlations among the variables. Below is the code for this.

nlschools$class<-NULL p.cor<-cor(nlschools[,-5]) corrplot.mixed(p.cor) No problems with collinearity. We will now setup are training and testing sets. ind<-sample(2,nrow(nlschools),replace=T,prob = c(0.7,0.3)) train<-nlschools[ind==1,] test<-nlschools[ind==2,] Remember that the ‘glmnet’ function does not like factor variables. So we need to convert our “COMB” variable to a dummy variable. In addition, “glmnet” function does not like data frames so we need to make two data frames. The first will include all the predictor variables and the second we include only the outcome variable. Below is the code train$COMB<-model.matrix( ~ COMB - 1, data=train ) #convert to dummy variable
test$COMB<-model.matrix( ~ COMB - 1, data=test ) predictor_variables<-as.matrix(train[,2:4]) language_score<-as.matrix(train$lang)

We can now run our model. We place both matrices inside the “glmnet” function. The family is set to “gaussian” because our outcome variable is continuous. The “alpha” is set to 1 as this indicates that we are using lasso regression.

lasso<-glmnet(predictor_variables,language_score,family="gaussian",alpha=1)

Now we need to look at the results using the “print” function. This function prints a lot of information as explained below.

• Df = number of variables including in the model (this is always the same number in a ridge model)
• %Dev = Percent of deviance explained. The higher the better
• Lambda = The lambda used to obtain the %Dev

When you use the “print” function for a lasso model it will print up to 100 different models. Fewer models are possible if the percent of deviance stops improving. 100 is the default stopping point. In the code below we will use the “print” function but, I only printed the first 5 and last 5 models in order to reduce the size of the printout. Fortunately, it only took 60 models to converge.

print(lasso)
##
## Call:  glmnet(x = predictor_variables, y = language_score, family = "gaussian",      alpha = 1)
##
##       Df    %Dev  Lambda
##  [1,]  0 0.00000 5.47100
##  [2,]  1 0.06194 4.98500
##  [3,]  1 0.11340 4.54200
##  [4,]  1 0.15610 4.13900
##  [5,]  1 0.19150 3.77100
............................
## [55,]  3 0.39890 0.03599
## [56,]  3 0.39900 0.03280
## [57,]  3 0.39900 0.02988
## [58,]  3 0.39900 0.02723
## [59,]  3 0.39900 0.02481
## [60,]  3 0.39900 0.02261

The results from the “print” function will allow us to set the lambda for the “test” dataset. Based on the results we can set the lambda at 0.02 because this explains the highest amount of deviance at .39.

The plot below shows us lambda on the x-axis and the coefficients of the predictor variables on the y-axis. The numbers next to the coefficient lines refers to the actual coefficient of a particular variable as it changes from using different lambda values. Each number corresponds to a variable going from left to right in a dataframe/matrix using the “View” function. For example, 1 in the plot refers to “IQ” 2 refers to “GS” etc.

plot(lasso,xvar="lambda",label=T)

As you can see, as lambda increase the coefficient decrease in value. This is how regularized regression works. However, unlike ridge regression which never reduces a coefficient to zero, lasso regression does reduce a coefficient to zero. For example, coefficient 3 (SES variable) and coefficient 2 (GS variable) are reduced to zero when lambda is near 1.

You can also look at the coefficient values at a specific lambda values. The values are unstandardized and are used to determine the final model selection. In the code below the lambda is set to .02 and we use the “coef” function to do see the results

lasso.coef<-coef(lasso,s=.02,exact = T)
lasso.coef
## 4 x 1 sparse Matrix of class "dgCMatrix"
##                       1
## (Intercept)  9.35736325
## IQ           2.34973922
## GS          -0.02766978
## SES          0.16150542

Results indicate that for a 1 unit increase in IQ there is a 2.41 point increase in language score. When GS (class size) goes up 1 unit there is a .03 point decrease in language score. Finally, when SES (socioeconomic status) increase  1 unit language score improves .13 point.

The second plot shows us the deviance explained on the x-axis. On the y-axis is the coefficients of the predictor variables. Below is the code

plot(lasso,xvar='dev',label=T)

If you look carefully, you can see that the two plots are completely opposite to each other. increasing lambda cause a decrease in the coefficients. Furthermore, increasing the fraction of deviance explained leads to an increase in the coefficient. You may remember seeing this when we used the “print”” function. As lambda became smaller there was an increase in the deviance explained.

Now, we will assess our model using the test data. We need to convert the test dataset to a matrix. Then we will use the “predict”” function while setting our lambda to .02. Lastly, we will plot the results. Below is the code.

test.matrix<-as.matrix(test[,2:4])
lasso.y<-predict(lasso,newx = test.matrix,type = 'response',s=.02)
plot(lasso.y,test$lang) The visual looks promising. The last thing we need to do is calculated the mean squared error. By its self this number does not mean much. However, it provides a benchmark for comparing our current model with any other models that we may develop. Below is the code lasso.resid<-lasso.y-test$lang
mean(lasso.resid^2)
## [1] 46.74314

Knowing this number, we can, if we wanted, develop other models using other methods of analysis to try to reduce it. Generally, the lower the error the better while keeping in mind the complexity of the model.

# Ridge Regression in R

In this post, we will conduct an analysis using ridge regression. Ridge regression is a type of regularized regression. By applying a shrinkage penalty, we are able to reduce the coefficients of many variables almost to zero while still retaining them in the model. This allows us to develop models that have many more variables in them compared to models using best subset or stepwise regression.

In the example used in this post, we will use the “SAheart” dataset from the “ElemStatLearn” package. We want to predict systolic blood pressure (sbp) using all of the other variables available as predictors. Below is some initial code that we need to begin.

library(ElemStatLearn);library(car);library(corrplot)
library(leaps);library(glmnet);library(caret)
data(SAheart)
str(SAheart)
## 'data.frame':    462 obs. of  10 variables:
##  $sbp : int 160 144 118 170 134 132 142 114 114 132 ... ##$ tobacco  : num  12 0.01 0.08 7.5 13.6 6.2 4.05 4.08 0 0 ...
##  $ldl : num 5.73 4.41 3.48 6.41 3.5 6.47 3.38 4.59 3.83 5.8 ... ##$ adiposity: num  23.1 28.6 32.3 38 27.8 ...
##  $famhist : Factor w/ 2 levels "Absent","Present": 2 1 2 2 2 2 1 2 2 2 ... ##$ typea    : int  49 55 52 51 60 62 59 62 49 69 ...
##  $obesity : num 25.3 28.9 29.1 32 26 ... ##$ alcohol  : num  97.2 2.06 3.81 24.26 57.34 ...
##  $age : int 52 63 46 58 49 45 38 58 29 53 ... ##$ chd      : int  1 1 0 1 1 0 0 1 0 1 ...

A look at the object using the “str” function indicates that one variable “famhist” is a factor variable. The “glmnet” function that does the ridge regression analysis cannot handle factors so we need to converts this to a dummy variable. However, there are two things we need to do before this. First, we need to check the correlations to make sure there are no major issues with multi-collinearity Second, we need to create our training and testing data sets. Below is the code for the correlation plot.

p.cor<-cor(SAheart[,-5])
corrplot.mixed(p.cor)

First we created a variable called “p.cor” the -5 in brackets means we removed the 5th column from the “SAheart” data set which is the factor variable “Famhist”. The correlation plot indicates that there is one strong relationship between adiposity and obesity. However, one common cut-off for collinearity is 0.8 and this value is 0.72 which is not a problem.

We will now create are training and testing sets and convert “famhist” to a dummy variable.

ind<-sample(2,nrow(SAheart),replace=T,prob = c(0.7,0.3))
train<-SAheart[ind==1,]
test<-SAheart[ind==2,]
train$famhist<-model.matrix( ~ famhist - 1, data=train ) #convert to dummy variable test$famhist<-model.matrix( ~ famhist - 1, data=test )

We are still not done preparing our data yet. “glmnet” cannot use data frames, instead, it can only use matrices. Therefore, we now need to convert our data frames to matrices. We have to create two matrices, one with all of the predictor variables and a second with the outcome variable of blood pressure. Below is the code

predictor_variables<-as.matrix(train[,2:10])
blood_pressure<-as.matrix(train$sbp) We are now ready to create our model. We use the “glmnet” function and insert our two matrices. The family is set to Gaussian because “blood pressure” is a continuous variable. Alpha is set to 0 as this indicates ridge regression. Below is the code ridge<-glmnet(predictor_variables,blood_pressure,family = 'gaussian',alpha = 0) Now we need to look at the results using the “print” function. This function prints a lot of information as explained below. • Df = number of variables including in the model (this is always the same number in a ridge model) • %Dev = Percent of deviance explained. The higher the better • Lambda = The lambda used to attain the %Dev When you use the “print” function for a ridge model it will print up to 100 different models. Fewer models are possible if the percent of deviance stops improving. 100 is the default stopping point. In the code below we have the “print” function. However, I have only printed the first 5 and last 5 models in order to save space. print(ridge) ## ## Call: glmnet(x = predictor_variables, y = blood_pressure, family = "gaussian", alpha = 0) ## ## Df %Dev Lambda ## [1,] 10 7.622e-37 7716.0000 ## [2,] 10 2.135e-03 7030.0000 ## [3,] 10 2.341e-03 6406.0000 ## [4,] 10 2.566e-03 5837.0000 ## [5,] 10 2.812e-03 5318.0000 ................................ ## [95,] 10 1.690e-01 1.2290 ## [96,] 10 1.691e-01 1.1190 ## [97,] 10 1.692e-01 1.0200 ## [98,] 10 1.693e-01 0.9293 ## [99,] 10 1.693e-01 0.8468 ## [100,] 10 1.694e-01 0.7716 The results from the “print” function are useful in setting the lambda for the “test” dataset. Based on the results we can set the lambda at 0.83 because this explains the highest amount of deviance at .20. The plot below shows us lambda on the x-axis and the coefficients of the predictor variables on the y-axis. The numbers refer to the actual coefficient of a particular variable. Inside the plot, each number corresponds to a variable going from left to right in a data-frame/matrix using the “View” function. For example, 1 in the plot refers to “tobacco” 2 refers to “ldl” etc. Across the top of the plot is the number of variables used in the model. Remember this number never changes when doing ridge regression. plot(ridge,xvar="lambda",label=T) As you can see, as lambda increase the coefficient decrease in value. This is how ridge regression works yet no coefficient ever goes to absolute 0. You can also look at the coefficient values at a specific lambda value. The values are unstandardized but they provide a useful insight when determining final model selection. In the code below the lambda is set to .83 and we use the “coef” function to do this ridge.coef<-coef(ridge,s=.83,exact = T) ridge.coef ## 11 x 1 sparse Matrix of class "dgCMatrix" ## 1 ## (Intercept) 105.69379942 ## tobacco -0.25990747 ## ldl -0.13075557 ## adiposity 0.29515034 ## famhist.famhistAbsent 0.42532887 ## famhist.famhistPresent -0.40000846 ## typea -0.01799031 ## obesity 0.29899976 ## alcohol 0.03648850 ## age 0.43555450 ## chd -0.26539180 The second plot shows us the deviance explained on the x-axis and the coefficients of the predictor variables on the y-axis. Below is the code plot(ridge,xvar='dev',label=T) The two plots are completely opposite to each other. Increasing lambda cause a decrease in the coefficients while increasing the fraction of deviance explained leads to an increase in the coefficient. You can also see this when we used the “print” function. As lambda became smaller there was an increase in the deviance explained. We now can begin testing our model on the test data set. We need to convert the test dataset to a matrix and then we will use the predict function while setting our lambda to .83 (remember a lambda of .83 explained the most of the deviance). Lastly, we will plot the results. Below is the code. test.matrix<-as.matrix(test[,2:10]) ridge.y<-predict(ridge,newx = test.matrix,type = 'response',s=.83) plot(ridge.y,test$sbp)

The last thing we need to do is calculated the mean squared error. By it’s self this number is useless. However, it provides a benchmark for comparing the current model with any other models you may develop. Below is the code

ridge.resid<-ridge.y-test$sbp mean(ridge.resid^2) ## [1] 372.4431 Knowing this number, we can develop other models using other methods of analysis to try to reduce it as much as possible. # Linear VS Quadratic Discriminant Analysis in R In this post we will look at linear discriminant analysis (LDA) and quadratic discriminant analysis (QDA). Discriminant analysis is used when the dependent variable is categorical. Another commonly used option is logistic regression but there are differences between logistic regression and discriminant analysis. Both LDA and QDA are used in situations in which there is a clear separation between the classes you want to predict. If the categories are fuzzier logistic regression is often the better choice. For our example, we will use the “Mathlevel” dataset found in the “Ecdat” package. Our goal will be to predict the sex of a respondent based on SAT math score, major, foreign language proficiency, as well as the number of math, physic, and chemistry classes a respondent took. Below is some initial code to start our analysis. library(MASS);library(Ecdat) data("Mathlevel") The first thing we need to do is clean up the data set. We have to remove any missing data in order to run our model. We will create a dataset called “math” that has the “Mathlevel” dataset but with the “NA”s removed use the “na.omit” function. After this, we need to set our seed for the purpose of reproducibility using the “set.seed” function. Lastly, we will split the data using the “sample” function using a 70/30 split. The training dataset will be called “math.train” and the testing dataset will be called “math.test”. Below is the code math<-na.omit(Mathlevel) set.seed(123) math.ind<-sample(2,nrow(math),replace=T,prob = c(0.7,0.3)) math.train<-math[math.ind==1,] math.test<-math[math.ind==2,] Now we will make our model and it is called “lda.math” and it will include all available variables in the “math.train” dataset. Next we will check the results by calling the modle. Finally, we will examine the plot to see how our model is doing. Below is the code. lda.math<-lda(sex~.,math.train) lda.math ## Call: ## lda(sex ~ ., data = math.train) ## ## Prior probabilities of groups: ## male female ## 0.5986079 0.4013921 ## ## Group means: ## mathlevel.L mathlevel.Q mathlevel.C mathlevel^4 mathlevel^5 ## male -0.10767593 0.01141838 -0.05854724 0.2070778 0.05032544 ## female -0.05571153 0.05360844 -0.08967303 0.2030860 -0.01072169 ## mathlevel^6 sat languageyes majoreco majoross majorns ## male -0.2214849 632.9457 0.07751938 0.3914729 0.1472868 0.1782946 ## female -0.2226767 613.6416 0.19653179 0.2601156 0.1907514 0.2485549 ## majorhum mathcourse physiccourse chemistcourse ## male 0.05426357 1.441860 0.7441860 1.046512 ## female 0.07514451 1.421965 0.6531792 1.040462 ## ## Coefficients of linear discriminants: ## LD1 ## mathlevel.L 1.38456344 ## mathlevel.Q 0.24285832 ## mathlevel.C -0.53326543 ## mathlevel^4 0.11292817 ## mathlevel^5 -1.24162715 ## mathlevel^6 -0.06374548 ## sat -0.01043648 ## languageyes 1.50558721 ## majoreco -0.54528930 ## majoross 0.61129797 ## majorns 0.41574298 ## majorhum 0.33469586 ## mathcourse -0.07973960 ## physiccourse -0.53174168 ## chemistcourse 0.16124610 plot(lda.math,type='both') Calling “lda.math” gives us the details of our model. It starts be indicating the prior probabilities of someone being male or female. Next is the means for each variable by sex. The last part is the coefficients of the linear discriminants. Each of these values is used to determine the probability that a particular example is male or female. This is similar to a regression equation. The plot provides us with densities of the discriminant scores for males and then for females. The output indicates a problem. There is a great deal of overlap between male and females in the model. What this indicates is that there is a lot of misclassification going on as the two groups are not clearly separated. Furthermore, this means that logistic regression is probably a better choice for distinguishing between male and females. However, since this is for demonstrating purposes we will not worry about this. We will now use the “predict” function on the training set data to see how well our model classifies the respondents by gender. We will then compare the prediction of the model with thee actual classification. Below is the code. math.lda.predict<-predict(lda.math) math.train$lda<-math.lda.predict$class table(math.train$lda,math.train$sex) ## ## male female ## male 219 100 ## female 39 73 mean(math.train$lda==math.train$sex) ## [1] 0.6774942 As you can see, we have a lot of misclassification happening. A large amount of false negatives which is a lot of males being classified as female. The overall accuracy us only 59% which is not much better than chance. We will now conduct the same analysis on the test data set. Below is the code. lda.math.test<-predict(lda.math,math.test) math.test$lda<-lda.math.test$class table(math.test$lda,math.test$sex) ## ## male female ## male 92 43 ## female 23 20 mean(math.test$lda==math.test$sex) ## [1] 0.6292135 As you can see the results are similar. To put it simply, our model is terrible. The main reason is that there is little distinction between males and females as shown in the plot. However, we can see if perhaps a quadratic discriminant analysis will do better QDA allows for each class in the dependent variable to have it’s own covariance rather than a shared covariance as in LDA. This allows for quadratic terms in the development of the model. To complete a QDA we need to use the “qda” function from the “MASS” package. Below is the code for the training data set. math.qda.fit<-qda(sex~.,math.train) math.qda.predict<-predict(math.qda.fit) math.train$qda<-math.qda.predict$class table(math.train$qda,math.train$sex) ## ## male female ## male 215 84 ## female 43 89 mean(math.train$qda==math.train$sex) ## [1] 0.7053364 You can see there is almost no difference. Below is the code for the test data. math.qda.test<-predict(math.qda.fit,math.test) math.test$qda<-math.qda.test$class table(math.test$qda,math.test$sex) ## ## male female ## male 91 43 ## female 24 20 mean(math.test$qda==math.test$sex) ## [1] 0.6235955 Still disappointing. However, in this post we reviewed linear discriminant analysis as well as learned about the use of quadratic linear discriminant analysis. Both of these statistical tools are used for predicting categorical dependent variables. LDA assumes shared covariance in the dependent variable categories will QDA allows for each category in the dependent variable to have it’s own variance. # Validating a Logistic Model in R In this post, we are going to continue are analysis of the logistic regression model from the the post on logistic regression in R. We need to rerun all of the code from the last post to be ready to continue. As such the code form the last post is all below library(MASS);library(bestglm);library(reshape2);library(corrplot); library(ggplot2);library(ROCR) data(survey) survey$Clap<-NULL
survey$W.Hnd<-NULL survey$Fold<-NULL
survey$Exer<-NULL survey$Smoke<-NULL
survey$M.I<-NULL survey<-na.omit(survey) pm<-melt(survey, id.var="Sex") ggplot(pm,aes(Sex,value))+geom_boxplot()+facet_wrap(~variable,ncol = 3) pc<-cor(survey[,2:5]) corrplot.mixed(pc) set.seed(123) ind<-sample(2,nrow(survey),replace=T,prob = c(0.7,0.3)) train<-survey[ind==1,] test<-survey[ind==2,] fit<-glm(Sex~.,binomial,train) exp(coef(fit)) train$probs<-predict(fit, type = 'response')
train$predict<-rep('Female',123) train$predict[train$probs>0.5]<-"Male" table(train$predict,train$Sex) mean(train$predict==train$Sex) test$prob<-predict(fit,newdata = test, type = 'response')
test$predict<-rep('Female',46) test$predict[test$prob>0.5]<-"Male" table(test$predict,test$Sex) mean(test$predict==test$Sex) Model Validation We will now do a K-fold cross validation in order to further see how our model is doing. We cannot use the factor variable “Sex” with the K-fold code so we need to create a dummy variable. First, we create a variable called “y” that has 123 spaces, which is the same size as the “train” dataset. Second, we fill “y” with 1 in every example that is coded “male” in the “Sex” variable. In addition, we also need to create a new dataset and remove some variables from our prior analysis otherwise we will confuse the functions that we are going to use. We will remove “predict”, “Sex”, and “probs” train$y<-rep(0,123)
train$y[train$Sex=="Male"]=1
my.cv<-train[,-8]
my.cv$Sex<-NULL my.cv$probs<-NULL

We now can do our K-fold analysis. The code is complicated so you can trust it and double check on your own.

bestglm(Xy=my.cv,IC="CV",CVArgs = list(Method="HTF",K=10,REP=1),family = binomial)
## Morgan-Tatar search since family is non-gaussian.
## CV(K = 10, REP = 1)
## BICq equivalent for q in (6.66133814775094e-16, 0.0328567092272112)
## Best Model:
##                Estimate Std. Error   z value     Pr(>|z|)
## (Intercept) -45.2329733 7.80146036 -5.798014 6.710501e-09
## Height        0.2615027 0.04534919  5.766425 8.097067e-09

The results confirm what we alreaedy knew that only the “Height” variable is valuable in predicting Sex. We will now create our new model using only the recommendation of the kfold validation analysis. Then we check the new model against the train dataset and with the test dataset. The code below is a repeat of prior code but based on the cross-validation

reduce.fit<-glm(Sex~Height, family=binomial,train)
train$cv.probs<-predict(reduce.fit,type='response') train$cv.predict<-rep('Female',123)
train$cv.predict[train$cv.probs>0.5]='Male'
table(train$cv.predict,train$Sex)
##
##          Female Male
##   Female     61   11
##   Male        7   44
mean(train$cv.predict==train$Sex)
## [1] 0.8536585
test$cv.probs<-predict(reduce.fit,test,type = 'response') test$cv.predict<-rep('Female',46)
test$cv.predict[test$cv.probs>0.5]='Male'
table(test$cv.predict,test$Sex)
##
##          Female Male
##   Female     16    7
##   Male        1   22
mean(test$cv.predict==test$Sex)
## [1] 0.826087

The results are consistent for both the train and test dataset. We are now going to create the ROC curve. This will provide a visual and the AUC number to further help us to assess our model. However, a model is only good when it is compared to another model. Therefore, we will create a really bad model in order to compare it to the original model, and the cross validated model. We will first make a bad model and store the probabilities in the “test” dataset. The bad model will use “age” to predict “Sex” which doesn’t make any sense at all. Below is the code followed by the ROC curve of the bad model.

bad.fit<-glm(Sex~Age,family = binomial,test)
test$bad.probs<-predict(bad.fit,type='response') pred.bad<-prediction(test$bad.probs,test$Sex) perf.bad<-performance(pred.bad,'tpr','fpr') plot(perf.bad,col=1) The more of a diagonal the line is the worst it is. As, we can see the bad model is really bad. What we just did with the bad model we will now repeat for the full model and the cross-validated model. As before, we need to store the prediction in a way that the ROCR package can use them. We will create a variable called “pred.full” to begin the process of graphing the the original full model from the last blog post. Then we will use the “prediction” function. Next, we will create the “perf.full” variable to store the performance of the model. Notice, the arguments ‘tpr’ and ‘fpr’ for true positive rate and false positive rate. Lastly, we plot the results pred.full<-prediction(test$prob,test$Sex) perf.full<-performance(pred.full,'tpr','fpr') plot(perf.full, col=2) We repeat this process for the cross-validated model pred.cv<-prediction(test$cv.probs,test$Sex) perf.cv<-performance(pred.cv,'tpr','fpr') plot(perf.cv,col=3) Now let’s put all the different models on one plot plot(perf.bad,col=1) plot(perf.full, col=2, add=T) plot(perf.cv,col=3,add=T) legend(.7,.4,c("BAD","FULL","CV"), 1:3) Finally, we can calculate the AUC for each model auc.bad<-performance(pred.bad,'auc') auc.bad@y.values ## [[1]] ## [1] 0.4766734 auc.full<-performance(pred.full,"auc") auc.full@y.values ## [[1]] ## [1] 0.959432 auc.cv<-performance(pred.cv,'auc') auc.cv@y.values ## [[1]] ## [1] 0.9107505 The higher the AUC the better. As such, the full model with all variables is superior to the cross-validated or bad model. This is despite the fact that there are many high correlations in the full model as well. Another point to consider is that the cross-validated model is simpler so this may be a reason to pick it over the full model. As such, the statistics provide support for choosing a model but the do not trump the ability of the research to pick based on factors beyond just numbers. # Logistic Regression in R In this post, we will conduct a logistic regression analysis. Logistic regression is used when you want to predict a categorical dependent variable using continuous or categorical dependent variables. In our example, we want to predict Sex (male or female) when using several continuous variables from the “survey” dataset in the “MASS” package. library(MASS);library(bestglm);library(reshape2);library(corrplot) data(survey) ?MASS::survey #explains the variables in the study The first thing we need to do is remove the independent factor variables from our dataset. The reason for this is that the function that we will use for the cross-validation does not accept factors. We will first use the “str” function to identify factor variables and then remove them from the dataset. We also need to remove in examples that are missing data so we use the “na.omit” function for this. Below is the code survey$Clap<-NULL
survey$W.Hnd<-NULL survey$Fold<-NULL
survey$Exer<-NULL survey$Smoke<-NULL
survey$M.I<-NULL survey<-na.omit(survey) We now need to check for collinearity using the “corrplot.mixed” function form the “corrplot” package. pc<-cor(survey[,2:5]) corrplot.mixed(pc) corrplot.mixed(pc) We have extreme correlation between “We.Hnd” and “NW.Hnd” this makes sense because people’s hands are normally the same size. Since this blog post is a demonstration of logistic regression we will not worry about this too much. We now need to divide our dataset into a train and a test set. We set the seed for. First we need to make a variable that we call “ind” that is randomly assigns 70% of the number of rows of survey 1 and 30% 2. We then subset the “train” dataset by taking all rows that are 1’s based on the “ind” variable and we create the “test” dataset for all the rows that line up with 2 in the “ind” variable. This means our data split is 70% train and 30% test. Below is the code set.seed(123) ind<-sample(2,nrow(survey),replace=T,prob = c(0.7,0.3)) train<-survey[ind==1,] test<-survey[ind==2,] We now make our model. We use the “glm” function for logistic regression. We set the family argument to “binomial”. Next, we look at the results as well as the odds ratios. fit<-glm(Sex~.,family=binomial,train) summary(fit) ## ## Call: ## glm(formula = Sex ~ ., family = binomial, data = train) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -1.9875 -0.5466 -0.1395 0.3834 3.4443 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -46.42175 8.74961 -5.306 1.12e-07 *** ## Wr.Hnd -0.43499 0.66357 -0.656 0.512 ## NW.Hnd 1.05633 0.70034 1.508 0.131 ## Pulse -0.02406 0.02356 -1.021 0.307 ## Height 0.21062 0.05208 4.044 5.26e-05 *** ## Age 0.00894 0.05368 0.167 0.868 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 169.14 on 122 degrees of freedom ## Residual deviance: 81.15 on 117 degrees of freedom ## AIC: 93.15 ## ## Number of Fisher Scoring iterations: 6 exp(coef(fit)) ## (Intercept) Wr.Hnd NW.Hnd Pulse Height ## 6.907034e-21 6.472741e-01 2.875803e+00 9.762315e-01 1.234447e+00 ## Age ## 1.008980e+00 The results indicate that only height is useful in predicting if someone is a male or female. The second piece of code shares the odds ratios. The odds ratio tell how a one unit increase in the independent variable leads to an increase in the odds of being male in our model. For example, for every one unit increase in height there is a 1.23 increase in the odds of a particular example being male. We now need to see how well our model does on the train and test dataset. We first capture the probabilities and save them to the train dataset as “probs”. Next we create a “predict” variable and place the string “Female” in the same number of rows as are in the “train” dataset. Then we rewrite the “predict” variable by changing any example that has a probability above 0.5 as “Male”. Then we make a table of our results to see the number correct, false positives/negatives. Lastly, we calculate the accuracy rate. Below is the code. train$probs<-predict(fit, type = 'response')
train$predict<-rep('Female',123) train$predict[train$probs>0.5]<-"Male" table(train$predict,train$Sex) ## ## Female Male ## Female 61 7 ## Male 7 48 mean(train$predict==train$Sex) ## [1] 0.8861789 Despite the weaknesses of the model with so many insignificant variables it is surprisingly accurate at 88.6%. Let’s see how well we do on the “test” dataset. test$prob<-predict(fit,newdata = test, type = 'response')
test$predict<-rep('Female',46) test$predict[test$prob>0.5]<-"Male" table(test$predict,test$Sex) ## ## Female Male ## Female 17 3 ## Male 0 26 mean(test$predict==test$Sex) ## [1] 0.9347826 As you can see, we do even better on the test set with an accuracy of 93.4%. Our model is looking pretty good and height is an excellent predictor of sex which makes complete sense. However, in the next post we will use cross-validation and the ROC plot to further assess the quality of it. # Best Subset Regression in R In this post, we will take a look at best subset regression. Best subset regression fits a model for all possible feature or variable combinations and the decision for the most appropriate model is made by the analyst based on judgment or some statistical criteria. Best subset regression is an alternative to both Forward and Backward stepwise regression. Forward stepwise selection adds one variable at a time based on the lowest residual sum of squares until no more variables continues to lower the residual sum of squares. Backward stepwise regression starts with all variables in the model and removes variables one at a time. The concern with stepwise methods is they can produce biased regression coefficients, conflicting models, and inaccurate confidence intervals. Best subset regression bypasses these weaknesses of stepwise models by creating all models possible and then allowing you to assess which variables should be including in your final model. The one drawback to best subset is that a large number of variables means a large number of potential models, which can make it difficult to make a decision among several choices. In this post, we will use the “Fair” dataset from the “Ecdat” package to predict marital satisfaction based on age, Sex, the presence of children, years married, religiosity, education, occupation, and number of affairs in the past year. Below is some initial code. library(leaps);library(Ecdat);library(car);library(lmtest) data(Fair) We begin our analysis by building the initial model with all variables in it. Below is the code fit<-lm(rate~.,Fair) summary(fit) ## ## Call: ## lm(formula = rate ~ ., data = Fair) ## ## Residuals: ## Min 1Q Median 3Q Max ## -3.2049 -0.6661 0.2298 0.7705 2.2292 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3.522875 0.358793 9.819 < 2e-16 *** ## sexmale -0.062281 0.099952 -0.623 0.53346 ## age -0.009683 0.007548 -1.283 0.20005 ## ym -0.019978 0.013887 -1.439 0.15079 ## childyes -0.206976 0.116227 -1.781 0.07546 . ## religious 0.042142 0.037705 1.118 0.26416 ## education 0.068874 0.021153 3.256 0.00119 ** ## occupation -0.015606 0.029602 -0.527 0.59825 ## nbaffairs -0.078812 0.013286 -5.932 5.09e-09 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.03 on 592 degrees of freedom ## Multiple R-squared: 0.1405, Adjusted R-squared: 0.1289 ## F-statistic: 12.1 on 8 and 592 DF, p-value: 4.487e-16 The initial results are already interesting even though the r-square is low. When couples have children the have less martial satisfaction than couples without children when controlling for the other factors and this is the strongest regression weight. In addition, the more education a person has there is an increase in marital satisfaction. Lastly, as the number of affairs increases there is also a decrease in martial satisfaction. Keep in mind that the “rate” variable goes from 1 to 5 with one meaning a terrible marriage to five being a great one. The mean marital satisfaction was 3.52 when controlling for the other variables. We will now create our subset models. Below is the code. sub.fit<-regsubsets(rate~.,Fair) best.summary<-summary(sub.fit) In the code above we create the sub models using the “regsubsets” function from the “leaps” package and saved it in the variable called “sub.fit”. We then saved the summary of “sub.fit” in the variable “best.summary”. We will use the “best.summary” “sub.fit variables several times to determine which model to use. There are many different ways to assess the model. We will use the following statistical methods that come with the results from the “regsubset” function. • Mallow’ Cp • Bayesian Information Criteria We will make two charts for each of the criteria above. The plot to the left will explain how many features to include in the model. The plot to the right will tell you which variables to include. It is important to note that for both of these methods, the lower the score the better the model. Below is the code for Mallow’s Cp. par(mfrow=c(1,2)) plot(best.summary$cp)
plot(sub.fit,scale = "Cp")

The plot on the left suggests that a four feature model is the most appropriate. However, this chart does not tell me which four features. The chart on the right is read in reverse order. The high numbers are at the bottom and the low numbers are at the top when looking at the y-axis. Knowing this, we can conclude that the most appropriate variables to include in the model are age, children presence, education, and number of affairs. Below are the results using the Bayesian Information Criterion

par(mfrow=c(1,2))
plot(best.summary$bic) plot(sub.fit,scale = "bic") These results indicate that a three feature model is appropriate. The variables or features are years married, education, and number of affairs. Presence of children was not considered beneficial. Since our original model and Mallow’s Cp indicated that presence of children was significant we will include it for now. Below is the code for the model based on the subset regression. fit2<-lm(rate~age+child+education+nbaffairs,Fair) summary(fit2) ## ## Call: ## lm(formula = rate ~ age + child + education + nbaffairs, data = Fair) ## ## Residuals: ## Min 1Q Median 3Q Max ## -3.2172 -0.7256 0.1675 0.7856 2.2713 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 3.861154 0.307280 12.566 < 2e-16 *** ## age -0.017440 0.005057 -3.449 0.000603 *** ## childyes -0.261398 0.103155 -2.534 0.011531 * ## education 0.058637 0.017697 3.313 0.000978 *** ## nbaffairs -0.084973 0.012830 -6.623 7.87e-11 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.029 on 596 degrees of freedom ## Multiple R-squared: 0.1352, Adjusted R-squared: 0.1294 ## F-statistic: 23.29 on 4 and 596 DF, p-value: < 2.2e-16 The results look ok. The older a person is the less satisfied they are with their marriage. If children are present the marriage is less satisfying. The more educated the more satisfied they are. Lastly, the higher the number of affairs indicate less marital satisfaction. However, before we get excited we need to check for collinearity and homoscedasticity. Below is the code vif(fit2) ## age child education nbaffairs ## 1.249430 1.228733 1.023722 1.014338 No issues with collinearity.For vif values above 5 or 10 indicate a problem. Let’s check for homoscedasticity par(mfrow=c(2,2)) plot(fit2) The normal qqplot and residuals vs leverage plot can be used for locating outliers. The residual vs fitted and the scale-location plot do not look good as there appears to be a pattern in the dispersion which indicates homoscedasticity. To confirm this we will use Breusch-Pagan test from the “lmtest” package. Below is the code bptest(fit2) ## ## studentized Breusch-Pagan test ## ## data: fit2 ## BP = 16.238, df = 4, p-value = 0.002716 There you have it. Our model violates the assumption of homoscedasticity. However, this model was developed for demonstration purpose to provide an example of subset regression. # Data Wrangling in R Collecting and preparing data for analysis is the primary job of a data scientist. This experience is called data wrangling. In this post, we will look at an example of data wrangling using a simple artificial data set. You can create the table below in r or excel. If you created it in excel just save it as a csv and load it into r. Below is the initial code library(readr) apple <- read_csv("~/Desktop/apple.csv") ## # A tibble: 10 × 2 ## weight location ## <chr> <chr> ## 1 3.2 Europe ## 2 4.2kg europee ## 3 1.3 kg U.S. ## 4 7200 grams USA ## 5 42 United States ## 6 2.3 europee ## 7 2.1kg Europe ## 8 3.1kg USA ## 9 2700 grams U.S. ## 10 24 United States This a small dataset with the columns of “weight” and “location”. Here are some of the problems • Weights are in different units • Weights are written in different ways • Location is not consistent In order to have any success with data wrangling you need to state specifically what it is you want to do. Here are our goals for this project • Convert the “Weight variable” to a numerical variable instead of character • Remove the text and have only numbers in the “weight variable” • Change weights in grams to kilograms • Convert the “location” variable to a factor variable instead of character • Have consistent spelling for Europe and United States in the “location” variable We will begin with the “weight” variable. We want to convert it to a numerical variable and remove any non-numerical text. Below is the code for this corrected.weight<-as.numeric(gsub(pattern = "[[:alpha:]]","",apple$weight))
corrected.weight
##  [1]    3.2    4.2    1.3 7200.0   42.0    2.3    2.1    3.1 2700.0   24.0

Here is what we did.

1. We created a variable called “corrected.weight”
2. We use the function “as.numeric” this makes whatever results inside it to be a numerical variable
3. Inside “as.numeric” we used the “gsub” function which allows us to substitute one value for another.
4. Inside “gsub” we used the argument pattern and set it to “[[alpha:]]” and “” this told r to look for any lower or uppercase letters and replace with nothing or remove it. This all pertains to the “weight” variable in the apple dataframe.

We now need to convert the weights in grams to kilograms so that everything is the same unit. Below is the code

gram.error<-grep(pattern = "[[:digit:]]{4}",apple$weight) corrected.weight[gram.error]<-corrected.weight[gram.error]/1000 corrected.weight ## [1] 3.2 4.2 1.3 7.2 42.0 2.3 2.1 3.1 2.7 24.0 Here is what we did 1. We created a variable called “gram.error” 2. We used the grep function to search are the “weight” variable in the apple data frame for input that is a digit and is 4 digits in length this is what the “[[:digit:]]{4}” argument means. We do not change any values yet we just store them in “gram.error” 3. Once this information is stored in “gram.error” we use it as a subset for the “corrected.weight” variable. 4. We tell r to save into the “corrected.weight” variable any value that is changeable according to the criteria set in “gram.error” and to divided it by 1000. Dividing it by 1000 converts the value from grams to kilograms. We have completed the transformation of the “weight” and will move to dealing with the problems with the “location” variable in the “apple” dataframe. To do this we will first deal with the issues related to the values that relate to Europe and then we will deal with values related to United States. Below is the code. europe<-agrep(pattern = "europe",apple$location,ignore.case = T,max.distance = list(insertion=c(1),deletions=c(2)))
america<-agrep(pattern = "us",apple$location,ignore.case = T,max.distance = list(insertion=c(0),deletions=c(2),substitutions=0)) corrected.location<-apple$location
corrected.location[europe]<-"europe"
corrected.location[america]<-"US"
corrected.location<-gsub(pattern = "United States","US",corrected.location)
corrected.location
##  [1] "europe" "europe" "US"     "US"     "US"     "europe" "europe"
##  [8] "US"     "US"     "US"

The code is a little complicated to explain but in short We used the “agrep” function to tell r to search the “location” to look for values similar to our term “europe”. The other arguments provide some exceptions that r should change because the exceptions are close to the term europe. This process is repeated for the term “us”. We then store are the location variable from the “apple” dataframe in a new variable called “corrected.location” We then apply the two objects we made called “europe” and “america” to the “corrected.location” variable. Next we have to make some code to deal with “United States” and apply this using the “gsub” function.

We are almost done, now we combine are two variables “corrected.weight” and “corrected.location” into a new data.frame. The code is below

cleaned.apple<-data.frame(corrected.weight,corrected.location)
names(cleaned.apple)<-c('weight','location')
cleaned.apple
##    weight location
## 1     3.2   europe
## 2     4.2   europe
## 3     1.3       US
## 4     7.2       US
## 5    42.0       US
## 6     2.3   europe
## 7     2.1   europe
## 8     3.1       US
## 9     2.7       US
## 10   24.0       US

If you use the “str” function on the “cleaned.apple” dataframe you will see that “location” was automatically converted to a factor.

This looks much better especially if you compare it to the original dataframe that is printed at the top of this post.

# Principal Component Analysis in R

This post will demonstrate the use of principal component analysis (PCA). PCA is useful for several reasons. One it allows you place your examples into groups similar to linear discriminant analysis but you do not need to know beforehand what the groups are. Second, PCA is used for the purpose of dimension reduction. For example, if you have 50 variables PCA can allow you to reduce this while retaining a certain threshold of variance. If you are working with a large dataset this can greatly reduce the computational time and general complexity of your models.

Keep in mind that there really is not a dependent variable as this is unsupervised learning. What you are trying to see is how different examples can be mapped in space based on whatever independent variables are used. For our example, we will use the “Carseats” dataset form the “ISLR”. Our goal is to understanding the relationship among the variables when examining the shelve location of the car seat. Below is the initial code to begin the analysis

library(ggplot2)
library(ISLR)
data("Carseats")

We first need to rearrange the data and remove the variables we are not going to use in the analysis. Below is the code.

Carseats1<-Carseats
Carseats1<-Carseats1[,c(1,2,3,4,5,6,8,9,7,10,11)]
Carseats1$Urban<-NULL Carseats1$US<-NULL

Here is what we did 1. We made a copy of the “Carseats” data called “Careseats1” 2. We rearranged the order of the variables so that the factor variables are at the end. This will make sense later 3.We removed the “Urban” and “US” variables from the table as they will not be a part of our analysis

We will now do the PCA. We need to scale and center our data otherwise the larger numbers will have a much stronger influence on the results than smaller numbers. Fortunately, the “prcomp” function has a “scale” and a “center” argument. We will also use only the first 7 columns for the analysis  as “sheveLoc” is not useful for this analysis. If we hadn’t moved “shelveLoc” to the end of the dataframe it would cause some headache. Below is the code.

Carseats.pca<-prcomp(Carseats1[,1:7],scale. = T,center = T)
summary(Carseats.pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5     PC6     PC7
## Standard deviation     1.3315 1.1907 1.0743 0.9893 0.9260 0.80506 0.41320
## Proportion of Variance 0.2533 0.2026 0.1649 0.1398 0.1225 0.09259 0.02439
## Cumulative Proportion  0.2533 0.4558 0.6207 0.7605 0.8830 0.97561 1.00000

The summary of “Carseats.pca” Tells us how much of the variance each component explains. Keep in mind that number of components is equal to the number of variables. The “proportion of variance” tells us the contribution each component makes and the “cumulative proportion”.

If your goal is dimension reduction than the number of components to keep depends on the threshold you set. For example, if you need around 90% of the variance you would keep the first 5 components. If you need 95% or more of the variance you would keep the first six. To actually use the components you would take the “Carseats.pca$x” data and move it to your data frame. Keep in mind that the actual components have no conceptual meaning but is a numerical representation of a combination of several variables that were reduce using PCA to fewer variables such as going form 7 variables to 5 variables. This means that PCA is great for reducing variables for prediction purpose but is much harder for explanatory studies unless you can explain what the new components represent. For our purposes, we will keep 5 components. This means that we have reduce our dimensions from 7 to 5 while still keeping almost 90% of the variance. Graphing our results is tricky because we have 5 dimensions but the human mind can only conceptualize 3 at the best and normally 2. As such we will plot the first two components and label them by shelf location using ggplot2. Below is the code scores<-as.data.frame(Carseats.pca$x)
pcaplot<-ggplot(scores,(aes(PC1,PC2,color=Carseats1$ShelveLoc)))+geom_point() pcaplot From the plot you can see there is little separation when using the first two components of the PCA analysis. This makes sense as we can only graph to components so we are missing a lot of the variance. However for demonstration purposes the analysis is complete. # Linear Discriminant Analysis in R In this post we will look at an example of linear discriminant analysis (LDA). LDA is used to develop a statistical model that classifies examples in a dataset. In the example in this post, we will use the “Star” dataset from the “Ecdat” package. What we will do is try to predict the type of class the students learned in (regular, small, regular with aide) using their math scores, reading scores, and the teaching experience of the teacher. Below is the initial code library(Ecdat) library(MASS) data(Star) We first need to examine the data by using the “str” function str(Star) ## 'data.frame': 5748 obs. of 8 variables: ##$ tmathssk: int  473 536 463 559 489 454 423 500 439 528 ...
##  $treadssk: int 447 450 439 448 447 431 395 451 478 455 ... ##$ classk  : Factor w/ 3 levels "regular","small.class",..: 2 2 3 1 2 1 3 1 2 2 ...
##  $totexpk : int 7 21 0 16 5 8 17 3 11 10 ... ##$ sex     : Factor w/ 2 levels "girl","boy": 1 1 2 2 2 2 1 1 1 1 ...
##  $freelunk: Factor w/ 2 levels "no","yes": 1 1 2 1 2 2 2 1 1 1 ... ##$ race    : Factor w/ 3 levels "white","black",..: 1 2 2 1 1 1 2 1 2 1 ...
##  $schidkn : int 63 20 19 69 79 5 16 56 11 66 ... ## - attr(*, "na.action")=Class 'omit' Named int [1:5850] 1 4 6 7 8 9 10 15 16 17 ... ## .. ..- attr(*, "names")= chr [1:5850] "1" "4" "6" "7" ... We will use the following variables • dependent variable = classk (class type) • independent variable = tmathssk (Math score) • independent variable = treadssk (Reading score) • independent variable = totexpk (Teaching experience) We now need to examine the data visually by looking at histograms for our independent variables and a table for our dependent variable hist(Star$tmathssk)

hist(Star$treadssk) hist(Star$totexpk)

prop.table(table(Star$classk)) ## ## regular small.class regular.with.aide ## 0.3479471 0.3014962 0.3505567 The data mostly looks good. The results of the “prop.table” function will help us when we develop are training and testing datasets. The only problem is with the “totexpk” variable. IT is not anywhere near to be normally distributed. TO deal with this we will use the square root for teaching experience. Below is the code star.sqrt<-Star star.sqrt$totexpk.sqrt<-sqrt(star.sqrt$totexpk) hist(sqrt(star.sqrt$totexpk))

Much better. We now need to check the correlation among the variables as well and we will use the code below.

cor.star<-data.frame(star.sqrt$tmathssk,star.sqrt$treadssk,star.sqrt$totexpk.sqrt) cor(cor.star) ## star.sqrt.tmathssk star.sqrt.treadssk ## star.sqrt.tmathssk 1.00000000 0.7135489 ## star.sqrt.treadssk 0.71354889 1.0000000 ## star.sqrt.totexpk.sqrt 0.08647957 0.1045353 ## star.sqrt.totexpk.sqrt ## star.sqrt.tmathssk 0.08647957 ## star.sqrt.treadssk 0.10453533 ## star.sqrt.totexpk.sqrt 1.00000000 None of the correlations are too bad. We can now develop our model using linear discriminant analysis. First, we need to scale are scores because the test scores and the teaching experience are measured differently. Then, we need to divide our data into a train and test set as this will allow us to determine the accuracy of the model. Below is the code. star.sqrt$tmathssk<-scale(star.sqrt$tmathssk) star.sqrt$treadssk<-scale(star.sqrt$treadssk) star.sqrt$totexpk.sqrt<-scale(star.sqrt$totexpk.sqrt) train.star<-star.sqrt[1:4000,] test.star<-star.sqrt[4001:5748,] Now we develop our model. In the code before the “prior” argument indicates what we expect the probabilities to be. In our data the distribution of the the three class types is about the same which means that the apriori probability is 1/3 for each class type. train.lda<-lda(classk~tmathssk+treadssk+totexpk.sqrt, data = train.star,prior=c(1,1,1)/3) train.lda ## Call: ## lda(classk ~ tmathssk + treadssk + totexpk.sqrt, data = train.star, ## prior = c(1, 1, 1)/3) ## ## Prior probabilities of groups: ## regular small.class regular.with.aide ## 0.3333333 0.3333333 0.3333333 ## ## Group means: ## tmathssk treadssk totexpk.sqrt ## regular -0.04237438 -0.05258944 -0.05082862 ## small.class 0.13465218 0.11021666 -0.02100859 ## regular.with.aide -0.05129083 -0.01665593 0.09068835 ## ## Coefficients of linear discriminants: ## LD1 LD2 ## tmathssk 0.89656393 -0.04972956 ## treadssk 0.04337953 0.56721196 ## totexpk.sqrt -0.49061950 0.80051026 ## ## Proportion of trace: ## LD1 LD2 ## 0.7261 0.2739 The printout is mostly readable. At the top is the actual code used to develop the model followed by the probabilities of each group. The next section shares the means of the groups. The coefficients of linear discriminants are the values used to classify each example. The coefficients are similar to regression coefficients. The computer places each example in both equations and probabilities are calculated. Whichever class has the highest probability is the winner. In addition, the higher the coefficient the more weight it has. For example, “tmathssk” is the most influential on LD1 with a coefficient of 0.89. The proportion of trace is similar to principal component analysis Now we will take the trained model and see how it does with the test set. We create a new model called “predict.lda” and use are “train.lda” model and the test data called “test.star” predict.lda<-predict(train.lda,newdata = test.star) We can use the “table” function to see how well are model has done. We can do this because we actually know what class our data is beforehand because we divided the dataset. What we need to do is compare this to what our model predicted. Therefore, we compare the “classk” variable of our “test.star” dataset with the “class” predicted by the “predict.lda” model. table(test.star$classk,predict.lda$class) ## ## regular small.class regular.with.aide ## regular 155 182 249 ## small.class 145 198 174 ## regular.with.aide 172 204 269 The results are pretty bad. For example, in the first row called “regular” we have 155 examples that were classified as “regular” and predicted as “regular” by the model. In rhe next column, 182 examples that were classified as “regular” but predicted as “small.class”, etc. To find out how well are model did you add together the examples across the diagonal from left to right and divide by the total number of examples. Below is the code (155+198+269)/1748 ## [1] 0.3558352 Only 36% accurate, terrible but ok for a demonstration of linear discriminant analysis. Since we only have two-functions or two-dimensions we can plot our model. Below I provide a visual of the first 50 examples classified by the predict.lda model. plot(predict.lda$x[1:50])
text(predict.lda$x[1:50],as.character(predict.lda$class[1:50]),col=as.numeric(predict.lda$class[1:100])) abline(h=0,col="blue") abline(v=0,col="blue") The first function, which is the vertical line, doesn’t seem to discriminant anything as it off to the side and not separating any of the data. However, the second function, which is the horizontal one, does a good of dividing the “regular.with.aide” from the “small.class”. Yet, there are problems with distinguishing the class “regular” from either of the other two groups. In order improve our model we need additional independent variables to help to distinguish the groups in the dependent variable. # Generalized Additive Models in R In this post, we will learn how to create a generalized additive model (GAM). GAMs are non-parametric generalized linear models. This means that linear predictor of the model uses smooth functions on the predictor variables. As such, you do not need to specific the functional relationship between the response and continuous variables. This allows you to explore the data for potential relationships that can be more rigorously tested with other statistical models In our example, we will use the “Auto” dataset from the “ISLR” package and use the variables “mpg”,“displacement”,“horsepower”,and “weight” to predict “acceleration”. We will also use the “mgcv” package. Below is some initial code to begin the analysis library(mgcv) library(ISLR) data(Auto) We will now make the model we want to understand the response of “accleration” to the explanatory variables of “mpg”,“displacement”,“horsepower”,and “weight”. After setting the model we will examine the summary. Below is the code model1<-gam(acceleration~s(mpg)+s(displacement)+s(horsepower)+s(weight),data=Auto) summary(model1) ## ## Family: gaussian ## Link function: identity ## ## Formula: ## acceleration ~ s(mpg) + s(displacement) + s(horsepower) + s(weight) ## ## Parametric coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 15.54133 0.07205 215.7 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df F p-value ## s(mpg) 6.382 7.515 3.479 0.00101 ** ## s(displacement) 1.000 1.000 36.055 4.35e-09 *** ## s(horsepower) 4.883 6.006 70.187 < 2e-16 *** ## s(weight) 3.785 4.800 41.135 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## R-sq.(adj) = 0.733 Deviance explained = 74.4% ## GCV = 2.1276 Scale est. = 2.0351 n = 392 All of the explanatory variables are significant and the adjust r-squared is .73 which is excellent. edf stands for “effective degrees of freedom”. This modified version of the degree of freedoms is due to the smoothing process in the model. GCV stands for generalized cross validation and this number is useful when comparing models. The model with the lowest number is the better model. We can also examine the model visually by using the “plot” function. This will allow us to examine if the curvature fitted by the smoothing process was useful or not for each variable. Below is the code. plot(model1) We can also look at a 3d graph that includes the linear predictor as well as the two strongest predictors. This is done with the “vis.gam” function. Below is the code vis.gam(model1) If multiple models are developed. You can compare the GCV values to determine which model is the best. In addition, another way to compare models is with the “AIC” function. In the code below, we will create an additional model that includes “year” compare the GCV scores and calculate the AIC. Below is the code. model2<-gam(acceleration~s(mpg)+s(displacement)+s(horsepower)+s(weight)+s(year),data=Auto) summary(model2) ## ## Family: gaussian ## Link function: identity ## ## Formula: ## acceleration ~ s(mpg) + s(displacement) + s(horsepower) + s(weight) + ## s(year) ## ## Parametric coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 15.54133 0.07203 215.8 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df F p-value ## s(mpg) 5.578 6.726 2.749 0.0106 * ## s(displacement) 2.251 2.870 13.757 3.5e-08 *** ## s(horsepower) 4.936 6.054 66.476 < 2e-16 *** ## s(weight) 3.444 4.397 34.441 < 2e-16 *** ## s(year) 1.682 2.096 0.543 0.6064 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## R-sq.(adj) = 0.733 Deviance explained = 74.5% ## GCV = 2.1368 Scale est. = 2.0338 n = 392 #model1 GCV model1$gcv.ubre
##   GCV.Cp
## 2.127589
#model2 GCV
model2$gcv.ubre ## GCV.Cp ## 2.136797 As you can see, the second model has a higher GCV score when compared to the first model. This indicates that the first model is a better choice. This makes sense because in the second model the variable “year” is not significant. To confirm this we will calculate the AIC scores using the AIC function. AIC(model1,model2) ## df AIC ## model1 18.04952 1409.640 ## model2 19.89068 1411.156 Again, you can see that model1 s better due to its fewer degrees of freedom and slightly lower AIC score. Conclusion Using GAMs is most common for exploring potential relationships in your data. This is stated because they are difficult to interpret and to try and summarize. Therefore, it is normally better to develop a generalized linear model over a GAM due to the difficulty in understanding what the data is trying to tell you when using GAMs. # Generalized Models in R Generalized linear models are another way to approach linear regression. The advantage of of GLM is that allows the error to follow many different distributions rather than only the normal distribution which is an assumption of traditional linear regression. Often GLM is used for response or dependent variables that are binary or represent count data. THis post will provide a brief explanation of GLM as well as provide an example. Key Information There are three important components to a GLM and they are • Error structure • Linear predictor • Link function The error structure is the type of distribution you will use in generating the model. There are many different distributions in statistical modeling such as binomial, gaussian, poission, etc. Each distribution comes with certain assumptions that govern their use. The linear predictor is the sum of the effects of the independent variables. Lastly, the link function determines the relationship between the linear predictor and the mean of the dependent variable. There are many different link functions and the best link function is the one that reduces the residual deviances the most. In our example, we will try to predict if a house will have air conditioning based on the interactioon between number of bedrooms and bathrooms, number of stories, and the price of the house. To do this, we will use the “Housing” dataset from the “Ecdat” package. Below is some initial code to get started. library(Ecdat) data("Housing") The dependent variable “airco” in the “Housing” dataset is binary. This calls for us to use a GLM. To do this we will use the “glm” function in R. Furthermore, in our example, we want to determine if there is an interaction between number of bedrooms and bathrooms. Interaction means that the two independent variables (bathrooms and bedrooms) influence on the dependent variable (aircon) is not additive, which means that the combined effect of the independnet variables is different than if you just added them together. Below is the code for the model followed by a summary of the results model<-glm(Housing$airco ~ Housing$bedrooms * Housing$bathrms + Housing$stories + Housing$price, family=binomial)
summary(model)
##
## Call:
## glm(formula = Housing$airco ~ Housing$bedrooms * Housing$bathrms + ## Housing$stories + Housing$price, family = binomial) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.7069 -0.7540 -0.5321 0.8073 2.4217 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -6.441e+00 1.391e+00 -4.632 3.63e-06 ## Housing$bedrooms                  8.041e-01  4.353e-01   1.847   0.0647
## Housing$bathrms 1.753e+00 1.040e+00 1.685 0.0919 ## Housing$stories                   3.209e-01  1.344e-01   2.388   0.0170
## Housing$price 4.268e-05 5.567e-06 7.667 1.76e-14 ## Housing$bedrooms:Housing$bathrms -6.585e-01 3.031e-01 -2.173 0.0298 ## ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 681.92 on 545 degrees of freedom ## Residual deviance: 549.75 on 540 degrees of freedom ## AIC: 561.75 ## ## Number of Fisher Scoring iterations: 4 To check how good are model is we need to check for overdispersion as well as compared this model to other potential models. Overdispersion is a measure to determine if there is too much variablity in the model. It is calcualted by dividing the residual deviance by the degrees of freedom. Below is the solution for this 549.75/540 ## [1] 1.018056 Our answer is 1.01, which is pretty good because the cutoff point is 1, so we are really close. Now we will make several models and we will compare the results of them Model 2 #add recroom and garagepl model2<-glm(Housing$airco ~ Housing$bedrooms * Housing$bathrms + Housing$stories + Housing$price + Housing$recroom + Housing$garagepl, family=binomial)
summary(model2)
##
## Call:
## glm(formula = Housing$airco ~ Housing$bedrooms * Housing$bathrms + ## Housing$stories + Housing$price + Housing$recroom + Housing$garagepl, ## family = binomial) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.6733 -0.7522 -0.5287 0.8035 2.4239 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -6.369e+00 1.401e+00 -4.545 5.51e-06 ## Housing$bedrooms                  7.830e-01  4.391e-01   1.783   0.0745
## Housing$bathrms 1.702e+00 1.047e+00 1.626 0.1039 ## Housing$stories                   3.286e-01  1.378e-01   2.384   0.0171
## Housing$price 4.204e-05 6.015e-06 6.989 2.77e-12 ## Housing$recroomyes                1.229e-01  2.683e-01   0.458   0.6470
## Housing$garagepl 2.555e-03 1.308e-01 0.020 0.9844 ## Housing$bedrooms:Housing$bathrms -6.430e-01 3.054e-01 -2.106 0.0352 ## ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 681.92 on 545 degrees of freedom ## Residual deviance: 549.54 on 538 degrees of freedom ## AIC: 565.54 ## ## Number of Fisher Scoring iterations: 4 #overdispersion calculation 549.54/538 ## [1] 1.02145 Model 3 model3<-glm(Housing$airco ~ Housing$bedrooms * Housing$bathrms + Housing$stories + Housing$price + Housing$recroom + Housing$fullbase + Housing$garagepl, family=binomial) summary(model3) ## ## Call: ## glm(formula = Housing$airco ~ Housing$bedrooms * Housing$bathrms +
##     Housing$stories + Housing$price + Housing$recroom + Housing$fullbase +
##     Housing$garagepl, family = binomial) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -2.6629 -0.7436 -0.5295 0.8056 2.4477 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -6.424e+00 1.409e+00 -4.559 5.14e-06 ## Housing$bedrooms                  8.131e-01  4.462e-01   1.822   0.0684
## Housing$bathrms 1.764e+00 1.061e+00 1.662 0.0965 ## Housing$stories                   3.083e-01  1.481e-01   2.082   0.0374
## Housing$price 4.241e-05 6.106e-06 6.945 3.78e-12 ## Housing$recroomyes                1.592e-01  2.860e-01   0.557   0.5778
## Housing$fullbaseyes -9.523e-02 2.545e-01 -0.374 0.7083 ## Housing$garagepl                 -1.394e-03  1.313e-01  -0.011   0.9915
## Housing$bedrooms:Housing$bathrms -6.611e-01  3.095e-01  -2.136   0.0327
##
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
##     Null deviance: 681.92  on 545  degrees of freedom
## Residual deviance: 549.40  on 537  degrees of freedom
## AIC: 567.4
##
## Number of Fisher Scoring iterations: 4
#overdispersion calculation
549.4/537
## [1] 1.023091

Now we can assess the models by using the “anova” function with the “test” argument set to “Chi” for the chi-square test.

anova(model, model2, model3, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: Housing$airco ~ Housing$bedrooms * Housing$bathrms + Housing$stories +
##     Housing$price ## Model 2: Housing$airco ~ Housing$bedrooms * Housing$bathrms + Housing$stories + ## Housing$price + Housing$recroom + Housing$garagepl
## Model 3: Housing$airco ~ Housing$bedrooms * Housing$bathrms + Housing$stories +
##     Housing$price + Housing$recroom + Housing$fullbase + Housing$garagepl
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1       540     549.75
## 2       538     549.54  2  0.20917   0.9007
## 3       537     549.40  1  0.14064   0.7076

The results of the anova indicate that the models are all essentially the same as there is no statistical difference. The only criteria on which to select a model is the measure of overdispersion. The first model has the lowest rate of overdispersion and so is the best when using this criteria. Therefore, determining if a hous has air conditioning depends on examining number of bedrooms and bathrooms simultenously as well as the number of stories and the price of the house.

Conclusion

The post explained how to use and interpret GLM in R. GLM can be used primarilyy for fitting data to disrtibutions that are not normal.

# Proportion Test in R

Proportions are are a fraction or “portion” of a total amount. For example, if there are ten men and ten women in a room the proportion of men in the room is 50% (5 / 10). There are times when doing an analysis that you want to evaluate proportions in our data rather than individual measurements of mean, correlation, standard deviation etc.

In this post we will learn how to do a test of proportions using R. We will use the dataset “Default” which is found in the “ISLR” pacakage. We will compare the proportion of those who are students in the dataset to a theoretical value. We will calculate the results using the z-test and the binomial exact test. Below is some initial code to get started.

library(ISLR)
data("Default")

We first need to determine the actual number of students that are in the sample. This is calculated below using the “table” function.

table(Default$student) ## ## No Yes ## 7056 2944 We have 2944 students in the sample and 7056 people who are not students. We now need to determine how many people are in the sample. If we sum the results from the table below is the code. sum(table(Default$student))
## [1] 10000

There are 10000 people in the sample. To determine the proprtion of students we take the number 2944 / 10000 which equals 29.44 or 29.44%. Below is the code to calculate this

table(Default$student) / sum(table(Default$student))
##
##     No    Yes
## 0.7056 0.2944

The proportion test is used to compare a particular value with a theoretical value. For our example, the particular value we have is 29.44% of the people were students. We want to compare this value with a theoretical value of 50%. Before we do so it is better to state specificallt what are hypotheses are. NULL = The value of 29.44% of the sample being students is the same as 50% found in the population ALTERNATIVE = The value of 29.44% of the sample being students is NOT the same as 50% found in the population.

Below is the code to complete the z-test.

prop.test(2944,n = 10000, p = 0.5, alternative = "two.sided", correct = FALSE)
##
##  1-sample proportions test without continuity correction
##
## data:  2944 out of 10000, null probability 0.5
## X-squared = 1690.9, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.2855473 0.3034106
## sample estimates:
##      p
## 0.2944

Here is what the code means. 1. prop.test is the function used 2. The first value of 2944 is the total number of students in the sample 3. n = is the sample size 4. p= 0.5 is the theoretical proportion 5. alternative =“two.sided” means we want a two-tail test 6. correct = FALSE means we do not want a correction applied to the z-test. This is useful for small sample sizes but not for our sample of 10000

The p-value is essentially zero. This means that we reject the null hypothesis and conclude that the proprtion of students in our sample is different from a theortical proprition of 50% in the population.

Below is the same analysis using the binomial exact test.

binom.test(2944, n = 10000, p = 0.5)
##
##  Exact binomial test
##
## data:  2944 and 10000
## number of successes = 2944, number of trials = 10000, p-value <
## 2.2e-16
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.2854779 0.3034419
## sample estimates:
## probability of success
##                 0.2944

The results are the same. Whether to use the “prop.test”” or “binom.test” is a major argument among statisticians. The purpose here was to provide an example of the use of both

# Theoretical Distribution and R

This post will explore an example of testing if a dataset fits a specific theoretical distribution. This is a very important aspect of statistical modeling as it allows to understand the normality of the data and the appropriate steps needed to take to prepare for analysis.

In our example, we will use the “Auto” dataset from the “ISLR” package. We will check if the horsepower of the cars in the dataset is normally distributed or not. Below is some initial code to begin the process.

library(ISLR)
library(nortest)
library(fBasics)
data("Auto")

Determining if a dataset is normally distributed is simple in R. This is normally done visually through making a Quantile-Quantile plot (Q-Q plot). It involves using two functions the “qnorm” and the “qqline”. Below is the code for the Q-Q plot

qqnorm(Auto$horsepower) We now need to add the Q-Q line to see how are distribution lines up with the theoretical normal one. Below is the code. Note that we have to repeat the code above in order to get the completed plot. qqnorm(Auto$horsepower)
qqline(Auto$horsepower, distribution = qnorm, probs=c(.25,.75)) The “qqline” function needs the data you want to test as well as the distribution and probability. The distribution we wanted is normal and is indicated by the argument “qnorm”. The probs argument means probability. The default values are .25 and .75. The resulting graph indicates that the distribution of “horsepower”, in the “Auto” dataset is not normally distributed. That are particular problems with the lower and upper values. We can confirm our suspicion by running a statistical test. The Anderson-Darling test from the “nortest” package will allow us to test whether our data is normally distributed or not. The code is below ad.test(Auto$horsepower)
##  Anderson-Darling normality test
##
## data:  Auto$horsepower ## A = 12.675, p-value < 2.2e-16 From the results, we can conclude that the data is not normally distributed. This could mean that we may need to use non-parametric tools for statistical analysis. We can further explore our distribution in terms of its skew and kurtosis. Skew measures how far to the left or right the data leans and kurtosis measures how peaked or flat the data is. This is done with the “fBasics” package and the functions “skewness” and “kurtosis”. First we will deal with skewness. Below is the code for calculating skewness. horsepowerSkew<-skewness(Auto$horsepower)
horsepowerSkew
## [1] 1.079019
## attr(,"method")
## [1] "moment"

We now need to determine if this value of skewness is significantly different from zero. This is done with a simple t-test. We must calculate the t-value before calculating the probability. The standard error of the skew is defined as the square root of six divided by the total number of samples. The code is below

stdErrorHorsepower<-horsepowerSkew/(sqrt(6/length(Auto$horsepower))) stdErrorHorsepower ## [1] 8.721607 ## attr(,"method") ## [1] "moment" Now we take the standard error of Horsepower and plug this into the “pt” function (t probability) with the degrees of freedom (sample size – 1 = 391) we also put in the number 1 and subtract all of this information. Below is the code 1-pt(stdErrorHorsepower,391) ## [1] 0 ## attr(,"method") ## [1] "moment" The value zero means that we reject the null hypothesis that the skew is not significantly different form zero and conclude that the skew is different form zero. However, the value of the skew was only 1.1 which is not that non-normal. We will now repeat this process for the kurtosis. The only difference is that instead of taking the square root divided by six we divided by 24 in the example below. horsepowerKurt<-kurtosis(Auto$horsepower)
horsepowerKurt
## [1] 0.6541069
## attr(,"method")
## [1] "excess"
stdErrorHorsepowerKurt<-horsepowerKurt/(sqrt(24/length(Auto$horsepower))) stdErrorHorsepowerKurt ## [1] 2.643542 ## attr(,"method") ## [1] "excess" 1-pt(stdErrorHorsepowerKurt,391) ## [1] 0.004267199 ## attr(,"method") ## [1] "excess" Again the pvalue is essentially zero, which means that the kurtosis is significantly different from zero. With a value of 2.64 this is not that bad. However, when both skew and kurtosis are non-normally it explains why our overall distributions was not normal either. Conclusion This post provided insights into assessing the normality of a dataset. Visually inspection can take place using Q-Q plots. Statistical inspection can be done through hypothesis testing along with checking skew and kurtosis. # Probability Distribution and Graphs in R In this post, we will use probability distributions and ggplot2 in R to solve a hypothetical example. This provides a practical example of the use of R in everyday life through the integration of several statistical and coding skills. Below is the scenario. At a busing company the average number of stops for a bus is 81 with a standard deviation of 7.9. The data is normally distributed. Knowing this complete the following. • Calculate the interval value to use using the 68-95-99.7 rule • Calculate the density curve • Graph the normal curve • Evaluate the probability of a bus having less then 65 stops • Evaluate the probability of a bus having more than 93 stops Calculate the Interval Value Our first step is to calculate the interval value. This is the range in which 99.7% of the values falls within. Doing this requires knowing the mean and the standard deviation and subtracting/adding the standard deviation as it is multiplied by three from the mean. Below is the code for this. busStopMean<-81 busStopSD<-7.9 busStopMean+3*busStopSD ## [1] 104.7 busStopMean-3*busStopSD ## [1] 57.3 The values above mean that we can set are interval between 55 and 110 with 100 buses in the data. Below is the code to set the interval. interval<-seq(55,110, length=100) #length here represents 100 fictitious buses Density Curve The next step is to calculate the density curve. This is done with our knowledge of the interval, mean, and standard deviation. We also need to use the “dnorm” function. Below is the code for this. densityCurve<-dnorm(interval,mean=81,sd=7.9) We will now plot the normal curve of our data using ggplot. Before we need to put our “interval” and “densityCurve” variables in a dataframe. We will call the dataframe “normal” and then we will create the plot. Below is the code. library(ggplot2) normal<-data.frame(interval, densityCurve) ggplot(normal, aes(interval, densityCurve))+geom_line()+ggtitle("Number of Stops for Buses") Probability Calculation We now want to determine what is the provability of a bus having less than 65 stops. To do this we use the “pnorm” function in R and include the value 65, along with the mean, standard deviation, and tell R we want the lower tail only. Below is the code for completing this. pnorm(65,mean = 81,sd=7.9,lower.tail = TRUE) ## [1] 0.02141744 As you can see, at 2% it would be unusually to. We can also plot this using ggplot. First, we need to set a different density curve using the “pnorm” function. Combine this with our “interval” variable in a dataframe and then use this information to make a plot in ggplot2. Below is the code. CumulativeProb<-pnorm(interval, mean=81,sd=7.9,lower.tail = TRUE) pnormal<-data.frame(interval, CumulativeProb) ggplot(pnormal, aes(interval, CumulativeProb))+geom_line()+ggtitle("Cumulative Density of Stops for Buses") Second Probability Problem We will now calculate the probability of a bus have 93 or more stops. To make it more interesting we will create a plot that shades the area under the curve for 93 or more stops. The code is a little to complex to explain so just enjoy the visual. pnorm(93,mean=81,sd=7.9,lower.tail = FALSE) ## [1] 0.06438284 x<-interval ytop<-dnorm(93,81,7.9) MyDF<-data.frame(x=x,y=densityCurve) p<-ggplot(MyDF,aes(x,y))+geom_line()+scale_x_continuous(limits = c(50, 110)) +ggtitle("Probabilty of 93 Stops or More is 6.4%") shade <- rbind(c(93,0), subset(MyDF, x > 93), c(MyDF[nrow(MyDF), "X"], 0)) p + geom_segment(aes(x=93,y=0,xend=93,yend=ytop)) + geom_polygon(data = shade, aes(x, y)) Conclusion A lot of work was done but all in a practical manner. Looking at realistic problem. We were able to calculate several different probabilities and graph them accordingly. # Using Maps in ggplot2 It seems as though there are no limits to what can be done with ggplot2. Another example of this is the use of maps in presenting data. If you are trying to share information that depends on location then this is an important feature to understand. This post will provide some basic explanation for understanding how to use maps with ggplot2. The Maps Package One of several packages available for using maps with ggplot2 is the “maps” package. This package contains a limited number of maps along with several databases that contain information that can be used to create data-filled maps. The “maps” package cooperates with ggplot2 through the use of the “borders” function and plotting the plot using lattitude and longitude for the “aes” function. After you have installed the “maps” package you can run the example code below. library(ggplot2);library(maps) ggplot(us.cities,aes(long,lat))+geom_point()+borders("state") In the code above we told R to use the data from “us.cities” which comes with the “maps” package. We then told R to graph the latitude and longitude and to do this by placing a point for each city. Lastly, the “borders” function was use to place this information on the state map of the US. There are several points way off of the map. These represents datapoints for cities in Alaska and Hawaii. Below is an example that is limited to one state in America. To do this we first must subset the data to only include one state. tx_cities<-subset(us.cities,country.etc=="TX") ggplot(tx_cities,aes(long,lat))+geom_point()+borders(database = "state",regions = "texas") The map shows all the cities in the state of Texas that are pulled form the “us.cities” dataset. We can also play with the colors of the maps just like any other ggplot2 output. Below is an example. data("world.cities") Thai_cities<-subset(world.cities, country.etc=="Thailand") ggplot(Thai_cities,aes(long,lat))+borders("world","Thailand", fill="light blue",col="dark blue")+geom_point(aes(size=pop),col="dark red") In the example above, we took all of the cities in Thailand and saved them into the variable “Thai_cities”. We then made a plot of Thailand but we played with the color and fill features. Lastly, we plotted the population be location and we indicated that the size of the data point should depend on the size. In this example, all the data points were the same size which means that all the cities in Thailand in the dataset are about the same size. We can also add text to maps. In the example below, we will use a subset of the data from Thailand and add the names of cities to the map. Big_Thai_cities<-subset(Thai_cities, pop>100000) ggplot(Big_Thai_cities,aes(long,lat))+borders("world","Thailand", fill="light blue",col="dark blue")+geom_point(aes(size=pop),col="dark red")+geom_text(aes(long,lat,label=name),hjust=-.2,size=3) In this plot there is a messy part in the middle where Bangkok is a long with several other large cities. However, you can see the flexiability in the plot by adding the “geom_text” function which has been discussed previously. In the “geom_text” function we added some aesthetics as well add the “name” of the city. Conclusion In this post, we look at some of the basic was of using maps with ggplot2. There are many more ways and features that can be explored in future post. # Axis and Title Modifications in ggplot2 This post will provide explanation on how to customize the axis and title of a plot that utilizes ggplot2. We will use the “Computer” dataset from the “Ecdat” package looking specifically at the difference in price of computers based on the inclusion of a cd-rom. Below is some code needed to be prepared for the examples along with a printout of our initial boxplot. library(ggplot2);library(grid);library("Ecdat") data("Computers") theBoxplot<-ggplot(Computers,aes(cd, price, fill=cd))+geom_boxplot() theBoxplot In the example below, we change the color of the tick marks to purple and we bold them. This all involves the use of the “axis.text” argument in the “theme” function. theBoxplot + theme(axis.text=element_text(color="purple",face="bold")) In the example below, the y label “price” is rotated 90 degrees to be in line with text. This is accomplished using the “axis.title.y” argument along with additional code. theBoxplot+theme(axis.title.y=element_text(size=rel(1.5),angle = 0)) Below is an example that includes a title with a change to the default size and color theBoxplot+labs(title="The Example")+theme(plot.title=element_text(size=rel(1.5),color="orange")) You can also remove the axis label. IN the example below, we remove the x axis along with its tick marks. theBoxplot+theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.title.x=element_blank()) It is also possible to modify the plot background axis as well. In the example below, we change the background color to blue, the color of the lines to green, and yellow. This is not an attractive plot but it does provide an example of the various options available in ggplot2 theBoxplot+theme(panel.background=element_rect(fill="blue"), panel.grid.major=element_line(color="green", size = 3),panel.grid.minor=element_line(color="yellow",linetype="solid",size=2)) All of the tricks we have discussed so far can also apply when faceting data. Below we make a scatterplot using the same background as before but comparing trend and price. theScatter<-ggplot(Computers,aes(trend, price, color=cd))+geom_point() theScatter1<-theScatter+facet_grid(.~cd)+theme(panel.background=element_rect(fill="blue"), panel.grid.major=element_line(color="green", size = 3),panel.grid.minor=element_line(color="yellow",linetype="solid",size=2)) theScatter1 Right now the plots are too close to each other. We can account for this by modifying the panel margins. theScatter1 +theme(panel.margin=unit(2,"cm")) Conclusion These examples provide further evidence of the endless variety that is available when using ggplot2. Whatever are your purposes, it is highly probably that ggplot2 has some sort of a data visualization answer. # Modifying Legends in ggplot2 This post will provide information on fine tuning the legend of a graph using ggplot2. We will be using the “Wage” dataset from the “ISLR” package. Below is some initial code that is needed to complete the examples. The initial plot is saved as a variable to save time and avoid repeating the same code. library(ggplot2);library(ISLR); library(grid) myBoxplot<-ggplot(Wage, aes(education, wage,fill=education))+geom_boxplot() myBoxplot The default ggplot has a grey background with grey text. By adding the “theme_bw” function to a plot you can create a plot that has a white background with black text. The code is below. myBoxplot+theme_bw() If you desire, you can also add a rectangle around the legend with the “legend.baclground” argument You can even specify the color of the rectangle as shown below. myBoxplot+theme(legend.background=element_rect(color="blue")) It is also possible to add a highlighting color to the keys in the legend. In the code below we highlight the keys with the color red using the “legend.key” argument myBoxplot+theme(legend.key=element_rect(fill="red")) The code below provides an example of how to change the size of a plot. myBoxplot+theme(legend.margin= unit(2, "cm")) This example demonstrate how to modify the text in a legend. This requires the use of the “legend.text”, along with several other arguments and functions. The code below does the following. • Size 15 font • Dark red font color • Text at 35 degree angle • Italic font myBoxplot + theme(legend.text = element_text(size = 15,color="dark red",angle= 35, face="italic")) Lastly, you can even move the legend around the plot. The first example moves the legend to the top of the plot using “legend.position” argument. The second example moves the legend based on numerical input. The first number moves the plot from left to right or from 0 being left to 1 being all the way to the right. The second number moves the text from bottom to top with 0 being the bottom and 1 being the top. myBoxplot+theme(legend.position="top") myBoxplot+theme(legend.position=c(.6,.7)) Conclusion The examples provided here show how much control over plots is possible when using ggplot2. In many ways this is just an introduction into the nuance controlled that is available # Axis and Labels in ggplot2 In this post, we will look at how to manipulate the labels and positioning of the data when using ggplot2. We will use the “Wage” data from the “ISLR” package. Below is initial code needed to begin. library(ggplot2);library(ISLR) data("Wage") Manipulating Labels Our first example involves adding labels for the x, y axis as well as a title. To do this we will create a histgram of the wage variable and save it as a variable in R. By saving the histogram as a variable it saves time as we do not have to recreate all of the code but only add the additional information. After creating the histogram and saving it to a variable we will add the code for creating the labels. Below is the code myHistogram<-ggplot(Wage, aes(wage, fill=..count..))+geom_histogram() myHistogram+labs(title="This is My Histogram", x="Salary as a Wage", y="Number") By using the “labs” function you can add a title and information for the x and y axis. If your title is really long you can use the code “” to break the information into separate lines as shown below. myHistogram+labs(title="This is the Longest Title for a Histogram \n that I have ever Seen in My Entire Life", x="Salary as a Wage", y="Number") Discrete Axis Scale We will now turn our attention to working with discrete scales. Discrete scales deal with categorical data such as boxplots and bar charts. First, we will store a boxplot of the wages subsetted by level of education in a variable and we will display it. myBoxplot<-ggplot(Wage, aes(education, wage,fill=education))+geom_boxplot() myBoxplot Now, by using the “scale_x_discrete” function along with the “limits” argument we are able to change the order of the gorups as shown below myBoxplot+scale_x_discrete(limits=c("5. Advanced Degree","2. HS Grad","1. < HS Grad","4. College Grad","3. Some College")) Continuous Scale The most common modification to a continuous scale is to modify the range. In the code below, we change the default range of “myBoxplot” to something that is larger. myBoxplot+scale_y_continuous(limits=c(0,400)) Conclusion This post provided some basic insights into modifiying plots using ggplot2. # Pie Charts and More Using ggplot2 This post will explain several types of visuals that can be developed in using ggplot2. In particular, we are going to make three specific types of charts and they are… • Pie chart • Bullseye chart • Coxcomb diagram To complete this ask, we will use the “Wage” dataset from the “ISLR” package. We will nbe using the “education” variable which has five factors in it. Below is the initial code to get started. library(ggplot2);library(ISLR) data("Wage") Pie Chart In order to make a pie chart, we first need to make a bar chart and add several pieces of code to change it into a pie chart. Below is the code for making a regular bar plot. ggplot(Wage, aes(education, fill=education))+geom_bar() We will now modify two parts of the code. First, we do not want separate bars. Instead we want one bar. The reason being is that we only want one pie chart so before that we need one bar. Therefore, for the x value in the “aes” function we will use the argument “factor(1)” which tells R to force the data as one factor on the chart thus making one bar. We also need to add the “width=1” inside the “geom_bar” function. This helps with spacing. Below is the code for this ggplot(Wage, aes(factor(1), fill=education))+geom_bar(width=1) To make the pie chart, we need to add the “coord_polar” function to the code which adjusts the mapping. We will include the argument “theta=y” which tells R that the size of the pie a factor gets depends on the number of people in that factor. Below is the code for the pie chart. ggplot(Wage, aes(factor(1), fill=education))+ geom_bar(width=1)+coord_polar(theta="y") By changing the “width” argument you can place a circle in the middle of the chart as shown below. ggplot(Wage, aes(factor(1), fill=education))+ geom_bar(width=.5)+coord_polar(theta="y") Bullseye Chart A bullseye chart is a pie chart that share the information in a concentric way. The coding is mostly the same except that you remove the “theta” argument from the “coord_polar” function. The thicker the circle the more respondents within it. Below is the code ggplot(Wage, aes(factor(1), fill=education))+ geom_bar(width=1)+coord_polar() Coxcomb Diagram The Coxcomb Diagram is similiar to the pie chart but the data is not normalized to fit the entire area of the circle. To make this plot we have to modify the code to make the by removing the “factor(1)” argument and replacing it with the name of the variable and be readding the “coord_polor” function. Below is the code ggplot(Wage, aes(education, fill=education))+ geom_bar(width=1)+coord_polar() Conclusion These are just some of the many forms of visualizations available using ggplot2. Which to use depends on many factors from personal preference to the needs of the audience. # Histograms and Colors with ggplot2 In this post, we will look at how ggplot2 is able to create variables for the purpose of providing aesthetic information for a histogram. Specifically, we will look at how ggplot2 calculates the bin sizes and then assigns colors to each bin depending on the count or density of that particular bin. To do this we will use dataset called “Star” from the “Edat” package. From the dataset, we will look at total math score and make several different histograms. Below is the initial code you need to begin. library(ggplot2);library(Ecdat) data(Star) We will now create are initial histogram. What is new in the code below is the “..count..” for the “fill” argument. This information tells are to fill the bins based on their count or the number of data points that fall in this bin. By doing this, we get a gradation of colors with darker colors indicating more data points and lighter colors indicating fewer data points. The code is as follows. ggplot(Star, aes(tmathssk, fill=..count..))+geom_histogram() As you can see, we have a nice histogram that uses color to indicate how common data in a specific bin is. We can also make a histogram that has a line that indicates the density of the data using the kernal function. This is similar to adding a LOESS line on a plot. The code is below. ggplot(Star, aes(tmathssk)) + geom_histogram(aes(y=..density.., fill=..density..))+geom_density() The code is mostly the same but we moved the “fill” argument inside “geom_histogram” function and added a second “aes” function. We also included a y argument inside the second “aes” function. Instead of using the “..count..” information we used “..density..” as this is needed to create the line. Lastly, we added the “geom_density” function. The chart below uses the “alpha” argument to add transparency to the histogram. This allows us to communicate additional information. In the histogram below we can see visually information about gender and the how common a particular gender and bin are in the data. ggplot(Star, aes(tmathssk, col=sex, fill=sex, alpha=..count..))+geom_histogram() Conclusion What we have learned in this post is some of the basic features of ggplot2 for creating various histograms. Through the use of colors a researcher is able to display useful information in an interesting way. # Linear Regression Lines and Facets in ggplot2 In this post, we will look at how to add a regression line to a plot using the “ggplot2” package. This is mostly a review of what we learned in the post on adding a LOESS line to a plot. The main difference is that a regression line is a straight line that represents the relationship between the x and y variable while a LOESS line is used mostly to identify trends in the data. One new wrinkle we will add to this discussion is the use of faceting when developing plots. Faceting is the development of multiple plots simultaneously with each sharing different information about the data. The data we will use is the “Housing” dataset from the “Ecdat” package. We will examine how lotsize affects housing price when also considering whether the house has central air conditioning or not. Below is the initial code in order to be prepared for analysis library(ggplot2);library(Ecdat) ## Loading required package: Ecfun ## ## Attaching package: 'Ecdat' ## ## The following object is masked from 'package:datasets': ## ## Orange data("Housing") The first plot we will make is the basic plot of lotsize and price with the data being distinguished by having central air or not, without a regression line. The code is as follows ggplot(data=Housing, aes(x=lotsize, y=price, col=airco))+geom_point() We will now add the regression line to the plot. We will make a new plot with an additional piece of code. ggplot(data=Housing, aes(x=lotsize, y=price, col=airco))+geom_point()+stat_smooth(method='lm') As you can see we get two lines for the two conditions of the data. If we want to see the overall regression line we use the code that is shown below. ggplot()+geom_point(data=Housing, aes(x=lotsize, y=price, col=airco))+stat_smooth(data=Housing, aes(x=lotsize, y=price ),method='lm') We will now experiment with a technique called faceting. Faceting allows you to split the data by various subgroups and display the result via plot simultaneously. For example, below is the code for splitting the data by central air for examining the relationship between lot size and price. ggplot(data=Housing, aes(lotsize, price, col=airco))+geom_point()+stat_smooth(method='lm')+facet_grid(.~airco) By adding the “facet_grid” function we can subset the data by the categorical variable “airco”. In the code below we have three plots. The first two show the relationship between lotsize and price based on central air and the last plot shows the overall relationship. ggplot(data=Housing, aes(lotsize, price, col=airco))+geom_point()+stat_smooth(method="lm")+facet_grid(.~airco, margins = TRUE) By adding the argument “margins” and setting it to true we are able to add the third plot that shows the overall results. So far all of are facetted plots have had the same statistical transformation of the use of a regression. However, we can actually mix the type of transformations that happen when facetting the results. This is shown below. ggplot(data=Housing, aes(lotsize, price, col=airco))+geom_point()+stat_smooth(data=subset(Housing, airco=="yes"))+stat_smooth(data=subset(Housing, airco=="no"), method="lm")+facet_grid(.~airco) In the code we needed to use two functions of “stat_smooth” and indicate the information to transform inside the function. The plot to the left is a regression line with houses without central air and the plot to the right is a LOESS line with houses that have central air. Conclusion In this post we explored the use of regression lines and advance faceting techniques. Communicating data with ggplot2 is one of many ways in which a data analyst can portray valuable information. # Adding LOESS Lines to Plots in R A common goal of statistics is to try and identify trends in the data as well as to predict what may happen. Both of these goals can be partially achieved through the development of graphs and or charts. In this post, we will look at adding a smooth line to a scatterplot using the “ggplot2” package. To accomplish this, we will use the “Carseats” dataset from the “ISLR” package. We will explore the relationship between the price of carseats with actual sales along with whether the carseat was purchase in an urban location or not. Below is some initial code to prepare for the analysis. library(ggplot2);library(ISLR) data("Carseats") We are going to us a layering approach in this example. This means we will add one piece of code at a time until we have the complete plot.We are now going to plot the initial scatterplot. We simply want a scatterplot depicting the relationship between Price and Sales of carseats. ggplot(data=Carseats, aes(x=Price, y=Sales, col=Urban))+geom_point() The general trend appears to be negative. As price increases sales decrease regardless if carseat was purchase in an urban setting or not. We will now add ar LOESS line to the graph. LOESS stands for “Locally weighted smoothing” this is a commonly used tool in regression analysis. The addition of a LOESS line allows in identifying trends visually much easily. Below is the code ggplot(data=Carseats, aes(x=Price, y=Sales, col=Urban))+geom_point()+ stat_smooth() Unlike a regression line which is strictly straight, a LOESS line curves with the data. As you look at the graph the LOESS line is mostly straight with curves at the extremes and for small rise in fall in the middle for carseats purchased in urban areas. So far we have created LOESS lines by the categorical variable Urban. We can actually make a graph with three LOESS lines. One for Yes urban, another for No Urban, and a last one that is an overall line that does not take into account the Urban variable. Below is the code. ggplot()+ geom_point(data=Carseats, aes(x=Price, y=Sales, col=Urban))+ stat_smooth(data=Carseats, aes(x=Price, y=Sales))+stat_smooth(data=Carseats, aes(x=Price, y=Sales, col=Urban)) Notice that the code is slightly different with the information being mostly outside of the “ggplot” function. You can barely see the third line in the graph but if you look closely you will see a new blue line that was not there previously. This is the overall trend line. If you want you can see the overall trend line with the code below. ggplot()+ geom_point(data=Carseats, aes(x=Price, y=Sales, col=Urban))+ stat_smooth(data=Carseats, aes(x=Price, y=Sales)) The very first graph we generated in this post only contained points. This is because we used the “geom_point” function. Any of the graphs we created could be generated with points by removing the “geom_point” function and only using the “stat_smooth” function as shown below. ggplot(data=Carseats, aes(x=Price, y=Sales, col=Urban))+ stat_smooth() Conclusion This post provided an introduction to adding LOESS lines to a graph using ggplot2. For presenting data in a visually appealing way, adding lines can help in identifying key characteristics in the data. # Intro to the Grammar of Graphics In developing graphs, there are certain core principles that need to be addressed in order to provide a graph that communicates meaning clearly to an audience. Many of these core principles are addressed in the book “The Grammar of Graphics” by Leland Wilkinson. The concepts of Wilkinson’s book were used to create the “ggplot2” r package by Hadley Wickham. This post will explain some of the core principles needed in developing high quality visualizations. In particular we will look at the following. • Aesthetic attributes • Geometric objects • Statistical transformations • Scales • Coordinates • Faceting One important point to mention is that when using ggplot2 not all of these concepts have to be addressed in the code as R will auto-select certain features if you do not specify them. Aesthetic Attributes and Geometric Objects Aesthetic attributes is about how the data is perceived. This general involves arguments in the “ggplot” relating to the x/y coordinates as well as the actual data that is being used. Aesthetic atrributes is mandatory information for making a graph. Geometric objects determines what type of plot is generated. There are many different examples such as bar, point, boxplot, and histogram. To use the “ggplot” function you must provide the aesthetic and geometric object informatio to generate a plot. Below is coding containing only this information. library(ggplot2) ggplot(Orange, aes(circumference))+geom_histogram() ## stat_bin() using bins = 30. Pick better value with binwidth. The code is broken down as follows ggplot(data, aesthetic attribute(x-axis data at least)+geometric object()) Statistical Transformation Statistical transformation involves combining the data in one way or the other to get a general sense of the data. Examples of statistical transformation includes adding a smooth line, a regression line, or even binning the data for histograms. This feature is optional but can provide additional explanation of the data. Below we are going to look at two variables on one plot. For this we will need a different geomtric object as we will use points instead of a histogram. We will also use a statisitcal transformation. In particular, the statistical transformation is regression line. The code is as follows ggplot(Orange, aes(circumference, age))+geom_point()+stat_smooth(method="lm") The code is broken down as follows ggplot(data, aesthetic attribute(x-axis data at least)+geometric object()+ statistical transformation(type of transformation)) Scales Scales is a rather complicated feature. For simplicity, scales have to do with labelling the title, x and y-axis, creating a legend, as well as the coloring of data points. This use of this feature is optional. Below is a simple example using the “labs” function in the plot we develop in the previous example. ggplot(Orange, aes(circumference,age))+geom_point()+stat_smooth(method="lm") + labs(title="Example Plot", x="circumference of the tree", y="age of the tree") The plot now has a title and clearly labelled x and y axises Coordinates Coordinates is another complex feature. This feature allows for the adjusting of the mapping of the data. Two common mappin features are cartesian and polar. Cartesian is commonly used for plots in 2D while polor is often used for pie charts. In the example below, we will use the same data but this time use a polor mapping approach. The plot doesn’t make much sense but is really just an example of using this feature. This feature is also optional. ggplot(Orange, aes(circumference, age))+geom_point()+stat_smooth(method="lm")+labs(title="Example Plot",x="circumference of the tree", y="age of the tree")+coord_polar() The last feature is faceting. Faceting allows you to group data in subsets. This allows you to look at your data from the perspective of various subgroups in the sample population. In the example below, we will look at the relationship between circumference and age by tree type. ggplot(Orange, aes(circumference, age))+geom_point()+stat_smooth(method="lm")+labs(title="Example Plot",x="circumference of the tree", y="age of the tree")+facet_grid(Tree~.) Now we can see the relationship between the two variables based on the type of tree. One important thing to note about the “facet_grid” function is the use of the “.~” If this symbol “~.” is placed behind the categorical variable the charts will be stacked on top of each other is in the previous example. However, if the symbol is written differently “.~” and placed in front of the categorical variable the plots will be placed next to each other as in the example below ggplot(Orange, aes(circumference, age))+geom_point()+stat_smooth(method="lm")+labs(title="Example Plot",x="circumference of the tree", y="age of the tree")+facet_grid(.~Tree) Conclusion This post provided an introduction to the grammar of graphics. In order to appreciate the art of data visualization it requires understanding how the different principles of graphics work together to communicate information in a visually manner with an audience. # Using Qplots for Graphs in R In this post, we will explore the use of the “qplot” function from the “ggplot2” package. One of the major advantages of “ggplot” when compared to the base graphics package in R is that the “ggplot2” plots are much more visually appealing. This will make more sense when we explore the grammar of graphics. for now we will just make plots to get use to using the “qplot” function. We are going to use the “Carseats” dataset from the “ISLR” package in the examples. This dataset has data about the purchase of carseats for babies. Below is the initial code you will need to make the various plots. library(ggplot2);library(ISLR) data("Carseats") In the first scatterplot, we are going to compare the price of a carseat with the volumn of sales. Below is the code qplot(Price, Sales,data=Carseats) Most of this coding format you are familiar. “Price” is the x variable. “Sales” is the y variable and the data used is “Carseats. From the plot, we can see that as the price of the carseat increases there is normally a decline in the number of sales. For our next plot, we will compare sales based on shelf location. This requires the use of a boxplot. Below is the code qplot(ShelveLoc, Sales, data=Carseats, geom="boxplot") The new argument in the code is the “geom” argument. This argument indicates what type of plot is drawn. The boxplot appears to indicate that a “good” shelf location has the best sales. However, this would need to be confirmed with a statistical test. Perhaps you are wondering how many of the Carseats where in the bad, medium, and good shelf locations. To find out, we will make a barplot as shown in the code below qplot(ShelveLoc, data=Carseats, geom="bar") The most common location was medium with bad and good be almost equal. Lastly, we will now create a histogram using the “qplot” function. We want to see the distribution of “Sales”. Below is the code qplot(Sales, data=Carseats, geom="histogram") ## stat_bin() using bins = 30. Pick better value with binwidth. The distribution appears to be normal but again to know for certain requires a statistical test. For one last, trick we will add the median to the plot by using the following code qplot(Sales, data=Carseats, geom="histogram") + geom_vline(xintercept = median(Carseats$Sales), colour="blue")
## stat_bin() using bins = 30. Pick better value with binwidth.

To add the median all we needed to do was add an additional argument called “geom_vline” which adds a line to a plot. Inside this argument we had to indicate what to add by indicating the median of “Sales” from the “Carseats” package.

Conclusion

This post provided an introduction to the use of the “qplot” function in the “ggplot2” package. Understanding the basics of “qplor” is beneficial in providing visually appealing graphics

# Making Graphics in R

Data visualization is a critical component of communicate results with an audience. Fortunately, R provides many different ways to present numerical data it a clear way visually. This post will look specifically at making data visualizations with the base r package “graphics”.

Generally, functions available in the “graphics” package can be either high-level functions or low-level functions. High-level functions actually make the plot. Examples of high-level functions includes are the “hist” (histogram), “boxplot” (boxplot), and “barplot” (bar plot).

Low-level functions are used to add additional information to a plot. Some commonly used low-level functions includes “legend” (add legend) and “text” (add text). When coding we allows call high-level functions before low-level functions as the other way is not accepted by R.

We are going to begin with a simple graph. We are going to use the“Caschool” dataset from the “Ecdat” package. For now, we are going to plot the average expenditures per student by the average number of computers per student. Keep in mind that we are only plotting the data so we are only using a high-level function (plot). Below is the code

library(Ecdat)
data("Caschool")
plot(compstu~expnstu, data=Caschool)

The plot is not “beautiful” but it is a start in plotting data. Next, we are going to add a low-level function to our code. In particular, we will add a regression line to try and see the diretion of the relationship between the two variables via a straight line. In addition, we will use the “loess.smooth” function. This function will allow us to see the general shape of the data. The regression line is green and the loess smooth line is blue. The coding is mostly familiy but the “lwd” argument allows us to make the line thicker.

plot(compstu~expnstu, data=Caschool)
abline(lm(compstu~expnstu, data=Caschool), col="green", lwd=5)
lines(loess.smooth(Caschool$expnstu, Caschool$compstu), col="blue", lwd=5)

Boxplots allow you to show data that has been subsetted in some way. This allows for the comparisions of groups. In addition, one or more boxplots can be used to identify outliers.

In the plot below, the student-to-teacher ratio of k-6 and k-8 grades are displayed.

boxplot(str~grspan, data=Caschool)

As you look at the data you can see there is very little difference. However, one major differnce is that the K-8 group has much more extreme values than K-6.

Histograms are an excellent way to display information about one continuous variable. In the plot below, we can see the spread of the expenditure per student.

hist(Caschool$expnstu) We will now add median to the plot by calling the low-level function “abline”. Below is the code. hist(Caschool$expnstu)
abline(v=median(Caschool$expnstu), col="green", lwd=5) Conclusion In this post, we learned some of the basic structures of creating plots using the “graphics” package. All plots in include both low and high-level functions that work together to draw and provide additional information for communicating data in a visual manner # Bagging in R In this post, we will explore the potential of bagging. Bagging is a process in which the original data is boostrapped to make several different datasets. Each of these datasets are used to generate a model and voting is used to classify an example or averaging is used for numeric prediction. Bagging is especially useful for unstable learners. These are algorithms who generate models that can change a great deal when the data is modified a small amount. In order to complete this example, you will need to load the following packages, set the seed, as well as load the dataset “Wages1”. We will be using a decision tree that was developed in an earlier post. Below is the initial code library(caret); library(Ecdat);library(ipred);library(vcd) set.seed(1) data(Wages1) We will now use the “bagging” function from the “ipred” package to create our model as well as tell R how many bags to make. theBag<-bagging(sex~.,data=Wages1,nbagg=25) Next, we will make our predictions. Then we will check the accuracy of the model looking at a confusion matrix and the kappa statistic. The “kappa” function comes from the “vcd” package. bagPred<-predict(theBag, Wages1) keep<-table(bagPred, Wages1$sex)
keep
##
## bagPred  female male
##   female   1518   52
##   male       51 1673
Kappa(keep)
##             value      ASE     z Pr(>|z|)
## Unweighted 0.9373 0.006078 154.2        0
## Weighted   0.9373 0.006078 154.2        0

The results appearing exciting with almost 97% accuracy. In addition, the Kappa was almost 0.94 indicating a well-fitted model. However, in order to further confirm this we can cross-validate the model instead of using bootstrap aggregating as bagging does. Therefore we will do a 10-fold cross-validation using the functions from the “caret” package. Below is the code.

ctrl<-trainControl(method="cv", number=10)
trainModel<-train(sex~.,data=Wages1, method="treebag",trControl=ctrl)
trainModel
## Bagged CART
##
## 3294 samples
##    3 predictors
##    2 classes: 'female', 'male'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 2965, 2965, 2965, 2964, 2964, 2964, ...
## Resampling results
##
##   Accuracy   Kappa       Accuracy SD  Kappa SD
##   0.5504128  0.09712194  0.02580514   0.05233441
##
## 

Now the results are not so impressive. In many ways the model is terrible. The accuracy has fallen significantly and the kappa is almost 0. Remeber that cross-validation is an indicator of future performance. This means that our current model would not generalize well to other datasets.

Bagging is not limited to decision trees and can be used for all machine learning models. The example used in this post was one that required the least time to run. For real datasets, the processing time can be quite long for the average laptop.

# Developing a Customize Tuning Process in R

In this post, we will learn how to develop customize criteria for tuning a machine learning model using the “caret” package. There are two things that need to be done in order to complete assess a model using customized features. These two steps are…

• Determine the model evaluation criteria
• Create a grid of parameters to optimize

The model we are going to tune is the decision tree model made in a previous post with the C5.0 algorithm. Below is code for loading some prior information.

library(caret); library(Ecdat)
data(Wages1)

DETERMINE the MODEL EVALUATION CRITERIA

We are going to begin by using the “trainControl” function to indicate to R what re-sampling method we want to use, the number of folds in the sample, and the method for determining the best model. Remember, that there are many more options but these are the onese we will use. All this information must be saved into a variable using the “trainControl” function. Later, the information we place into the variable will be used when we rerun our model.

For our example, we are going to code the following information into a variable we will call “chck” for re sampling we will use k-fold cross-validation. The number of folds will be set to 10. The criteria for selecting the best model will be the through the use of the “oneSE” method. The “oneSE” method selects the simplest model within one standard error of the best performance. Below is the code for our variable “chck”

chck<-trainControl(method = "cv",number = 10, selectionFunction = "oneSE")

For now this information is stored to be used later

CREATE GRID OF PARAMETERS TO OPTIMIZE

We now need to create a grid of parameters. The grid is essential the characteristics of each model. For the C5.0 model we need to optimize the model, number of trials, and if winnowing was used. Therefore we will do the following.

• For model, we want decision trees only
• Trials will go from 1-35 by increments of 5
• For winnowing, we do not want any winnowing to take place.

In all we are developing 8 models. We know this based on the trial parameter which is set to 1, 5, 10, 15, 20, 25, 30, 35. To make the grid we use the “expand.grid” function. Below is the code.

modelGrid<-expand.grid(.model ="tree", .trials= c(1,5,10,15,20,25,30,35), .winnow="FALSE")

CREATE THE MODEL

We are now ready to generate our model. We will use the kappa statistic to evaluate each model’s performance

set.seed(1)
customModel<- train(sex ~., data=Wages1, method="C5.0", metric="Kappa", trControl=chck, tuneGrid=modelGrid)
customModel
## C5.0
##
## 3294 samples
##    3 predictors
##    2 classes: 'female', 'male'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 2966, 2965, 2964, 2964, 2965, 2964, ...
## Resampling results across tuning parameters:
##
##   trials  Accuracy   Kappa      Accuracy SD  Kappa SD
##    1      0.5922991  0.1792161  0.03328514   0.06411924
##    5      0.6147547  0.2255819  0.03394219   0.06703475
##   10      0.6077693  0.2129932  0.03113617   0.06103682
##   15      0.6077693  0.2129932  0.03113617   0.06103682
##   20      0.6077693  0.2129932  0.03113617   0.06103682
##   25      0.6077693  0.2129932  0.03113617   0.06103682
##   30      0.6077693  0.2129932  0.03113617   0.06103682
##   35      0.6077693  0.2129932  0.03113617   0.06103682
##
## Tuning parameter 'model' was held constant at a value of tree
##
## Tuning parameter 'winnow' was held constant at a value of FALSE
## Kappa was used to select the optimal model using  the one SE rule.
## The final values used for the model were trials = 5, model = tree
##  and winnow = FALSE.

The actually output is similar to the model that “caret” can automatically create. The difference here is that the criteria was set by us rather than automatically. A close look reveals that all of the models perform poorly but that there is no change in performance after ten trials.

CONCLUSION

This post provided a brief explanation of developing a customize way of assessing a models performance. To complete this, you need configure your options as well as setup your grid in order to assess a model. Understanding the customization process for evaluating machine learning models is one of the strongest ways to develop supremely accurate models that retain generalizability.

# Developing an Automatically Tuned Model in R

In this post, we are going to learn how to use the “caret” package to automatically tune a machine learning model. This is perhaps the simplest way to evaluate the performance of several models. In a later post, we will explore how to perform custom tuning to a model.

The model we are trying to tune is the decision tree we made using the C5.0 algorithm in a previous post. Specifically we were trying to predict sex based on the variables available in the “Wages1” dataset in the “Ecdat” package.

In order to accomplish our goal we will need to load the “caret” and “Ecdat”package, load the “Wages1” dataset as well as set the seed. Setting the seed will allow us to reproduce our results. Below is the code for these steps.

library(caret); library(Ecdat)
data(Wages1)
set.seed(1)

We will now build and display our model using the code below.

tuned_model<-train(sex ~., data=Wages1, method="C5.0")

tuned_model
## C5.0
##
## 3294 samples
##    3 predictors
##    2 classes: 'female', 'male'
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 3294, 3294, 3294, 3294, 3294, 3294, ...
## Resampling results across tuning parameters:
##
##   model  winnow  trials  Accuracy   Kappa      Accuracy SD  Kappa SD
##   rules  FALSE    1      0.5892713  0.1740587  0.01262945   0.02526656
##   rules  FALSE   10      0.5938071  0.1861964  0.01510209   0.03000961
##   rules  FALSE   20      0.5938071  0.1861964  0.01510209   0.03000961
##   rules   TRUE    1      0.5892713  0.1740587  0.01262945   0.02526656
##   rules   TRUE   10      0.5938071  0.1861964  0.01510209   0.03000961
##   rules   TRUE   20      0.5938071  0.1861964  0.01510209   0.03000961
##   tree   FALSE    1      0.5841768  0.1646881  0.01255853   0.02634012
##   tree   FALSE   10      0.5930511  0.1855230  0.01637060   0.03177075
##   tree   FALSE   20      0.5930511  0.1855230  0.01637060   0.03177075
##   tree    TRUE    1      0.5841768  0.1646881  0.01255853   0.02634012
##   tree    TRUE   10      0.5930511  0.1855230  0.01637060   0.03177075
##   tree    TRUE   20      0.5930511  0.1855230  0.01637060   0.03177075
##
## Accuracy was used to select the optimal model using  the largest value.
## The final values used for the model were trials = 10, model = rules
##  and winnow = TRUE.

There is a lot of information that is printed out. The first column is the type of model developed. Two types of models were developed either a rules-based classification tree or a normal decision tree. Next, is the winnow column. This column indicates if a winnowing process was used to remove poor predictor variables.

The next two columns are accuracy and kappa which have been explained previously. The last two columns are the standard deviations of accuarcy and kappa. None of the models are that good but the purpose here is for teaching.

At the bottom of the printout, r tells you which model was the best. For us, the best model was the fifth model from the top which was a rule-based, 10 trial model with winnow set to “TRUE”.

We will now use the best model (the caret package automatically picks it) to make predictions on the training data. We will also look at the confusion matrix of the correct classification followed by there proportions. Below is the code.

predict_model<-predict(tuned_model, Wages1)
table(predict_model, Wages1$sex) ## ## predict_model female male ## female 936 590 ## male 633 1135 prop.table(table(predict_model, Wages1$sex))
##
## predict_model    female      male
##        female 0.2841530 0.1791135
##        male   0.1921676 0.3445659

In term of prediction, the model was correct 62% of the time (.28 + .34 = .62). If we want to know, can also see the probabilities for each example using the following code.

probTable<-(predict(tuned_model, Wages1, type="prob"))
head(probTable)
##      female       male
## 1 0.6191287 0.38087132
## 2 0.2776770 0.72232303
## 3 0.2975327 0.70246734
## 4 0.7195866 0.28041344
## 5 1.0000000 0.00000000
## 6 0.9092993 0.09070072


Conclusion

In this post, we looked at an automated way to determine the best model among many using the “caret” package. Understanding how to improve the performance of a model is critical skill in machine learning.

# Improving the Performance of Machine Learning Model

For many, especially beginners, making a machine learning model is difficult enough. Trying to understand what to do, how to specify the model, among other things is confusing in itself. However, after developing a model it is necessary to assess ways in which to improve performance.

This post will serve as an introduction to understanding how to improving model performance. In particular, we will look at the following

• When it is necessary to improve performance
• Parameter tuning

When to Improve

It is not always necessary to try and improve the performance of a model. There are times when a model does well and you know this through the evaluating it. If the commonly used measures are adequate there is no cause for concern.

However, there are times when improvement is necessary. Complex problems, noisy data, and trying to look for subtle/unclear relationships can make improvement necessary. Normally, real-world data has the problems so model improvement is usually necessary.

Model improvement requires the application of scientific means in an artistic manner. It requires a sense of intuition at times and also brute trial-and-error effort as well. The point is that there is no singular agreed upon way to improve a model. It is better to focus on explaining how you did it if necessary.

Parameter Tuning

Parameter tuning is the actual adjustment of model fit options. Different machine learning models have different options that can be adjusted. Often, this process can be automated in r through the use of the “caret” package.

When trying to decide what to do when tuning parameters it is important to remember the following.

• What machine learning model and algorithm you are using for your data.
• Which parameters you can adjust.
• What criteria you are using to evaluate the model

Naturally, you need to know what kind of model and algorithm you are using in order to improve the model. There are three types of models in machine learning, those that classify, those that employ regression, and those that can do both. Understanding this helps you to make decision about what you are trying to do.

Next, you need to understand what exactly you or r are adjusting when analyzing the model. For example, for C5.0 decision trees “trials” is one parameter you can adjust. If you don’t know this, you will not know how the model was improved.

Lastly, it is important to know what criteria you are using to compare the various models. For classifying models you can look at the kappa and the various information derived from the confusion matrix. For regression based models you may look at the r-square, the RMSE (Root mean squared error), or the ROC curve.

Conclusion

As you can perhaps tell there is an incredible amount of choice and options in trying to improve a model. As such, model improvement requires a clearly developed strategy that allows for clear decision-making.

In a future post, we will look at an example of model improvement.

# K-Fold Cross-Validation

In this post, we are going to look at k-fold cross-validation and its use in evaluating models in machine learning.

K-fold cross-validation is use for determining the performance of statistical models. How it works is the data is divided into a pre-determined number of folds (called ‘k’). One fold is used to determine the model estimates and the other folds are used for evaluating. This is done k times and the results are average based on a statistic such as kappa to see how the model performs.

In our example we are going to review a model we made using the C5.0 algorithm. In that post, we were trying to predict gender based on several other features.

First, we need to load several packages into R as well as the dataset we are going to use. All of this is shared in the code below

library(caret);library(C50);library(irr);library(Ecdat)
data("Wages1")

We now need to set the seed. This is important for allowing us to reproduce the results. Every time a k-fold is perform the results can be slightly different but setting the seed prevents this. The code is as follows

set.seed(1)

We will now create are folds. How many folds to create is up to the researcher. For us, we are going to create ten folds. What this means is that R will divide are sample into ten equal parts. To do this we use the “createFolds” function from the “caret” package. After creating the folds, we will view the results using the “str” function which will tell us how many examples are in each fold. Below is the code to complete this.

folds<-createFolds(Wages1$sex, k=10) str(folds) ## List of 10 ##$ Fold01: int [1:328] 8 13 18 37 39 57 61 67 78 90 ...
##  $Fold02: int [1:329] 5 27 47 48 62 72 76 79 85 93 ... ##$ Fold03: int [1:330] 2 10 11 31 32 34 36 64 65 77 ...
##  $Fold04: int [1:330] 19 24 40 45 55 58 81 96 99 102 ... ##$ Fold05: int [1:329] 6 14 28 30 33 84 88 91 95 97 ...
##  $Fold06: int [1:330] 4 15 16 38 43 52 53 54 56 63 ... ##$ Fold07: int [1:330] 1 3 12 22 29 50 66 73 75 82 ...
##  $Fold08: int [1:330] 7 21 23 25 26 46 51 59 60 83 ... ##$ Fold09: int [1:329] 9 20 35 44 49 68 74 94 100 105 ...
##  $Fold10: int [1:329] 17 41 42 71 101 107 117 165 183 184 ... As you can see, we normally have about 330 examples per fold. In order to get the results that we need. We have to take fold 1 to make the model and fold 2-10 to evaluate it. We repeat this process until every combination possible is used. First fold 1 is used and 2-10 are the test data, then fold 2 is used and then folds 1, 3-10 are the test data etc. Manually coding this would take great deal of time. To get around this we will use the “lapply” function. Using “lapply” we will create a function that takes “x” (one of our folds) and makes it the “training set” shown here as “Wages1_train”. Next we assigned the rest of the folds to be the “test” (Wages1_test) set as depicted with the “-x”. The next two lines of code should look familiar as it is the code for developing a decision tree. The “Wages_actual” are the actual labels for gender in the “Wages_1” testing set. The “kappa2” function is new and it comes form the “irr” package. The kappa statistic is a measurement of accuracy of a model while taking into account chance. The closer the value is to 1 the better. Below is the code for what has been discussed. results<-lapply(folds, function(x) { Wages1_train<-Wages1[x, ] Wages1_test<-Wages1[-x, ] Wages1_model<-C5.0(sex~.,data=Wages1) Wages1_pred<-predict(Wages1_model, Wages1_test) Wages1_actual<-Wages1_test$sex
Wages1_kappa<-kappa2(data.frame(Wages1_actual, Wages1_pred))$value return(Wages1_kappa) }) To get our results, we will use the “str” function again to display them. This will tell us the kappa for each fold. To really see how are model does we need to calculate the mean kappa of the ten models. This is done with the “unlist” and “mean” function as shown below str(results) ## List of 10 ##$ Fold01: num 0.205
##  $Fold02: num 0.186 ##$ Fold03: num 0.19
##  $Fold04: num 0.193 ##$ Fold05: num 0.202
##  $Fold06: num 0.208 ##$ Fold07: num 0.196
##  $Fold08: num 0.202 ##$ Fold09: num 0.194
##  $Fold10: num 0.204 mean(unlist(results)) ## [1] 0.1978915 The final mean kappa was 0.19 which is really poor. It indicates that the model is no better at predicting then chance alone. However, for illustrative purposes we now understand how to perfrom a k-fold cross-validation # Receiver Operating Characteristic Curve The receiver operating characteristic curve (ROC curve) is a tool used in statistical research to assess the trade-off of detecting true positives and true negatives. The origins of this tool goes all the way back to WWII when engineers were trying to distinguish between true and false alarms. Now this technique is used in machine learning This post will explain the ROC curve and provide and example using R. Below is a diagram of an ROC curve On the X axis we have the false positive rate. As you move to the right the false positive rate increases which is bad. We want to be as close to zero as possible. On the y axis we have the true positive rate. Unlike the x axis we want the true positive rate to be as close to 100 as possible. In general we want a low value on the x-axis and a high value on the y-axis. In the diagram above, the diagonal line called “Test without diagnostic benefit” represents a model that cannot tell the difference between true and false positives. Therefore, it is not useful for our purpose. The L-shaped curve call “Good diagnostic test” is an example of an excellent model. This is because all the true positives are detected . Lastly, the curved-line called “Medium diagonistic test” represents an actually model. This model is a balance between the perfect L-shaped model and the useless straight-line model. The curved-line model is able to moderately distinguish between false and true positives. Area Under the ROC Curve The area under an ROC curve is literally called the “Area Under the Curve” (AUC). This area is calculated with a standardized value ranging from 0 – 1. The closer to 1 the better the model We will now look at an analysis of a model using the ROC curve and AUC. This is based on the results of a post using the KNN algorithm for nearest neighbor classification. Below is the code predCollege <- ifelse(College_test_pred=="Yes", 1, 0) realCollege <- ifelse(College_test_labels=="Yes", 1, 0) pr <- prediction(predCollege, realCollege) collegeResults <- performance(pr, "tpr", "fpr") plot(collegeResults, main="ROC Curve for KNN Model", col="dark green", lwd=5) abline(a=0,b=1, lwd=1, lty=2) aucOfModel<-performance(pr, measure="auc") unlist(aucOfModel@y.values) 1. The first to variables (predCollege & realCollege) is just for converting the values of the prediction of the model and the actual results to numeric variables 2. The “pr” variable is for storing the actual values to be used for the ROC curve. The “prediction” function comes from the “ROCR” package 3. With the information information of the “pr” variable we can now analyze the true and false positives, which are stored in the “collegeResults” variable. The “performance” function also comes from the “ROCR” package. 4. The next two lines of code are for plot the ROC curve. You can see the results below 6. The curve looks pretty good. To confirm this we use the last two lines of code to calculate the actually AUC. The actual AUC is 0.88 which is excellent. In other words, the model developed does an excellent job of discerning between true and false positives. Conclusion The ROC curve provides one of many ways in which to assess the appropriateness of a model. As such, it is yet another tool available for a person who is trying to test models. # Using Probability of the Prediction to Evaluate a Machine Learning Model In evaluating a model when employing machine learning techniques, there are three common types of data used for evaluation. • The actual classification values • The predicted classification values • The estimated probability of the prediction The first two types of data (actual and predicted) are used for assessing the accuracy of a model in several different ways such as error rate, sensitivity, specificity, etc. The benefit of the probabilities of prediction is that it is a measure of a model’s confidence in its prediction. If you need to compare to models and one is more confident in it’s prediction of its classification of examples, the more confident model is the better learner. In this post, we will look at examples of the probability predictions of several models that have been used in this blog in the past. Prediction Probabilities for Decision Trees Our first example come from the decision tree we made using the C5.0 algorithm. Below is the code for calculating the probability of the correct classification of each example in the model followed by an output of the first Wage_pred_prob<-predict(Wage_model, Wage_test, type="prob")   head(Wage_pred_prob) female male 497 0.2853016 0.7146984 1323 0.2410568 0.7589432 1008 0.5770177 0.4229823 947 0.6834378 0.3165622 695 0.5871323 0.4128677 1368 0.4303364 0.5696636 The argument “type” is added to the “predict” function so that R calculates the probability that the example is classified correctly. A close look at the results using the “head” function provides a list of 6 examples from the model. • For example 497, there is a 28.5% probability that this example is female and a 71.5% probability that this example is male. Therefore, the model predicts that this example is male. • For example 1322, there is a 24% probability that this example is female and a 76% probability that this example is male. Therefore, the model predicts that this example is male. • etc. Prediction Probabilities for KNN Nearest Neighbor Below is the code for finding the probilities for KNN algorithm. College_test_pred_prob<-knn(train=College_train, test=College_test, + cl=College_train_labels, k=27, prob=TRUE)  College_test_pred_prob The print for this is rather long. However, you can match the predict level with the actual probability by looking carefully at the data. • For example 1, there is a 77% probability that this example is a yes and a 23% probability that this example is a no. Therefore, the model predicts that this example as yes. • For example 2, there is a 71% probability that this example is no and a 29% probability that this example is yes. Therefore, the model predicts that this example is a no. Conclusion One of the primary purposes of the probabilities option is in comparing various models that are derived from the same data. This information combined with other techniques for evaluating models can help a researcher in determining the most appropriate model of analysis. # Kmeans Analysis in R In this post, we will conduct a kmeans analysis on some data on student alcohol consumption. The data is available at the UC Irvine Machine Learning Repository and is available at https://archive.ics.uci.edu/ml/datasets/STUDENT+ALCOHOL+CONSUMPTION We want to see what segments are within the sample of students who participated in this study on several factors in addition to alcohol consumption. Understanding the characteristics that groups of students have in common could be useful in reaching out to them for various purposes. We will begin by loading the “stats” package for the kmeans analysis. Then we will combine the data at it is in two different files and we will explore the data using tables and histograms. I will not print the results of the exploration of the data here as there are too many variables. Lastly, we need to set the seed in order to get the same results each time. The code is still found below. library(stats) student.mat <- read.csv("~/Documents/R working directory/student-mat.csv", sep=";") student.por <- read.csv("~/Documents/R working directory/student-por.csv", sep=";") student_alcohol <- rbind(student.mat, student.por) set.seed(123) options(digits = 2) • str(student_alcohol) • hist(student_alcoholage) • table(studentalcoholage) • table(studentalcoholaddress) • table(student_alcoholfamsize) • table(studentalcoholfamsize) • table(studentalcoholPstatus) • hist(student_alcoholMedu) • hist(studentalcoholMedu) • hist(studentalcoholFedu) • hist(student_alcoholtraveltime) • hist(studentalcoholtraveltime) • hist(studentalcoholstudytime) • hist(student_alcoholfailures) • table(studentalcoholfailures) • table(studentalcoholschoolsup) • table(student_alcoholfamsup) • table(studentalcoholfamsup) • table(studentalcoholpaid) • table(student_alcoholactivities) • table(studentalcoholactivities) • table(studentalcoholnursery) • table(student_alcoholhigher) • table(studentalcoholhigher) • table(studentalcoholinternet) • hist(student_alcoholfamrel) • hist(studentalcoholfamrel) • hist(studentalcoholfreetime) • hist(student_alcoholgoout) • hist(studentalcoholgoout) • hist(studentalcoholDalc) • hist(student_alcoholWalc) • hist(studentalcoholWalc) • hist(studentalcoholhealth) • hist(student_alcohol$absences)

The details about the variables can be found at the website link in the first paragraph of this post. The study look at students alcohol use and other factors related to school and family life.

Before we do the actual kmeans clustering we need to normalize the variables. This is because are variables are measured using different scales. Some are Likert with 5 steps, while others are numeric going from 0 to over 300. The different ranges have an influence on the results. To deal with this problem we will use the “scale” for the variables that will be included in the analysis. Below is the code

student_alcohol_clean<-as.data.frame(student_alcohol)
student_alcohol_clean$age<-scale(student_alcohol$age)
student_alcohol_clean$address<-scale(as.numeric(student_alcohol$address))
student_alcohol_clean$famsize<-scale(as.numeric(student_alcohol$famsize))
student_alcohol_clean$Pstatus<-scale(as.numeric(student_alcohol$Pstatus))
student_alcohol_clean$Medu<-scale(student_alcohol$Medu)
student_alcohol_clean$Fedu<-scale(student_alcohol$Fedu)
student_alcohol_clean$traveltime<-scale(student_alcohol$traveltime)
student_alcohol_clean$studytime<-scale(student_alcohol$studytime)
student_alcohol_clean$failures<-scale(student_alcohol$failures)
student_alcohol_clean$schoolsup<-scale(as.numeric(student_alcohol$schoolsup))
student_alcohol_clean$famsup<-scale(as.numeric(student_alcohol$famsup))
student_alcohol_clean$paid<-scale(as.numeric(student_alcohol$paid))
student_alcohol_clean$activities<-scale(as.numeric(student_alcohol$activities))
student_alcohol_clean$internet<-scale(as.numeric(student_alcohol$internet))
student_alcohol_clean$famrel<-scale(student_alcohol$famrel)
student_alcohol_clean$freetime<-scale(student_alcohol$freetime)
student_alcohol_clean$goout<-scale(student_alcohol$goout)
student_alcohol_clean$Dalc<-scale(student_alcohol$Dalc)
student_alcohol_clean$Walc<-scale(student_alcohol$Walc)
student_alcohol_clean$health<-scale(student_alcohol$health)
student_alcohol_clean$absences<-scale(student_alcohol$absences)
student_alcohol_clean$G1<-scale(student_alcohol$G1)
student_alcohol_clean$G2<-scale(student_alcohol$G2)
student_alcohol_clean$G3<-scale(student_alcohol$G3)

We also need to create a matrix in order to deal with the factor variables. All factor variables need to be converted so that they have dummy variables for the analysis. To do this we use the “matrix” function as shown in the code below.

student_alcohol_clean_matrix<-(model.matrix(~.+0, data=student_alcohol_clean))

We are now ready to conduct our kmeans cluster analysis using the “kmeans” function. We have to determine how many clusters to develop before the analysis. There are statistical ways to do this but another method is domain knowledge. Since we are dealing with teenagers, it is probably that there will be about four distinct groups because of how high school is structured. Therefore, we will use four segments for our analysis. Our code is below.

alcohol_cluster<-kmeans(student_alcohol_clean_matrix, 4)

To view the results we need to view two variables in are “alcohol_cluster” list. The “size” variable will tell us how many people are in each cluster and the “centers” variable describes a clusters characteristics on a particular variable. Below is the code

alcohol_cluster$size # size of clusters ## [1] 191 381 325 147 alcohol_cluster$centers #center of clusters
##   schoolGP schoolMS sexM   age address famsize Pstatus  Medu  Fedu
## 1     0.74     0.26 0.70  0.06  0.0017   0.265  -0.030  0.25  0.29
## 2     0.69     0.31 0.27 -0.15 -0.1059  -0.056  -0.031 -0.53 -0.45
## 3     0.88     0.12 0.45 -0.13  0.3363   0.005   0.016  0.73  0.58
## 4     0.56     0.44 0.48  0.59 -0.4712  -0.210   0.086 -0.55 -0.51
##   Mjobhealth Mjobother Mjobservices Mjobteacher Fjobhealth Fjobother
## 1      0.079      0.30         0.27       0.199     0.0471      0.54
## 2      0.031      0.50         0.17       0.024     0.0210      0.62
## 3      0.154      0.26         0.28       0.237     0.0708      0.48
## 4      0.034      0.46         0.22       0.041     0.0068      0.60
##   Fjobservices Fjobteacher reasonhome reasonother reasonreputation
## 1         0.33       0.042       0.27       0.120             0.20
## 2         0.25       0.031       0.27       0.113             0.20
## 3         0.28       0.123       0.24       0.077             0.34
## 4         0.29       0.034       0.17       0.116             0.15
##   guardianmother guardianother traveltime studytime failures schoolsup
## 1           0.69         0.079       0.17     -0.32    -0.12    -0.079
## 2           0.70         0.052       0.10      0.10    -0.26     0.269
## 3           0.73         0.040      -0.37      0.29    -0.35    -0.213
## 4           0.65         0.170       0.33     -0.49     1.60    -0.123
##   famsup   paid activities nurseryyes higheryes internet romanticyes
## 1 -0.033  0.253     0.1319       0.79      0.92    0.228        0.34
## 2 -0.095 -0.098    -0.2587       0.76      0.93   -0.387        0.35
## 3  0.156  0.079     0.2237       0.88      1.00    0.360        0.31
## 4 -0.057 -0.250     0.0047       0.73      0.70   -0.091        0.49
##   famrel freetime goout  Dalc  Walc health absences    G1    G2     G3
## 1 -0.184     0.43  0.76  1.34  1.29  0.273    0.429 -0.23 -0.17 -0.129
## 2 -0.038    -0.31 -0.35 -0.37 -0.40 -0.123   -0.042 -0.17 -0.12 -0.053
## 3  0.178    -0.01 -0.14 -0.41 -0.35 -0.055   -0.184  0.90  0.87  0.825
## 4 -0.055     0.25  0.22  0.11  0.14  0.087   -0.043 -1.24 -1.40 -1.518

The size of the each cluster is about the same which indicates reasonable segmentation of the sample. The output for “centers” tells us how much above or below the mean a particular cluster is. For example, for the variable “age” we see the following

age
1    0.06
2   -0.15
3   -0.13
4    0.59

What this means is that people in cluster one have an average age 0.06 standard deviations above the mean, cluster two is -0.14 standard deviations below the mean etc. To give our clusters meaning we have to look at the variables and see which one the clusters are extremely above or below the mean. Below is my interpretation of the clusters. The words in parenthesis is the variable from which I made my conclusion

Cluster 1 doesn’t study much (studytime), lives in the biggest families (famsize), requires litle school support (schoolsup), has a lot of free time (freetime), and consumes the most alcohol (Dalc, Walc), lives in an urban area (address), loves to go out the most (goout), and has the most absences. This is the underachieving party alcoholics of the sample

Cluster 2 have parents that are much less educated (Medu, Fedu), requires the most school support (schoolsup), while receiving the less family support (famsup), have the least involvement in extra-curricular activities (activities), has the least internet access at home (internet), socialize the least (goout), and lowest alcohol consumption. (Dalc, Walc) This cluster is the unsupported non-alcoholic loners of the sample

Cluster 3 has the most educated parents (Medu, Fedu), live in an urban setting (address) choose their school based on reputation (reasonreputation), have the lowest travel time to school (traveltime), study the most (studytime), rarely fail a course (failures), have the lowest support from the school while having the highest family support and family relationship (schoolsup, famsup, famrel), most involve in extra-curricular activities (activities), best internet acess at home (internet), least amount of free time (freetime) low alcohol consumption (Dalc, Walc). This cluster represents the affluent non-alcoholic high achievers.

CLuster 4 is the oldest (age), live in rural setting (address), has the smallest families (famsize), the least educated parents (Medu, Fedu), spends the most time traveling to school (traveltime), doesnt study much (studytime), has the highest failure rate (failures), never pays for extra classes (paid), most likely to be in a relationship (romanticyes), consumes alcohol moderately (Dalc, Walc), does poorly in school (G3). These students are the students in the greatest academic danger.

To get better insights, we can add the cluster results to our original dataset that was not normalize we can then identify what cluster each student belongs to ad calculate unstandardized means if we wanted.

student_alcohol$cluster<-alcohol_cluster$cluster # add clusters back to original data normalize does not mean much
View(student_alcohol)
aggregate(data=student_alcohol, G3~cluster, mean)
##   cluster   G3
## 1       1 10.8
## 2       2 11.1
## 3       3 14.5
## 4       4  5.5

The “aggregate function” tells us the average for each cluster on this question. We could do this for all of our variables to learn more about our clusters.

Kmeans provides a researcher with an understanding of the homogeneous characteristics of individuals within a sample. This information can be used to develop intervention plans in education or marketing plans in business. As such, kmeans is another powerful tool in machine learning

# K-Means Clustering

There are times in research when you neither want to predict nor classify examples. Rather, you want to take a dataset and segment the examples within the dataset so that examples with similar characteristics are gather together.

When your goal is to great several homogeneous groups within a sample it involves clustering. One type of clustering used in machine learning is k-means clustering. In this post we will learn the following about k-means clustering.

• The purpose of k-means
• Pros and cons of k-means

The Purpose of K-means

K-means has a number of applications. In the business setting, k-means has been used to segment customers. Businesses use this information to adjusts various aspects of their strategy for reaching their customers.

Another purpose for k-means is in simplifying large datasets so that you have several smaller datasets with similar characteristics. This subsetting could be useful in finding distinct patterns

K-means is a form of unsupervised classification. This means that the results label examples that the researcher must give meaning too. When R gives the results of an analysis it just labels the clusters as 1,2,3 etc. It is the researchers job to look at the clusters and give a qualitative meaning to them.

Pros and Cons of K-Means

The pros of k-means is that it is simple, highly flexible, and efficient. The simplicity of k-means makes it easy to explain the results in contrast to artificial neural networks or support vector machines. The flexibility of k-means allows for easy adjust if there are problems. Lastly, the efficiency of k-means implies that the algorithm is good at segmenting a dataset.

Some drawbacks to k-means is that it does not allows develop the most optimal set of clusters and that the number of clusters to make must be decided before the analysis. When doing the analysis, the k-means algorithm will randomly selecting several different places from which to develop clusters. This can be good or bad depending on where the algorithm chooses to begin at. From there, the center of the clusters is recalculated until an adequate “center” is found for the number of clusters requested.

How many clusters to include is left at the discretion of the researcher. This involves a combination of common sense, domain knowledge, and statistical tools. Too many clusters tells you nothing because the groups becoming very small and there are too many of them.

There are statistical tools that measure within group homogeneity and with group heterogeneity. In addition, there is a technique called a dendrogram. The results of a dendrogram analysis provides a recommendation of how many clusters to use. However, calculating a dendrogram for a large dataset could potential crash a computer due to the computational load and the limits of RAM.

Conclusion.

K-means is an analytical tool that helps to separate apples from oranges to give you one example. If you are in need of labeling examples based on the features in the dataset this method can be useful.

# Market Basket Analysis in R

In this post, we will conduct a market basket analysis on the shopping habits of people at a grocery store. Remember that a market basket analysis provides insights through indicating relationships among items that are commonly purchased together.

The first thing we need to do is load the package that makes association rules, which is the “arules” package. Next, we need to load our dataset groceries. This dataset is commonly used as a demonstration for market basket analysis.

However, you don’t won’t to load this dataset as dataframe because it leads to several technical issues during the analysis. Rather you want to load it as a sparse matrix. The function for this is “read.transactions” and is available in the “arules” pacakge

library(arules)
## Loading required package: Matrix
##
## Attaching package: 'arules'
##
## The following objects are masked from 'package:base':
##
##     abbreviate, write
#make sparse matrix
groceries<-read.transactions("/home/darrin/Documents/R working directory/Machine-Learning-with-R-datasets-master/groceries.csv", sep = ",")

Please keep in mind that the location of the file on your computer will be different from my hard drive.

We will now explore the data set by using several different functions. First, we will use the “summary” function as indicated below.

summary(groceries)
## transactions as itemMatrix in sparse format with
##  9835 rows (elements/itemsets/transactions) and
##  169 columns (items) and a density of 0.02609146
##
## most frequent items:
##       whole milk other vegetables       rolls/buns             soda
##             2513             1903             1809             1715
##           yogurt          (Other)
##             1372            34055
##
## element (itemset/transaction) length distribution:
## sizes
##    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
## 2159 1643 1299 1005  855  645  545  438  350  246  182  117   78   77   55
##   16   17   18   19   20   21   22   23   24   26   27   28   29   32
##   46   29   14   14    9   11    4    6    1    1    1    1    3    1
##
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   1.000   2.000   3.000   4.409   6.000  32.000
##
## includes extended item information - examples:
##             labels
## 1 abrasive cleaner
## 2 artif. sweetener
## 3   baby cosmetics

The output tells us the number of rows in our dataset (9835) columns (169) as well as the density, which is the percentage of columns that are not empty (2.6%). This may seem small but remember that the number of purchases varies from person to person so this affects how many empty columns there are.

Next, we have the most commonly purchased items. Milk and other vegetables were the two most common followed by other foods. After the most frequent items we have the size of each transaction. For example, 2159 people purchased one item during a transaction. While one person purchased 32 items in a transaction.

Lastly, we summary statistics about transactions. On average, a person would purchased 4.4 items per transaction.

We will now look at the support of different items. Remember, that the support is the frequency of an item in the dataset. We will use the “itemFrequencyPlot” function to do this and we will add the argument “topN” to sort the items from most common to less for the 15 most frequent transactions. Below is the code

itemFrequencyPlot(groceries, topN=15)

The plot that is produce gives you an idea of what people were purchasing. We will now attempt to develop association rules using the “apriori” function.

For now we will use the default settings for support and confidence (confidence is the proportion of transactions that have they same item(s)). The default for support is 0.1 and for confidence it is 0.8. Below is the code.

apriori(groceries)
## Apriori
##
## Parameter specification:
##  confidence minval smax arem  aval originalSupport support minlen maxlen
##         0.8    0.1    1 none FALSE            TRUE     0.1      1     10
##  target   ext
##   rules FALSE
##
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
##
## Absolute minimum support count: 983
##
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[169 item(s), 9835 transaction(s)] done [0.00s].
## sorting and recoding items ... [8 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 done [0.00s].
## writing ... [0 rule(s)] done [0.00s].
## creating S4 object  ... done [0.00s].
## set of 0 rules

As you can see from the printout, nothing meets the criteria of a support of 0.1 and confidence of 0.8. How to play with these numbers is a matter of experience as there are few strong rules for this matter. Below, I set the support to 0.006, confidence to 0.25, and the minimum number of rules items to 2. The support of 0.006 means that this item must have been purchased at least 60 times out of 9835 items and the confidence of 0.25 means that rule needs to be accurate 25% of the time. Lastly, I want at least two items in each rule that is produce as indicated by minlen = 2. Below is the code with the “summary” as well.

groceriesrules<-apriori(groceries, parameter = list(support=0.006, confidence = 0.25, minlen=2))
## Apriori
##
## Parameter specification:
##  confidence minval smax arem  aval originalSupport support minlen maxlen
##        0.25    0.1    1 none FALSE            TRUE   0.006      2     10
##  target   ext
##   rules FALSE
##
## Algorithmic control:
##  filter tree heap memopt load sort verbose
##     0.1 TRUE TRUE  FALSE TRUE    2    TRUE
##
## Absolute minimum support count: 59
##
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[169 item(s), 9835 transaction(s)] done [0.00s].
## sorting and recoding items ... [109 item(s)] done [0.00s].
## creating transaction tree ... done [0.00s].
## checking subsets of size 1 2 3 4 done [0.00s].
## writing ... [463 rule(s)] done [0.00s].
## creating S4 object  ... done [0.00s].
summary(groceriesrules)
## set of 463 rules
##
## rule length distribution (lhs + rhs):sizes
##   2   3   4
## 150 297  16
##
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   2.000   2.000   3.000   2.711   3.000   4.000
##
## summary of quality measures:
##     support           confidence          lift
##  Min.   :0.006101   Min.   :0.2500   Min.   :0.9932
##  1st Qu.:0.007117   1st Qu.:0.2971   1st Qu.:1.6229
##  Median :0.008744   Median :0.3554   Median :1.9332
##  Mean   :0.011539   Mean   :0.3786   Mean   :2.0351
##  3rd Qu.:0.012303   3rd Qu.:0.4495   3rd Qu.:2.3565
##  Max.   :0.074835   Max.   :0.6600   Max.   :3.9565
##
## mining info:
##       data ntransactions support confidence
##  groceries          9835   0.006       0.25

Are current analysis has 463 rules. This is a major improvement from 0. We can also see how many rules contain 2 (150), 3 (297), and 4 items (16). We also have a summary of of the average number of items per rule. The next is descriptive stats on the support and confidence of the rules generated.

Something that is new for us is the “lift” column. Lift is a measure how much more likely an item is to be purchased above its number rate. Anything above 1 means that the likelihood of purchase is higher than chance.

We are now going to look for useful rules from our dataset. We are looking for rules that we can use to make decisions. It takes industry experiences in the field of your data to really glean useful rules. For now, we can only determine this statistically by sorting the rules by their lift. Below is the code for this

inspect(sort(groceriesrules, by="lift")[1:7])
##   lhs             rhs                     support    confidence  lift
## 1 {herbs}         => {root vegetables}    0.007015760  0.4312500 3.956
## 2 {berries}       => {whipped/sour cream} 0.009049314  0.2721713 3.796
## 3 {other vegetables,                                                  ##    tropical fruit,
##    whole milk}    => {root vegetables}    0.007015760  0.4107143 3.768
## 4 {beef,
##    other vegetables} => {root vegetables} 0.007930859  0.4020619 3.688
## 5 {other vegetables,
##    tropical fruit} => {pip fruit}         0.009456024  0.2634561 3.482
## 6 {beef,
##    whole milk}    => {root vegetables}    0.008032537  0.3779904 3.467
## 7 {other vegetables,
##    pip fruit}     => {tropical fruit}     0.009456024  0.3618677 3.448

The first three rules rules are translated into simply English below as ordered by lift.

1. If herbs are purchased then root vegetables are purchased
2. If berries are purchased then whipped sour cream is purchased
3. If other vegetables, tropical fruit and whole milk are purchased then root vegetables are purchased

Conclusion

Since we are making no predictions, there is no way to really objectively improve the model. This is normal when the learning is unsupervised. If we had to make a recommendation based on the results we could say that the store should place all vegetables near each other.

The power of market basket analysis is allowing the researcher to identify relationships that may not have been noticed any other way. Naturally, insights gained from this approach must be use for practical actions in the setting in which they apply.

# Support Vector Machines in R

In this post, we will use support vector machine analysis to look at some data available on kaggle. In particular, we will predict what number a person wrote by analyze the pixels that were used to make the number. The file for this example is available at https://www.kaggle.com/c/digit-recognizer/data

To do this analysis you will need to use the ‘kernlab’ package. While playing with this dataset I noticed a major problem, doing the analysis with the full data set of 42000 examples took forever. To alleviate this problem. We are going to practice with a training set of 7000 examples and a test set of 3000. Below is the code for the first few sets. Remember that the dataset was download separately

#load packages
library(kernlab)
#split data
digitRedux<-digitTrain[1:7000,]
digitReduxTest<-digitTrain[7001:10000,]
#explore data
str(digitRedux)
## 'data.frame':    7000 obs. of  785 variables:
##  $label : int 1 0 1 4 0 0 7 3 5 3 ... ##$ pixel0  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $pixel1 : int 0 0 0 0 0 0 0 0 0 0 ... ##$ pixel2  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $pixel3 : int 0 0 0 0 0 0 0 0 0 0 ... ##$ pixel4  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $pixel5 : int 0 0 0 0 0 0 0 0 0 0 ... ##$ pixel6  : int  0 0 0 0 0 0 0 0 0 0 ...
##   [list output truncated]

From the “str” function you can tell we have a lot of variables (785). This is what slowed the analysis down so much when I tried to run the full 42000 examples in the original dataset.

SVM need a factor variable as the predictor if possible. We are trying to predict the “label” variable so we are going to change this to a factor variable because that is what it really is. Below is the code

#convert label variable to factor
digitRedux$label<-as.factor(digitRedux$label)
digitReduxTest$label<-as.factor(digitReduxTest$label)

Before we continue with the analysis we need to scale are variables. This makes all variables to be within the same given range which helps to equalizes the influence of them. However, we do not want to change our “label” variable as this is the predictor variable and scaling it would make the results hard to understand. Therefore, we are going to temporarily remove the “label” variable from both of our data sets and save them in a temporary data frame. The code is below.

#temporary dataframe for the label results
keep<-as.data.frame(digitRedux$label) keeptest<-as.data.frame(digitReduxTest$label)
#null label variable in both datasets
digitRedux$label<-NA digitReduxTest$label<-NA

Next, we scale the remaining variable and reinsert the label variables for each data set as show in our code below.

digitRedux<-as.data.frame(scale(digitRedux))
digitRedux[is.na(digitRedux)]<- 0 #replace NA with 0
digitReduxTest<-as.data.frame(scale(digitReduxTest))
digitReduxTest[is.na(digitReduxTest)]<- 0
digitRedux$label<-keep$digitRedux$label digitReduxTest$label<-keeptest$digitReduxTest$label

Now we make our model using the “ksvm” function in the “kernlab” package. We set the kernel to “vanilladot” which is a linear kernel. We will aslo print the results. However, the results do not make any sense on their own and the model can only be assess througn other means. Below is the code. If you get a warning message about scaling do not worry about this as we scaled the data ourselves.

#make the model
number_classify<-ksvm(label~.,  data=digitRedux,
kernel="vanilladot")
##  Setting default kernel parameters
## Warning in .local(x, ...): Variable(s) ' constant. Cannot scale data.
#look at the results
number_classify
## Support Vector Machine object of class "ksvm"
##
## SV type: C-svc  (classification)
##  parameter : cost C = 1
##
## Linear (vanilla) kernel function.
##
## Number of Support Vectors : 2218
##
## Objective Function Value : -0.0623 -0.207 -0.1771 -0.0893 -0.3207 -0.4304 -0.0764 -0.2719 -0.2125 -0.3575 -0.2776 -0.1618 -0.3408 -0.1108 -0.2766 -1.0657 -0.3201 -1.0509 -0.2679 -0.4565 -0.2846 -0.4274 -0.8681 -0.3253 -0.1571 -2.1586 -0.1488 -0.2464 -2.9248 -0.5689 -0.2753 -0.2939 -0.4997 -0.2429 -2.336 -0.8108 -0.1701 -2.4031 -0.5086 -0.0794 -0.2749 -0.1162 -0.3249 -5.0495 -0.8051
## Training error : 0

We now need to use the “predict” function so that we can determine the accuracy of our model. Remember that for predicting, we use the answers in the test data and compare them to what our model would guess based on what it knows.

number_predict<-predict(number_classify, digitReduxTest)
table(number_predict, digitReduxTest$label) ## ## number_predict 0 1 2 3 4 5 6 7 8 9 ## 0 297 0 3 3 1 4 6 1 1 1 ## 1 0 307 4 1 0 4 2 5 11 1 ## 2 0 2 268 10 5 1 3 10 12 3 ## 3 0 1 7 291 1 11 0 1 8 3 ## 4 0 1 3 0 278 4 3 2 0 9 ## 5 2 0 1 10 1 238 4 1 11 1 ## 6 2 1 1 0 2 1 287 1 0 0 ## 7 0 1 1 0 1 0 0 268 3 10 ## 8 1 3 4 10 0 11 1 0 236 2 ## 9 0 0 2 2 9 2 0 14 2 264 accuracy<-number_predict == digitReduxTest$label
prop.table(table(accuracy))
## accuracy
##      FALSE       TRUE
## 0.08866667 0.91133333

The table allows you to see how many were classified correctly and how they were misclassified. The prop.table allows you to see an overall percentage. This particular model was highly accurate at 91%. It would be difficult to improve further. Below is code for a model that is using a different kernel with results that are barely better. However, if you ever enter a data science competition any improve ususally helps even if it is not practical for everyday use.

number_classify_rbf<-ksvm(label~.,  data=digitRedux,
kernel="rbfdot")
## Warning in .local(x, ...): Variable(s) ' constant. Cannot scale data.
#evaluate improved model
number_predict_rbf<-predict(number_classify_rbf, digitReduxTest)
table(number_predict_rbf, digitReduxTest$label) ## ## number_predict_rbf 0 1 2 3 4 5 6 7 8 9 ## 0 294 0 2 3 1 1 2 1 1 0 ## 1 0 309 1 0 1 0 1 3 4 2 ## 2 4 2 277 12 4 4 8 9 9 8 ## 3 0 1 3 297 1 3 0 1 3 3 ## 4 0 1 3 0 278 4 1 5 0 6 ## 5 0 0 1 6 1 254 4 0 9 1 ## 6 2 1 0 0 2 6 289 0 2 0 ## 7 0 0 3 2 3 0 0 277 2 13 ## 8 2 2 4 3 0 3 1 0 253 3 ## 9 0 0 0 4 7 1 0 7 1 258 accuracy_rbf<-number_predict_rbf == digitReduxTest$label
prop.table(table(accuracy_rbf))
## accuracy_rbf
##      FALSE       TRUE
## 0.07133333 0.92866667


Conclusion

From this demonstration we can see the power of support vector machines with numerical data. This type of analysis can be used for things beyond the conventional analysis and can be used to predict things such as hand written numbers. As such, SVM is yet another tool available for the data scientist.

# Developing an Artificial Neural Network in R

In this post, we are going make an artificial neural network (ANN) by analyzing some data about computers. Specifically, we are going to make an ANN to predict the price of computers.

We will be using the “ecdat” package and the data set “Computers” from this pacakge. In addition, we are going to use the “neuralnet” pacakge to conduct the ANN analysis. Below is the code for the packages and dataset we are using

library(Ecdat);library(neuralnet)
## Loading required package: Ecfun
##
## Attaching package: 'Ecdat'
##
## The following object is masked from 'package:datasets':
##
##     Orange
##
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:Ecdat':
##
##     SP500
#load data set
data("Computers")

Explore the Data

The first step is always data exploration. We will first look at nature of the data using the “str” function and then used the “summary” function. Below is the code.

str(Computers)
## 'data.frame':    6259 obs. of  10 variables:
##  $price : num 1499 1795 1595 1849 3295 ... ##$ speed  : num  25 33 25 25 33 66 25 50 50 50 ...
##  $hd : num 80 85 170 170 340 340 170 85 210 210 ... ##$ ram    : num  4 2 4 8 16 16 4 2 8 4 ...
##  $screen : num 14 14 15 14 14 14 14 14 14 15 ... ##$ cd     : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 2 1 1 1 ...
##  $multi : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ... ##$ premium: Factor w/ 2 levels "no","yes": 2 2 2 1 2 2 2 2 2 2 ...
##  $ads : num 94 94 94 94 94 94 94 94 94 94 ... ##$ trend  : num  1 1 1 1 1 1 1 1 1 1 ...
lapply(Computers, summary)
## $price ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 949 1794 2144 2220 2595 5399 ## ##$speed
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   25.00   33.00   50.00   52.01   66.00  100.00
##
## $hd ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 80.0 214.0 340.0 416.6 528.0 2100.0 ## ##$ram
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   2.000   4.000   8.000   8.287   8.000  32.000
##
## $screen ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 14.00 14.00 14.00 14.61 15.00 17.00 ## ##$cd
##   no  yes
## 3351 2908
##
## $multi ## no yes ## 5386 873 ## ##$premium
##   no  yes
##  612 5647
##
## $ads ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 39.0 162.5 246.0 221.3 275.0 339.0 ## ##$trend
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##    1.00   10.00   16.00   15.93   21.50   35.00

ANN is primarily for numerical data and not categorical or factors. As such, we will remove the factor variables cd,multi, and premium from further analysis. Below are is the code and histograms of the remaining variables.

hist(Computers$price);hist(Computers$speed);hist(Computers$hd);hist(Computers$ram)

hist(Computers$screen);hist(Computers$ads);hist(Computers$trend) Clean the Data Looking at the summary combined with the histograms indicates that we need to normalize are data as this is a requirement for ANN. We want all of the variables to have an equal influence initially. Below is the code for the function we wil use to normalize are variables. normalize<-function(x) { return((x-min(x)) / (max(x))) } Explore the Data Again We now need to make a new dataframe that has only the variables we are going to use for the analysis. Then we will use our “normalize” function to scale and center the variables appropriately. Lastly, we will re-explore the data to make sure it is ok using the “str” “summary” and “hist” functions. Below is the code. #make dataframe without factor variables Computers_no_factors<-data.frame(Computers$price,Computers$speed,Computers$hd,
Computers$ram,Computers$screen,Computers$ad, Computers$trend)
#make a normalize dataframe of the data
Computers_norm<-as.data.frame(lapply(Computers_no_factors, normalize))
#reexamine the normalized data
lapply(Computers_norm, summary)
## $Computers.price ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.0000 0.1565 0.2213 0.2353 0.3049 0.8242 ## ##$Computers.speed
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##  0.0000  0.0800  0.2500  0.2701  0.4100  0.7500
##
## $Computers.hd ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.00000 0.06381 0.12380 0.16030 0.21330 0.96190 ## ##$Computers.ram
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##  0.0000  0.0625  0.1875  0.1965  0.1875  0.9375
##
## $Computers.screen ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.00000 0.00000 0.00000 0.03581 0.05882 0.17650 ## ##$Computers.ad
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##  0.0000  0.3643  0.6106  0.5378  0.6962  0.8850
##
## $Computers.trend ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.0000 0.2571 0.4286 0.4265 0.5857 0.9714 hist(Computers_norm$Computers.price);hist(Computers_norm$Computers.speed); hist(Computers_norm$Computers.hd);hist(Computers_norm$Computers.ram); hist(Computers_norm$Computers.screen);hist(Computers_norm$Computers.ad); hist(Computers_norm$Computers.trend)

Develop the Model

Everything looks good, so we will now split are data into a training and testing set and develop the ANN model that will predict computer price. Below is the code for this.

#split data into training and testing
Computer_train<-Computers_norm[1:4694,]
Computers_test<-Computers_norm[4695:6259,]
#model
computer_model<-neuralnet(Computers.price~Computers.speed+Computers.hd+
Computers.ram+Computers.screen+
Computers.ad+Computers.trend, Computer_train)

Our intial model is a simple feedforward networm with a single hidden node. You can visualize the model using the “plot” function as shown in the code below.

plot(computer_model)

We now need to evaluate our model’s performance. We will use the “compute” function to generate predictions. The predictions generated by the “compute” function will then be compared to the actual prices in the test data set using a pearson correlation. Since we are not classifying we can’t measure accuracy with a confusion matrix but rather with correlation. Below is the code followed by the results

evaluate_model<-compute(computer_model, Computers_test[2:7])
predicted_price<-evaluate_model$net.result cor(predicted_price, Computers_test$Computers.price)
##              [,1]
## [1,] 0.8809571295

The correlation between the predict results and the actual results is 0.88 which is a strong relationship. THis indicates that our model does an excellent job in predicting the price of a computer based on ram, screen size, speed, hard drive size, advertising, and trends.

Develop Refined Model

Just for fun we are going to make a more complex model with three hidden nodes and see how the results change below is the code.

computer_model2<-neuralnet(Computers.price~Computers.speed+Computers.hd+
Computers.ram+Computers.screen+
Computer_train, hidden =3)
plot(computer_model2)
evaluate_model2<-compute(computer_model2, Computers_test[2:7])
predicted_price2<-evaluate_model2$net.result cor(predicted_price2, Computers_test$Computers.price)


##              [,1]
## [1,] 0.8963994092

The correlation improves to 0.89. As such, the increased complexity did not yield much of an improvement in the overall correlation. Therefore, a single node model is more appropriate.

Conclusion

In this post we explored an application of artificial neural networks. This black box method is useful for making powerful predictions in highly complex data.

# Making Regression and Modal Trees in R

In this post, we will look at an example of regression trees. Regression trees use decision tree-like approach to develop predicition models involving numerical data. In our example we will be trying to predict how many kids a person has based on several independent variables in the “PSID” data set in the “Ecdat” package.

Lets begin by loading the necessary packages and data set. The code is below

library(Ecdat);library(rpart);library(rpart.plot)
library(RWeka)
data(PSID)
str(PSID)
## 'data.frame':    4856 obs. of  8 variables:
##  $intnum : int 4 4 4 4 5 6 6 7 7 7 ... ##$ persnum : int  4 6 7 173 2 4 172 4 170 171 ...
##  $age : int 39 35 33 39 47 44 38 38 39 37 ... ##$ educatn : int  12 12 12 10 9 12 16 9 12 11 ...
##  $earnings: int 77250 12000 8000 15000 6500 6500 7000 5000 21000 0 ... ##$ hours   : int  2940 2040 693 1904 1683 2024 1144 2080 2575 0 ...
##  $kids : int 2 2 1 2 5 2 3 4 3 5 ... ##$ married : Factor w/ 7 levels "married","never married",..: 1 4 1 1 1 1 1 4 1 1 ...
summary(PSID)
##      intnum        persnum            age           educatn
##  Min.   :   4   Min.   :  1.00   Min.   :30.00   Min.   : 0.00
##  1st Qu.:1905   1st Qu.:  2.00   1st Qu.:34.00   1st Qu.:12.00
##  Median :5464   Median :  4.00   Median :38.00   Median :12.00
##  Mean   :4598   Mean   : 59.21   Mean   :38.46   Mean   :16.38
##  3rd Qu.:6655   3rd Qu.:170.00   3rd Qu.:43.00   3rd Qu.:14.00
##  Max.   :9306   Max.   :205.00   Max.   :50.00   Max.   :99.00
##                                                  NA's   :1
##     earnings          hours           kids                 married
##  Min.   :     0   Min.   :   0   Min.   : 0.000   married      :3071
##  1st Qu.:    85   1st Qu.:  32   1st Qu.: 1.000   never married: 681
##  Median : 11000   Median :1517   Median : 2.000   widowed      :  90
##  Mean   : 14245   Mean   :1235   Mean   : 4.481   divorced     : 645
##  3rd Qu.: 22000   3rd Qu.:2000   3rd Qu.: 3.000   separated    : 317
##  Max.   :240000   Max.   :5160   Max.   :99.000   NA/DF        :   9
##                                                   no histories :  43

The variables “intnum” and “persnum” are for identification and are useless for our analysis. We will now explore our dataset with the following code.

hist(PSID$age) hist(PSID$educatn)

hist(PSID$earnings) hist(PSID$hours)

hist(PSID$kids) table(PSID$married)
##
##       married never married       widowed      divorced     separated
##          3071           681            90           645           317
##         NA/DF  no histories
##             9            43

Almost all of the variables are non-normal. However, this is not a problem when using regression trees. There are some major problems with the “kids” and “educatn” variables. Each of these variables have values at 98 and 99. When the data for this survey was collected 98 meant the respondent did not know the answer and a 99 means they did not want to say. Since both of these variables are numerical we have to do something with them so they do not ruin our analysis.

We are going to recode all values equal to or greater than 98 as 3 for the “kids” variable. The number 3 means they have 3 kids. This number was picked because it was the most common response for the other respondents. For the “educatn” variable all values equal to or greater than 98 are recoded as 12, which means that they completed 12th grade. Again this was the most frequent response. Below is the code.

PSID$kids[PSID$kids >= 98] <- 3
PSID$educatn[PSID$educatn >= 98] <- 12

Another peek at the histograms for these two variables and things look much better.

hist(PSID$kids) hist(PSID$educatn)

Make Model and Visualization

Now that everything is cleaned up we now need to make are training and testing data sets as seen in the code below.

PSID_train<-PSID[1:3642,]
PSID_test<-PSID[3643:4856,]

We will now make our model and also create a visual of it. Our goal is to predict the number of children a person has based on their age, education, earnings, hours worked, marital status. Below is the code

#make model
PSID_Model<-rpart(kids~age+educatn+earnings+hours+married, PSID_train)
#make visualization
rpart.plot(PSID_Model, digits=3, fallen.leaves = TRUE,type = 3, extra=101)

The first split on the tree is by income. On the left we have those who make more than 20k and on the right those who make less than 20k. On the left the next split is by marriage, those who are never married or not applicable have on average 0.74 kids. Those who are married, widowed, divorced, separated, or have no history have on average 1.72.

The left side of the tree is much more complicated and I will not explain all of it. The after making less than 20k the next split is by marriage. Those who are married, widowed, divorced, separated, or no history with less than 13.5 years of education have 2.46 on average.

Make Prediction Model and Conduct Evaluation

Our next task is to make the prediction model. We will do this with the folllowing code

PSID_pred<-predict(PSID_Model, PSID_test)

We will now evaluate the model. We will do this three different ways. The first involves looking at the summary statistics of the prediction model and the testing data. The numbers should be about the same. After that we will calculate the correlation between the prediction model and the testing data. Lastly, we will use a technique called the mean absolute error. Below is the code for the summary statistics and correlation.

summary(PSID_pred)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   0.735   2.041   2.463   2.226   2.463   2.699
summary(PSID_test$kids) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.000 2.000 2.000 2.494 3.000 10.000 cor(PSID_pred, PSID_test$kids)
## [1] 0.308116

Looking at the summary stats are model has a hard time predicting extreme values because the max value of the two models are vey far a part. However, how often do people have ten kids? As such, this is not a major concern.

A look at the correlation finds that it is pretty low (0.30) this means that the two models have little in common and this means we need to make some changes. The mean absolute error is a measure of the difference between the predicted and actual values in a model. We need to make a function first before we analyze our model.

MAE<-function(actual, predicted){
mean(abs(actual-predicted))
}

We now assess the model with the code below

MAE(PSID_pred, PSID_test$kids) ## [1] 1.134968 The results indicate that on average the difference between our model’s prediction of the number of kids and the actual number of kids was 1.13 on a scale of 0 – 10. That’s a lot of error. However, we need to compare this number to how well the mean does to give us a benchmark. The code is below. ave_kids<-mean(PSID_train$kids)
MAE(ave_kids, PSID_test$kids) ## [1] 1.178909 Model Tree Our model with a score of 1.13 is slightly better than using the mean which is 1.17. We will try to improve our model by switching from a regression tree to a model tree which uses a slightly different approach for prediction. In a model tree each node in the tree ends in a linear regression model. Below is the code. PSIDM5<- M5P(kids~age+educatn+earnings+hours+married, PSID_train) PSIDM5 ## M5 pruned model tree: ## (using smoothed linear models) ## ## earnings <= 20754 : ## | earnings <= 2272 : ## | | educatn <= 12.5 : LM1 (702/111.555%) ## | | educatn > 12.5 : LM2 (283/92%) ## | earnings > 2272 : LM3 (1509/88.566%) ## earnings > 20754 : LM4 (1147/82.329%) ## ## LM num: 1 ## kids = ## 0.0385 * age ## + 0.0308 * educatn ## - 0 * earnings ## - 0 * hours ## + 0.0187 * married=married,divorced,widowed,separated,no histories ## + 0.2986 * married=divorced,widowed,separated,no histories ## + 0.0082 * married=widowed,separated,no histories ## + 0.0017 * married=separated,no histories ## + 0.7181 ## ## LM num: 2 ## kids = ## 0.002 * age ## - 0.0028 * educatn ## + 0.0002 * earnings ## - 0 * hours ## + 0.7854 * married=married,divorced,widowed,separated,no histories ## - 0.3437 * married=divorced,widowed,separated,no histories ## + 0.0154 * married=widowed,separated,no histories ## + 0.0017 * married=separated,no histories ## + 1.4075 ## ## LM num: 3 ## kids = ## 0.0305 * age ## - 0.1362 * educatn ## - 0 * earnings ## - 0 * hours ## + 0.9028 * married=married,divorced,widowed,separated,no histories ## + 0.2151 * married=widowed,separated,no histories ## + 0.0017 * married=separated,no histories ## + 2.0218 ## ## LM num: 4 ## kids = ## 0.0393 * age ## - 0.0658 * educatn ## - 0 * earnings ## - 0 * hours ## + 0.8845 * married=married,divorced,widowed,separated,no histories ## + 0.3666 * married=widowed,separated,no histories ## + 0.0037 * married=separated,no histories ## + 0.4712 ## ## Number of Rules : 4 It would take too much time to explain everything. You can read part of this model as follows earnings greater than 20754 use linear model 4earnings less than 20754 and less than 2272 and less than 12.5 years of education use linear model 1 earnings less than 20754 and less than 2272 and greater than 12.5 years of education use linear model 2 earnings less than 20754 and greater than 2272 linear model 3 The print out then shows each of the linear model. Lastly, we will evaluate our model tree with the following code PSIDM5_Pred<-predict(PSIDM5, PSID_test) summary(PSIDM5_Pred) ## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 0.3654 2.0490 2.3400 2.3370 2.6860 4.4220 cor(PSIDM5_Pred, PSID_test$kids)
## [1] 0.3486492
MAE(PSID_test\$kids, PSIDM5_Pred)
## [1] 1.088617

This model is slightly better. FOr example, it is better at predict extreme values at 4.4 compare to 2.69 for the regression tree model. The correlation is 0.34 which is better than 0.30 for the regression tree model. Lastly. the mean absolute error shows a slight improve to 1.08 compared to 1.13 in the regression tree model

Conclusion

This provide examples of the use of regression trees and model trees. Both of these models make prediction a key component of their analysis.