# Extra 4

## a) fit true model
n <- 30                     # sample size
beta0 <- 1.0                # true intercept value
sig <- 1.0                  # variance
beta1v <- seq(0,2,by=0.1)   # true slope value(s)
N <- 100                    # number of simulation runs
Pval <- matrix(nrow=length(beta1v),ncol=N) # empty matrix for storing p-values
for(i in 1:length(beta1v))
{
  for(j in 1:N)
  {
    beta1 <- beta1v[i]
    x <- rnorm(n)                          # draw n predictor values from normal dist
    y <- rnorm(n,beta0+beta1*x,sig)        # draw n outcome values from conditional normal dist
    fit <- lm(y~x)                         # fit linear regression model
    Pval[i,j] <- summary(fit)$coef["x",4]  # record p-value
  }
}

# compute rejection rate
alpha <- 0.05
rej.rate <- rowMeans(Pval<alpha)     # rate of significant coefficients
plot(beta1v,rej.rate,type="l", ylim = c(0, 1))


## c) extended model
Pval <- matrix(nrow=length(beta1v),ncol=N)
Pvalx <- matrix(nrow=length(beta1v),ncol=N)
Pvalz <- matrix(nrow=length(beta1v),ncol=N)
Pvalxz  <- matrix(nrow=length(beta1v),ncol=N)
coef <- matrix(nrow=length(beta1v),ncol=N)
coefx <-  matrix(nrow=length(beta1v),ncol=N)

rho = 0.99 # set correlation coefficient
for(i in 1:length(beta1v))
{
  for(j in 1:N)
  {
    beta1 <- beta1v[i]
    beta2 <- 0
    x <- rnorm(n)
    z <- rho*x+sqrt(1-rho^2)*rnorm(n)
    y <- rnorm(n,beta0+beta1*x+beta2*z,sig)

    # fit true model
    fit <- lm(y~x)
    Pval[i, j] <- summary(fit)$coef["x",4]
    coef[i, j] <- summary(fit)$coef["x", "Estimate"]


    # fit model including z
    fit <- lm(y~x+z)
    Pvalx[i,j] <- summary(fit)$coef["x",4]
    Pvalz[i,j] <- summary(fit)$coef["z",4]
    coefx[i, j] <- summary(fit)$coef["x", "Estimate"]
  }
}

# compute rejecetion rates
alpha <- 0.05
rej.rate  <- rowMeans(Pval<alpha)
rej.ratex <- rowMeans(Pvalx<alpha)
rej.ratez <- rowMeans(Pvalz<alpha)
rej.ratexz <- rowMeans(Pvalx<alpha | Pvalz<alpha)

matplot(cbind(beta1v,beta1v,beta1v,beta1v,beta1v), cbind(rej.rate,rej.ratex,rej.ratez,rej.ratexz),
        type = "l", lty = 1, col = 1:4, xlab = expression(paste(beta[1])), ylab = "rejection rate")
legend("topright", legend = c("x with only x","x","z","x or z"), lty = 1, col = 1:4, bg = "white")


# plot histogram of estimated coefficient on x
# - much more variation in the estimate when we include in the model,
#   as parts of the variation in y is assigned to z.
breaks <- seq(-5, 5, length.out = 100)
hist(coef[11, ], breaks = breaks)
hist(coefx[11, ], breaks = breaks, add = T, col=rgb(0,0,1,1/4))
abline(v = 1, col = 2)