###################################################
### chunk number 1: 
###################################################
options(width=65, prompt = "R> ", continue = "+  ", useFancyQuotes = FALSE)
library("flexmix")
ltheme <- canonical.theme("postscript", FALSE)
lattice.options(default.theme=ltheme)
data("NPreg")
data("dmft")
source("v28i04-concomitant.R")
source("v28i04-ziglm.R")


###################################################
### chunk number 2: 
###################################################
par(mfrow=c(1,2))
plot(yn~x, col=class, pch=class, data=NPreg)
plot(yp~x, col=class, pch=class, data=NPreg)


###################################################
### chunk number 3: 
###################################################
set.seed(1802)
library("flexmix")
data("NPreg")
Model_n <- FLXMRglm(yn ~ . + I(x^2))
Model_p <- FLXMRglm(yp ~ ., family = "poisson")
m1 <- flexmix(. ~ x, data = NPreg, k = 2, model = list(Model_n, Model_p),
  control = list(verbose = 10))


###################################################
### chunk number 4: 
###################################################
print(plot(m1))


###################################################
### chunk number 5: 
###################################################
m1.refit <- refit(m1)
summary(m1.refit, which = "model", model = 1)


###################################################
### chunk number 6: 
###################################################
print(plot(m1.refit, layout = c(1,3), bycluster = FALSE,
      main = expression(paste(yn *tilde(" ")* x + x^2))),
      split= c(1,1,2,1), more = TRUE)
print(plot(m1.refit, model = 2, 
           main = expression(paste(yp *tilde(" ")* x)), 
           layout = c(1,2), bycluster = FALSE), 
      split = c(2,1,2,1))


###################################################
### chunk number 7: 
###################################################
Model_n2 <- FLXMRglmfix(yn ~ . + 0, nested = list(k = c(1, 1), 
  formula = c(~ 1 + I(x^2), ~ 0)))
m2 <- flexmix(. ~ x, data = NPreg, cluster = posterior(m1), 
  model = list(Model_n2, Model_p))
m2


###################################################
### chunk number 8: 
###################################################
c(BIC(m1), BIC(m2))


###################################################
### chunk number 9: 
###################################################
data("betablocker")
betaGlm <- glm(cbind(Deaths, Total - Deaths) ~ Treatment, 
  family = "binomial", data = betablocker)
betaGlm


###################################################
### chunk number 10: 
###################################################
betaMixFix <- stepFlexmix(cbind(Deaths, Total - Deaths) ~ 1 | Center,
  model = FLXMRglmfix(family = "binomial", fixed = ~ Treatment), 
  k = 2:4, nrep = 5, data = betablocker)


###################################################
### chunk number 11: 
###################################################
betaMixFix


###################################################
### chunk number 12: 
###################################################
betaMixFix_3 <- getModel(betaMixFix, which = "BIC")


###################################################
### chunk number 13: 
###################################################
parameters(betaMixFix_3)


###################################################
### chunk number 14: 
###################################################
library("grid")
betablocker$Center <- with(betablocker, factor(Center, levels = Center[order((Deaths/Total)[1:22])]))
clusters <- factor(clusters(betaMixFix_3), labels = paste("Cluster", 1:3))
print(dotplot(Deaths/Total ~ Center | clusters, groups  = Treatment, as.table = TRUE,
              data  = betablocker, xlab = "Center", layout = c(3, 1), 
              scales  = list(x = list(cex = 0.7, tck = c(1, 0))),
              key = simpleKey(levels(betablocker$Treatment), lines = TRUE, corner = c(1,0))))
betaMixFix.fitted <- fitted(betaMixFix_3)
for (i in 1:3) {
  seekViewport(trellis.vpname("panel", i, 1))
  grid.lines(unit(1:22, "native"), unit(betaMixFix.fitted[1:22, i], "native"), gp = gpar(lty = 1))
  grid.lines(unit(1:22, "native"), unit(betaMixFix.fitted[23:44, i], "native"), gp = gpar(lty = 2))
}


###################################################
### chunk number 15: 
###################################################
betaMix <- stepFlexmix(cbind(Deaths, Total - Deaths) ~ Treatment | Center,
  model = FLXMRglm(family = "binomial"), k = 3, nrep = 5, 
  data = betablocker)
parameters(betaMix)
c(BIC(betaMixFix_3), BIC(betaMix))


###################################################
### chunk number 16: 
###################################################
print(plot(betaMixFix_3, nint = 10, mark = 1, col = "grey", layout = c(3, 1)))


###################################################
### chunk number 17: 
###################################################
print(plot(betaMixFix_3, nint = 10, mark = 3, col = "grey", layout = c(3, 1)))


###################################################
### chunk number 18: 
###################################################
table(clusters(betaMix))


###################################################
### chunk number 19: 
###################################################
predict(betaMix, 
  newdata = data.frame(Treatment = c("Control", "Treated")))


###################################################
### chunk number 20: 
###################################################
betablocker[c(1, 23), ]
fitted(betaMix)[c(1, 23), ]


###################################################
### chunk number 21: 
###################################################
summary(refit(betaMix))


