The Goal

In our last class together, we started to explore splines. These are a special type of regression model that we use when the relationship between a numeric feature \(X\) and a numeric response variable \(Y\) is a curve which is not well matched to a polynomial.

Today, we are going to work on splines. Our goal is to increase our understanding of spline models, to practice building spline models in R, and to try assessing the predictive accuracy of splines models.

The Data

The data you need can be found on Canvas. When you load the data in, it should have 50598 rows and 79 columns. The data set contains information on different countries, and there are multiple rows of information for each country. These variables include things like population and GDP for a country for a given year. We are interested in \(X=\) the number of people living in the country and \(Y\) = the amount of carbon dioxide (CO2) the country produces.

We are going to focus on just the United States today, so we can subset our data down to only the rows relating to the United States.

dataUS <- subset(owid.co2.data, country == "United States")

When you run the code above, you should end up with a data set of \(n=222\) rows. This is the data set we will work with for today, so you can remove the original data set if you wish.

rm(owid.co2.data)

EDA

Now that we have our data, the next step is to visualize the relationship that we are interested in. Recall that we want to build a model for \(Y\) = the amount of carbon dioxide (CO2) in kilotons the country produces using \(X=\) the number of people living in the country as a feature.

table(cut_number(dataUS$population,2))

Building the Spline

With our knot chosen, we are ready to progress to building our spline. To create a cubic basis spline model, or just spline model for short, we need the splines package in R. You will need to install the package, and once you do you can run the library command.

library(splines)

To build a spline in R, the code is very similar to linear regression:

spline_1knot <- lm( co2 ~ bs( population, knots = ),data = dataUS)

We use bs to let R know we want a cubic basis spline. There is also an input called knots which allows us to specify where we want the knot(s) to go. You will notice that I have left knots = blank, so you will need to fill in the value of the knot you chose in Question 4.

Once we have run the code to create the spline_1knot object, our model is trained. Just like in linear regression, to view the actual coefficients for the trained model we have to use the summary command.

Now that we have trained the model, it can be helpful to visualize the spline that we have built to see how it matches the patterns in the data. To do this, we need to adapt our scatter plot from Question 5.

compute_RMSE <- function(truth, predictions){
  
  # Part 1
  part1 <- truth - predictions
  
  #Part 2
  part2 <- part1^2
  
  #Part 3: RSS 
  part3 <- sum(part2)
  
  #Part4: MSE
  part4 <- 1/length(predictions)*part3
  
  # Part 5:RMSE
  sqrt(part4)
}

# Option 1: Specify the knots

spline_2knots <- lm( co2 ~ bs( population, knots = ),data = dataUS)

# Option 2: Specify how many knots 

spline_2knots <- lm( co2 ~ bs( population, df =  ),data = dataUS)

Natural Splines

When you run the 10-fold CV code, we get warning messages. This message appears exactly twice: once for the smallest population value in the data set and once for largest population in the data set. In other words, we get the warnings with the the min and max values of the domain They define what are called boundary knots.

Okay, cool, but what do the warnings actually mean? What this warning message is telling us is that we will have trouble estimating \(\hat{y}_i\) at the boundary knots, i.e., min and max points in the domain. The reason is that there is no data beyond the boundary knots, and the spline models are so flexible that this can be a problem if we want to estimate beyond the boundary knots.

If we are interested in estimating at the boundary knots, or beyond the boundary knots, we generally make a very small adjustment to our spline. Instead of using a cubic basis spline, we use a natural cubic spline.

A natural cubic spline is an adaptation of the cubic spline that adjusts for the fact that cubic splines often have unstable estimates at the boundary knots. To correct this, natural splines require that the spline function be linear at the boundary, i.e., in the region smaller than the minimum \(x\) value and larger than the maximum \(x\) value. To do this, natural splines make the restriction that the second and third derivatives of the basis form of the spline model is 0 at the boundary knots.

We don’t really need to know this math, but we do need to know that we can use natural splines to stabilize estimation at the boundary knots.

To fit a natural spline in R, all we have to do is change the bs in our code to an ns.

# Option 1: Specify the knots

spline_2knots <- lm( co2 ~ ns( population, knots = ),data = dataUS)

# Option 2: Specify how many knots 

spline_2knots <- lm( co2 ~ ns( population, df =  ),data = dataUS)

Now, when we switch to a natural spline, we do unfortunately need to check again to see what the optimal number of knots is. Luckily, this is a very simple adaption to the code you have already done!