<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/machine-learning-with-r-caret-part-1/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://datascienceplus.com/machine-learning-with-r-caret-part-1/"> R Programming – DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This blog post series is on machine learning with R. We will use the Caret package in R. In this part, we will first perform exploratory Data Analysis (EDA) on a real-world dataset, and then apply non-regularized linear regression to solve a supervised regression problem on the dataset. We will predict power output given a set of environmental readings from various sensors in a natural gas-fired power generation plant. </p>
<p>In the second part of the post, we will work with regularized linear regression models (ridge, lasso and elasticnet). Next, we will see the other non-linear regression models. The real-world <a rel="nofollow" href="https://archive.ics.uci.edu/ml/datasets/Combined+Cycle+Power+Plant" target="_blank">data we are using</a> in this post consists of 9,568 data points, each with 4 environmental attributes collected from a Combined Cycle Power Plant over 6 years (2006-2011), and is provided by the University of California, Irvine at UCI Machine Learning Repository Combined Cycle Power Plant Data Set. You can find more details about the dataset on the UCI page. The task is a regression problem since the label (or target) we are trying to predict is numeric.</p>
<h2>Libraries</h2>
<p>The caret package (short for _C_lassification _A_nd _RE_gression _T_raining) is a set of functions that attempt to streamline the process for creating predictive models.</p>
<pre>library(readxl)
library(ggplot2)
library(corrplot)
library(tidyverse)</pre>
<h2>Load Data</h2>
<pre>power_plant = read_excel(&quot;Folds5x2_pp.xlsx&quot;)</pre>
<h2>Exploratory Data Analysis (EDA)</h2>
<p>This is a step that we should always perform before trying to fit a model to the data, as this step will often lead to important insights about our data.</p>
<pre>class(power_plant)
<em>'tbl_df' 'tbl' 'data.frame'</em></pre>
<pre>
# Caret faces problems working with tbl, 
# so let's change the data to simple data frame
power_plant = data.frame(power_plant)
message(&quot;The class is now &quot;, class(power_plant))
# See first few rows
head(power_plant)
<em>AT	V	AP	RH	PE
14.96	41.76	1024.07	73.17	463.26
25.18	62.96	1020.04	59.08	444.37
5.11	39.40	1012.16	92.14	488.56
20.86	57.32	1010.24	76.64	446.48
10.82	37.50	1009.23	96.62	473.90
26.27	59.44	1012.23	58.77	443.67</em></pre>
<p>The columns in the DataFrame are:</p>
<ul>AT = Atmospheric Temperature in C<br>
V = Exhaust Vacuum Speed<br>
AP = Atmospheric Pressure<br>
RH = Relative Humidity<br>
PE = Power Output</ul>
<p>Power Output is the value we are trying to predict given the measurements above.</p>
<pre># Size of  DataFrame
dim(power_plant)     
<em>9568 5</em></pre>
<pre>map(power_plant, class) 
<em># all columns are numeric
$AT
'numeric'
$V
'numeric'
$AP
'numeric'
$RH
'numeric'
$PE
'numeric'</em></pre>
<p>Are there any missing values in any of the columns?</p>
<pre>map(power_plant, ~sum(is.na(.)))  # using purrr package
<em># There is no missing data in all of the columns
$AT
0
$V
0
$AP
0
$RH
0
$PE
0</em></pre>
<h2>Visualize relationship between variables</h2>
<p>Before we perform any modeling, it is a good idea to explore correlations between the predictors and the predictand. This step can be important as it helps us to select appropriate models. If our features and the outcome are linearly related, we may start with linear regression models. However, if the relationships between the label and the features are non-linear, non-linear ensemble models such as random forest can be better.</p>
<p>Correlation between atmospheric temperature and power output</p>
<pre>power_plant %>% ggplot(aes(AT, PE)) +
    geom_point(color= &quot;blue&quot;, alpha = 0.3) +
    ggtitle(&quot;Temperature vs Power Output&quot;) +
    xlab(&quot;Atmospheric Temperature&quot;) +
    ylab(&quot;Power Output&quot;) +
    theme(plot.title = element_text(color=&quot;darkred&quot;,
                                    size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
                      axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))</pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot-caret.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot-caret-490x490.png?w=450&ssl=1" alt="" class="wp-image-10572" data-recalc-dims="1"></a></p>
