|Data Scientist|Business Analyst|Business Intelligence Specialist|Senior Project Management Officer|Senior Project Planning|
|EASA Part M CAMO 145 Technical Service Engineer|Aircraft Asset Management Specialist|
options(width=60, prompt = "R> ", continue = "+ ", useFancyQuotes = FALSE)
library("graphics")
library("stats")
library("flexmix")
library("lattice")
ltheme <- canonical.theme("postscript", FALSE)
lattice.options(default.theme=ltheme)
data("NPreg", package = "flexmix")
data("dmft", package = "flexmix")
source("C:/Users/Test/Documents/R/win-library/3.5/flexmix/doc/myConcomitant.R")
par(mfrow=c(1,2))
plot(yn~x, col=class, pch=class, data=NPreg)
plot(yp~x, col=class, pch=class, data=NPreg)
set.seed(1802)
library("flexmix")
data("NPreg", package = "flexmix")
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))
## Classification: weighted
## 10 Log-likelihood : -1044.7688
## 11 Log-likelihood : -1044.7678
## converged
m1.refit <- refit(m1)
summary(m1.refit, which = "model", model = 1)
## $Comp.1
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 14.58965 1.24635 11.706 < 2.2e-16 ***
## x 9.91572 0.55294 17.933 < 2.2e-16 ***
## I(x^2) -0.97578 0.05201 -18.762 < 2.2e-16 ***
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.2
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.140549 0.961868 -0.1461 0.8838
## x 4.732610 0.474428 9.9754 <2e-16 ***
## I(x^2) 0.042722 0.046890 0.9111 0.3622
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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))
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
##
## Call:
## flexmix(formula = . ~ x, data = NPreg, cluster = posterior(m1),
## model = list(Model_n2, Model_p))
##
## Cluster sizes:
## 1 2
## 96 104
##
## convergence after 3 iterations
## [1] 2158.414 2149.956
data("betablocker", package = "flexmix")
betaGlm <- glm(cbind(Deaths, Total - Deaths) ~ Treatment,
family = "binomial", data = betablocker)
betaGlm
##
## Call: glm(formula = cbind(Deaths, Total - Deaths) ~ Treatment, family = "binomial",
## data = betablocker)
##
## Coefficients:
## (Intercept) TreatmentTreated
## -2.1971 -0.2574
##
## Degrees of Freedom: 43 Total (i.e. Null); 42 Residual
## Null Deviance: 333
## Residual Deviance: 305.8 AIC: 527.2
betaMixFix <- stepFlexmix(cbind(Deaths, Total - Deaths) ~ 1 | Center,
model = FLXMRglmfix(family = "binomial", fixed = ~ Treatment),
k = 2:4, nrep = 5, data = betablocker)
## 2 : * * * * *
## 3 : * * * * *
## 4 : * * * * *
##
## Call:
## stepFlexmix(cbind(Deaths, Total - Deaths) ~ 1 |
## Center, model = FLXMRglmfix(family = "binomial",
## fixed = ~Treatment), data = betablocker, k = 2:4,
## nrep = 5)
##
## iter converged k k0 logLik AIC BIC ICL
## 2 12 TRUE 2 2 -181.3308 370.6617 377.7984 380.2105
## 3 9 TRUE 3 3 -159.3605 330.7210 341.4262 343.3257
## 4 13 TRUE 4 4 -158.2465 332.4929 346.7664 354.4909
betaMixFix_3 <- getModel(betaMixFix, which = "BIC")
betaMixFix_3 <- relabel(betaMixFix_3, "model", "Intercept")
## Comp.1 Comp.2 Comp.3
## coef.TreatmentTreated -0.2581849 -0.2581849 -0.2581849
## coef.(Intercept) -2.8336816 -2.2501814 -1.6097872
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))
}
betaMix <- stepFlexmix(cbind(Deaths, Total - Deaths) ~ Treatment | Center,
model = FLXMRglm(family = "binomial"), k = 3, nrep = 5,
data = betablocker)
## 3 : * * * * *
betaMix <- relabel(betaMix, "model", "Treatment")
parameters(betaMix)
## Comp.1 Comp.2 Comp.3
## coef.(Intercept) -1.5800031 -2.2476980 -2.91633722
## coef.TreatmentTreated -0.3248497 -0.2630017 -0.08047829
c(BIC(betaMixFix_3), BIC(betaMix))
## [1] 341.4262 346.8925
print(plot(betaMixFix_3, nint = 10, mark = 1, col = "grey", layout = c(3, 1)))
print(plot(betaMixFix_3, nint = 10, mark = 2, col = "grey", layout = c(3, 1)))
##
## 1 2 3
## 10 24 10
predict(betaMix,
newdata = data.frame(Treatment = c("Control", "Treated")))
## $Comp.1
## [,1]
## 1 0.1707950
## 2 0.1295602
##
## $Comp.2
## [,1]
## 1 0.09554822
## 2 0.07511149
##
## $Comp.3
## [,1]
## 1 0.05135184
## 2 0.04756995
## Deaths Total Center Treatment
## 1 3 39 1 Control
## 23 3 38 1 Treated
fitted(betaMix)[c(1, 23), ]
## Comp.1 Comp.2 Comp.3
## [1,] 0.1707950 0.09554822 0.05135184
## [2,] 0.1295602 0.07511149 0.04756995
## $Comp.1
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.579965 0.065997 -23.9401 < 2.2e-16
## TreatmentTreated -0.324833 0.092882 -3.4973 0.0004701
##
## (Intercept) ***
## TreatmentTreated ***
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.2
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.247678 0.045181 -49.7481 < 2.2e-16
## TreatmentTreated -0.262990 0.065598 -4.0091 6.095e-05
##
## (Intercept) ***
## TreatmentTreated ***
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.3
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.91634 0.09921 -29.3958 <2e-16 ***
## TreatmentTreated -0.08048 0.14104 -0.5706 0.5683
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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)
## $Comp.1
## coef.TreatmentTreated coef.(Intercept)
## -0.2837779 -1.5985089
##
## $Comp.2
## coef.TreatmentTreated coef.(Intercept)
## -0.2837779 -2.2379835
##
## $Comp.3
## coef.(Intercept)
## -2.956159
c(BIC(betaMix), BIC(betaMixNested), BIC(betaMixFix_3))
## [1] 346.8925 339.9429 341.4262
data("bioChemists", package = "flexmix")
data("bioChemists", package = "flexmix")
Model1 <- FLXMRglm(family = "poisson")
ff_1 <- stepFlexmix(art ~ ., data = bioChemists, k = 1:3, model = Model1)
## 1 : * * *
## 2 : * * *
## 3 : * * *
ff_1 <- getModel(ff_1, "BIC")
print(plot(refit(ff_1), bycluster = FALSE,
scales = list(x = list(relation = "free"))))
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))
## [1] 3212.991 3200.071
## $Comp.1
## Estimate Std. Error z value Pr(>|z|)
## kid5 -0.2070933 0.0521307 -3.9726 7.11e-05 ***
## marMarried 0.1646178 0.0768505 2.1421 0.0321892 *
## ment 0.0274302 0.0032899 8.3377 < 2.2e-16 ***
## (Intercept) -0.4782683 0.2203178 -2.1708 0.0299455 *
## femWomen 0.6045078 0.1822082 3.3177 0.0009077 ***
## phd 0.0775564 0.0408161 1.9001 0.0574148 .
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Comp.2
## Estimate Std. Error z value Pr(>|z|)
## kid5 -0.2070933 0.0521307 -3.9726 7.110e-05 ***
## marMarried 0.1646178 0.0768505 2.1421 0.03219 *
## ment 0.0274302 0.0032899 8.3377 < 2.2e-16 ***
## (Intercept) 1.3208461 0.2337508 5.6507 1.598e-08 ***
## femWomen -2.2054873 0.4313329 -5.1132 3.168e-07 ***
## phd -0.0818211 0.0554627 -1.4752 0.14015
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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))
## [1] 3200.071 3192.816
print(plot(refit(ff_3), bycluster = FALSE, scales = list(x = list(relation = "free"))))
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)
## Comp.1 Comp.2
## coef.kid5 -0.1819162 -0.1819162
## coef.marMarried 0.1883990 0.1883990
## coef.ment 0.0288501 0.0288501
## coef.(Intercept) -0.2583828 0.9950050
summary(refit(ff_4), which = "concomitant")
## $Comp.2
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.02262 0.28385 -3.6027 0.0003149 ***
## femWomen -0.61281 0.27280 -2.2464 0.0246782 *
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] 3182.328
Model5 <- FLXMRglmfix(family = "poisson", fixed = ~ kid5 + ment + fem)
ff_5 <- flexmix(art ~ 1, data = bioChemists, cluster = posterior(ff_2),
model = Model5)
BIC(ff_5)
## [1] 3174.266
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")
data("dmft", package = "flexmix")
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))
## $Comp.2
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1470652 0.0963031 -1.5271 0.126734
## log(Begin + 0.5) 0.7303431 0.0402403 18.1496 < 2.2e-16
## Gendermale 0.0068805 0.0550486 0.1250 0.900531
## Ethnicwhite 0.0500281 0.0592975 0.8437 0.398848
## Ethnicblack -0.0472578 0.0899984 -0.5251 0.599516
## Treatmenteduc -0.2371848 0.0905877 -2.6183 0.008837
## Treatmentall -0.3277747 0.1011637 -3.2400 0.001195
## Treatmentenrich 0.0172642 0.0838730 0.2058 0.836918
## Treatmentrinse -0.2414627 0.0871032 -2.7721 0.005569
## Treatmenthygiene -0.1026315 0.0916676 -1.1196 0.262882
##
## (Intercept)
## log(Begin + 0.5) ***
## Gendermale
## Ethnicwhite
## Ethnicblack
## Treatmenteduc **
## Treatmentall **
## Treatmentenrich
## Treatmentrinse **
## Treatmenthygiene
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## print(plot(refit(Fitted), components = 2, box.ratio = 3))
print(plot(refit(Fitted), components = 2, box.ratio = 3))
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)
##
## Call:
## flexmix(formula = . ~ x, data = NPreg, k = 2,
## model = list(Model_n, Model_p), concomitant = Concomitant)
##
## prior size post>0 ratio
## Comp.1 0.493 96 138 0.696
## Comp.2 0.507 104 137 0.759
##
## 'log Lik.' -1044.724 (df=14)
## AIC: 2117.448 BIC: 2163.625
##
## Call:
## flexmix(formula = . ~ x, data = NPreg, k = 2,
## cluster = posterior(m2), model = list(Model_n,
## Model_p), concomitant = MyConcomitant)
##
## prior size post>0 ratio
## Comp.1 0.493 96 138 0.696
## Comp.2 0.507 104 137 0.759
##
## 'log Lik.' -1044.724 (df=14)
## AIC: 2117.448 BIC: 2163.625
determinePrior <- function(object) {
object@concomitant@fit(object@concomitant@x,
posterior(object))[!duplicated(object@concomitant@x), ]
}
## [,1] [,2]
## 1 0.4830052 0.5169948
## 4 0.5050200 0.4949800
## [,1] [,2]
## 1 0.4830037 0.5169963
## 2 0.5050168 0.4949832
SI <- sessionInfo()
pkgs <- paste(sapply(c(SI$otherPkgs, SI$loadedOnly), function(x)
paste("\\\\pkg{", x$Package, "} ",
x$Version, sep = "")), collapse = ", ")