Visualizing Statistics and Regressions from a Spreadsheet using R

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

These notebooks are typically this is designed to create a pleasing viewing environment of data analysis that allows you to include figures, text, links, etc. so that your work is better understood and can be reproduced and used with confidence.

The source code for this R notebook (Rmd suffixed files), when stored as web pages (html files), can be downloaded by clicking the button at the top of the page.

If viewing the source code in R Studio, try executing each R “chunk” by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

Warning. Typos are Legion!

1. Introduction

When you’re in MATH 381 (Intro to Probability and Stats) you’ll get a taste of R. R is an open-source statistical package build off of an earlier generation of commercial.

The goal here is to demonstrate cracking open an excel spreadsheet in R and calculate some basic stats, create various plots to view the statistics, and finally, do some linear and multivariate regression

Another goal here is to show off some of R’s features. R is a very powerful tool. When translating “powerful” from computereese to any frustrated human dialect, that means “steep learning curve.” It’s also a community-supported environment. When translating “powerful” from computereese to any over-scheduled human dialect, that means “there are LOTS of people donating packages and libraries to R.” Some have evolved to be a standard in the community. Others are highly specialized for a given discipline (but have one or two items that people outside their user communities find handy.)

But don’t let that intimidate you. Once you learn one language you can slowly pick up more. Also with this demo we aren’t going to get to to be an R guru in a day.

If you want a good stepping off point to learn R I’d recommend some of the resources at Data Camp which have some free starter tutorials for R.

2. Loading the Libraries

To work with R we will first have to load some libraries. This is like in C where you have the #include statement to do things like raise things to powers and stuff like that.

Some of these libraries or “packages” come with R. Others will have to be installed. Here are the ones we are using for this exercise.

Also in this exercise, we’re going to use the tidyverse set of packages. Tidyverse is a set of co-developed tools for data science in R. This is the new big thing in R and is widely used so we are just going to jump in here. SD Mines has a course beyond Engineering Stats, MATH 443/543 (Data Analysis) that leverages this set of packages.

• Install Us First
  • tidyverse : Set of commonly-used Data Science packages for R that it can install and load all at once. In the long-run you probably also want to install the tidyverse package suite anyway. For this exercise this will include…
    • ggplot2 : Create Elegant Data Visualizations Using the Grammar of Graphics
    • tibble : Simple Data Frames
    • tidyr : Tools for shepherding data in data frames.
    • readr : Read Rectangular Text Data
    • purr : Functional Programming Tools
    • dplyr : A grammar of data manipulation
    • stringr : Simple, Consistent Wrappers for Common String Operations
    • forcats : Tools for Working with Categorical Variables (Factors)
  • readxl : also part of the tidyverse package suite for reading traditional excel spreadsheets.
    
  • moderndive : Tidyverse-Friendly Introductory Linear Regression
• This should come with R’s core install, if not install ’em.
  • MASS : Has a lot of resources for regression.
• This doesn’t come with R’s core install so install that one…
  • moments : This has a load of good stuff for data analysis and plotting, more than you will need here, but get it anyway.
• This is a nice contributed library that lets us make pretty statistics tables. It was written for ecological applications but it’s still pretty handy for looking at concrete
  • pastecs: Package for Analysis of Space-Time Ecological Series
• Another nice contributed library that makes matrices of correlation coefficients look pretty (and graphically informative).
  • corrplot Visualization of a Correlation Matrix

  # Tidyverse Handling Libraries

  library(package = "tidyverse")   # main tidyverse suite
  library(package = "readxl")     # Read Excel Files
  library(package = "moderndive") # regression support

  # Statistics Libraries

  library(package = "moments")   # Moments, cumulants, skewness, kurtosis and related tests
  library(package = "MASS", warn.conflicts=FALSE)      # Support Functions and Datasets for Venables & Ripley's MASS text

  # Extra Graphics Libraries

  library(package = "corrplot")  # Visualization of a Correlation Matrix


  # Data Processing Libraries

  library(package = "pastecs")   # Package for Analysis of Space-Time Ecological Series

3. Cracking a Spreadsheet

The spreadsheet example below is a more complicated than what you hopefully have.

The original data set is from a set of papers on Concrete by I-Cheng Yeh

• Yeh, I-Cheng, “Modeling slump of concrete with fly ash and superplasticizer,” Computers and Concrete, 5(6), 559-572, 2008. doi: 10.12989/cac.2008.5.6.559.

• Yeh, I-Cheng, “Simulation of concrete slump using neural networks,” Construction Materials, 162(1), 11-18, 2009. doi: 10.1680/coma.2009.162.1.11

• Yeh, I-Cheng, “Prediction of workability of concrete using design of experiments for mixtures,” Computers and Concrete, 5(1), 1-20, 2008. doi: 10.12989/cac.2008.5.1.001

• Yeh, I-Cheng, “Modeling slump flow of concrete using second-order regressions and artificial neural networks,” Cement and Concrete Composites, 29(6), 474-480, 2007. doi: 10.1016/j.cemconcomp.2007.02.001

• Yeh, I-Cheng, “Exploring concrete slump model using artificial neural networks,” ASCE J. of Computing in Civil Engineering, 20(3), 217-221, 2006. doi: 10.1061/(ASCE)0887-3801(2006)20:3(217)

