Data Science in the Cloud with Microsoft Azure Machine Learning and R

Mixing native modules and R in Azure ML

Azure ML provides a wide range of modules for data I/O, data transformation, predictive modeling, and model evaluation. Most native Azure ML modules are computationally efficient and scalable.

The deep and powerful R language and its packages can be used to meet the requirements of specific data science problems. For example, solution-specific data transformation and cleaning can be coded in R. R language scripts contained in Execute R Script modules can be run in-line with native Azure ML modules. Additionally, the R language gives Azure ML powerful data visualization capabilities. In other cases, data science problems that require specific models available in R can be integrated with Azure ML.

As we work through the examples in subsequent sections, you will see how to mix native Azure ML modules with Execute R Script modules.

Module I/O

In the AzureML Studio, input ports are located above module icons, and output ports are located below module icons.

For example, the Execute R Script module has five ports:

• The Dataset1 and Dataset2 ports are inputs for rectangular Azure data tables.

• The Script Bundle port accepts a zipped R script file (.R file) or R dataset file.

• The Result Dataset output port produces an Azure rectangular data table from a data frame.

• The R Device port produces output of text or graphics from R.

Workflows are created by connecting the appropriate ports between modules—output port to input port. Connections are made by dragging your mouse from the output port of one module to the input port of another module.

In Figure 1, you can see that the output of the data is connected to the Dataset1 input port of the Execute R Script module.

Model training workflow

Figure 2 shows a generalized workflow for training, scoring, and evaluating a model in Azure ML. This general workflow is the same for most regression and classification algorithms.

Key points on the model training workflow:

• Data input can come from a variety of data interfaces, including HTTP connections, SQLAzure, and Hive Query.

• For training and testing models, you will use a saved dataset.

• Transformations of the data can be performed using a combination of native Azure ML modules and the R language.

• A Model Definition module defines the model type and properties. On the lefthand pane of the Studio you will see numerous choices for models. The parameters of the model are set in the properties pane.

• The Training module trains the model. Training of the model is scored in the Score module and performance summary statistics are computed in the Evaluate module.

The following sections include specific examples of each of the steps illustrated in Figure 2.

Workflow for R model training

The Azure ML workflow changes slightly if you are using an R model. The generalized workflow for this case is shown in Figure 3.

In the R model workflow shown in Figure 3, the computation and prediction steps are in separate Execute R Script modules. The R model object is serialized, passed to the Prediction module, and unserialized. The model object is used to make predictions, and the Evaluate module measures the performance of the model.

Two advantages of separating the model computation step from the prediction step are:

• Predictions can be made rapidly on any number of new data, without recomputing the model.

• The Prediction module can be published as a web service.

Publishing a model as a web service

Once you have developed a satisfactory model you can publish it as a web service. You will need to create streamlined workflow for promotion to production. A generalized example is shown in Figure 4.

Key points on the workflow for publishing a web service:

• Data transformations are typically the same as those used to create the trained model.

• The product of the training processes (discussed above) is the trained model.

• You can apply transformations to results produced by the model. Examples of transformations include deleting unneeded columns, and converting units of numerical results.

A first set of transformations

For our first step, we’ll perform some transformations on the raw input data using the code shown below in an Azure ML Execute R Script module:

## This file contains the code for the transformation 
## of the raw bike rental data. It is intended to run in an 
## Azure ML Execute R Script module. By changing
## some comments you can test the code in RStudio 
## reading data from a .csv file. 

## The next lines are used for testing in RStudio only.
## These lines should be commented out and the following 
## line should be uncommented when running in Azure ML.
#BikeShare <- read.csv("BikeSharing.csv", sep = ",", 
#                      header = T, stringsAsFactors = F )
#BikeShare$dteday <- as.POSIXct(strptime(
#                         paste(BikeShare$dteday, " ", 
#                               "00:00:00", 
#                               sep = ""), 
#                         "%Y-%m-%d %H:%M:%S"))
BikeShare <- maml.mapInputPort(1)

## Select the columns we need
BikeShare <- BikeShare[, c(2, 5, 6, 7, 9, 10, 
                           11, 13, 14, 15, 16, 17)] 

## Normalize the numeric perdictors
BikeShare[, 6:9] <- scale(BikeShare[, 6:9])  

