FLXMRnegbin {countreg} | R Documentation |
FlexMix driver for fitting of negative binomial regression models.
FLXMRnegbin(formula = . ~ ., theta = NULL, offset = NULL, control = list(reltol = .Machine$double.eps^(1/1.5), maxit = 500))
formula |
formula. This is interpreted relative to the formula
specified in the call to |
theta |
numeric or |
offset |
numeric. Optional offset vector for the linear predictor. |
control |
list with control parameters passed to
|
The driver function FLXMRnegbin
enables estimation of finite mixtures
of negative binomial regression models via flexmix
or
stepFlexmix
. The driver is modeled after FLXMRglm
and supports both fixed and unknown theta. In the M-step for fixed theta,
glm.fit
is employed along with the negative.binomial
family. If the theta
is unknown and has be estimated
along with the regression coefficients, direct optimization using
optim
with analytical gradients is employed.
An object of class FLXMRglm
.
flexmix
,
stepFlexmix
,
FLXMRglm
,
negative.binomial
## artificial data from a two-component mixture of geometric regressions set.seed(1) d <- data.frame(x = runif(500, -1, 1)) d$cluster <- rep(1:2, each = 250) d$y <- rnbinom(500, mu = exp(c(1, -1)[d$cluster] + c(0, 3)[d$cluster] * d$x), size = 1) if(require("flexmix")) { ## fit mixture models with known correct theta and unknown theta fm1 <- flexmix(y ~ x, data = d, k = 2, model = FLXMRnegbin(theta = 1)) fm0 <- flexmix(y ~ x, data = d, k = 2, model = FLXMRnegbin()) ## parameter recovery parameters(fm1) parameters(fm0) ## refit to obtain joint summary summary(refit(fm1, gradient = NULL)) summary(refit(fm0, gradient = NULL)) ## refitting both components manually for rootograms rf1 <- lapply(1:2, function(i) glm(y ~ x, data = d, family = negative.binomial(1), weights = posterior(fm1)[,i])) rf0 <- lapply(1:2, function(i) glm.nb(y ~ x, data = d, weights = posterior(fm0)[,i])) ## Rootograms par(mfrow = c(1, 2)) r11 <- rootogram(rf1[[1]]) r12 <- rootogram(rf1[[2]]) r01 <- rootogram(rf0[[1]]) r02 <- rootogram(rf0[[2]]) rootogram(glm.nb(y ~ x, data = d)) plot(r01 + r02) } ## Not run: ## two-component mixture model fro NMES1988 physician office visits ## (fitting takes some time...) if(require("flexmix") & require("AER")) { ## data from AER data("NMES1988", package = "AER") nmes <- NMES1988[, c(1, 7:8, 13, 15, 18:19)] ## single-component model nmes_nb <- glm.nb(visits ~ ., data = nmes) ## two-component model set.seed(1090) nmes_fnb <- stepFlexmix(visits ~ ., data = nmes, k = 2, model = FLXMRnegbin()) ## refit to obtain summary with estimate of joint covariance matrix summary(refit(nmes_fnb, gradient = NULL)) ## refit individual models manually for rootograms nmes_fnb_rf <- lapply(1:2, function(i) glm.nb(visits ~ ., data = nmes, weights = posterior(nmes_fnb)[,i])) r1 <- rootogram(nmes_fnb_rf[[1]], max = 50, plot = FALSE) r2 <- rootogram(nmes_fnb_rf[[2]], max = 50, plot = FALSE) par(mfrow = c(2, 2)) rootogram(nmes_nb, max = 50, main = "Negative Binomial", ylim = c(-1, 25)) plot(r1 + r2, xlab = "visits", main = "Mixture Negative Binomial", ylim = c(-1, 25)) plot(r1, main = "Mixture Negative Binomial (Component 1)", ylim = c(-1, 25)) plot(r2, main = "Mixture Negative Binomial (Component 2)", ylim = c(-1, 25)) } ## End(Not run)