# STAT3530 -- randomization and bootstrap inference
# Savings rates for 50 countries. Two predictors describe the age structure of
# the population (pop15, pop75) and two describe the economy (dpi, ddpi).
# Run a line at a time.

site <- 'https://tdruiz-stat3530.share.connect.posit.cloud/_data/'

load(url(paste0(site, 'savings.RData')))

m1 <- lm(sr ~ ., data = savings)               # full model
m0 <- lm(sr ~ pop15 + pop75, data = savings)   # model under H0: no economic effect

anova(m0, m1)                                  # the test we are reproducing

# ===========================================================================
# THE TWO FUNCTIONS
# ===========================================================================
# sim.data() makes a new dataset. estimand() computes the statistic of interest

sim.data <- function(fit) fitted(fit) + sample(resid(fit))

estimand <- function(y) anova(lm(y ~ pop15 + pop75, data = savings),
                              lm(y ~ pop15 + pop75 + dpi + ddpi, data = savings))$F[2]

estimand(savings$sr)        # on the real data

estimand(sim.data(m0))      # on one simulated dataset

# run the simulation
set.seed(1)
null <- replicate(2000, {
  yb <- sim.data(m0)
  Fb <- estimand(yb)
  Fb
})

# visualize
hist(null, breaks = 40, freq = FALSE, col = 'grey90', border = 'grey60',
     xlab = 'F', main = '')
abline(v = estimand(savings$sr), col = 'red', lwd = 2, lty = 2)

# compute the p-value
mean(null >= estimand(savings$sr))

# ===========================================================================
# EXERCISE 1: change the simulation method
# ===========================================================================
# Rewrite sim.data() to draw fresh errors instead of reshuffling the old ones:
#
#     fitted(fit) + rnorm(nobs(fit), 0, sigma(fit))
#
# Rebuild null, take the p-value again, and compare the three.

# replace the ...
sim.data <- function(fit) ...

# no change here, same as above
set.seed(1)
null <- replicate(2000, {
  yb <- sim.data(m0)
  Fb <- estimand(yb)
  Fb
})

# visualize
hist(null, breaks = 40, freq = FALSE, col = 'grey90', border = 'grey60',
     xlab = 'F', main = '')
abline(v = estimand(savings$sr), col = 'red', lwd = 2, lty = 2)

# compute p-value: you supply code


# ===========================================================================
# EXERCISE 2: change the estimand
# ===========================================================================
# Rewrite estimand() to return R^2 instead of the F statistic:
#
#     summary(lm(y ~ pop15 + pop75 + dpi + ddpi, data = savings))$r.squared
#
# Simulate from m1 rather than m0, and report a 95% percentile interval.


# replace the ...
estimand <- function(y) ...

# now you write the simulation code
set.seed(1)
null <- replicate(2000, {
  ...
})

# visualize: you supply code


# confidence interval