extrapolate-basis-coverage.R

library(mgcv)

#creating data for the smooth matrix and for predictions
inner_range = c(1900, 2000)
outer_range = c(1850, 2050)

dat = data.frame(year = seq(inner_range[1],inner_range[2],length=500))
pred = data.frame(year = seq(outer_range[1],outer_range[2],length=500))


k = 10

#Creating knots. 

#The 4-knot setup suggested for extrapolation in `?smooth.construct.bs.smooth.spec`.
outerknot = c(outer_range[1],inner_range[1],inner_range[2],outer_range[2])

#Creating evenly spaced knots across the range you want to extrapolate over. 
allknot_firstorder = seq(outer_range[1],outer_range[2], length=2*k+2)
allknot_secondorder = seq(outer_range[1],outer_range[2], length=2*k+3)

# smooth basis and predicted basis functions for 1st-order B-spline using the 4-knot approach
first_order_outerknot = smooth.construct(s(year,bs="bs",m=c(1,1),
                                           k=k), 
                                         data=dat,
                                         knots = list(year=outerknot))

# Basis functions for predictig the extrapolated data
first_order_outerknot_pred = Predict.matrix(first_order_outerknot, 
                                            data= pred)

# smooth basis and predicted basis functions for 1st-order B-spline using the multi-knot approach
first_order_allknot = smooth.construct(s(year,bs="bs",m=c(1,1),k=2*k), 
                                       data=dat,
                                       knots = list(year=allknot_firstorder))

first_order_allknot_pred = Predict.matrix(first_order_allknot, data= pred)



# smooth basis and predicted basis functions for 2nd-order B-spline using the 4-knot approach
second_order_outerknot = smooth.construct(s(year,bs="bs",m=c(2,1),k=k), 
                                          data=dat,
                                          knots = list(year=outerknot))

second_order_outerknot_pred = Predict.matrix(second_order_outerknot, 
                                             data= pred)

# smooth basis and predicted basis functions for 2nd-order B-spline using the multi-knot approach
second_order_allknot = smooth.construct(s(year,bs="bs",m=c(2,1),k=2*k), 
                                        data=dat,
                                        knots = list(year=allknot_secondorder))

second_order_allknot_pred = Predict.matrix(second_order_allknot, 
                                           data= pred)



#plotting all four sets of basis functions. Red dashed lines delineate the interior region where data points are present in the training data.
par(mfrow =c(2,2))

matplot(pred$year, first_order_outerknot_pred,
        type="l", 
        main = "first order, 4-knot extrapolation",
        ylim = c(0,1),
        col="black",
        lty=1)
abline(v = inner_range, lty=3, col="red")

matplot(pred$year, second_order_outerknot_pred,type="l", 
        main = "second order, 4-knot extrapolation",
        ylim = c(0,1),
        col="black",
        lty=1)
abline(v = inner_range, lty=3, col="red")

matplot(pred$year, first_order_allknot_pred,
        main = "first order, extra bases extrapolation",
        type="l",
        ylim = c(0,1),
        col="black",
        lty=1)
abline(v = inner_range, lty=3, col="red")

matplot(pred$year, second_order_allknot_pred,type="l", 
        main = "second order, extra bases extrapolation",
        ylim = c(0,1),
        col="black",
        lty=1)
abline(v = inner_range, lty=3, col="red")