1
$\begingroup$

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)

base ggplot

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?
$\endgroup$
2
$\begingroup$

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

enter image description here

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – PLY
    May 10 '20 at 11:32
  • $\begingroup$ 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. $\endgroup$
    – Edward
    May 10 '20 at 12:16
1
$\begingroup$

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

enter image description here

And to answer your questions:

  1. Why is there a difference in the two plots? Does geom_smooth() somehow force the intercept by any chance?

geom_smooth() takes the fitted values from scam.model. The fit is the same. The difference you see comes from plot.scam having a different choice of plotting.

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – PLY
    May 10 '20 at 11:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.