---
title: "CSSS/POLS 510 - Maximum Likelihood Estimation"
subtitle: "Lab 4. Prediction and Quantities of Interest"
author: "Minji Jeong"
fontfamily: accanthis
linkcolor: blue
output:
  html_document:
    toc: true
    toc_depth: 4
editor_options:
  chunk_output_type: console
fontsize: 11pt
header-includes:
 \usepackage{float}
 \usepackage{xcolor}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, error=FALSE, warning=FALSE, message = FALSE)

# clean environment

rm(list=ls())

# Load packages
library(MASS)
library(simcf)
library(tidyverse)    # For data wrangling


```


# Simulation of Quantities of Interest

**Simulation** is a powerful method for generating meaningful **quantities of interest (QI)** from statistical models. Instead of only reporting coefficients, researchers can simulate **expected values**, **expected probabilities**, **first differences**, and measures of **uncertainty** to give a clearer picture of the model results.

QI are especially valuable because they go beyond mere coefficients to provide a more intuitive understanding of the model’s output. While it is possible to compute expected probabilities using the model’s coefficients alone, simulation adds value by incorporating **confidence intervals**, thus addressing **uncertainty**.

In any statistical model, the coefficients are subject to **estimation uncertainty**—we estimate parameters from a sample but want to generalize to the broader population. If a different sample were collected, the estimated coefficients could vary, leading to uncertainty in any computed quantities from those coefficients. **Simulation quantifies this uncertainty** by drawing from the **sampling distribution** of the model’s estimates.

These draws reflect the variability of each coefficient and their covariances, resulting in a large number of plausible estimates. Rather than simply summarizing these draws, it is more informative to compute a **substantively meaningful QI** from the simulated estimates, such as a predicted or expected value of the dependent variable.

### Two Types of Uncertainty

When computing **predictions** and **simulations** of **quantities of interest (QoI)**, it is essential to recognize that the model we use to approximate the **data-generating process (DGP)** is subject to two types of **uncertainties**:


  - **Estimation uncertainty** arises because the model's parameter estimates $\hat{\beta}$ are **random variables** because they are a function of the sample/data. Hence, we must reflect the uncertainty of $\hat{\beta}$.
  - **Fundamental uncertainty** reflects the inherent randomness of the world by incorporating uncertainty from the stochastic component and applies when making predictions beyond the model's estimation.

Thus, the choice between **predicted** and **expected values** depends on the type of uncertainty the analyst wants to highlight—whether focusing on the uncertainty about the model itself (**estimation uncertainty**) or on the variability in the world (**fundamental uncertainty**).


In a simple linear regression model, we can highlight both sources of uncertainty. The equation is as follows:

$$
\hat{y_i} = \color{red}{\hat{\beta_0}} + \color{red}{\hat{\beta_1}} x_i + \color{blue}{\hat{\epsilon}_i}
$$




  - $\color{red}{\text{Estimation uncertainty}}$: Represents uncertainty in the estimation of the coefficients $\hat{\beta_0}$ and $\hat{\beta_1}$, which are estimates derived from the data.

  - $\color{blue}{\text{Fundamental uncertainty}}$: combines uncertainty from the 
$\color{red}{\text{systematic}}$ and $\color{blue}{\text{stochastic}}$ components, accounting for the randomness in the dependent variable that the model cannot explain.


Which **quantity of interest** should we report: **expected values** or **predicted values**? When analyzing **estimands** at the population or **data-generating process (DGP)** level, we focus on **expected values**. However, if the goal is to make specific predictions for real-world cases—such as in **policy analysis** for government agencies, private banks, or consulting firms—**predicted values** should be reported. 


In summary, **estimation uncertainty** relates to the estimates $\hat{\beta}$ in the model, while **fundamental uncertainty** incorporates the **residuals** $\hat{\epsilon}$ and reflects real-world randomness.


## Application: Ehlirhc (1973) Crime Data


The following application is based on [Carse and Harden (2014, Chapter 9)](https://methods.sagepub.com/book/monte-carlo-simulation-and-resampling-methods-for-social-science), and uses [Ehlirhc (1973)](https://www.journals.uchicago.edu/doi/abs/10.1086/260058) Crime Data.

The **Ehrlich (1973) Crime Data** was originally collected by Isaac Ehrlich from various U.S. government sources, including the **Uniform Crime Report** of the FBI. It has been used in multiple studies to analyze how different variables influence the rate of crime at the state level in the United States during 1960.

We will use the following variables:
  
  - **Ed**: Mean years of schooling of the population aged 25 years or over.
  - **Pop**: State population in 1960 (in hundreds of thousands).
  - **Wealth**: Median value of transferable assets or family income.
  - **Ineq**: Income inequality, defined as the percentage of families earning below half the median income.
  - **Prob**: Probability of imprisonment, calculated as the ratio of the number of commitments to the number of offenses.



### Simulation: Expected Values


If we are interested in analyzing the **marginal effects** of laws and socioeconomic characteristics at the **DGP/population level**, we aim to generalize results (external validity). In these cases, we primarily report **expected values**. Specifically:

$$
\hat{\mathbf{y}} = \mathbf{X} \color{red}{\hat{\boldsymbol{\beta}}} + \hat{\epsilon}
$$

We simulate expected values based on the estimates of the $\color{red}{\text{systematic component}}$ ($\hat{\boldsymbol{\beta}}$), reflecting **estimation uncertainty** from the sample.

To compute expected values $\mathbb{E}[y | X_c \tilde{\beta}] = \tilde{\mu_c}$ via simulation:

  1. Extract the model's $k$ point **estimates** $\hat{\beta}_k$ and their **variance-covariance matrix** $\hat{V}(\hat{\beta}_k)$. Use these to draw $n$ simulated random draws from the estimated coefficients: $\boldsymbol{\tilde{\beta}} \sim \mathcal{MVN}(\boldsymbol{\hat{\beta}},\hat{V}[\boldsymbol{\hat{\beta}}])$. Here, $\boldsymbol{\tilde{\beta}}$ is an $n \times k$ matrix of simulated coefficients.
      - This assumes the dependent variable's **data-generating process** follows a multivariate normal distribution $\boldsymbol{y} \sim \mathcal{MVN}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$.
      - The assumption is plausible because (i) the MLE estimates $\hat{\beta}$ are normally distributed and (ii) normality is guaranteed asymptotically by the **Central Limit Theorem**.
  
  2. With $n$ simulated draws of $\boldsymbol{\tilde{\beta}}$, create a $1 \times k$ transposed vector of predictor values $\boldsymbol{x'}$ for some scenario $c$ of interest, e.g., at their means $\boldsymbol{x_c'}$.
  
  3. Compute the expected values $\boldsymbol{\tilde{\mu_c}}$ using matrix algebra: $\boldsymbol{x_c'} \times \boldsymbol{\tilde{\beta}}$. Here, $\boldsymbol{\tilde{\mu_c}}$ is an $n \times 1$ vector of expected values.
  
  4. Visualize $\tilde{\mu_c}$ and/or report statistics from its simulated distribution, such as the mean or quantiles to represent the confidence intervals.

The code below demonstrates this using **Ehlirch's crime data**.



```{r}
crime <- read.csv(file="data/crime.csv")

