# Extrapolating from a smooth monotonic curve in R

I am trying to extrapolate a smooth monotonic curve made with scam(). Please see the following code and produced images:

library(scam)
library(mgcv)
library(nlme)
library(ggplot2)

DF <- structure(list(Date = structure(c(15548, 16149, 16513,
16858,  17211, 17617, 17985), class = "numeric"),
percentile = c(1.90982475506864, 2.94546018485894,
3.00516383648267, 2.96295221921635, 3.7070240276891,
6.52469881377116, 7.70470323658277)), row.names = c(NA,
-7L), class = "data.frame")

# Fit  monotonically increasing convex spline
scam.model <- scam(percentile ~ s(Date, k = 4, bs = 'micx'),
data = DF)

# Missing an intercept perhaps?
plot(scam.model); points(as.numeric(DF$$Date), DF$$percentile,
col = 'red')

# Fit as expected
ggplot(DF, mapping = aes(Date,percentile)) +
geom_point() +
geom_smooth(method = 'scam',
formula = y ~ s(x, k = 4, bs = 'micx'),
se = T)


My questions are the following:

1. Why is there a difference in the two plots? Does geom_smooth() somehow force the intercept by any chance?
2. Is it possible to determine the value of Date for which percentile would surpass a threshold? The other way around would be easier with some predict() function I assume, but for my use I am interested in this direction.
3. Would you recommend another techniques and packages for this kind of extrapolation?

I think I can only answer 2.

Create a new data frame containing many Date values within the range of your data and add the predicted values from your model.

newdata <- data.frame(Date=seq(min(DF$$Date), max(DF$$Date), by=1))
newdata$pr <- predict.scam(scam.model, newdata = newdata)  You seek the value of Date for which percentile would exceed a certain threshold. Let's say this threshold is 4. newdata[which.max(newdata$pr>4),]
#       Date      pr
# 1544 17091 4.00036


• I see, clever solution! However, not to undermine your answer, I was hoping to find a more explicit way of finding the threshold rather than predicting for several values for Date. This would still work in my case though, so thank you.
– PLY
May 10 '20 at 11:32
• Those gam/scam models always tie me up in knots; I can never figure out the closed-form formula. So I just let the computer figure it out. May 10 '20 at 12:16

I think I can only answer 1

When you do plot(scam.model), it is in fact plot.scam(scam.model) and in the source code of plot.scam, it chooses to omit the intercept and use only your other coefficients. This is the relevant lines:

if (m>0) for (i in 1:m) { ## work through smooth terms
first <- x$$smooth[[i]]$$first.para
last <- x$$smooth[[i]]$$last.para


I guess it is meant for comparisons with other models, and you might have to read it in more detail. So to get the values fitted, you can do:

plot(scam.model,shift=scam.model$$coefficients[1,]) points(as.numeric(DF$$Date), DF\$percentile, col = 'red')


• Thank you for clarifying! It still makes me wonder why their base plot() function would not simply plot the fit. For comparing with other models I would look at the coefficients explicitly. Regardless, thank you for taking the time to check in the code.