and is kept at the UC-Irvine Machine Learning Repository.

It can be found here at http://kyrill.ias.sdsmt.edu/cee_284/Base_Concrete_Slump_Test_for_R.xlsx

The relevant page and screenshot is below. For drama-free R import you are probably best off keeping a page on your spreadsheet file that is very simple, with numbers going down, and a single line for Row-1 with the headers of each column. If you want to get fancy on other pages that you’d turn in as tables in reports, you can do that on another spreadsheet page.

Concrete Spreadsheet Screenshot

To crack open the spreadsheet we will want to use the readxl::read_excel function.

You can read the spreadsheet from a local drive or from a website.

  # you will need the full path to the file you are using (either online or locally on your disk)

  # The if else block should query your machine to determine which operating system.
  #  if you are not bi-platform, you likely don't need this.

  if(.Platform$OS.type == "windows") {
    # Windows
    spreadsheet_name     = "%HOMEPATH%/Downloads/Base_Concrete_Slump_Test_for_R.xlsx"
  } else {
    # Unix (Linux, MacOS, Solaris)
    spreadsheet_name     = "~/Downloads/Base_Concrete_Slump_Test_for_R.xlsx"
  }


  # I am keeping a copy of these spreadsheet at the URL below.  It can be downloaded automatically
  #   and then loaded.  We can also discretely delete it when done.

      spreadsheet_url = "http://kyrill.ias.sdsmt.edu/wjc/eduresources/Base_Concrete_Slump_Test_for_R.xlsx"
   
      download.file(url      =   spreadsheet_url, # URL location
                    destfile = spreadsheet_name) # local downloaded location

trying URL 'http://kyrill.ias.sdsmt.edu/wjc/eduresources/Base_Concrete_Slump_Test_for_R.xlsx'
Content type 'application/vnd.openxmlformats-officedocument.spreadsheetml.sheet' length 22230 bytes (21 KB)
==================================================
downloaded 21 KB

      
      remove(spreadsheet_url) # clean up variables
  
  # this command will read the file

  concrete = read_excel(path      = spreadsheet_name,  # remove spreadsheet location
                        sheet     = "Data",            # page of spreadsheet
                        col_names = TRUE)              # first row are the column headers
  
  
  # clean up your hard drive!  Don't be like me!

  

concrete

With the data read in we can now look at the table of the data. This looks much nicer when working in R Notebooks instead of Plain Ordinary R.

  # Print data frame
  colnames(concrete)[1] = "Test_Number"
  print(concrete)

NA

Extra: Units: AVOID!

… Uh… Just don’t. Unlike Mathcad which is designed to work with units or python which claims to be able to work with units, R is awful with it. I’ve used units once or twice. But I’ve found it to be more work than it’s worth.

4. Some Basic Statistics and Traditional Single Variable Plots

Lets start with some basic statistics and plotting of them.

4.1. The “classic” stats

Let’s get the mom-and-apple-pie stats for Concrete That second argument allows you to deal with missing data.

  # statistics for cement


  print(str_c("    Mean Cement : ",
              mean(x     = concrete$Cement, # variable to crunch
                   na.rm =            TRUE) # ignore msissing data
              ))

[1] "    Mean Cement : 229.894174757282"

  print(str_c("   Stdev Cement : ",
              sd(x     = concrete$Cement, # variable to crunch
                 na.rm =            TRUE) # ignore msissing data
              ))

[1] "   Stdev Cement : 78.8772300268858"

  
  print(str_c("Skewness Cement : ",
              skewness(x     = concrete$Cement, # variable to crunch
                       na.rm =            TRUE) # ignore msissing data
              ))

[1] "Skewness Cement : 0.143018080025135"

  
  print(str_c("Kurtosis Cement : ",
              kurtosis(x     = concrete$Cement, # variable to crunch
                       na.rm =            TRUE) # ignore msissing data
              ))

[1] "Kurtosis Cement : 1.33448397363582"

OK this is a little clunky. It would be nice if someone somewhere made a support library for R that will make nice tables of statistics.

In this case Vive La France! A team from French Research Institute for Exploitation of the Sea thought the same question and as is often the case for the R community not only drafted a set of tools to do this, and made it public.

Here we ware using their stat.desc function.

This will hopefully give people wanting to make basic tables maximum satisfaction with minimal effort.

  # Plot a statistics table -- all the classics nice and handy and pretty.

  options(digits=2) # this simply set the decimal count in the table to be created below  
                    # this particular function creates the table in scientific notation
  
  concrete_statistics = stat.desc(x    = concrete,  # data frame
                                  basic =    TRUE,  # includes counts and extremes 
                                  desc =     TRUE,  # include classic stats (mean etc)
                                  norm =     TRUE,  # include normal dist stats (skewness etc)
                                  p    =     0.95)  # use 95% confidence limits


  print(concrete_statistics)

NA

4.2. Reorganizing Your Data to Handle Multiple Variables at Once

To leverage some of R’s more nifty features we will need to reorganize our data from a “spreadsheet style” format to what some people have called a “long form” table so that the column headers of our concrete traits become a single column with the values in the columns placed all into a single column similar to the graphics below.

Example of the “pivot_long” Function

Example of the “pivot_wide” Function

In R this is now done with an analog to a function in Python’s Pandas’s package (and excel), the “pivot”.