## Take the log of response variables. First we 
## must ensure there are no zero values. The difference 
## between 0 and 1 is inconsequential. 
BikeShare[, 10:12] <- lapply(BikeShare[, 10:12], 
                             function(x){ifelse(x == 0, 
                                                1,x)})
BikeShare[, 10:12] <- lapply(BikeShare[, 10:12], 
                             function(x){log(x)})

## Create a new variable to indicate workday
BikeShare$isWorking <- ifelse(BikeShare$workingday & 
                                !BikeShare$holiday, 1, 0)  ## Create a new variable to indicate workday

## Add a column of the count of months which could 
## help model trend. Next line is only needed running
## in Azure ML
Dteday <- strftime(BikeShare$dteday, 
                   format = "%Y-%m-%dT%H:%M:%S")
yearCount <- as.numeric(unlist(lapply(strsplit(
                                    Dteday, "-"), 
                                  function(x){x[1]}))) - 2011 
BikeShare$monthCount <- 12 * yearCount + BikeShare$mnth

## Create an ordered factor for the day of the week 
## starting with Monday. Note this factor is then 
## converted to an "ordered" numerical value to be
## compatible with Azure ML table data types.
BikeShare$dayWeek <- as.factor(weekdays(BikeShare$dteday))
BikeShare$dayWeek <- as.numeric(ordered(BikeShare$dayWeek, 
                                        levels = c("Monday", 
                                                   "Tuesday",
                                                   "Wednesday", 
                                                   "Thursday", 
                                                   "Friday", 
                                                   "Saturday", 
                                                   "Sunday")))

## Output the transformed data frame.
maml.mapOutputPort('BikeShare')

In this case, five basic types of transformations are being performed:

• A filter, to remove columns we will not be using.

• Transforming the values in some columns. The numeric predictor variables are being centered and scaled and we are taking the log of the response variables. Taking a log of a response variable is commonly done to transform variables with non-negative values to a more symmetric distribution.

• Creating a column indicating whether it’s a workday or not.

• Counting the months from the start of the series. This variable is used to model trend.

• Creating a variable indicating the day of the week.

Exploring the data

Let’s have a first look at the data by walking through a series of exploratory plots.

At this point, our Azure ML experiment looks like Figure 5. The first Execute R Script module, titled “Transform Data,” contains the code shown here.

The Execute R Script module shown at the bottom of Figure 5 runs code for exploring the data, using output from the Execute R Script module that transforms the data.

Our first step is to read the transformed data and create a correlation matrix using the following code:

## This code will create a series of data visualizations
## to explore the bike rental dataset. This code is 
## intended to run in an Azure ML Execute R 
## Script module. By changing some comments you can 
## test the code in RStudio.

## Source the zipped utility file
source("src/utilities.R")

## Read in the dataset. 
BikeShare <- maml.mapInputPort(1)

## Extract the date in character format 
BikeShare$dteday <- get.date(BikeShare$dteday)

## Look at the correlation between the predictors and 
## between predictors and quality. Use a linear 
## time series regression to detrend the demand.
Time <- POSIX.date(BikeShare$dteday, BikeShare$hr)
BikeShare$count <- BikeShare$cnt - fitted(
  lm(BikeShare$cnt ~ Time, data = BikeShare))
cor.BikeShare.all <- cor(BikeShare[, c("mnth", 
                                       "hr", 
                                       "weathersit", 
                                       "temp",
                                       "hum", 
                                       "windspeed",
                                       "isWorking", 
                                       "monthCount", 
                                       "dayWeek", 
                                       "count")])

diag(cor.BikeShare.all) <- 0.0 
cor.BikeShare.all
library(lattice)
plot( levelplot(cor.BikeShare.all, 
        main ="Correlation matrix for all bike users",
        scales=list(x=list(rot=90), cex=1.0)) )

We’ll use lm() to compute a linear model used for de-trending the response variable column in the data frame. De-trending removes a source of bias in the correlation estimates. We are particularly interested in the correlation of the predictor variables with this de-trended response.

This code requires a few functions, which are defined in the utilities.R file. This file is zipped and used as an input to the Execute R Script module on the Script Bundle port. The zipped file is read with the familiar source() function.