###################################################
### chunk number 22: 
###################################################
ModelNested <- FLXMRglmfix(family = "binomial", nested = list(k = c(2, 1),
  formula = c(~ Treatment, ~ 0)))
betaMixNested <- flexmix(cbind(Deaths, Total - Deaths) ~ 1 | Center,
  model = ModelNested, k = 3, data = betablocker, 
  cluster = posterior(betaMix))
parameters(betaMixNested)
c(BIC(betaMix), BIC(betaMixNested), BIC(betaMixFix_3))


###################################################
### chunk number 23: 
###################################################
data("bioChemists")


###################################################
### chunk number 24: 
###################################################
data("bioChemists")
Model1 <- FLXMRglm(family = "poisson")
ff_1 <- stepFlexmix(art ~ ., data = bioChemists, k = 1:3, model = Model1)
ff_1 <- getModel(ff_1, "BIC")


###################################################
### chunk number 25: 
###################################################
print(plot(refit(ff_1), bycluster = FALSE, 
           scales = list(x = list(relation = "free"))))


###################################################
### chunk number 26: 
###################################################
Model2 <- FLXMRglmfix(family = "poisson", fixed = ~ kid5 + mar + ment)
ff_2 <- flexmix(art ~ fem + phd, data = bioChemists, 
  cluster = posterior(ff_1), model = Model2)
c(BIC(ff_1), BIC(ff_2))


###################################################
### chunk number 27: 
###################################################
summary(refit(ff_2))


###################################################
### chunk number 28: 
###################################################
Model3 <- FLXMRglmfix(family = "poisson", fixed = ~ kid5 + mar + ment)
ff_3 <- flexmix(art ~ fem, data = bioChemists, cluster = posterior(ff_2),
  model = Model3)
c(BIC(ff_2), BIC(ff_3))


###################################################
### chunk number 29: 
###################################################
print(plot(refit(ff_3), bycluster = FALSE, scales = list(x = list(relation = "free"))))


###################################################
### chunk number 30: 
###################################################
Model4 <- FLXMRglmfix(family = "poisson", fixed = ~ kid5 + mar + ment)
ff_4 <- flexmix(art ~ 1, data = bioChemists, cluster = posterior(ff_2),
  concomitant = FLXPmultinom(~ fem), model = Model4)
parameters(ff_4)
summary(refit(ff_4), which = "concomitant")
BIC(ff_4)


###################################################
### chunk number 31: 
###################################################
Model5 <- FLXMRglmfix(family = "poisson", fixed = ~ kid5 + ment + fem)
ff_5 <- flexmix(art ~ 1, data = bioChemists, cluster = posterior(ff_2),
  model = Model5)
BIC(ff_5)


###################################################
### chunk number 32: 
###################################################
pp <- predict(ff_5, newdata = data.frame(kid5 = 0, 
  mar = factor("Married", levels = c("Single", "Married")),
  fem = c("Men", "Women"),  ment = mean(bioChemists$ment)))
matplot(0:12, sapply(unlist(pp), function(x) dpois(0:12, x)), 
        type = "b", lty = 1, xlab = "Number of articles", ylab = "Probability")
legend("topright", paste("Comp.", rep(1:2, each = 2), ":",
  c("Men", "Women")), lty = 1, col = 1:4, pch = paste(1:4), bty = "n")


###################################################
### chunk number 33: 
###################################################
data("dmft")
Model <- FLXMRziglm(family = "poisson")
Fitted <- flexmix(End ~ log(Begin + 0.5) + Gender + Ethnic + Treatment, 
  model = Model, k = 2 , data = dmft, control = list(minprior = 0.01))
summary(refit(Fitted))


###################################################
### chunk number 34: refit eval=FALSE
###################################################
## print(plot(refit(Fitted), components = 2, box.ratio = 3))


###################################################
### chunk number 35: 
###################################################
print(plot(refit(Fitted), components = 2, box.ratio = 3))


###################################################
### chunk number 36: 
###################################################
Concomitant <- FLXPmultinom(~ yb)
MyConcomitant <- myConcomitant(~ yb)

m2 <- flexmix(. ~ x, data = NPreg, k = 2, model = list(Model_n, Model_p), 
  concomitant = Concomitant)
m3 <- flexmix(. ~ x, data = NPreg, k = 2, model = list(Model_n, Model_p), 
  cluster = posterior(m2), concomitant = MyConcomitant)


###################################################
### chunk number 37: 
###################################################
summary(m2)
summary(m3)


###################################################
### chunk number 38: 
###################################################
determinePrior <- function(object) {
  object@concomitant@fit(object@concomitant@x, 
    posterior(object))[!duplicated(object@concomitant@x), ]
}


###################################################
### chunk number 39: 
###################################################
determinePrior(m2)
determinePrior(m3)


###################################################
### chunk number 40: 
###################################################
SI <- sessionInfo()
pkgs <- paste(sapply(c(SI$otherPkgs, SI$loadedOnly), function(x) 
                     paste("\\\\pkg{", x$Package, "} ", 
                           x$Version, sep = "")), collapse = ", ")