In R this is managed in the tidr package using tidyr::pivot_longer() to gathers a group of columns into a single longer column, and tidyr::pivot_wider() which spreads a column into a shorter-but-wider table.

The documentation doesn’t help me a lot here. I need cartoons for that!

  # Gathering our components into a single column.

  # the pivot_longer command will group everything in the column name group 
  
  concrete_tidy = concrete %>%    # pipe your data frame from "concrete"
    pivot_longer(cols      = "Cement":"Compressive_Strength_28dy", # the list for the columns to "gather"
                 names_to  = "Parameter",                          # column name for your former columns
                 values_to =     "Value")                          # column name for your data
  

  # this will let us sort future plots in the same order as our plots.  
  
  concrete_tidy$Parameter = factor(x      = concrete_tidy$Parameter,
                                   levels = unique(concrete_tidy$Parameter))
  
  # we can also split things between our dependant variables and independant variables.
  
  
  concrete_independent = subset(x      = concrete_tidy,
                                subset = (Parameter != "Slump") &
                                         (Parameter != "Flow")  &
                                         (Parameter != "Compressive_Strength_28dy")
                                ) 
    
    
  concrete_dependent = subset(x      = concrete_tidy,
                              subset = (Parameter == "Slump") |
                                       (Parameter == "Flow")  |
                                       (Parameter == "Compressive_Strength_28dy")
                              )

 
                       

  print(concrete_tidy)

  print(concrete_independent)

  print(concrete_dependent)

NA
NA

5. Plotting Graphics using Tidyverse Resources

R has a few ways to do the basic histograms, Boxplots and other distribution plots.

There are a number of spiffy ways to plot these statistical plots in R. We’re just using one here…

5.1. SLOOOOWWWWLLLLLYYY Making a Simple Plot (Histogram Edition)

Now I’m going to do this one tiny step at a time until we get to a viable product. (This is how I work through cryptic procedures so I can see what each little additional mystery thingie does.)

Graphing is invoked by the ggplot2 command.. which has a heluvalot under its hood! For me all that detail was what had me a little shy to adopt this way of printing data.