fact.conv <- function(inVec){
  ## Function gives the day variable meaningful 
  ## level names.
  outVec <- as.factor(inVec)
  levels(outVec) <- c("Monday", "Tuesday", "Wednesday", 
                      "Thursday", "Friday", "Saturday", 
                      "Sunday")
  outVec
}


get.date <- function(Date){
  ## Funciton returns the data as a character 
  ## string from a POSIXct datatime object. 
  strftime(Date, format = "%Y-%m-%d %H:%M:%S")
}


POSIX.date <- function(Date,Hour){
  ## Function returns POSIXct time series object 
  ## from date and hour arguments.
  as.POSIXct(strptime(paste(Date, " ", as.character(Hour), 
                            ":00:00", sep = ""), 
                        "%Y-%m-%d %H:%M:%S"))
}

Using the cor() function, we’ll compute the correlation matrix. This correlation matrix is displayed using the levelplot() function in the lattice package.

A plot of the correlation matrix showing the relationship between the predictors, and the predictors and the response variable, can be seen in Figure 6. If you run this code in an Azure ML Execute R Script, you can see the plots at the R Device port.

This plot is dominated by the strong correlation between dayWeek and isWorking—this is hardly surprising. It’s clear that we don’t need to include both of these variables in any model, as they are proxies for each other.

To get a better look at the correlations between other variables, see the second plot, in Figure 7, without the dayWeek variable.

In this plot we can see that a few of the predictor variables exhibit fairly strong correlation with the response. The hour (hr), temp, and month (mnth) are positively correlated, whereas humidity (hum) and the overall weather (weathersit) are negatively correlated. The variable windspeed is nearly uncorrelated. For this plot, the correlation of a variable with itself has been set to 0.0. Note that the scale is asymmetric.

We can also see that several of the predictor variables are highly correlated—for example, hum and weathersit or hr and hum. These correlated variables could cause problems for some types of predictive models.

Next, time series plots for selected hours of the day are created, using the following code:

## Make time series plots for certain hours of the day
times <- c(7, 9, 12, 15, 18, 20, 22)
lapply(times, function(x){
       plot(Time[BikeShare$hr == x], 
            BikeShare$cnt[BikeShare$hr == x], 
            type = "l", xlab = "Date", 
            ylab = "Number of bikes used",
            main = paste("Bike demand at ",
                      as.character(x), ":00", spe ="")) } )

Two examples of the time series plots for two specific hours of the day are shown in Figures 8 and 9.

Notice the differences in the shape of these curves at the two different hours. Also, note the outliers at the low side of demand.

Next, we’ll create a number of box plots for some of the factor variables using the following code:

## Convert dayWeek back to an ordered factor so the plot is in
## time order.
BikeShare$dayWeek <- fact.conv(BikeShare$dayWeek)

## This code gives a first look at the predictor values vs the demand for bikes.
library(ggplot2)
labels <- list("Box plots of hourly bike demand",
            "Box plots of monthly bike demand",
            "Box plots of bike demand by weather factor",
            "Box plots of bike demand by workday vs. holiday",
            "Box plots of bike demand by day of the week")
xAxis <- list("hr", "mnth", "weathersit", 
              "isWorking", "dayWeek")
capture.output( Map(function(X, label){ 
      ggplot(BikeShare, aes_string(x = X, 
                                   y = "cnt", 
                                  group = X)) + 
      geom_boxplot( ) + ggtitle(label) +
                           theme(text = 
                                   element_text(size=18)) },
    xAxis, labels),
  file = "NUL" )

If you are not familiar with using Map() this code may look a bit intimidating. When faced with functional code like this, always read from the inside out. On the inside, you can see the ggplot2 package functions. This code is wrapped in an anonymous function with two arguments. Map() iterates over the two argument lists to produce the series of plots.

Three of the resulting box plots are shown in Figures 10, 11, and 12.

From these plots you can see a significant difference in the likely predictive power of these three variables. Significant and complex variation in hourly bike demand can be seen in Figure 10. In contrast, it looks doubtful that weathersit is going to be very helpful in predicting bike demand, despite the relatively high (negative) correlation value observed.

The result shown in Figure 12 is surprising—we expected bike demand to depend on the day of the week.

Once again, the outliers at the low end of bike demand can be seen in the box plots.

Finally, we’ll create some plots to explore the continuous variables, using the following code:

## Look at the relationship between predictors and bike demand
labels <- c("Bike demand vs temperature",
            "Bike demand vs humidity",
            "Bike demand vs windspeed",
            "Bike demand vs hr")
xAxis <- c("temp", "hum", "windspeed", "hr")
capture.output( Map(function(X, label){ 
      ggplot(BikeShare, aes_string(x = X, y = "cnt")) + 
      geom_point(aes_string(colour = "cnt"), alpha = 0.1) + 
      scale_colour_gradient(low = "green", high = "blue") + 
      geom_smooth(method = "loess") + 
      ggtitle(label) +
      theme(text = element_text(size=20)) },
    xAxis, labels), 
  file = "NUL" )

This code is quite similar to the code used for the box plots. We have included a “loess” smoothed line on each of these plots. Also, note that we have added a color scale so we can get a feel for the number of overlapping data points. Examples of the resulting scatter plots are shown in Figures 13 and 14.

Figure 13 shows a clear trend of generally decreasing bike demand with increased humidity. However, at the low end of humidity, the data are sparse and the trend is less certain. We will need to proceed with care.

Figure 14 shows the scatter plot of bike demand by hour. Note that the “loess” smoother does not fit parts of these data very well. This is a warning that we may have trouble modeling this complex behavior.

Once again, in both scatter plots we can see the prevalence of outliers at the low end of bike demand.

Exploring a potential interaction

Perhaps there is an interaction between time of day and day of the week. A day of week effect is not apparent from Figure 12, but we may need to look in more detail. This idea is easy to explore. Adding the following code to the visualization Execute R Script module creates box plots for working and non-working days for peak demand hours:

## Explore the interaction between time of day
## and working or non-working days.
labels <- list("Box plots of bike demand at 0900 for \n working and non-working days",
               "Box plots of bike demand at 1800 for \n working and non-working days")
Times <- list(8, 17)
capture.output( Map(function(time, label){ 
      ggplot(BikeShare[BikeShare$hr == time, ], 
         aes(x = isWorking, y = cnt, group = isWorking)) + 
      geom_boxplot( ) + ggtitle(label) +
      theme(text = element_text(size=18)) },
    Times, labels),
  file = "NUL" )

The result of running this code can be seen in Figures 15 and 16.

Now we can clearly see that we are missing something important in the initial set of features. There is clearly a different demand between working and non-working days at peak demand hours.

Creating a new variable

We need a new variable that differentiates the time of the day by working and non-working days; to do this, we will add the following code to the transform Execute R Script module:

## Add a variable with unique values for time of day for working and non-working days.
BikeShare$workTime <- ifelse(BikeShare$isWorking, 
                             BikeShare$hr, 
                             BikeShare$hr + 24) 

Transformed time: Another new variable

As noted earlier, the complex hour-to-hour variation bike demand shown in Figures 10 and 14 may be difficult for some models to deal with. Perhaps, if we shift the time axis we will create a new variable where demand is closer to a simple hump shape. The following code shifts the time axis by five hours:

## Shift the order of the hour variable so that it is smoothly
## "humped over 24 hours.
BikeShare$xformHr <- ifelse(BikeShare$hr > 4, 
                            BikeShare$hr - 5, 
                            BikeShare$hr + 19)

We can add one more plot type to the scatter plots we created in the visualization model, with the following code:

## Look at the relationship between predictors and bike demand
labels <- c("Bike demand vs temperature",
            "Bike demand vs humidity",
            "Bike demand vs windspeed",
            "Bike demand vs hr",
            "Bike demand vs xformHr")
xAxis <- c("temp", "hum", "windspeed", "hr", "xformHr")
capture.output( Map(function(X, label){ 
      ggplot(BikeShare, aes_string(x = X, y = "cnt")) + 
      geom_point(aes_string(colour = "cnt"), alpha = 0.1) + 
      scale_colour_gradient(low = "green", high = "blue") + 
      geom_smooth(method = "loess") + 
      ggtitle(label) +
      theme(text = element_text(size=20)) },
    xAxis, labels), 
  file = "NUL" )

The resulting plot is shown in Figure 17.

The bike demand by transformed hour is definitely more of a hump shape. However, there is still a bit of residual structure at the lower end of the horizontal axis. The question is, will this new variable improve the performance of any of the models?