# Interpolation of errors from model predictions over time-series

I have a regression model:

# Create Dataset ####
# simulate exponential relationship
set.seed(123)
# generate random x values between 0 and 60
x <- runif(400, 0, 60)
y <- 1 - exp(-0.1 * x) * rnorm(400, 0.7, 0.1)

data = data.frame(sr= x, fipar = y)

# Fit a GLM model #
glm_mod <- glm(fipar ~ log(sr), data = data, family = binomial)


Let's say that I want to predict this model on new data that has been collected over time:


newdat <- data.frame(sample = c(rep('A',3),rep('B',3),rep('C',3)),
sr = c(1,20,55,4,11,12,2,45,20),
timepoint = c(1,2,3,1,2,3,1,2,3))

newpreds <- predict(glm_mod,
newdata = newdat,
type = "response", se.fit = T)

newdat$$predicted_fipar <- newpreds$$fit
newdat$$predicted_fipar_se <- newpreds$$se.fit



From this model, I get standard errors and can plot errorbars at the timepoint where the predictions were made, i.e

newdat %>%
ggplot(aes(x=timepoint,y=predicted_fipar,group=sample)) +
geom_point(aes(color=sample)) +
geom_line() +
geom_errorbar(aes(ymin=predicted_fipar-predicted_fipar_se,
ymax=predicted_fipar+predicted_fipar_se,
color=sample),
width=0.1)


But let's say I want to know the error at timepoint 1.5. If the interpolation between points is linear, I think that it's possible to also interpolate the error linearly (not sure if this overlooks any statistical rules), i.e:

newdat %>%
ggplot(aes(x=timepoint,y=predicted_fipar,group=sample)) +
geom_point(aes(color=sample)) +
geom_line() +
geom_ribbon(aes(ymin=predicted_fipar-predicted_fipar_se,
ymax=predicted_fipar+predicted_fipar_se,
fill=sample),
alpha=0.4) +
theme(legend.position = 'none')


But what if the interpolation between predicted points is non-linear? Are there other methods that allow this? I was thinking I should first interpolate the independent variable over time (in this case sr) and then make predictions from the model using the interpolated values. This way I can calculate the errors for the in-between points, i.e:

# Interpolate SR over time for each sample #
interpolated_sr <-
plyr::rbind.fill(
lapply(unique(newdat$sample),function(samp){ sampledat <- newdat %>% filter(sample == samp) sr_loess <- loess(sr ~ timepoint, data = sampledat) topred <- data.frame(sample = samp, timepoint = seq(1,3,0.1)) topred$sr <- predict(sr_loess, topred)

return(topred)

}))

interpolated_preds <-
predict(glm_mod,
newdata = interpolated_sr,
type = "response", se.fit = T)

interpolated_sr$$predicted_fipar <- interpolated_preds$$fit
interpolated_sr$$predicted_fipar_se <- interpolated_preds$$se.fit

interpolated_sr %>%
ggplot(aes(x=timepoint,y=predicted_fipar,group=sample)) +
geom_point(aes(color=sample)) +
geom_line() +
geom_ribbon(aes(ymin=predicted_fipar-predicted_fipar_se,
ymax=predicted_fipar+predicted_fipar_se,
fill=sample),
alpha=0.4) +
theme(legend.position = 'none')


Is this approach valid? I understand that by interpolating the independent variable over time I am introducing some additional error from the loess model (which in this case will be very small, but could be very big depending on the fit). Is this where error propagation comes into play?

Your predictions and associated prediction intervals (not really errors as such) will only vary according to the value of the sr variable. So you cannot really get much sense of how these intervals will change between timepoints if you don't include sr values for any times between those timepoints. To illustrate, I below show how you can use plot_predictions() from the {marginaleffects} package to automatically compute and plot predictions and their associated uncertainties. Hopefully this makes it clear that you don't actually get any measure of uncertainty for the intervening timepoints because there is no sr value there to predict on. The plot simply draws a polygon that interpolates linearly, but this is merely a function of the plotting mechanism and not the underlying predictions:

# Simulate exponential relationship
set.seed(123)

# Generate random x values between 0 and 60
x <- runif(400, 0, 60)
y <- 1 - exp(-0.1 * x) * rnorm(400, 0.7, 0.1)

# y is a proportional variable here
hist(y)


data = data.frame(sr = x, fipar = y)

# Fit a GLM model using a Beta observation model, more appropriate
# for the proportional outcome that doesn't have a number of trials
library(mgcv)
glm_mod <- gam(fipar ~ log(sr), data = data, family = betar())

# Investigate predictions across an evenly-spaced grid of sr values
library(marginaleffects)
library(ggplot2)
plot_predictions(glm_mod, condition = 'sr') +
theme_classic()


# Predict over the 'newdat' grid
newdat <- data.frame(sample = c(rep('A',3),
rep('B',3),
rep('C',3)),
sr = c(1,20,55,4,11,12,2,45,20),
timepoint = c(1,2,3,1,2,3,1,2,3))
plot_predictions(glm_mod, newdata = newdat, by = 'timepoint') +
theme_classic()


# Look at the data and you'll see that there is no prediction for the
# intervening timepoint; the intervals you see are merely a byproduct of the
# plotting mechanism
plot_predictions(glm_mod, newdata = newdat, by = 'timepoint',
draw = FALSE)
#>   timepoint  estimate   std.error statistic       p.value  s.value  conf.low
#> 1         1 0.3012604 0.011258944  26.75743 1.011827e-157 521.5257 0.2791933
#> 2         2 0.9015958 0.002435880 370.13147  0.000000e+00      Inf 0.8968215
#> 3         3 0.9099221 0.002274289 400.09078  0.000000e+00      Inf 0.9054645
#>   conf.high
#> 1 0.3233276
#> 2 0.9063700
#> 3 0.9143796


Created on 2024-04-04 with reprex v2.1.0

• Thank you, I understand that the prediction intervals are only computed for the timepoints where the independent variable was observed. For certain situations, it would be reasonable to interpolate these to produce 'new' observations, as that variable might follow a predictable pattern over time. Is it possible to use the interpolated values of that variable to then predict and compute intervals like I show in the third graph? Apr 4 at 2:28
• Ok I see. Yes I believe that is reasonable, though you could generate many replicate data grids that each take a different possible 'draw' from the LOESS and make predictions for those. To get better prediction intervals, I'd first generate 'draws' of model coefficients (i.e. coefs = MASS::mvrnorm(1000, coef(glm_mod), vcov(glm_mod))) and make predictions for each data grid using each 'draw' of coefficients. You can then summarise predictions using quantile() to get more realistic intervals that incorporate the uncertainty in the LOESS smooths Apr 4 at 2:46