names(crime)

# OLS model form the crime data

ols.crime <- lm(Crime ~ Prob + Ed + Wealth + Ineq + Pop,
                data= crime)

summary(ols.crime)

# Required library to draw from a multivariate distribution
library(MASS)

set.seed(139742)

### Simulation EV: step 1

# Number of simulated sets of estimates
sims <- 1000

# Create empty vector to store expected values
crime.ev <- numeric(length=sims)

# Extract coefficient estimates
(ols.pe <- ols.crime$coef)

# Extract variance-covariance matrix
(ols.vcv <- vcov(ols.crime))


# create matrix of random draws from point estimates and variance-covariance
crime.sim <- MASS::mvrnorm(n = sims,
                           mu = ols.pe,
                           Sigma = ols.vcv)

# Observe the output
head(crime.sim)

# Check out the mean
apply(crime.sim, 2, mean)

### Simulation EV: step 2

# Set hypothetical values for the independent variables (at their means)
crime.x <- tibble(intercept = 1,
                  Prob = mean(crime$Prob),
                  Ed = mean(crime$Ed),
                  Wealth = mean(crime$Wealth),
                  Ineq = mean(crime$Ineq),
                  Pop = mean(crime$Pop))


### Simulation EV: step 3

# Compute one expected value to get the intuition

(ev1 <- as.matrix(crime.x) %*% crime.sim[1,]) # arguments must be conformable!

# Once you get the desired output, program the loop

for (i in 1:sims) {
  crime.ev[i] <- as.matrix(crime.x) %*% crime.sim[i,]
}

# Browse expected values
head(crime.ev)


### Simulation EV: step 4

(mean_ev <- mean(crime.ev))

# We can now represent the expected value and its uncertainty