Tidyverse uses what is sometimes called the “grammar of graphics” method… to make a long story longer, the GoG presents separate commands to do separate things rather bundle stuff in a single graphing function. Sometimes it makes a lot of sense… other times it may be confusion. (Hence me demonstrating making a graph this one tiny step at a time!

First thing we are going to do is open a plotting space with the command ggplot2::ggplot()

# invoke the ggplot plotting environmnent.

ggplot() 

Wow. We have a… big square of… grey. All it’s doing is setting up our plot environment… so let’s do some more…

If we want to do a histogram we are going to have to tell it what we want to print and where to get the stuff

When we add things to a plot command in Tidyverse we “add” to the steps incrementally.

This involves a “mapping” function called “ggplot2::aes” (short for aesthetics)

here, we are working with the data frame “concrete” and are working on the variable Cement which we are tossing onto the x axis because that’s where the bins of cement go!

ggplot(data = concrete) +   # EDIT:  invoke graphics environment using a given dataframe
  
  aes(x    = Cement)        # NEW: select variable to print... You can get really fancy here later

OK now we have something that looks like we may have the making of the graph. If you don’t like grey outlines and white grids, no worries, we can change that shortly.

OK.. we are now ready to make a histogram…

Here we will use one of the gglot2’s “geom_*” (draw stuff) resources. The default should work for us here.

ggplot(data = concrete) +   # invoke graphics environment using a given dataframe
  
  aes(x = Cement)   +       # select variable to print... You can get really fancy here later

  geom_histogram()          # NEW: insert histogram

(you may have gotten a warning about using the bin=X, you can adjust it.)

Now quickly before moving on… I am not keen on the grey background with white lines.

There are a number of out-of-the-box “themes” for ggplot2.

I’m partial to theme_bw() and theme_light() but try the ones that you prefer or stick with the default, theme_gray().

These plots shown here are mine. You should fidget about so they are yours and so you can adapt to this new way of working with data.

ggplot(data = concrete) + # invoke graphics environment using a given dataframe
  
  theme_bw() +            # NEW: changing the plotting theme
  
  aes(x = Cement) +       # select variable to print... You can get really fancy here later

  geom_histogram()        # insert histogram (including controlling number of bins)

My OCD hates axes where the labels don’t envelop all of the data…

We can fix that with ggplot2::xlim() or ggplot2::ylim()

ggplot(data = concrete) +     # invoke graphics environment using a given dataframe
  
  theme_bw() +                # changing the plotting theme
  
  aes(x = Cement) +           # select variable to print... You can get really fancy here later
  
  xlim( 100, 400 ) +          # NEW: adding x-axis limits

  geom_histogram()            # insert histogram

How about changing the color of the fill in the bars…

You really don’t want to know about all the colors you can use.

ggplot(data = concrete) +     # invoke graphics environment using a given dataframe
  
  theme_bw() +                # changing the plotting theme
  
  aes(x = Cement) +           # select variable to print... You can get really fancy here later
  
  xlim( 100, 400 ) +          # NEW: adding x-axis limits

  geom_histogram(fill="gray") # EDIT: insert histogram (with a single chosen color)

Want to customize the labels and titles so we can have units?

You can add custom labels and titles! (https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/colorPaletteCheatsheet.pdf)

For the superscripting in the x-axis label, I am using the expression() tool in R.

ggplot(data = concrete) +     # invoke graphics environment using a given dataframe
  
  theme_bw() +                # changing the plotting theme
  
  aes(x = Cement) +           # select variable to print... You can get really fancy here later
  
  xlim( 100, 400 ) +          # adding x-axis limits

  ggtitle("Yeh Superplasticizer Tests") +          # NEW : Custom Title
  
  xlab(expression('Cement Amount (kg m'^-3*")")) + # NEW : Custom Axis Label

  geom_histogram(fill="gray") # insert histogram (with a single chosen color)

And I could keep tweaking this graph all day, but good enough is good enough so this is a good place to stop…

We also can plot a few other fields with some trial and error..

# Histogram of Water

ggplot(data = concrete) +     # invoke graphics environment using a given dataframe
  
  theme_bw() +                # changing the plotting theme
  
  aes(x = Water) +           # select variable to print... You can get really fancy here later
  
  xlim( 150, 250 ) +          # adding x-axis limits

  ggtitle("Yeh Superplasticizer Tests") + #Custom Title
  
  xlab(expression('Water Amount (kg m'^-3*")")) + # NEW : Custom Axis Label note use of superscripts from above

  geom_histogram(fill="blue") # insert histogram (with a single chosen color)

# Histogram of Strength

ggplot(data = concrete) +     # invoke graphics environment using a given dataframe
  
  theme_bw() +                # changing the plotting theme
  
  aes(x = Compressive_Strength_28dy) + # select variable to print... You can get really fancy here later
  
  xlim( 10, 60 ) +          # adding x-axis limits

  ggtitle("Yeh Superplasticizer Tests") + #Custom Title
  
  xlab("28-dy Compressive Strength (MPa)") + # NEW : Custom Axis Label

  geom_histogram(fill="red") # insert histogram (with a single chosen color)

(And from our Intro to Stats Lecture…)

# Histogram of Strength

ggplot(data = concrete) +     # invoke graphics environment using a given dataframe
  
  theme_bw() +                # changing the plotting theme
  
  aes(x = Slump) + # select variable to print... You can get really fancy here later
  
  xlim( 0, 30 ) +          # adding x-axis limits

  ggtitle("Yeh Superplasticizer Tests") + #Custom Title
  
  xlab("Slump (cm)") + # NEW : Custom Axis Label

  geom_histogram(fill="darkgreen") # insert histogram (with a single chosen color)

5.2 Distribution Plot [not so good an] Example

There are some other plots that we can use to describe our data.

Here to play with them we will take a quick step back and address that “tidy”’ed (should that say “tidied”?) dataframe “concrete_tidy”

We can now use all the parameters in the “tidy” (long) data frame to print by specific traits.

ggplot(data = concrete_tidy) +            # invoke graphics environment using a given dataframe
  
  theme_bw() +                            # changing the plotting theme
  
  aes(x      = Value,                     # map x-axis value
      color  = Parameter) +               # map colors for different quality
  
  ggtitle("Yeh Superplasticizer Tests") + # Custom Title
  
  xlab("Value") +                         #  Custom Axis Label

  geom_density()                          # insert crete a relative density plot 

In the past, I’ve gotten good results with this but in this case, I think it’s too messy in part due to the disparity in the dynamic range of our parameters.

5.3. Box-Whisker Plot Example

How about leveraging a box whisker? (I’m using only the independent variables this time.)

ggplot(data = concrete_independent) +      # EDIT Changing dataframe
  
  theme_bw( ) +                            # changing the plotting theme
  
  theme(axis.text.x = element_blank()) +   # adding an extra trait to the x-axis
                                           # to not print labels on the x-axis 
                                           # (the labels overlap and doesn't look
                                           # pretty...)
  
  aes(y      = Value,                     # map y-axis value
      x      = Parameter,                 # map x-axis value
      color  = Parameter) +               # map colors for different quality
  
  ggtitle(label    = "Yeh Superplasticizer Tests",
          subtitle = "Concrete Test Components") + # Custom Title
  
  ylab(expression('Amount (kg m'^-3*")")) + # EDIT : Changing Custom Axis Label

  geom_boxplot()                          # insert crete a relative density plot 

What about our dependent variables? We can start by changing the data frame…

ggplot(data = concrete_dependent) +      # EDIT Changing dataframe
  
  theme_bw( ) +                            # changing the plotting theme
  
  theme(axis.text.x = element_blank()) +   # adding an extra trait to the x-axis
                                           # to not print labels on the x-axis 
                                           # (the labels overlap and doesn't look
                                           # pretty...)
  
  aes(y      = Value,                     # map y-axis value
      x      = Parameter,                 # map x-axis value
      color  = Parameter) +               # map colors for different quality
  
  ggtitle(label    = "Yeh Superplasticizer Tests",
          subtitle = "Concrete Test Results") + # Custom Title
  
  ylab("Values") +

  geom_boxplot()                          # insert crete a relative density plot 

Want units? That’s a little tougher here since the units differ by parameter. We can force the values to into new names though.

ggplot(data = concrete_dependent) +      # EDIT Changing dataframe
  
  theme_bw( ) +                            # changing the plotting theme
  
  theme(axis.text.x = element_blank()) +   # adding an extra trait to the x-axis
                                           # to not print labels on the x-axis 
                                           # (the labels overlap and doesn't look
                                           # pretty...)
  
  aes(y      = Value,                     # map y-axis value
      x      = Parameter,                 # map x-axis value
      color  = Parameter) +               # map colors for different quality
  
  ggtitle(label    = "Yeh Superplasticizer Tests",
          subtitle = "Concrete Test Results") + # Custom Title
  
  ylab("Values") +

  # NEW: It says scale color but "color" is how we are distinguishing
  #      out boxplots (as seen in the mapping/aes command)
  #      we can then use the same plot order above to rewrite the labels
  #      (likewise we could change the plot order and of coruse the colors.)
  scale_color_discrete(labels = c("Slump (cm)",
                                  "Flow (cm)", 
                                  "28dy-Compresional Stress (mPa)")) + 
  
  geom_boxplot() # insert crete a relative density plot 

NA
NA

5.4. Violin Plot Example

How about leveraging a “violin” plot? A violin plot’s width swells in areas with more observations and contracts with sparser data so it is like looking at a probability distribution.

ggplot(data = concrete_independent) +      # EDIT Changing dataframe
  
  theme_bw( ) +                            # changing the plotting theme
  
  theme(axis.text.x = element_blank()) +   # adding an extra trait to the x-axis
                                           # to not print labels on the x-axis 
                                           # (the labels overlap and doesn't look
                                           # pretty...)
  
  aes(y      = Value,                     # map y-axis value
      x      = Parameter,                 # map x-axis value
      color  = Parameter) +               # map colors for different quality
  
  ggtitle(label    = "Yeh Superplasticizer Tests",
          subtitle = "Concrete Test Components") + # Custom Title
  
  ylab(expression('Amount (kg m'^-3*")")) + #  Changing Custom Axis Label

  geom_violin(scale="width") # EDIT: change to a violin plot 

                             #   the width argument 
                             # gives every plot the same width
  

and…

ggplot(data = concrete_dependent) +      # EDIT Changing dataframe
  
  theme_bw( ) +                            # changing the plotting theme
  
  theme(axis.text.x = element_blank()) +   # adding an extra trait to the x-axis
                                           # to not print labels on the x-axis 
                                           # (the labels overlap and doesn't look
                                           # pretty...)
  
  aes(y      = Value,                     # map y-axis value
      x      = Parameter,                 # map x-axis value
      color  = Parameter) +               # map colors for different quality
  
  ggtitle(label    = "Yeh Superplasticizer Tests",
          subtitle = "Concrete Test Results") + # Custom Title
  
  ylab("Values") +

  # NEW: It says scale color but "color" is how we are distinguishing
  #      out boxplots (as seen in the mapping/aes command)
  #      we can then use the same plot order above to rewrite the labels
  #      (likewise we could change the plot order and of coruse the colors.)
  scale_color_discrete(labels = c("Slump (cm)",
                                  "Flow (cm)", 
                                  "28dy-Compresional Stress (mPa)")) + 
  

  geom_violin(scale="width") # EDIT: change to a violin plot 

                             #   the width argument 
                             # gives every plot the same width  

This is basically the above “density” plot but “looking down” as with a box plot. Also here we are trimming the plot so that when we leave the range of any of the data points, the “violins” are truncated.

5.5. Stacked Column or Bar Plot Example

We also can do bar plots or stacked column plots. The one produced here shows the combined components by test unit.

ggplot(data = concrete_independent) +      # EDIT Changing dataframe
  
  theme_bw( ) +                            # changing the plotting theme

  
  aes(x     = Test_Number,
      y     = Value,
      fill  = Parameter) +               # map colors for different quality
  
  ggtitle(label    = "Yeh Superplasticizer Tests",
          subtitle = "Concrete Test Components") + # Custom Title
  
  ylab(expression('Amount (kg m'^-3*")")) + #  Changing Custom Axis Label

  geom_col(position = "stack",  # new, create a stacekd column graph 
           width    = 1.0    )  # with no space between columns

6. Correlation of Variables

6.1. Correlating and then Fitting Cement to Compressive Strength

Let’s start by doing a “simple”” plot . In this case since I already know the answer because the spreadsheet also has a table of how well our independent variables correlate against the dependent variables (e.g., Slump, Flow, or in our case Strength). The Cement correlates the best against Compressive Strength (OK, truth be told, it correlates the least badly).

We can actually do this with a correlate function, boot::cor()…

To grab a value in the table “concrete” we call the data frame (concrete) and the variable name (Cement or Water vs Compressive_Strength_28dy), separating the frame and variable names by a $ sign.

print("Cement vs Compressive Strength Correlation, r")

[1] "Cement vs Compressive Strength Correlation, r"

cor(x = concrete$Cement,                    # the x-value 
    y = concrete$Compressive_Strength_28dy, # the y-value
    method = "pearson"                      # method of correlation
    )

[1] 0.45

or if you like to do everything at once…

# calculate all correlation values against each other

correlation_matrix = cor(x      = concrete, # using our dataframe to correlate evything
                         method = "pearson" )

tbl_df(correlation_matrix)

Warning: `tbl_df()` was deprecated in dplyr 1.0.0.
Please use `tibble::as_tibble()` instead.

Lots of numbers… not all that insightful on their own…

You also can graph the look-n-feel of what all of the different correlations are… (it works best with a much smaller number of variables). This function is called corrplot::corrplot() and makes a nice presentation of the data.

  # draw a correlation graphic...

  corrplot(corr   = correlation_matrix,
           type   = "upper")

We can now see for example that cement, slag, and fly ash amounts have a nominal but not thrilling correlation to compression strength while water has a good correlation with the resulting slump values. One thing that this does not show is how well these parameters play with other parameters. As we’ll see when all of our independent values are working together we’ll discover that cement and water, followed by fly ash and coarse aggregates will, together, contribute the most of our independent parameters in calculating the compressive strength.

6.2. Scatter Plot Example

But for now, let’s plot plot the Cement amount against Compressive Strength

# Making a simple X-Y scatter plot.

ggplot(data = concrete) +                # invoke graphics environment using a given dataframe
  
  theme_bw( ) +                           # changing the plotting theme
  
  aes(x      = Cement,                       # x-value
      y      = Compressive_Strength_28dy) +  # y-value

  ggtitle("Yeh Superplasticizer Tests") +    # Custom Title
  
  xlab(expression('Cement Amount (kg m'^3*")")) +   # x-label
  ylab("28-dy Compressive Strength (MPa)")      +   # y-label

  geom_point(colour="grey")   # EDIT: plot points the color keyword part was

                              #       writen by an anglophile!

Here’s a cute trick: Could we color those dots by a variable?

Sure!

# Making a simple X-Y scatterplot now coloured by another parameter

ggplot(data = concrete) +                # invoke graphics environment using a given dataframe
  
  theme_bw( ) +                           # changing the plotting theme
  
  aes(x      = Cement,                       # x-value
      y      = Compressive_Strength_28dy,    # y-value
      color  = Superplasticizer)          +  # ADD: we can color by a variable!

  ggtitle("Yeh Superplasticizer Tests") +    # Custom Title
  
  xlab(expression('Cement Amount (kg m'^3*")")) +   # x-label
  ylab("28-dy Compressive Strength (MPa)")      +   # y-label

  geom_point() +  # plot points 
  scale_color_distiller(palette = "Spectral") # NEW: pick a custom "colour" palate.

Love overkill without any distinct numerical score and look at how everything in your data set correlates with every other variables…?

Try graphics::pairs()

(I like the corrplot function better!)

# way too many tiny plots!

pairs(x   = concrete, # do everything in the dataframe
      pch = ".")      # plot dots (the default is circles)

(Obviously the more variables in your dataframe the messier it gets!)

6.3. Creating our linear model and “calibrating” it

We weren’t all that thrilled with the correlation between these components and strength but let’s go ahead and demonstrate a regression.

But let’s move on and create a regression model from this.

Here we will use the stats::lm() (linear model) from the basic toolboxes that come with R.

For the regression formula

\(\widehat{y}(x) = {\alpha_0}+{\alpha_1}\ x\)

or

\(\widehat{Strength}(concrete) = {\alpha_0}+{\alpha_1}\ concrete\)

the “prototype” (formula) for the function is written as …

“Y ~ X” (with the y-intercept implicit in the formula… you don’t put it in but it’ll be there when you’re done.)

The above syntax is works like this….

Dependent Variable [~ is a function of ] Independent Variable [and any other parameter you need gets added with a plus]

If this were a \(\widehat{y}(x)={\alpha_0}+{\alpha_0}\ x^3\), then the prototype for the function would be y ~ x^3

This will hopefully make more sense as we continue!

(lm and similar linear regression functions don’t play well with units.)

linear_model.S_v_c =  lm(formula = Compressive_Strength_28dy ~ Cement, # your formula y ~ x
                         data    = concrete)                           # the data frame

Let’s see what we have… This summary command will provide the details of the lm() function’s important results

For us we want to see the Y-Intercept [the (Intercept) under “Estimate”] and the slope that goes with our independent value (“Concrete” under “Estimate”)

The Standard Error of the Estimate is there (Residual Standard Error) as is the Coefficient of Determination (Multiple R-squared)

We’ll talk about a few of the other features when we do the larger multivariate regression

 summary(object = linear_model.S_v_c)

Call:
lm(formula = Compressive_Strength_28dy ~ Cement, data = concrete)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.134  -5.313   0.832   5.155  17.968 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 25.85676    2.15022      12  < 2e-16 ***
Cement       0.04429    0.00885       5  2.4e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7 on 101 degrees of freedom
Multiple R-squared:  0.199, Adjusted R-squared:  0.191 
F-statistic:   25 on 1 and 101 DF,  p-value: 2.38e-06

In the above output, the asterisk identify the most significant independent variables. Here it’s trivial even though this is a terrible relationship between cement and strength. Later we will use all of our available independent variables and the use of these asterisks will become more important.

Want to plot it?

Good news?

Like Excel, you have some automated features to give you quick satisfaction and happiness. More still, it will give you confidence limits.

For this we use an extension to the graphics package called ggplot2::geom_smooth()

# Making a simple X-Y scatterplot and adding a regression to it

ggplot(data = concrete) +                # invoke graphics environment using a given dataframe
  
  theme_bw( ) +                           # changing the plotting theme
  
  aes(x      = Cement,                       # x-value
      y      = Compressive_Strength_28dy) +  # y-value

  ggtitle("Yeh Superplasticizer Tests") +    # Custom Title
  
  xlab(expression('Cement Amount (kg m'^-3*")")) +   # x-label
  ylab("28-dy Compressive Strength (MPa)")      +   # y-label

  geom_point(colour="darkgrey") +  # plot points
  geom_smooth(method  = "lm",    # use a simple linar model
              formula = y ~ x,   # lm-style formula
              se      = TRUE,    # splay Confidence Intervals
              level   = 0.95,    # Confidene Level to Map Out
              colour  = "black", # regression line color
              size    = 0.5)     # line thickness

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

The line here looks like a positive correlation between the cement amount and the resulting strength.

Let’s try water:

# getting the linear model


linear_model.S_v_w =  lm(formula = Compressive_Strength_28dy ~ Water, # your formula y ~ x
                         data    = concrete   )                           # the data frame

summary(linear_model.S_v_w)

Call:
lm(formula = Compressive_Strength_28dy ~ Water, data = concrete)

Residuals:
    Min      1Q  Median      3Q     Max 
-19.359  -5.451  -0.986   4.690  18.825 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  55.4824     7.3978    7.50  2.5e-11 ***
Water        -0.0986     0.0373   -2.64   0.0096 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.6 on 101 degrees of freedom
Multiple R-squared:  0.0646,    Adjusted R-squared:  0.0554 
F-statistic: 6.98 on 1 and 101 DF,  p-value: 0.00956

# Making a simple X-Y scatterplot and adding a regression to it

ggplot(data = concrete) +                # invoke graphics environment using a given dataframe
  
  theme_bw( ) +                           # changing the plotting theme
  
  aes(x      = Water,                      # x-value
      y      = Compressive_Strength_28dy) +  # y-value

  ggtitle("Yeh Superplasticizer Tests") +    # Custom Title
  
  xlab(expression('Water Amount (kg m'^-3*")")) +  # x-label
  ylab("28-dy Compressive Strength (MPa)")      +   # y-label

  geom_point(colour="darkblue") +  # plot points
  
  geom_smooth(method  = "lm",    # use a simple linar model
              formula = y ~ x,   # lm-style formula
              se      = TRUE,    # splay Confidence Intervals
              level   = 0.95,    # Confidene Level to Map Out
              colour  = "blue",  # regression line color
              fill    = "cyan",  # NEW: fill for confidence limits
              size    = 0.5)     # line thickness

Looking up back the tables none of the variables

7. Multivariate Linear Regression

And now we’re going to do something about that!

We’re now going to use not just one independent variable… but all 7 of them!

The good news is that it follows the same form as the simple linear regression. This time we string along all of our independent variables with in our formula prototype.

Our formula now has multiple independent values but still follows the same style of solution…

\(\widehat{y}(\mathbf{x}) = {\alpha_0}+{\alpha_1} x_1 + {\alpha_2} x_2 + {\alpha_2} x_3 + ... +{\alpha_n} x_n\)

linear_model.S_v_all <- lm(data    = concrete,                             # your data frame
                           formula = Compressive_Strength_28dy ~ Cement +  # your formula
                                                                 Slag +
                                                                 Fly_Ash +
                                                                 Water +
                                                                 Superplasticizer +
                                                                 Fine_Aggregates +
                                                                 Coarse_Aggregates)  

And here are these results…

summary(object = linear_model.S_v_all)

Call:
lm(formula = Compressive_Strength_28dy ~ Cement + Slag + Fly_Ash + 
    Water + Superplasticizer + Fine_Aggregates + Coarse_Aggregates, 
    data = concrete)

Residuals:
   Min     1Q Median     3Q    Max 
-5.841 -1.706 -0.283  1.299  7.942 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)   
(Intercept)       139.7815    71.1013    1.97   0.0522 . 
Cement              0.0614     0.0228    2.69   0.0084 **
Slag               -0.0297     0.0318   -0.94   0.3520   
Fly_Ash             0.0505     0.0232    2.18   0.0316 * 
Water              -0.2327     0.0717   -3.25   0.0016 **
Superplasticizer    0.1031     0.1346    0.77   0.4453   
Fine_Aggregates    -0.0391     0.0288   -1.36   0.1783   
Coarse_Aggregates  -0.0556     0.0274   -2.03   0.0455 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.6 on 95 degrees of freedom
Multiple R-squared:  0.897, Adjusted R-squared:  0.889 
F-statistic:  118 on 7 and 95 DF,  p-value: <2e-16