<p>As shown in the above figure, there is strong linear correlation between Atmospheric Temperature and Power Output.</p>
<p>Correlation between Exhaust Vacuum Speed and power output </p>
<pre>power_plant %>% ggplot(aes(V, PE)) +
    geom_point(color= &quot;darkgreen&quot;, alpha = 0.3) +
    ggtitle(&quot;Exhaust Vacuum Speed vs Power Output&quot;) +
    xlab(&quot;Exhaust Vacuum Speed&quot;) +
    ylab(&quot;Power Output&quot;) +
    theme(plot.title = element_text(color=&quot;darkred&quot;,size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
                      axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))</pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot2caret.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot2caret-490x490.png?w=450&ssl=1" alt="" class="wp-image-10574" data-recalc-dims="1"></a></p>
<p>Correlation between atmospheric pressure and power output </p>
<pre>power_plant %>% ggplot(aes(AP, PE)) +
    geom_point(color= &quot;red&quot;, alpha = 0.3) +
    ggtitle(&quot;Atmospheric Pressure vs Power Output&quot;) +
    xlab(&quot;Atmospheric Pressure&quot;) +
    ylab(&quot;Power Output&quot;) +
    theme(plot.title = element_text(color=&quot;darkred&quot;,size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
                      axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))</pre>
<p><a rel="nofollow" href="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot3caret.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot3caret-490x490.png?w=450&ssl=1" alt="" class="wp-image-10575" data-recalc-dims="1"></a></p>
<p>Correlation between relative humidity and power output </p>
<pre>power_plant %>% ggplot(aes(RH, PE)) +
    geom_point(color= &quot;magenta&quot;, alpha = 0.3) +
    ggtitle(&quot;Relative Humidity vs Power Output&quot;) +
    xlab(&quot;Relative Humidity&quot;) +
    ylab(&quot;Power Output&quot;) +
    theme(plot.title = element_text(color=&quot;darkred&quot;,size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
                      axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))</pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot4caret.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot4caret-490x490.png?w=450&ssl=1" alt="" class="wp-image-10576" data-recalc-dims="1"></a></p>
<h3>Correlation heatmap </h3>
<pre>correlations = cor(power_plant)
corrplot(correlations, method=&quot;color&quot;)</pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot5caret.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot5caret-490x490.png?w=450&ssl=1" alt="" class="wp-image-10577" data-recalc-dims="1"></a></p>
<p>As shown in the correlation heatmap above, the target is correlated with the features. However, we also observe correlation among the features, hence we have multi-collinearity problem. He will use regularization to check if the collinearity we observe has a significant impact on the performance of linear regression model.</p>
<h2>Data Modeling</h2>
<p>All the columns are numeric and there are no missing values, which makes our modeling task straightforward. Now, let’s model our data to predict what the power output will be given a set of sensor readings. Our first model will be based on simple linear regression since we saw some linear patterns in our data based on the scatter plots and correlation heatmap during the exploration stage. We need a way of evaluating how well our linear regression model predicts power output as a function of input parameters. We can do this by splitting up our initial data set into a Training Set, used to train our model and a Test Set, used to evaluate the model’s performance in giving predictions.</p>
<pre>library(caret)
# Caret streamlines the process for creating predictive models.
# The package contains tools for:
#  data splitting
#  pre-processing
#  feature selection
#  model tuning using resampling
#  variable importance estimation
#  as well as other functionality.</pre>
<p>The Caret <code>createDataPartition</code> function creates indices for training and test portion of the data for us.</p>
<pre># create training set indices with 80% of data
set.seed(100)  # For reproducibility
# Create index for testing and training data
inTrain &lt;- createDataPartition(y = power_plant$PE, 
                               p = 0.8, list = FALSE)
# subset power_plant data to training
training &lt;- power_plant[inTrain,]
# subset the rest to test
 testing &lt;- power_plant[-inTrain,]

# Size ratio of training and test dataset
message(&quot;As shown below, the training set is about 80%  and the test set is about 20% of the original data&quot;)

rbind(&quot;Training set&quot; = nrow(training)/nrow(power_plant),
      &quot;Testing set&quot; = nrow(testing)/nrow(power_plant)) %>% 
       round(2)*100
<em>Training set	80
Testing set	20</em></pre>
<p>As shown above, the training set is about 80% and the test set is about 20% of the original data</p>
<h2>Linear Regression</h2>
<pre># Fit linear regression model
# put the predictors on the same scale: mean of zero and unit variance
my_lm = train(training[,1:4], training[,5],
               method = &quot;lm&quot;,
               preProc = c(&quot;center&quot;, &quot;scale&quot;)
              )
message(&quot;Linear Regression: Model performance on \n the training set&quot;)
my_lm$results[c(&quot;RMSE&quot;,&quot;Rsquared&quot;)] %>%
        round(2)