# CIs at a 0.05 significance level
(lwr <- quantile(crime.ev, probs = 0.025))
(upr <- quantile(crime.ev, probs = 0.975))
```


```{r, echo=FALSE, fig.align='center'}
# Visualize distribution of expected values
tibble(pe=crime.ev) |> 
  ggplot(aes(x=pe)) +
  geom_density(color="dodgerblue4",
               fill="skyblue",
               alpha=0.3) +
  geom_rug(color="darkslateblue") +
  geom_vline(aes(xintercept = mean(crime.ev)),
             linetype = "dashed", 
             color = "red", linewidth = 1)  +
  geom_vline(aes(xintercept = lwr),
             linetype = "dotted", 
             color = "red", linewidth = 1)  +
  geom_vline(aes(xintercept = upr),
             linetype = "dotted", 
             color = "red", linewidth = 1)  +
  labs(title = "Simulated Expected Values From the Ehrlich (1973) Crime Data",
       subtitle = "Simulated Expected Values, Mean EV = 909",
       x = "Expected Values", 
       y = "Density") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        text = element_text(size = 12),
        plot.caption = element_text(size = 8, face = "italic"))
```


#### Simulating QoI: First Differences

Although displaying a distribution of **expected values** is an effective way to communicate model results, by itself it doesn't add much beyond simply reporting the average predicted expectation.

The real advantage of reporting **estimands** or **quantities of interest (QoI)** via simulation is to provide and contrast **meaningful scenarios** that address your research questions. For example, suppose we hypothesize that **economic inequality** moderates the **crime rate** at the state level. We can compute expected values for two scenarios by adding the following steps:

  1. **Low inequality**: Compute the expected value at the minimum observed level of inequality.
  2. **High inequality**: Compute the expected value at the maximum observed level of inequality.
  3. **First Difference (QoI)**: Compute the difference between the expected values of high and low inequality and determine if this difference is significantly different from zero.



```{r}
# Set hypothetical values for the independent variables

(crime.x2 <- tibble(intercept = 1,
                    Prob = mean(crime$Prob),
                    Ed = mean(crime$Ed),
                    Wealth = mean(crime$Wealth),
                    # notice values of inequality
                    Ineq = c(min(crime$Ineq),
                             max(crime$Ineq)),
                    Pop = mean(crime$Pop)))


# Matrix multiply both sets of hypothetical values 
# each set of a simulated coefficients
# Make 3 columns - 1st: expected values of low inequality; 2st: exp. val. of high inequality; 3rd: compue the difference

crime.ev2 <- matrix(NA,
                    nrow = sims,
                    ncol = 3)

# Because we have two scenarios, we need parallel simulations

for (i in 1:sims) {
  # two scenarios: low and high inequality
  crime.ev2[i, 1] <- as.matrix(crime.x2[1, ]) %*% crime.sim[i, ]
  crime.ev2[i, 2] <- as.matrix(crime.x2[2, ]) %*% crime.sim[i, ]
  
  # estimand - quantity of interest: first difference
  # Think low inequality = treatment, high inequality = control, difference = association; causal effect
  crime.ev2[i, 3] <-  crime.ev2[i, 2] - crime.ev2[i, 1]
}


## Prepare the data for visualization (long format for ggplot)

data_vis <- tibble(
  pe=c(crime.ev2[,1], # low inequality estimates
       crime.ev2[,2], # high inequality estimates
       crime.ev2[,3]), # high - low = difference
  
  # labels for scenarios
  scenario = c(rep("Low Ineq",sims),
               rep("High Ineq",sims),
               rep("First-Diff",sims)),
  
  # labels of quantities
  QIs = c(rep("Expected Values",sims*2),
          rep("First Differences",sims))
)
```


```{r, echo=FALSE, fig.align='center'}
## Visualize Quantities of interest

data_vis |> 
    ggplot(aes(x=pe,
             group=scenario,
             color=scenario,
             fill=scenario)) +
  geom_density(alpha=0.2) +
  geom_rug() +
  scale_color_manual(values = c("Low Ineq" = "#85DD41",
                                "High Ineq" = "#F9994A",
                                "First-Diff"= "#633EAA")) +
  scale_fill_manual(values = c("Low Ineq" = "#85DD41",
                               "High Ineq" = "#F9994A",
                                "First-Diff"= "#633EAA")) +
  geom_vline(aes(xintercept = 0),
             linetype = "dashed", 
             color = "black") +
  labs(title = "Simulated Expected Values From the Ehrlich (1973) Crime Data",
       subtitle = "QoIs: Marginal Effect of Inequality on Crime (from low to high)",
       x = "Estimates", 
       y = "Density",
       fill = "", color = "") +
  theme_minimal() +
  facet_wrap(~QIs, scales = "free") +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        legend.position = "bottom",
        text = element_text(size = 12),
        legend.text = element_text(size = 10),
        plot.caption = element_text(size = 8, face = "italic"))

# the first difference distribution shows relative marginal effect of increasing inequality (from low to high inequality) 