Our regression coefficients are still here under the “Estimate” column as are our Standard Error of our Estimate and our Coeff of Determination.

Also we can now take a good look at those asterisks at the end of line with the parameter coefficients. These can explain which independent variables do the heaviest lifting in our regression. The more asterisks, the more important the dependent variable is to the larger multivariate regression. Here, we can see that the Cement and Water are doing most of the “work” in fitting our suite of independent variables to our dependent variable of Compressive Strength.

Finally there is the P parameter for which the smaller it is, the better we can say that the relationship that we’ve made with our regression represents our dependent variable.

Now… on to looking at our results.

Here is where viewing the results of the regression is tricky.

We have 7 independent variables but we’d like to see the impact of the fit if all 7 variables on our strength

When I do this I like to plot the true y value against my regression y(x1,x2,x3,..)

So to do this I will take the fitted values of y and plot them against the original values of y

Getting the fitted values is easy.

I’m using the get_regression_points function which adds the modeled “y-hat” value to the dataframe of all of the other values stats::get_regression_points() function.

The fitted version is the dependent variable w/ a “_hat”” at the end

fitted.S_v_all = get_regression_points(model = linear_model.S_v_all)

print(fitted.S_v_all)

NA

And finally we can plot our actual vs modeled values. (I’m adding a trend line)

