--- 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 (QoI)** 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. QoI 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")) ``` 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). #### 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",2*sims), 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: 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, # which creates heteroskedasticity mu <- x %*% beta # Create systemtic component for the mean sigma2 <- exp(z %*% gamma) # Create systematic component for the variance # 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") # Linear regression fro comparison 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(par = stval, # Initial guesses fn = 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 # Negative sign # bc optim is set as minimizer 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 # only w1 changes (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} w1mean <- mean(w1) # Set up mean for w1 w2range <- seq(from = 0, to = 1, by = 0.05) # Set up hypothetical values 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() ```