```


Simulation results show that the **simulated coefficients** average out to be very close to the actual estimated coefficients. Importantly, you could increase the precision of these simulations by increasing the number of iterations (e.g., the `sims` object).


### Simulation: Predicted Values


When we are interested to provide predictions for real, concrete cases (e.g. forecast of crime in Washington state), then we need to compute predicted values:

$$
\hat{\mathbf{y}} = \mathbf{X} \color{red}{\hat{\boldsymbol{\beta}}} + \color{blue}{\boldsymbol{\hat{\epsilon}}}
$$

To compute predicted values, we need to simulate both the estimates from the $\color{red}{\text{systematic component}}$ ($\hat{\boldsymbol{\beta}}$) **and** the $\color{blue}{\text{stochastic component}}$ ($\hat{\boldsymbol{\epsilon}}$)


To generate **predicted values** of the dependent variable, we first compute the **expected values**, and then add random draws from the appropriate distribution. For **OLS**, we assume a **normal distribution** for the dependent variable, with the independent variables affecting the parameter $\mu$.

However, we must also account for the normal distribution's other parameter, **$\sigma$**. The `lm()` function in R estimates $\hat{\sigma}$ using the formula:

$$
\hat{\sigma} = \sqrt{\frac{\sum_{i=1}^n \left( Y_i - X_i \hat{\beta} \right)^2}{n - p}}
$$

Where:

  - $Y$ is the dependent variable,
  - $X$ is the matrix of independent variables,
  - $\hat{\beta}$ is the vector of estimated coefficients,
  - $n$ is the sample size,
  - $p$ is the number of estimated coefficients.

The numerator represents the **sum of squared residuals (SSR)**, while the denominator adjusts for **degrees of freedom**.

#### Simulating $\hat{\sigma}$

We cannot simply use the model's $\hat{\sigma}$ for generating predicted values because it, too, is an estimate. To generate **predicted values**, we must simulate estimates of both the **coefficients** and **variance**. We simulate $\hat{\sigma}$ by applying the formula above to each set of simulated coefficients.



```{r}

# model residual degrees of freedom
crime.DoF <- ols.crime

# Empty vector
crime.sigma <- numeric(sims)

# Loop thorugh sets of simulated coefficients

for (i in 1:sims) {
  # create the matrix of independent variables
  x <- as.matrix(cbind(1, ols.crime$model[ , -1]))
  
  # Compute the SSR for each set of simulated coefficeints
  ssr <- sum((crime$Crime - (x %*% crime.sim[i, ]))^2)
  
  # Estimate sigma
  crime.sigma[i] <- sqrt(ssr/ols.crime$df.residual)
}



```


To generate a single predicted value, we provide two quantities to `rnorm()`:

  1. **Mean**: The expected value calculated from the simulated estimates of $\hat{\beta}$.
  2. **Standard deviation**: The value of $\hat{\sigma}$ calculated from the simulated estimates.

By doing this for each set of simulated estimates, we incorporate both **estimation uncertainty** (from the coefficients $\hat{\beta}$) and **fundamental uncertainty** (from the error term $\hat{\sigma}$) into the predicted values.

The code example below shows how to do this when all predictors are set to their means:


```{r}
## Example: one predicted value with both estimation and fundamental uncertainty

(mu <-  as.matrix(crime.x) %*% crime.sim[1, ])

(sigma <- crime.sigma[1])

(example.pred <- rnorm(1, mean=mu, sd=sigma))

# Once you get the intuition, and the desired output, you can program the loop.

# Matrix multiply x predictors and each set of simulated coefficients to create mu
# then take a random draw from the normal distribution with mean equal to mu and
# sd equal to sigma
# as we introduce uncertainty, the confidence interval will become larger

crime.pred <- numeric(sims)

for (i in 1:sims) {
  mu <- as.matrix(crime.x) %*% crime.sim[i, ]
  sigma <- crime.sigma[i] 
  crime.pred[i] <- rnorm(1, mean=mu, sd=sigma)
}


# Browse expected values
head(crime.pred)

(mean_pred <- mean(crime.pred))

# We can now represent the expected value and its uncertainty

# CIs at a 0.05 significance level
(lwr <- quantile(crime.pred, probs = 0.025))
(upr <- quantile(crime.pred, probs = 0.975))
```



```{r, echo=FALSE, fig.align='center'}
# Visualize distribution of expected values
tibble(pe=crime.pred) |> 
  ggplot(aes(x=pe)) +
  geom_density(color="dodgerblue4",
               fill="skyblue",
               alpha=0.3) +
  geom_rug(color="darkslateblue") +
  geom_vline(aes(xintercept = mean_pred),
             linetype = "dashed", 
             color = "red", linewidth = 1)  +
  geom_vline(aes(xintercept = lwr),
             linetype = "dotted", 
             color = "red", linewidth = 1)  +
  geom_vline(aes(xintercept = upr),
             linetype = "dotted", 
             color = "red", linewidth = 1)  +
  labs(title = "Simulated Predicted Values From the Ehrlich (1973) Crime Data",
       subtitle = "Simulated Predicted Values, Mean EV = 909",
       x = "Predicted Values", 
       y = "Density") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        plot.subtitle = element_text(hjust = 0.5),
        text = element_text(size = 12),
        plot.caption = element_text(size = 8, face = "italic"))