# Making a simple X-Y scatterplot and adding a regression to it

ggplot(data = fitted.S_v_all) +           # invoke graphics environment using a given dataframe
  
  theme_bw( ) +                           # changing the plotting theme
  
  aes(x      = Compressive_Strength_28dy,    # x-value
      y      = Compressive_Strength_28dy_hat) +  # y-value

  ggtitle("Yeh Superplasticizer Tests",
          subtitle = "28-dy Compressive Strength (MPa)") +    # EDITED: Custom Title now with a subtitle
  
  ylab("Modelled")     + # y-label
  xlab("Observed")     + # x-label

  geom_point(colour="darkred") +  # plot points
  
  geom_smooth(method  = "lm",      # use a simple linar model
              formula = y ~ x,     # lm-style formula
              se      = TRUE,      # display Confidence Intervals
              level   = 0.95,      # Confidene Level to Map Out
              colour  = "red",     # regression line color
              fill    = "magenta", # fill for confidence limits
              size    = 0.5)  +    # line thickness
  
  geom_abline(slope     = 1,       # NEW: add a very simple line
              intercept = 0,       #  (for a 1:1 reference)
              color     = "grey",
              linetype  = "dashed") +

  coord_fixed(ratio = 1)           # NEW: make the aspect ratio 

                                   #   (I like my plots square)

