This is an R example of our model.

simple-example.R

### Let's simulate some data very simply

## Load libraries
library(ggplot2)
library(dplyr)
library(gam)
library(patchwork)

# Set parameters
n_sample = 100
a = 1
b = 2
sig = 2

## Create two covariates
dat = tibble(x1 = rnorm(n_sample*3),
             x2 = rnorm(n_sample*3),
             eps = rnorm(n_sample*3,sd=sig),
             split = rep(c("tr","te","va"),each=n_sample))

dat = dat %>% 
  mutate(y = a*x1 + b*x2^2 + eps )



## Look at the data
p1 = dat %>% ggplot(aes(x1,y)) + geom_point() + theme_minimal()
p2 = dat %>% ggplot(aes(x2,y)) + geom_point() + theme_minimal()
p1 + p2

## Make prediction

g1 = gam(y ~ s(x1) + s(x2), data=dat %>% filter(split=="tr"))

## Add predictions
dat$pred = predict(g1,newdata=dat)

## Plot predictions
dat %>%
  filter(split=="te") %>%
  ggplot(aes(y,pred)) + geom_point() + theme_minimal()


## Fit model and calculate standard errors

rel_mod = lm(y ~ pred,data=dat %>% filter(split=="te"))
sigma(rel_mod)

## Get standard error and compare to unadjusted

inf_pred = lm(pred ~ x1,data=dat %>% filter(split=="va"))

inf_real = lm(y ~ x1,data=dat %>% filter(split=="va"))

inf_pred %>% tidy()
inf_real %>% tidy()


## Get adjustment
mod_matrix = model.matrix(~ x1,data=dat %>% filter(split=="va"))

inf_factor = sigma(rel_mod)^2 + rel_mod$coefficients[2]^2*sigma(inf_pred)^2

adj_se = sqrt(solve(t(mod_matrix) %*% mod_matrix)[2,2]*inf_factor)
adj_se