```


This process ensures our predictions account for **fundamental uncertainty**.





# Recap: Heteroskedastic Normal Example


In this section, we will revisit the Heteroskedastic Normal model from last week. We will later compute quantities of interest using `simcf()`.


```{r, warning=F, message=F}
rm(list=ls())       # Clear memory

set.seed(123456)    # For reproducible random numbers
library(MASS)       # Load packages
library(simcf) 
library(tidyverse)

n <- 1500           # Generate 1500 observations

w0 <- rep(1, n)     # Create the constant
w1 <- runif(n)      # Create two covariates
w2 <- runif(n)

x <- cbind(w0, w1, w2)       # Create a matrix of the covariates
z <- x                       # i.e., same covariates affect mu and sigma

beta <- c(0, 5, 15) # Set a parameter vector for the mean
                    # One for constant, one for covariate 1, 
                    # one for covariate 2.

gamma <- c(1, 0, 3) # Set a parameter vector for the variance 
                    # Gamma estimate for covariate 2 is set to be 3,
                    # whic creates heteroskedasticity 

mu <- x %*% beta             # Create systemtic component for the mean

sigma2 <- exp(z %*% gamma)   # Create systematic component for the variance
                             # Since ith row of sigma2 = exp(1 + 0 + w2_i * 3)
                             # i.e., it is a function of w2 (heteroskedastic)

y <- rnorm(n = n,            # Create the outcome variable
           mean = mu,        # Think about the stochastic component!
           sd = sqrt(sigma2)
           ) 

data <- cbind(y, w1, w2)     # Save the data to a data frame
data <- as.data.frame(data)
names(data) <- c("y", "w1", "w2")


ls.result <- lm(y ~ w1 + w2, data = data) # Fit a linear model
                                          # using simulated data

ls.aic <- AIC(ls.result) # Calculate and print the AIC 
                         # i.e. Akaike Information Criterion
                         # to assess goodness of fit; lower AIC is better


# A likelihood function for ML heteroskedastic Normal
llk.hetnorm <- function(param, y, x, z) {
  x <- as.matrix(x)		                 # x (some covariates) as a matrix
  z <- as.matrix(z)		                 # z (some covariates) as a matrix
  os <- rep(1, nrow(x)) 	             # Set the intercept as 1 (constant)
  x <- cbind(os, x)		                 # Add intercept to covariates x
  z <- cbind(os, z)                    # Add intercept to covariates z
  
  b <- param[1 : ncol(x)]                         # Parameters for x
  g <- param[(ncol(x) + 1) : (ncol(x) + ncol(z))] # Parameters for z
  
  xb <- x %*% b                          # Systematic components for mean
  s2 <- exp(z %*% g)                     # Systematic components for variance
  
  sum(0.5 * (log(s2) + (y - xb)^2 / s2)) # Likelihood we want to maximize 
                                         # optim is a minimizer by default
                                         # To maximize lnL is to minimize -lnL
                                         # so the +/- signs are reversed
}


# Create input matrices
xcovariates <- cbind(w1, w2)
zcovariates <- cbind(w1, w2)

# Initial guesses of beta0, beta1, ..., gamma0, gamma1, ...
# We need one entry per parameter, in order!
# Note: also include beta and gamma estiamtes for constants
stval <- c(0, 0, 0, 0, 0, 0)  

# Run ML, and get the output we need
hetnorm.result <- optim(stval,           # Initial guesses
                        llk.hetnorm,     # Likelihood function
                        method = "BFGS", # Gradient method
                        hessian = TRUE,  # Return Hessian matrix
                        y = y,           # Outcome variable
                        x = xcovariates, # Covariates x (w/o constant)
                        z = zcovariates  # Covariates z (w/o constant)
                        ) 
```

Note: By default, `optim()` performs a minimizing procedure. You can make optim a maximizer by adding `control = list(fnscale = -1)`


```{r}
pe <- hetnorm.result$par            # Point estimates
round(pe, 3)

vc <- solve(hetnorm.result$hessian) # (inverse of) Var-cov matrix (for computing s.e.)
round(vc, 5)

se <- sqrt(diag(vc))                # To compute standard errors (s.e.)
                                    # take the diagonal of the Hessian; 
                                    # then take square root
round(se, 3)