And here we have a nice plot showing our true vs predicted values.

8. Regression Quality Metrics

And to close things off, we can do some general error metrics that may be useful..

First, the Mean Squared Error (MSE) or Bias… (if we are too high or too low)

\(BIAS = MSE = \frac{1}{N} \sum_{i=1}^{n} [\widehat{y}(\overrightarrow{x_i})-y_i] = \overline{[\widehat{y}(\overrightarrow{x_i})-y_i]}\)

  # Calculate Bias (MSE)

  bias = mean(fitted.S_v_all$Compressive_Strength_28dy_hat - 
                 fitted.S_v_all$Compressive_Strength_28dy)
  
  print(str_c(" Mean Squared Error (MSE) or Bias: ", bias))

[1] " Mean Squared Error (MSE) or Bias: 2.91262135922143e-05"

For a linear or multivariate regression the average of our residuals (the difference between each observation and prediction) should be zero.

The root mean squared error (RMSE) is shown here. It shouldn’t be zero since the residuals are squared before summing them up. We technically should use the standard error of the estimate, but RMSE remains a common error metric. We can always do both. The standard error of the estimate takes into account the degrees of freedom which which now includes all of the independent variables (p). We can get the standard error of the estimate from our

\(RMSE = \sqrt{ \frac{1}{N} \sum_{i=1}^{n} [\widehat{y}(\overrightarrow{x_i})-y_i]^2 } = \sqrt{\overline{[\widehat{y}(\overrightarrow{x_i})-y_i]^2} }\)

