Modeling in R with Log Likelihood Function
[This article was first published on Yet Another Blog in Statistical Computing » S+/R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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
To leave a comment for the author, please follow the link and comment on their blog: Yet Another Blog in Statistical Computing » S+/R.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.