mle.result <- cbind(pe[1:3], se[1:3])         # Report ML results 
colnames(mle.result) <- c("Estimate", "Std.Error")
rownames(mle.result) <- c("(Intercept)", "w1", "w2")
round(mle.result, 3)

round(coef(summary(ls.result))[, c(1, 2)], 3) # Compare with lm() results

ll <- -hetnorm.result$value                   # Report likelihood at maximum 
                                              # No need to have negative sign 
                                              # if optim is set as maximizer
ll


# AIC is 2 * number of parameters - 2 * ll (i.e. likelihood at maximum)
hetnorm.aic <- 2 * length(stval) - 2 * ll  

# Compare AIC from ML and from lm(); lower is better
print(hetnorm.aic)
print(ls.aic)
```


## Prediction of QoI using R's `predict` function.

**Motivation**: We want to study how the change in a particular explanatory variable affects the outcome variable, *all else being equal*

### Scenario 1: Vary covariate 1; hold covariate 2 constant

  1. Create a data frame with a set of hypothetical scenarios for covariate 1, while keeping covariate 2 at its mean
      + What is the sensible range of some hypothetical scenarios for covariate 1? Consider the original range of `w1`.

  2. Calculate the predicted values using the `predict()` function 
      + Hint: you need at least the following arguments: `predict(object = ... , newdata = ... , interval = ... , level = ...)`

  3. Plot the prediction intervals

  4. Similarly, calculate the confidence intervals using the `predict()` function

  5. Plot the confidence intervals; compare them with the predictive intervals


```{r}
# Set up

# Set up hypothetical values for w1 and w2
(w1range <- seq(from = 0, to = 1, by = 0.05))
(w2mean <- mean(w2))                            

# Create a new dataset using crossing()
(xhypo <- crossing(w1 = w1range, w2 = w2mean))

# Note: if you are more comfortable, you can directly create a data frame
# with the hypothetical values

# Let's predict values using base R function

# Calculate predicted values
simPI.w1 <- predict(ls.result,                 # A model object
                    newdata = xhypo,           # New dataset 
                    interval = "prediction",   # What kind of intervals?
                    level = 0.95)              # What levels of confidence?

simPI.w1 <- simPI.w1 |>
  as_tibble() |>         # Coerce it into a tibble
  bind_cols(w1 = w1range) # Combine hypo w1 with predicted y

# Calculate confidence intervals using predict()
simCI.w1 <- predict(ls.result,
                    newdata = xhypo,
                    interval = "confidence",
                    level = 0.95)

simCI.w1 <- simCI.w1 |>
  as_tibble() |>
  bind_cols(w1 = w1range)


# For plot
simALL.w1 <- bind_rows(
  simPI.w1 |> mutate(type = "PI"),
  simCI.w1 |> mutate(type = "CI")
  )

#print(simALL.w2, n=Inf)
```


```{r, echo=FALSE, fig.align='center'}
# Plot confidence intervals and predictive intervals side by side
# positive slope

ggplot(simALL.w1, aes(x = w1, y = fit, ymax = upr, ymin = lwr)) +
  geom_line() +
  geom_ribbon(alpha = 0.1) +
  labs(y = "Predicted Y", x = "Covariate 1") +
  facet_grid(~ type)+
  theme_classic()
```


The Prediction interval (fundamental uncertainty) has a wide interval and also presents a lot of heteroskedasticity. The model is not good when we have high heteroskedasticity.

### Scenario 2: Vary covariate 2; hold covariate 1 constant

Now, instead of providing a range of values from `w1`, we are going to hold constant this variable and instead we are going to compute the marginal effect of change in `w2` by simulating a range of values.

```{r}
w2range <- seq(from = 0, to = 1, by = 0.05)   # Set up hypothetical values for w1
w1mean <- mean(w1)                            # Set up mean for w2
xhypo <- crossing(w1 = w1mean,               # Create a new dataset 
                  w2 = w2range)                # crossing() is from tidyr package
# Calculate predicted values
simPI.w2 <- predict(ls.result,                 # A model object
                    newdata = xhypo,           # New dataset 
                    interval = "prediction",   # What kind of intervals?
                    level = 0.95)              # What levels of confidence?

simPI.w2 <- simPI.w2 |>
  as_tibble() |>
  bind_cols(w2 = w2range)


# Calculate confidence intervals using predict()
simCI.w2 <- predict(ls.result,
                    newdata = xhypo,
                    interval = "confidence",
                    level = 0.95)

simCI.w2 <- simCI.w2 |>
  as_tibble() |>
  bind_cols(w2 = w2range)


# For plot

simALL.w2 <- bind_rows(
  simPI.w2 |> mutate(type = "PI"),
  simCI.w2 |> mutate(type = "CI")
  )