summary(my_lm)
<em>Linear Regression: Model performance on the training set
RMSE	Rsquared
4.61	0.93
Call:
lm(formula = .outcome ~ ., data = dat)
Residuals:
Min 1Q Median 3Q Max -43.281 -3.164 -0.117 3.239 17.704 Coefficients:
Estimate Std. Error t value Pr(>|t|) (Intercept) 454.39620 0.05233 8683.725 &lt; 2e-16 ***
AT -14.62996 0.12728 -114.943 &lt; 2e-16 ***
V -3.03266 0.10347 -29.309 &lt; 2e-16 ***
AP 0.39110 0.06320 6.189 6.38e-10 ***
RH -2.23520 0.06804 -32.849 &lt; 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.579 on 7651 degrees of freedom
Multiple R-squared: 0.9281,	Adjusted R-squared: 0.9281 F-statistic: 2.469e+04 on 4 and 7651 DF, p-value: &lt; 2.2e-16</em></pre>
<p>Predict using the test data</p>
<pre>pred = predict(my_lm, testing[, 1:4])
SSE = sum((testing[,5] -pred)^2)    # sum of squared errors
SST = sum((testing[,5] - mean(training[,5]))^2) # total sum of squares, remember to use training data here
R_square = 1 - SSE/SST
message('R_squared on the test data:')
round(R_square, 2)
SSE = sum((testing[,5] - pred)^2)
RMSE = sqrt(SSE/length(pred))
message(&quot;Root mean square error on the test data: &quot;)
round(RMSE, 2)
<em>R_squared on the test data:
0.93
Root mean square error on the test data: 4.48</em></pre>
<p>Preparing data for ploting</p>
<pre>my_data = as.data.frame(cbind(predicted = pred,
                            observed = testing$PE))
# Plot predictions vs test data
ggplot(my_data,aes(predicted, observed)) +
      geom_point(color = &quot;darkred&quot;, alpha = 0.5) + 
      geom_smooth(method=lm)+ ggtitle('Linear Regression ') +
      ggtitle(&quot;Linear Regression: Prediction vs Test Data&quot;) +
      xlab(&quot;Predecited Power Output &quot;) +
      ylab(&quot;Observed Power Output&quot;) +
      theme(plot.title = element_text(color=&quot;darkgreen&quot;,size=18,hjust = 0.5),
                     axis.text.y = element_text(size=12),
          
            axis.text.x = element_text(size=12,hjust=.5),
                      axis.title.x = element_text(size=14),
                      axis.title.y = element_text(size=14))</pre>
<p>Gives this plot:<br>
<a rel="nofollow" href="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot6caret.png?ssl=1" target="_blank"><img sizes="(max-width: 490px) 100vw, 490px" src="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2018/01/plot6caret-490x490.png?w=450&ssl=1" alt="" class="wp-image-10578" data-recalc-dims="1"></a></p>
<p>This blog post series is on machine learning with R. In this first post, we saw non-regularized multivariate linear regression. In the second part of the post, we will work with regularized linear regression models. Next, we will see the other non-linear regression models. See you in the next post.</p>
<ul></ul>
<p><strong>
<p>Related Post</p>
<p></p></strong>
<ol>
<li><a rel="nofollow" href="https://datascienceplus.com/how-to-perform-logistic-regression-lda-qda-in-r/" title="Permanent Link toHow to perform Logistic Regression, LDA,  &amp; QDA in R" target="_blank">How to perform Logistic Regression, LDA, &amp; QDA in R</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/how-to-apply-linear-regression-in-r/" title="Permanent Link toHow to apply Linear Regression in R" target="_blank">How to apply Linear Regression in R</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/linear-regression-in-python-predict-the-bay-areas-home-prices/" title="Permanent Link toLinear Regression in Python; Predict The Bay Area’s Home Prices" target="_blank">Linear Regression in Python; Predict The Bay Area’s Home Prices</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/building-a-logistic-regression-in-python-step-by-step/" title="Permanent Link toBuilding A Logistic Regression in Python, Step by Step" target="_blank">Building A Logistic Regression in Python, Step by Step</a></li>
<li><a rel="nofollow" href="https://datascienceplus.com/multicollinearity-in-r/" title="Permanent Link toMulticollinearity in R" target="_blank">Multicollinearity in R</a></li>
</ol>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://datascienceplus.com/machine-learning-with-r-caret-part-1/"> R Programming – DataScience+</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/8HCMaDSKcg0" height="1"></a></div>