\(s_{e}\) or \(s_{y/x} = \sqrt{ \frac{1}{N-p-1} \sum_{i=1}^{n} [\widehat{y}(\overrightarrow{x_i})-y_i]^2 }\)

  # Calculate RMSE

  rmse = sqrt(mean( (fitted.S_v_all$Compressive_Strength_28dy_hat -
                       fitted.S_v_all$Compressive_Strength_28dy)^2)  )
  
  print(str_c("     Root Mean Squared Error (RMSE): ",  
              rmse))

[1] "     Root Mean Squared Error (RMSE): 2.50527978593714"

  print(str_c("Standard Error of the Estimate (se): ", 
              summary(linear_model.S_v_all)$sigma))  # you have to dig for this one!

[1] "Standard Error of the Estimate (se): 2.6086576339523"

And finally our correlation coefficient (which is basically our coefficient of determination before the “R” is “squared”)

  # Get The Unadjusted Correlation Coefficient

  r = cor(x = fitted.S_v_all$Compressive_Strength_28dy,     # the x-value 
          y = fitted.S_v_all$Compressive_Strength_28dy_hat, # the y-value
          method = "pearson"                                # method of correlation
          )
  
  print(str_c("                        correlation coefficient (r): ", r))

[1] "                        correlation coefficient (r): 0.94701611900088"

  print(str_c("                  coefficient of determination (r²): ", r^2, 
                                                                 " ", 
                                 summary(linear_model.S_v_all)$r.squared))

[1] "                  coefficient of determination (r²): 0.89683952964749 0.896837609814009"

  print(str_c("adjusted coefficient of determination (Adjusted r²): ", 
               summary(linear_model.S_v_all)$adj.r.squared))

[1] "adjusted coefficient of determination (Adjusted r²): 0.889236170537146"

9. Closing

And with that, we’re done… Once again, this exercise demonstrates a lot of tricks just to show how you can use R for various statistics. You may not use all of them in your encounters with R for linear or multivariate regression or even at all, but you may be able to cannibalize some of the tricks here for other applications.