Modeling in R with Log Likelihood Function
This article is originally published at https://statcompute.wordpress.com
Similar to NLMIXED procedure in SAS, optim() in R provides the functionality to estimate a model by specifying the log likelihood function explicitly. Below is a demo showing how to estimate a Poisson model by optim() and its comparison with glm() result.
> df <- read.csv('credit_count.csv')
> # ESTIMATE A POISSON MODEL WITH GLM()
> mdl <- glm(MAJORDRG ~ AGE + ACADMOS + ADEPCNT + MINORDRG, family = poisson, data = df)
> summary(mdl)
Call:
glm(formula = MAJORDRG ~ AGE + ACADMOS + ADEPCNT + MINORDRG,
family = poisson, data = df)
Deviance Residuals:
Min 1Q Median 3Q Max
-9.9940 -0.8907 -0.8079 -0.7633 11.6866
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.3813012 0.0450281 -30.676 < 2e-16 ***
AGE 0.0056126 0.0013616 4.122 3.76e-05 ***
ACADMOS 0.0013437 0.0001975 6.803 1.03e-11 ***
ADEPCNT 0.0803056 0.0093378 8.600 < 2e-16 ***
MINORDRG 0.4499422 0.0068969 65.238 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 24954 on 13443 degrees of freedom
Residual deviance: 22026 on 13439 degrees of freedom
AIC: 28520
Number of Fisher Scoring iterations: 6
> # ESTIMATE A POISSON MODEL WITH OPTIM()
> log.like <- function(par) {
+ xb <- par[1] + par[2] * df$AGE + par[3] * df$ACADMOS + par[4] * df$ADEPCNT + par[5] * df$MINORDRG
+ mu <- exp(xb)
+ ll <- sum(log(exp(-mu) * (mu ^ df$MAJORDRG) / factorial(df$MAJORDRG)))
+ return(-ll)
+ }
> result <- optim(c(0, 0, 0, 0, 0), log.like, hessian = TRUE, method = "BFGS")
> stder <- sqrt(diag(solve(result$hessian)))
> estimate <- data.frame(beta = result$par, stder = stder, z_values = result$par / stder)
> print(estimate)
beta stder z_values
1 -1.380911081 0.0450398804 -30.659741
2 0.005656423 0.0013611828 4.155520
3 0.001298029 0.0001956315 6.635072
4 0.080171673 0.0093427325 8.581180
5 0.450468859 0.0068922289 65.358952
Thanks for visiting r-craft.org
This article is originally published at https://statcompute.wordpress.com
Please visit source website for post related comments.
