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


outerknot = c(outer_range[1],inner_range[1],inner_range[2],outer_range[2])
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)

first_order_outerknot = smooth.construct(s(year,bs="bs",m=c(1,1),
                                           k=k), 
                                         data=dat,
                                         knots = list(year=outerknot))

first_order_outerknot_pred = Predict.matrix(first_order_outerknot, 
                                            data= pred)


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)




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)


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)



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")