maxLik: An Exercise in Max Likelihood Estimation

---

Using obsevations of income and education in the turnout dataset from the Zelig package, we are going to perform Max Likelihood Estimation to generate the parameters for a normal distribution under which our observed values have the maximum likelihood of occuring.

---

library(Zelig)
library(maxLik)
library(dplyr)
library(knitr)
library(texreg)
library(ggplot2)
data(turnout)

---

Income x Education

Is there a relatinoship between educaiton & income?

---

Let’s explore by plotting the average income for each level of educational attainment.

turnout%>%
  group_by(educate)%>%
  summarize(mean_income=mean(income))%>%
  ggplot()+
  geom_col(aes(x=educate, y=mean_income, fill=mean_income))+
  theme(legend.position="none")

Writing the likelihood function in R.

• sigma tells us the standard deviation of income
• beta1 is the intercept.
• beta 2 is the slope for education
• beta 3 is the slope for age

It appears that education is a statistically singificant predictor of income (p<.001), such that the model estimates a .37 unit increase in income for every additional unit of education. The standard deviation is estimated to remain constant at a value of 2.52.

Age is not statistically significant as a predictor of income in this model.

ols.lf <- function(param) {
  beta <- param[-1] #Regression Coefficients
  sigma <- param[1] #Standard Deviation
  y <- as.vector(turnout$income) #DV
  x <- cbind(1, turnout$educate, turnout$age) #IV
  mu <- x%*%beta #multiply matrices
  sum(dnorm(y, mu, sigma, log = TRUE)) #normal distribution(vector of observations, mean, sd)
  }    
mle_ols <- maxLik(logLik = ols.lf, start = c(sigma = 1, beta1 = 1, beta2 = 1, beta3=1))
summary(mle_ols)

--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 16 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -4690.815 
4  free parameters
Estimates:
       Estimate Std. error t value Pr(> t)    
sigma  2.525576   0.039919  63.268  <2e-16 ***
beta1 -0.446047   0.300583  -1.484   0.138    
beta2  0.371011   0.017493  21.209  <2e-16 ***
beta3 -0.003184   0.003373  -0.944   0.345    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

The values from our maxLik output should be identical to the output from a linear model function in R.

#Checking against linear regression
htmlreg(lm(income~educate+age, data=turnout))

Statistical models

	Model 1
(Intercept)	-0.45
	(0.30)
educate	0.37***
	(0.02)
age	-0.00
	(0.00)
R2	0.20
Adj. R2	0.20
Num. obs.	2000
RMSE	2.53
p < 0.001, p < 0.01, p < 0.05

---

What does this regression look like?

visreg::visreg(lm(income~educate, data=turnout))

---

Income Dispersion x Education

Is there more income inequality among people with higher education?

---

Let’s explore by plotting the standard deviation for income for each level of educational attainment.

turnout%>%
  group_by(educate)%>%
  summarize(income_sd=sd(income))%>%
  ggplot()+
  geom_point(aes(x=educate,y=income_sd))

We can modify the variables that we define in our Likelihood function to engineer our max likelihood estimation procedure to estimate the standard deviation of income rather than the mean value of income.

Writing the likelihood function in R.

• mu tells us the mean income
• theta1 is the intercept.
• theta2 is the slope of education
• theta3 is the slope of age

It appears that education is a statistically singificant predictor of income dispersion (p<.001), such that while holding the mean constant at 3.55, the model estimates a .13 unit increase in the deviation of income for every additional unit of education.

Age is also a statistically singificant predictor of income dispersion according to our model, such that while holding the mean constant at 3.55, the model estimates a .02 unit increase in the deviation of income for every 1 unit increase in age.

ols.lf2 <- function(param) {
  mu <- param[1]
  theta <- param[-1]
  y <- as.vector(turnout$income) #DV
  x <- cbind(1, turnout$educate, turnout$age) #IV
  sigma <- x%*%theta   #multiply matrices
  sum(dnorm(y, mu, sigma, log = TRUE)) #normal distribution(vector of observations, mean, sd)
  }
mle_ols2 <- maxLik(logLik = ols.lf2, start = c(mu = 1, theta1 = 1, theta2 = 1, theta3 = 1), method="BFGS")
summary(mle_ols2)

--------------------------------------------
Maximum Likelihood estimation
BFGS maximization, 150 iterations
Return code 0: successful convergence 
Log-Likelihood: -4843.15 
4  free parameters
Estimates:
       Estimate Std. error t value  Pr(> t)    
mu     3.555011   0.069193  51.378  < 2e-16 ***
theta1 0.362114   0.204550   1.770   0.0767 .  
theta2 0.133349   0.010756  12.398  < 2e-16 ***
theta3 0.017507   0.002852   6.139 8.32e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------