#print(simALL.w2, n=Inf)
```


```{r, echo=FALSE, fig.align='center'}
# Plot confidence intervals and predictive intervals side by side
# w2 has larger marginal effect
# but still we have large predictive intervals

ggplot(simALL.w2, aes(x = w2, y = fit, ymax = upr, ymin = lwr)) +
  geom_line() +
  geom_ribbon(alpha = 0.1) +
  labs(y = "Predicted Y", x = "Covariate 2") +
  facet_grid(~ type) +
  theme_classic()
```





## Prediction of QoI using Chris's `simcf` package

**Motivation**: Can we use simulation methods to produce the same prediction and confidence intervals? Recall the lecture: we can draw a bunch of $\tilde{\beta}$ from a multivariate normal distribution

### Scenario 1: Vary covariate 1; hold covariate 2 constant

  1. Create a data frame with a set of hypothetical scenarios for covariate 1 while keeping covariate 2 at its mean 

  2. Simulate the predicted values using `MASS` and `simcf`

  3. Plot the results



In order to use `simcf` to generate quantities of interest, the following steps are needed:

  1. Estimate: MLE $\hat{\beta}$ and its variance $\hat{V}(\hat{\beta})$\
$\color{red}{\rightarrow  \texttt{optim(), glm()}}$



  2. Simulate estimation uncertainty from a multivariate normal distribution:\
      Draw $\tilde{\beta} \sim MVN \big[\hat{\beta}, \hat{V}(\hat{\beta})\big]$\
$\color{red}{\rightarrow  \texttt{MASS::mvrnorm()}}$ 



  3. Create hypothetical scenarios of your substantive interest:\
      Choose valuese of X: $X_c$ 
$\color{red}{\rightarrow  \texttt{simcf::cfmake(), cfchange()} \dots}$ 


  4. Calculate expected values:\
      $\tilde{\mu_c} = g(X_c, \tilde{\beta})$ \ 



  5. Simulate fundamental uncertainty:\
      $\tilde{y_c} \sim f(\tilde{\mu_c}, \tilde{\alpha})$\
$\color{red}{\rightarrow  \texttt{simcf::hetnormsimpv()} \dots}$\
      
      or


  6. Compute EVs, First Differences or Relative Risks\
      EV: $\mathbb{E}(y|X_{c1})$\
      $\color{red}{\rightarrow  \texttt{simcf::logitsimev()} \dots}$\
      FD: $\mathbb{E}(y|X_{c2}) - \mathbb{E}(y|X_{c1})$\
      $\color{red}{\rightarrow  \texttt{simcf::logitsimfd()} \dots}$\
      RR: $\frac{\mathbb{E}(y|X_{c2})}{\mathbb{E}(y|X_{c1})}$\
      $\color{red}{\rightarrow  \texttt{simcf::logitsimrr()} \dots}$\




Recall our likelihood function is:


$$
\begin{aligned} \mathrm{P}\left(\mathbf{y} | \boldsymbol{\mu}, \boldsymbol{\sigma}^{2}\right) &=\prod_{i=1}^{n} f_{\mathcal{N}}\left(y_{i} | \mu_{i}, \sigma_{i}^{2}\right) & \tiny[\text{Joint probability}] \\
\mathrm{P}\left(\mathbf{y} | \boldsymbol{\mu}, \boldsymbol{\sigma}^{2}\right) &=\prod_{i=1}^{n}\left(2 \pi \sigma_{i}^{2}\right)^{-1 / 2} \exp \left[\frac{-\left(y_{i}-\mu_{i}\right)^{2}}{2 \sigma_{i}^{2}}\right] & \tiny[\text{Expressed in Normal distribution}]\\ 
\log \mathcal{L}\left(\mathbf{\mu}, \boldsymbol{\sigma}^{2} | \mathbf{y}\right) & \propto-\frac{1}{2} \sum_{i=1}^{n} \log \sigma_{i}^{2}-\frac{1}{2} \sum_{i=1}^{n} \frac{\left(y_{i}-\mathbf{\mu}_{i} \right)^{2}}{\sigma_{i}^{2}} & \tiny[\text{Converted to log likelihood; simplify}] \\ 
\log \mathcal{L}(\boldsymbol{\beta}, \boldsymbol{\gamma} | \mathbf{y}) & \propto-\frac{1}{2} \sum_{i=1}^{n} \log \mathbf{z}_{i} \gamma-\frac{1}{2} \sum_{i=1}^{n} \frac{\left(y_{i}-\mathbf{x}_{i} \boldsymbol{\beta}\right)^{2}}{\exp \left(\mathbf{z}_{i} \gamma\right)} & \tiny[\text{Substitute in systematic components}]
\end{aligned}
$$


```{r}
# Draw parameters from the model predictive distribution
sims <- 10000
simparam <- mvrnorm(n = sims,
                    mu = pe,    # M.N.D. with population mean = point estimates (from ML);
                    Sigma = vc) # population var-covar matrix = vc (from ML)
