--- title: "Lab 3: Practice Key" subtitle: "Ordinary Least Squares and Maximum Likelihood Estimation" author: "Minji Jeong" output: pdf_document editor_options: chunk_output_type: console --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE) # Good practice to remove all objects from the workspace rm(list = ls()) # Use library() for packages you need, or source() for other R files. library(tidyverse) ``` # Optimization algorithms The purpose of this short coding exercise is to start practicing with the `optim()` function. This function allows us to estimate likelihood parameters from any Probability Density Function (PDF) using Maximum Likelihood Estimation (MLE). Recall that `optim()` has several key arguments: - **par**: The initial guess or starting values for the parameters to be optimized. You must provide as many starting values as there are parameters to estimate. - **fn**: The objective function to be optimized. For OLS, we can provide the Sum of Squared Residuals (SSR). In contrast, for MLE, we should provide a function that characterizes the log-likelihood of the data-generating process (DGP) of interest. - **IMPORTANT**: `optim()` is a minimizer by default. Therefore, when estimating likelihoods, you should multiply the log-likelihood function by $-1$, so you are maximizing the likelihood instead of minimizing it. - **method**: Specifies the optimization algorithm to be used (e.g., BFGS, Nelder-Mead, etc.). - **hessian**: whether the Hessian matrix should be returned. This argument is required to compute the standard errors for ML estimates, but not for OLS. - **Other arguments**: Additional arguments for the objective function can be passed to `optim()`. For example, you might pass the predictor variables ($X_1$, $X_2$) and response ($Y$) as additional arguments. # Data ## Simulation example Let's estimate a linear model with three parameters: $\beta_0$, $\beta_1$, and $\beta_2$, by minimizing the sum of squared residuals (SSR) using the `optim()` function. The model you want to estimate is: $$ Y \sim \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \varepsilon $$ where: - $Y$ is normally distributed, and $\varepsilon \sim \mathcal{N}(0,1)$. - $X_1 \sim \mathcal{N}(10,2)$ and $X_2 \sim \mathcal{N}(5,1)$. - The true estimands are $\beta_0 = 5$, $\beta_1 = 2$, and $\beta_2 = -1.5$. ```{r, fig.align='center', fig.width=5, fig.height=3.5} # Set simulation parameters set.seed(123) # for reproducibility n <- 100 # Sample size b0 <- 5 # True intercept b1 <- 2 # True coefficient for X1 b2 <- -1.5 # True coefficient for X2 # Simulate random variables ## Note: the square root in the sd argument of rnorm() ## Predictors x1 <- rnorm(n, mean = 10, sd = sqrt(2)) x2 <- rnorm(n, mean = 5, sd = sqrt(1)) ## Residuals epsilon <- rnorm(n, mean = 0, sd = sqrt(1)) # DGP of the outcome variable Y <- b0 + b1 * x1 + b2 * x2 + epsilon # Visualize the distribution of Y hist(Y, main = "Histogram of Y", xlab = "Y") # model <- lm(Y ~ x1 + x2) # coef(model) ``` ## Salaries dataset **Description:** This dataset contains information on the 2008-09 nine-month academic salaries for Assistant Professors, Associate Professors, and Professors at a U.S. college. The data were collected as part of the college’s ongoing efforts to monitor salary differences between male and female faculty members. Below, I am cleaning the data for you. ```{r} # install.packages("carData") # Load the data data("Salaries", package = "carData") # Prepare data df <- Salaries |> # optional: turn it into a tibble as_tibble() |> # re-scale salary mutate(salary=salary/1000) |> # select predictors for analysis select(salary,sex,yrs.since.phd) ``` # Ordinary Least Squares In this section, we will use the built-in dataset *Salaries* from the `car` / `carData` package, that we used in the previous lab code practice. However, before we will fit a toy example with simulated data. The goal is to fit least squares estimates using `optim()` instead of `lm()`. ## 1. Ordinary Least Squares estimation using OLS Create a function called `SSR` that 1- Maps a parameter argument `params=` with each parameter to be estimated ($\beta_0, \beta_1$, and $\beta_2$). 2- Provides an objective function to estimate the linear projection of $Y$. That is $\hat{Y} = \hat{\beta_0} + \hat{\beta_1} x1 + \hat{\beta_2} X_2$. Do not forget to include an intercept. 3- An objective function to estimate the Sum of Squared Residuals (SSR) which will return an output. This output is the one `optim()` will optimize until finding the global minimum. Once you have created `SSR()`, use `optim()` to estimate $\hat{\beta_0}$, $\hat{\beta_1}$, and $\hat{\beta_2}$. Double-check your results with `lm()`. ```{r} # Define the residual sum of squares (SSR) function SSR <- function(param, y, x1, x2) { # 1- Map a parameter arguments b0 <- param[1] # Intercept b1 <- param[2] # Coefficient for X1 b2 <- param[3] # Coefficient for X2 # 2- Linear projection model of y (aka model to predict y) # rep(1,n) is constant y_hat <- b0 * rep(1, n) + b1 * x1 + b2 * x2 # Objective funtion to estimate SSR SSR <- sum((y - y_hat) ^ 2) return(SSR) } ``` ### 2. Use `SSR()` and `optim()` to estimate the following model: $Y_i = \alpha + \beta_1 X_{1i} + \beta_2 X_{2i} + \epsilon$ using the *Salaries* dataset. Compare `optim()` results with `lm()`. where: - **salary** ($Y$): Outcome variable representing the nine-month salary in dollars. - **sex** ($x1$): Binary indicator for the respondent’s reported sex (Female or Male). - **yrs.since.phd** ($X_2$): Count variable indicating the number of years since the respondent received their PhD. - Note that input the variable *sex* as **numeric** class. ```{r} # Set vectors/matrix of covariates y <- df$salary x1 <- as.numeric(df$sex) x2 <- df$yrs.since.phd # adjust n for the intercept n <- length(y) # Initial parameter guesses for b0, b1, and b2 # optim requrires initial parameters stval <- c(b0 = 0, b1 = 0, b2 = 0) # Use optim to minimize SSR optim_ols <- optim(par = stval, fn = SSR, method = "BFGS", y = y, x1 = x1, x2 = x2) # Compare with lm() lm_ols <- lm(y ~ x1 + x2) # Display results: cat("Minimum SSR:", optim_ols$value, "\n") cat("Estimated parameters with optim (b0, b1, b2):", coef(lm_ols), "\n") cat("Estimated parameters with lm (b0, b1, b2):", optim_ols$par, "\n") ``` ```{r} ### Alternative version # SSR <- function(param, y, x) { # # # 1- Map a parameter arguments, do not forget to include an intercept! # x <- as.matrix(x) # x (covariates) as a matrix for mean # # os <- rep(1, nrow(x)) # Set the intercept as 1 (constant) # # x <- cbind(os, x) # Add intercept to covariates x # # # 2- Define vector of parameters (b) # b <- param[1 : ncol(x)] # Parameters for x # # # 3- Define systematic components # # Systematic components for mean # y_hat <- x %*% b # # # Objective function to estimate SSR # SSR <- sum((y - y_hat) ^ 2) # # return(SSR) # } # # y <- df$salary # x1 <- as.numeric(df$sex) # x2 <- df$yrs.since.phd # # x <- cbind(x1,x2) # stval <- c(0,0,0) # # # optim_ols <- optim(par = stval, # fn = SSR, # method = "BFGS", # y = y, # x = x) # optim_ols$par ``` # Maximum Likelihood Estimation From `lab2_practice.Rmd`, we know that model 3 (the one estimated in the previous section) violates the assumption of homoskedasticity. Therefore, in this section, we will re-estimate the model using a heteroskedastic normal model via Maximum Likelihood Estimation (MLE). ## 3. Estimating a Heteroskedastic Normal Model Using MLE To estimate the heteroskedastic normal model, you will first need to create a function that connects the parameters to the objective function for estimation. - You can modify the function we covered during the lab session or create your own. This time, ensure that the covariates/predictors are set into matrices for each parameter of the normal distribution that we want to model ($\beta$ and $\gamma$). ```{r} # Set vectors/matrix of covariates ## vector of outcomes y <- df$salary ## matrix of covariates for mean and variance parameters xcov <- cbind(as.numeric(df$sex),df$yrs.since.phd) zcov <- cbind(as.numeric(df$sex),df$yrs.since.phd) llk.hetnorm <- function(param, y, x, z) { # 1- Map a parameter arguments, do not forget to include an intercept! x <- as.matrix(x) # x (some covariates) as a matrix for mean z <- as.matrix(z) # z (some covariates) as a matrix for variance 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 # 2- Define vector of parameters (b and g) b <- param[1 : ncol(x)] # Parameters for x g <- param[(ncol(x) + 1) : (ncol(x) + ncol(z))] # Parameters for z # 3- Define systematic components # Systematic components for mean xb <- x %*% b # Systematic components for variance s2 <- exp(z %*% g) # 4- Find Objective function sum(0.5 * (log(s2) + (y - xb)^2 / s2)) # Likelihood we want to maximize } # Initial guess, how many parameters we are going to estimate? stval <- c(0, 0, 0, 0, 0, 0) # inverse the "value" (multiply -1) # Run ML, and get the output we need opmtim_mle <- optim(par = stval, # Initial guesses fn= llk.hetnorm, # Likelihood function method = "BFGS", # Gradient method hessian = TRUE, # Return Hessian matrix y = y, # Outcome variable x = xcov, # Covariates x (w/o constant) z = zcov # Covariates z (w/o constant) ) # Extract point estimates and standard errors pe <- opmtim_mle$par # Point estimates vc <- solve(opmtim_mle$hessian) # Var-cov matrix (for computing s.e.) se <- sqrt(diag(vc)) # Standard errors (s.e.) # Compare with mle and lm() results #options(tibble.width = Inf) tibble( "Variables" = c("(Intercept)", "x1", "x2"), "betas pe" = pe[1:3], "betas se" = se[1:3], "gamma pe" = pe[4:6], "gamma se" = se[4:6] ) round(coef(summary(lm_ols))[, 1:2],3) # MLE-heteroskedastic normal is better estimator when we have heteroscedasticity??? ``` # Optional: Monte Carlo Experiment ```{r} ### Monte Carlo Experiment OLS vs. MLE:--------------------------- set.seed(123456) # sample size n <- 600 # systematic component b0 <- 0.2 # True value for the intercept b1 <- 0.5 # True value for the slope X <- rnorm(n, 0, 2) # independent variable # stochastic component e <- rnorm(n, 0, 1) # error term # Heteroskedasticity parameter gamma <- 0.35 # The DGP: Heteroskedasticity Y <- b0 + b1*X + rnorm(n, 0, exp(X*gamma)) # Simulate repeated samples sims <- 500 ## Simulation for heteroskedasticity par.est.1 <- tibble( b0_hat = rep(NA,sims), b1_hat = rep(NA,sims), model = "OLS" ) par.est.2 <- tibble( b0_hat = rep(NA,sims), b1_hat = rep(NA,sims), model = "MLE" ) ## Likelihood function llk.hetnormlin <- 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 } ## starting values for searching algorithm stval <- c(0, 0, 0, 0) for (i in 1:sims) { # Start the loop # The DGP: Heteroskedasticity Y <- b0 + b1*X + rnorm(n, 0, exp(X*gamma)) # Estimate OLS model out.ols <- lm(Y ~ X) # Save OLS estimates par.est.1[i, 1] <- out.ols$coef[1] par.est.1[i, 2] <- out.ols$coef[2] # Estimate MLE-Normal Heteroskedastic out.mle <- optim(par = stval, # Initial guesses fn= llk.hetnormlin, # Likelihood function method = "BFGS", # Gradient method hessian = TRUE, # Return Hessian matrix y = Y, # Outcome variable x = X, # Covariates x (w/o constant) z = X # Covariates z (w/o constant) ) # Save MLE estimates par.est.2[i, 1] <- out.mle$par[1] par.est.2[i, 2] <- out.mle$par[2] } # visualize rbind(par.est.1, par.est.2) |> pivot_longer(cols = 1:2, names_to = "par", values_to = "estimate")|> ggplot(aes(x=estimate, group=model, color=model, fill=model)) + geom_density(alpha=0.2) + scale_y_continuous(breaks = seq(0,100,5)) + geom_rug() + scale_color_manual(values = c("OLS" = "#D24306", "MLE" = "#0D5088")) + scale_fill_manual(values = c("OLS" = "#D24306", "MLE" = "#0D5088")) + geom_vline(aes(xintercept = ifelse(par == "b0_hat", 0.2, 0.5)), linetype = "dashed", color = "black") + labs(title = "Heteroskedasticity: OLS vs MLE-Heteroskedastic Normal", subtitle = "Simulation: 500 regression estimate samples (n=600)", x = "Estimates", y = "Density", fill = "", color = "", caption = "Note: Each facet's scales are set free.") + theme_minimal() + facet_wrap(~par, scales = "free") + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), legend.position = "bottom", text = element_text(size = 14), legend.text = element_text(size = 10), plot.caption = element_text(size = 8, face = "italic")) ## MLE estimate is very accurate -> why we use MLE ```