Comparing results reported on maximizing log likelihood function and MLE results

logLikFun <- function(param) {
  mu <- param[1]
  sigma <- param[2]
  sum(dnorm(x, mean = mu, sd = sigma, log = TRUE))
}

Maximizing the Log Likelihood Function

The Effects of Education on Income From the “Turnout” dataset

library(Zelig)
library(maxLik)
data(turnout)
mle_ols <- maxLik(logLik = ols.lf, start = c(sigma = 1, beta1 = 1, beta2 = 1))
summary(mle_ols)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 12 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -4691.256 
3  free parameters
Estimates:
      Estimate Std. error t value Pr(> t)    
sigma  2.52613    0.03989  63.326 < 2e-16 ***
beta1 -0.65207    0.20827  -3.131 0.00174 ** 
beta2  0.37613    0.01663  22.612 < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

In the above results from slide 4.14 Education is our independent variable and income is our dependent variable. In order to understand the results shown above we must first understand our three parameters that are being estimated which are sigma, beta1 and beta2. In this model Sigma is our standard deviation which shows that there is variation (2.526). Beta1 represents the y-intercept of what would be the best fit line. Beta1 means that people who have no education would have a average income of -.065207. Beta2 represents the slope which shows the effect of education has on a voters income. We see that one unit increase in education has a .38 increase on a person’s income. The results show that there is a positive correlation between education and income and this correlation is statistically significant at a 95% confidence level. Below we check the results using a regression.

Checking and Replicating the Results Above Using a Regression Model

summary(lm(income ~ educate, data = turnout))

Call:
lm(formula = income ~ educate, data = turnout)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2028 -1.7363 -0.4273  1.3150 11.0632 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.65207    0.21016  -3.103  0.00194 ** 
educate      0.37613    0.01677  22.422  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.527 on 1998 degrees of freedom
Multiple R-squared:  0.201, Adjusted R-squared:  0.2006 
F-statistic: 502.8 on 1 and 1998 DF,  p-value: < 2.2e-16

As we can see in this regression the intercept and educate estimate are the same as the beta1 and beta2 from the maxLik function we ran above. Shown below the coefficients is the residual standard error of 2.527 which is shown as the sigma in the MaxLik model.

Creating the function

ols.lf2 <- function(param) {
  mu <- param[1]
  theta <- param[-1]
  y <- as.vector(turnout$income)
  x <- cbind(1, turnout$educate)
  sigma <- x%*%theta
  sum(dnorm(y, mu, sigma, log = TRUE))
}    

Maximum Likelihood Estimation Results

library(maxLik)
mle_ols2 <- maxLik(logLik = ols.lf2, start = c(mu = 1, theta1 = 1, theta2 = 1))
summary(mle_ols2)
--------------------------------------------
Maximum Likelihood estimation
Newton-Raphson maximisation, 9 iterations
Return code 2: successive function values within tolerance limit
Log-Likelihood: -4861.964 
3  free parameters
Estimates:
       Estimate Std. error t value Pr(> t)    
mu     3.516764   0.070320   50.01  <2e-16 ***
theta1 1.461011   0.106745   13.69  <2e-16 ***
theta2 0.109081   0.009185   11.88  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
--------------------------------------------

Here we have a very similar model to the one above except for some slight differences. In this model you can see that the three parameters are now mu, theta1 and theta 2 as opposed to sigma and beta. The new parameters shown signify that you are now looking for standard deviation across different education levels. Our variabes remain the same in this model however we are no longer looking for the mean income levels but rather the effect education has on the variation of income. Now getting to the results of this model(slide 4.18), once again theta1 represents the y-intercept but this time it means that the voters who have no education have a 1.461 variation on income. Theta2 which is once again the slope however in this model it represents the effect education has on the SD on income. Meaning for every one unit increase in education this is an increase of .11 on the variation of income. Both theta1 and theta2 are statistically significant at a 99% confidence level.

Effect of Age on Income

  1. In both models, if you introduce a second independent variable age, what do you think the relationship between age and income might be?
summary(lm(income ~ educate + age, data = turnout))

Call:
lm(formula = income ~ educate + age, data = turnout)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2128 -1.7471 -0.4217  1.3042 11.1256 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.446084   0.303955  -1.468    0.142    
educate      0.371013   0.017641  21.031   <2e-16 ***
age         -0.003183   0.003394  -0.938    0.348    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.527 on 1997 degrees of freedom
Multiple R-squared:  0.2014,    Adjusted R-squared:  0.2006 
F-statistic: 251.8 on 2 and 1997 DF,  p-value: < 2.2e-16