head(simparam)                  # Inspect

# Separate into the simulated betas and simulated gammas
simbetas <- simparam[ , 1:(ncol(xcovariates) + 1)]
simgammas <- simparam[ , (ncol(simbetas) + 1):ncol(simparam)]

# Specify model formulas
varmodel <- (y ~ w1 + w2)

# Create hypothetical scenarios for covariate 1 
w1range <- seq(from = 0, to = 1, by = 0.2) 


# Use `cfMake` to create a baseline dataset
xhypo <- cfMake(formula = varmodel, 
                data = data, 
                nscen = length(w1range)) 

xhypo


# Use `cfChange` and loop function to loop through each hypothetical x values 
for (i in 1:length(w1range)) {
  xhypo <- cfChange(xscen = xhypo, 
                    covname = "w1", 
                    x = w1range[i], 
                    scen = i) 
}

xhypo

# Repeat the same procedures for z 
zhypo <- cfMake(varmodel, data, nscen = length(w1range))

for (i in 1:length(w1range)) {
  zhypo <- cfChange(zhypo, "w1", x = w1range[i], scen = i) 
  }

zhypo


# Simulate predictive intervals for heteroskedastic linear models
# `hetnormsimpv()` is from simcf package
simRES.w1 <- hetnormsimpv(x = xhypo,
                          b = simbetas,
                          z = zhypo,
                          g = simgammas,
                          ci = 0.95,
                          constant = 1, 
                          varconstant = 1)

simRES.w1
```


```{r, fig.asp = 0.618, out.width = '60%', fig.align = "center"}
simRES.w1 <- simRES.w1 |>
  bind_rows() |>         # Collaspe the list into d.f.
  bind_cols(w1 = w1range) # Combine with hypo. w1 values

simRES.w1                 # Inspect
```


```{r, echo=FALSE, fig.align='center'}
# ggplot2
ggplot(simRES.w1, aes(x = w1, y = pe, ymax = upper, ymin = lower)) +
  geom_line() +
  geom_ribbon(alpha = 0.1) +
  labs(y = "Predicted Y", x = "Covariate 1") + 
  theme_classic()
```


### Scenario 2: Vary covariate 2; hold covariate 1 constant

1. Create a data frame with a set of hypothetical scenarios for covariate 2 while keeping covariate 1 at its mean 

2. Simulate the predicted values using `MASS` and `simcf`

3. Plot the results


```{r}
# Create hypothetical scenarios for w2
w2range <- seq(0, 1, by = 0.05)

xhypo <- cfMake(varmodel, data, nscen = length(w2range))

for (i in 1:length(w2range)) {
  xhypo <- cfChange(xhypo, "w2", x = w2range[i], scen = i) 
}

zhypo <- cfMake(varmodel, data, nscen = length(w2range))

for (i in 1:length(w2range)) {
  zhypo <- cfChange(zhypo, "w2", x = w2range[i], scen = i) 
}


# Simulate the predicted Y's and PI's
simRES.w2 <- hetnormsimpv(xhypo, simbetas,
                          zhypo, simgammas,
                          ci = 0.95,
                          constant = 1, 
                          varconstant = 1)

simRES.w2 <- simRES.w2 |>
  bind_rows() |>
  bind_cols(w2 = w2range)
```


```{r, fig.asp = 0.618, out.width = '70%', fig.align = "center"}
# Plot the predictive intervals for hypothetical w2
# confidence interval increases as we increase w2 
ggplot(simRES.w2, aes(x = w2, y = pe, ymax = upper, ymin = lower)) +
  geom_line() +
  geom_ribbon(alpha = 0.1) +
  labs(y = "Predicted Y", x = "Covariate 2") +
  theme_classic()
```
Can we compare them with the prediction intervals generated by `predict()`?


```{r}
# Combine two dataframes
simPI.w2 <- simPI.w2 |> 
  rename(pe = fit, lower = lwr, upper = upr)

simALL.w2 <- bind_rows(
  simRES.w2 |> mutate(method = "sim"),
  simPI.w2 |> mutate(method = "lm")
)
```


```{r, fig.asp = 0.618, out.width = '70%', fig.align = "center"}
# Plot the predictive intervals for hypothetical w2 (compare lm and sim)
ggplot(simALL.w2, aes(x = w2, y = pe, ymax = upper, ymin = lower)) +
  geom_line() +
  geom_ribbon(alpha = 0.1) +
  labs(y = "Predicted Y", x = "Covariate 2") +
  facet_grid(~ method) +
  theme_classic()
```


