[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.
In my post on 06/05/2013 (http://statcompute.wordpress.com/2013/06/05/estimating-composite-models-for-count-outcomes-with-fmm-procedure), I’ve shown how to estimate finite mixture models, e.g. zero-inflated Poisson and 2-class finite mixture Poisson models, with FMM and NLMIXED procedure in SAS. Today, I am going to demonstrate how to achieve the same results with flexmix package in R.
R Code
library(flexmix) data <- read.table('../data/credit_count.txt', header = TRUE, sep = ',') set.seed(2013) class <- FLXPmultinom(~ AGE + ACADMOS + MINORDRG + LOGSPEND) formula <- "MAJORDRG ~ AGE + ACADMOS + MINORDRG + LOGSPEND" control <- list(verbose = 10, iter.max = 500, minprior = 0.1, tol = 0.01) cat("\n### TWO-CLASS FINITE MIXTURE POSSION MODEL ###\n") mdl1 <- FLXMRglm(family = "poisson") fit1 <- flexmix(as.formula(formula), data = data[data$CARDHLDR == 1, ], k = 2, model = mdl1, concomitant = class, control = control) refit1 <- refit(fit1, method = 'optim') cat("\n=== MODEL THE RESPONSE ===\n") summary(refit1, which = 'model') cat("\n=== MODEL THE MIXTURE DISTRIBUTION ===\n") summary(refit1, which = 'concomitant') cat("\n### ZERO-INFLATED POSSION MODEL ###\n") mdl2 <- FLXMRziglm(family = "poisson") fit <- flexmix(as.formula(formula), data = data[data$CARDHLDR == 1, ], k = 2 , model = mdl2, concomitant = class, control = control) refit2 <- refit(fit, method = 'optim') cat("\n=== MODEL THE RESPONSE ===\n") summary(refit2, which = 'model') cat("\n=== MODEL THE MIXTURE DISTRIBUTION ===\n") summary(refit2, which = 'concomitant')
R Output for 2-Class Finite Mixture Poisson
### TWO-CLASS FINITE MIXTURE POSSION MODEL ### Classification: weighted 2 Log-likelihood : -4303.5967 converged === MODEL THE RESPONSE === $Comp.1 Estimate Std. Error z value Pr(>|z|) (Intercept) -8.0940843 1.6102947 -5.0265 4.996e-07 *** AGE 0.0114988 0.0129009 0.8913 0.372759 ACADMOS 0.0045677 0.0020299 2.2502 0.024438 * MINORDRG 0.2641256 0.6769000 0.3902 0.696390 LOGSPEND 0.6826690 0.2224763 3.0685 0.002151 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 $Comp.2 Estimate Std. Error z value Pr(>|z|) (Intercept) -2.44490331 0.34951344 -6.9952 2.650e-12 *** AGE 0.02214164 0.00662479 3.3422 0.0008311 *** ACADMOS 0.00052922 0.00077757 0.6806 0.4961209 MINORDRG 0.05054178 0.04048630 1.2484 0.2118965 LOGSPEND 0.21398000 0.04127917 5.1837 2.175e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 === MODEL THE MIXTURE DISTRIBUTION === $Comp.2 Estimate Std. Error z value Pr(>|z|) (Intercept) -1.4274523 0.5275366 -2.7059 0.006812 ** AGE -0.0027648 0.0100981 -0.2738 0.784246 ACADMOS 0.0016143 0.0014455 1.1168 0.264096 MINORDRG 1.5865202 0.1791074 8.8579 < 2.2e-16 *** LOGSPEND -0.0695020 0.0745171 -0.9327 0.350976 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R Output for Zero-Inflated Poisson
### ZERO-INFLATED POSSION MODEL ### Classification: weighted 3 Log-likelihood : -4175.7431 converged === MODEL THE RESPONSE === $Comp.2 Estimate Std. Error z value Pr(>|z|) (Intercept) -2.27782892 0.30030188 -7.5851 3.322e-14 *** AGE 0.01955236 0.00602142 3.2471 0.001166 ** ACADMOS 0.00024907 0.00067394 0.3696 0.711698 MINORDRG 0.11764062 0.02711666 4.3383 1.436e-05 *** LOGSPEND 0.16441222 0.03531969 4.6550 3.240e-06 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 === MODEL THE MIXTURE DISTRIBUTION === $Comp.2 Estimate Std. Error z value Pr(>|z|) (Intercept) -1.91132158 0.41706481 -4.5828 4.588e-06 *** AGE -0.00081716 0.00841240 -0.0971 0.922617 ACADMOS 0.00293407 0.00109997 2.6674 0.007644 ** MINORDRG 1.44243620 0.13613625 10.5955 < 2.2e-16 *** LOGSPEND 0.09561048 0.05080464 1.8819 0.059846 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
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.