Extending the Maximum Likelihood Estimation to include age as a second independent variable.

ols.lf3 <- function(param) {
  mu <- param[1]
  theta <- param[-1]
  y <- as.vector(turnout$income)
  x <- cbind(1, turnout$educate, turnout$age)
  sigma <- x%*%theta
  sum(dnorm(y, mu, sigma, log = TRUE))
} 

New MLE results

mle_ols3 <- maxLik(logLik = ols.lf3, start = c(mu = 1, theta1 = 1, theta2 = 1, theta3 = 1), method = "BFGS")
summary(mle_ols3)
--------------------------------------------
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
--------------------------------------------
  1. In both models, if you introduce a second independent variable age, what do you think the relationship between age and income might be?

So as you can see we now have produced the two same models discussed above but have introduced age as a second independent variable. In our first model we see that the RSE or sigma remained almost the same suggesting that age might have or explain any difference in education. When you look at the coefficient for age you can actually see that the estimation given is slightly negative meaning that for every 1 unit increase in age your income will descrease slightly. However this finding is not statistically significant. Education did however remain statistically significant after adding age as an IV.

For our last model we see that both age and education impact the variance in average income. Before we introduced age as an IV the MU was 2.52 however after we introduced age as an IV the MU increased to 3.56 You can also see that for every year of education(theta2), income increases by.13 and every year of age(theta3), income increases by .02. These results are both signnificant at a 99% confidence level.

LS0tCnRpdGxlOiAiKkhvbWV3b3JrM19tYXhMaWsqIgphdXRob3I6ICJSb2JlcnQgUGVyZXoiCmRhdGU6ICJGZWJydWFyeSAyNCwgMjAxOCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKIyMjI0NvbXBhcmluZyByZXN1bHRzIHJlcG9ydGVkIG9uIG1heGltaXppbmcgbG9nIGxpa2VsaWhvb2QgZnVuY3Rpb24gYW5kIE1MRSByZXN1bHRzCmBgYHtyfQpsb2dMaWtGdW4gPC0gZnVuY3Rpb24ocGFyYW0pIHsKICBtdSA8LSBwYXJhbVsxXQogIHNpZ21hIDwtIHBhcmFtWzJdCiAgc3VtKGRub3JtKHgsIG1lYW4gPSBtdSwgc2QgPSBzaWdtYSwgbG9nID0gVFJVRSkpCn0KYGBgCgojI01heGltaXppbmcgdGhlIExvZyBMaWtlbGlob29kIEZ1bmN0aW9uIAojIyMjVGhlIEVmZmVjdHMgb2YgRWR1Y2F0aW9uIG9uIEluY29tZSBGcm9tIHRoZSAiVHVybm91dCIgZGF0YXNldCAKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmxpYnJhcnkoWmVsaWcpCmxpYnJhcnkobWF4TGlrKQpkYXRhKHR1cm5vdXQpCm1sZV9vbHMgPC0gbWF4TGlrKGxvZ0xpayA9IG9scy5sZiwgc3RhcnQgPSBjKHNpZ21hID0gMSwgYmV0YTEgPSAxLCBiZXRhMiA9IDEpKQpzdW1tYXJ5KG1sZV9vbHMpCmBgYApJbiB0aGUgYWJvdmUgcmVzdWx0cyBmcm9tIHNsaWRlIDQuMTQgRWR1Y2F0aW9uIGlzIG91ciBpbmRlcGVuZGVudCB2YXJpYWJsZSBhbmQgaW5jb21lIGlzIG91ciBkZXBlbmRlbnQgdmFyaWFibGUuIEluIG9yZGVyIHRvIHVuZGVyc3RhbmQgdGhlIHJlc3VsdHMgc2hvd24gYWJvdmUgd2UgbXVzdCBmaXJzdCB1bmRlcnN0YW5kIG91ciB0aHJlZSBwYXJhbWV0ZXJzIHRoYXQgYXJlIGJlaW5nIGVzdGltYXRlZCB3aGljaCBhcmUgc2lnbWEsIGJldGExIGFuZCBiZXRhMi4gSW4gdGhpcyBtb2RlbCBTaWdtYSBpcyBvdXIgc3RhbmRhcmQgZGV2aWF0aW9uIHdoaWNoIHNob3dzIHRoYXQgdGhlcmUgaXMgdmFyaWF0aW9uICgyLjUyNikuIEJldGExIHJlcHJlc2VudHMgdGhlIHktaW50ZXJjZXB0IG9mIHdoYXQgd291bGQgYmUgdGhlIGJlc3QgZml0IGxpbmUuIEJldGExIG1lYW5zIHRoYXQgcGVvcGxlIHdobyBoYXZlIG5vIGVkdWNhdGlvbiB3b3VsZCBoYXZlIGEgYXZlcmFnZSBpbmNvbWUgb2YgLS4wNjUyMDcuIEJldGEyIHJlcHJlc2VudHMgdGhlIHNsb3BlIHdoaWNoIHNob3dzIHRoZSBlZmZlY3Qgb2YgZWR1Y2F0aW9uIGhhcyBvbiBhIHZvdGVycyBpbmNvbWUuIFdlIHNlZSB0aGF0IG9uZSB1bml0IGluY3JlYXNlIGluIGVkdWNhdGlvbiBoYXMgYSAuMzggaW5jcmVhc2Ugb24gYSBwZXJzb24ncyBpbmNvbWUuIFRoZSByZXN1bHRzIHNob3cgdGhhdCB0aGVyZSBpcyBhIHBvc2l0aXZlIGNvcnJlbGF0aW9uIGJldHdlZW4gZWR1Y2F0aW9uIGFuZCBpbmNvbWUgYW5kIHRoaXMgY29ycmVsYXRpb24gaXMgc3RhdGlzdGljYWxseSBzaWduaWZpY2FudCBhdCBhIDk1JSBjb25maWRlbmNlIGxldmVsLiBCZWxvdyB3ZSBjaGVjayB0aGUgcmVzdWx0cyB1c2luZyBhIHJlZ3Jlc3Npb24uIAoKIyNDaGVja2luZyBhbmQgUmVwbGljYXRpbmcgdGhlIFJlc3VsdHMgQWJvdmUgVXNpbmcgYSBSZWdyZXNzaW9uIE1vZGVsCmBgYHtyfQpzdW1tYXJ5KGxtKGluY29tZSB+IGVkdWNhdGUsIGRhdGEgPSB0dXJub3V0KSkKYGBgCkFzIHdlIGNhbiBzZWUgaW4gdGhpcyByZWdyZXNzaW9uIHRoZSBpbnRlcmNlcHQgYW5kIGVkdWNhdGUgZXN0aW1hdGUgYXJlIHRoZSBzYW1lIGFzIHRoZSBiZXRhMSBhbmQgYmV0YTIgZnJvbSB0aGUgbWF4TGlrIGZ1bmN0aW9uIHdlIHJhbiBhYm92ZS4gU2hvd24gYmVsb3cgdGhlIGNvZWZmaWNpZW50cyBpcyB0aGUgcmVzaWR1YWwgc3RhbmRhcmQgZXJyb3Igb2YgMi41Mjcgd2hpY2ggaXMgc2hvd24gYXMgdGhlIHNpZ21hIGluIHRoZSBNYXhMaWsgbW9kZWwuIAoKIyNDcmVhdGluZyB0aGUgZnVuY3Rpb24gCmBgYHtyfQpvbHMubGYyIDwtIGZ1bmN0aW9uKHBhcmFtKSB7CiAgbXUgPC0gcGFyYW1bMV0KICB0aGV0YSA8LSBwYXJhbVstMV0KICB5IDwtIGFzLnZlY3Rvcih0dXJub3V0JGluY29tZSkKICB4IDwtIGNiaW5kKDEsIHR1cm5vdXQkZWR1Y2F0ZSkKICBzaWdtYSA8LSB4JSoldGhldGEKICBzdW0oZG5vcm0oeSwgbXUsIHNpZ21hLCBsb2cgPSBUUlVFKSkKfSAgICAKYGBgCgojI01heGltdW0gTGlrZWxpaG9vZCBFc3RpbWF0aW9uIFJlc3VsdHMKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9CmxpYnJhcnkobWF4TGlrKQptbGVfb2xzMiA8LSBtYXhMaWsobG9nTGlrID0gb2xzLmxmMiwgc3RhcnQgPSBjKG11ID0gMSwgdGhldGExID0gMSwgdGhldGEyID0gMSkpCnN1bW1hcnkobWxlX29sczIpCmBgYApIZXJlIHdlIGhhdmUgYSB2ZXJ5IHNpbWlsYXIgbW9kZWwgdG8gdGhlIG9uZSBhYm92ZSBleGNlcHQgZm9yIHNvbWUgc2xpZ2h0IGRpZmZlcmVuY2VzLiBJbiB0aGlzIG1vZGVsIHlvdSBjYW4gc2VlIHRoYXQgdGhlIHRocmVlIHBhcmFtZXRlcnMgYXJlIG5vdyBtdSwgdGhldGExIGFuZCB0aGV0YSAyIGFzIG9wcG9zZWQgdG8gc2lnbWEgYW5kIGJldGEuIFRoZSBuZXcgcGFyYW1ldGVycyBzaG93biBzaWduaWZ5IHRoYXQgeW91IGFyZSBub3cgbG9va2luZyBmb3Igc3RhbmRhcmQgZGV2aWF0aW9uIGFjcm9zcyBkaWZmZXJlbnQgZWR1Y2F0aW9uIGxldmVscy4gT3VyIHZhcmlhYmVzIHJlbWFpbiB0aGUgc2FtZSBpbiB0aGlzIG1vZGVsIGhvd2V2ZXIgd2UgYXJlIG5vIGxvbmdlciBsb29raW5nIGZvciB0aGUgbWVhbiBpbmNvbWUgbGV2ZWxzIGJ1dCByYXRoZXIgdGhlIGVmZmVjdCBlZHVjYXRpb24gaGFzIG9uIHRoZSB2YXJpYXRpb24gb2YgaW5jb21lLiBOb3cgZ2V0dGluZyB0byB0aGUgcmVzdWx0cyBvZiB0aGlzIG1vZGVsKHNsaWRlIDQuMTgpLCBvbmNlIGFnYWluIHRoZXRhMSByZXByZXNlbnRzIHRoZSB5LWludGVyY2VwdCBidXQgdGhpcyB0aW1lIGl0IG1lYW5zIHRoYXQgdGhlIHZvdGVycyB3aG8gaGF2ZSBubyBlZHVjYXRpb24gaGF2ZSBhIDEuNDYxIHZhcmlhdGlvbiBvbiBpbmNvbWUuIFRoZXRhMiB3aGljaCBpcyBvbmNlIGFnYWluIHRoZSBzbG9wZSBob3dldmVyIGluIHRoaXMgbW9kZWwgaXQgcmVwcmVzZW50cyB0aGUgZWZmZWN0IGVkdWNhdGlvbiBoYXMgb24gdGhlIFNEIG9uIGluY29tZS4gTWVhbmluZyBmb3IgZXZlcnkgb25lIHVuaXQgaW5jcmVhc2UgaW4gZWR1Y2F0aW9uIHRoaXMgaXMgYW4gaW5jcmVhc2Ugb2YgLjExIG9uIHRoZSB2YXJpYXRpb24gb2YgaW5jb21lLiBCb3RoIHRoZXRhMSBhbmQgdGhldGEyIGFyZSBzdGF0aXN0aWNhbGx5IHNpZ25pZmljYW50IGF0IGEgOTklIGNvbmZpZGVuY2UgbGV2ZWwuIAoKIyNFZmZlY3Qgb2YgQWdlIG9uIEluY29tZSAKMikgSW4gYm90aCBtb2RlbHMsIGlmIHlvdSBpbnRyb2R1Y2UgYSBzZWNvbmQgaW5kZXBlbmRlbnQgdmFyaWFibGUgYWdlLCB3aGF0IGRvIHlvdSB0aGluayB0aGUgcmVsYXRpb25zaGlwIGJldHdlZW4gYWdlIGFuZCBpbmNvbWUgbWlnaHQgYmU/CmBgYHtyfQpzdW1tYXJ5KGxtKGluY29tZSB+IGVkdWNhdGUgKyBhZ2UsIGRhdGEgPSB0dXJub3V0KSkKYGBgCgojIyNFeHRlbmRpbmcgdGhlIE1heGltdW0gTGlrZWxpaG9vZCBFc3RpbWF0aW9uIHRvIGluY2x1ZGUgYWdlIGFzIGEgc2Vjb25kIGluZGVwZW5kZW50IHZhcmlhYmxlLiAKYGBge3J9Cm9scy5sZjMgPC0gZnVuY3Rpb24ocGFyYW0pIHsKICBtdSA8LSBwYXJhbVsxXQogIHRoZXRhIDwtIHBhcmFtWy0xXQogIHkgPC0gYXMudmVjdG9yKHR1cm5vdXQkaW5jb21lKQogIHggPC0gY2JpbmQoMSwgdHVybm91dCRlZHVjYXRlLCB0dXJub3V0JGFnZSkKICBzaWdtYSA8LSB4JSoldGhldGEKICBzdW0oZG5vcm0oeSwgbXUsIHNpZ21hLCBsb2cgPSBUUlVFKSkKfSAKYGBgCgojIyNOZXcgTUxFIHJlc3VsdHMKYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9Cm1sZV9vbHMzIDwtIG1heExpayhsb2dMaWsgPSBvbHMubGYzLCBzdGFydCA9IGMobXUgPSAxLCB0aGV0YTEgPSAxLCB0aGV0YTIgPSAxLCB0aGV0YTMgPSAxKSwgbWV0aG9kID0gIkJGR1MiKQpzdW1tYXJ5KG1sZV9vbHMzKQpgYGAKMikgSW4gYm90aCBtb2RlbHMsIGlmIHlvdSBpbnRyb2R1Y2UgYSBzZWNvbmQgaW5kZXBlbmRlbnQgdmFyaWFibGUgYWdlLCB3aGF0IGRvIHlvdSB0aGluayB0aGUgcmVsYXRpb25zaGlwIGJldHdlZW4gYWdlIGFuZCBpbmNvbWUgbWlnaHQgYmU/CgpTbyBhcyB5b3UgY2FuIHNlZSB3ZSBub3cgaGF2ZSBwcm9kdWNlZCB0aGUgdHdvIHNhbWUgbW9kZWxzIGRpc2N1c3NlZCBhYm92ZSBidXQgaGF2ZSBpbnRyb2R1Y2VkIGFnZSBhcyBhIHNlY29uZCBpbmRlcGVuZGVudCB2YXJpYWJsZS4gSW4gb3VyIGZpcnN0IG1vZGVsIHdlIHNlZSB0aGF0IHRoZSBSU0Ugb3Igc2lnbWEgcmVtYWluZWQgYWxtb3N0IHRoZSBzYW1lIHN1Z2dlc3RpbmcgdGhhdCBhZ2UgbWlnaHQgaGF2ZSBvciBleHBsYWluIGFueSBkaWZmZXJlbmNlIGluIGVkdWNhdGlvbi4gV2hlbiB5b3UgbG9vayBhdCB0aGUgY29lZmZpY2llbnQgZm9yIGFnZSB5b3UgY2FuIGFjdHVhbGx5IHNlZSB0aGF0IHRoZSBlc3RpbWF0aW9uIGdpdmVuIGlzIHNsaWdodGx5IG5lZ2F0aXZlIG1lYW5pbmcgdGhhdCBmb3IgZXZlcnkgMSB1bml0IGluY3JlYXNlIGluIGFnZSB5b3VyIGluY29tZSB3aWxsIGRlc2NyZWFzZSBzbGlnaHRseS4gSG93ZXZlciB0aGlzIGZpbmRpbmcgaXMgbm90IHN0YXRpc3RpY2FsbHkgc2lnbmlmaWNhbnQuIEVkdWNhdGlvbiBkaWQgaG93ZXZlciByZW1haW4gc3RhdGlzdGljYWxseSBzaWduaWZpY2FudCBhZnRlciBhZGRpbmcgYWdlIGFzIGFuIElWLiAKCkZvciBvdXIgbGFzdCBtb2RlbCAgd2Ugc2VlIHRoYXQgYm90aCBhZ2UgYW5kIGVkdWNhdGlvbiBpbXBhY3QgdGhlIHZhcmlhbmNlIGluIGF2ZXJhZ2UgaW5jb21lLiBCZWZvcmUgd2UgaW50cm9kdWNlZCBhZ2UgYXMgYW4gSVYgdGhlIE1VIHdhcyAyLjUyIGhvd2V2ZXIgYWZ0ZXIgd2UgaW50cm9kdWNlZCBhZ2UgYXMgYW4gSVYgdGhlIE1VIGluY3JlYXNlZCB0byAzLjU2IFlvdSBjYW4gYWxzbyBzZWUgdGhhdCBmb3IgZXZlcnkgeWVhciBvZiBlZHVjYXRpb24odGhldGEyKSwgaW5jb21lIGluY3JlYXNlcyBieS4xMyBhbmQgZXZlcnkgeWVhciBvZiBhZ2UodGhldGEzKSwgaW5jb21lIGluY3JlYXNlcyBieSAuMDIuIFRoZXNlIHJlc3VsdHMgYXJlIGJvdGggc2lnbm5pZmljYW50IGF0IGEgOTklIGNvbmZpZGVuY2UgbGV2ZWwuIAoKCgoKCgoKCgoKCgoKCgoKCgoKCg==