Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
This is an R example of our model.
### 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment