R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("flexmix-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('flexmix')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "ExNPreg"
> 
> ### * ExNPreg
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ExNPreg
> ### Title: Artificial Example for Normal and Poisson Regression
> ### Aliases: ExNPreg NPreg
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(NPreg)
> plot(yn~x, data=NPreg, col=class)
> plot(yp~x, data=NPreg, col=class)
> 
> 
> 
> cleanEx(); ..nameEx <- "ExNclus"
> 
> ### * ExNclus
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ExNclus
> ### Title: Artificial Example for Normal Clustering
> ### Aliases: ExNclus Nclus
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(Nclus)
> require("MASS")
Loading required package: MASS
[1] TRUE
> eqscplot(Nclus, col=rep(1:4, c(100, 100, 150, 200)))
> 
> 
> 
> cleanEx(); ..nameEx <- "FLXcontrol-class"
> 
> ### * FLXcontrol-class
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: FLXcontrol-class
> ### Title: Class "FLXcontrol"
> ### Aliases: FLXcontrol-class coerce,list,FLXcontrol-method
> ###   coerce,NULL,FLXcontrol-method
> ### Keywords: classes
> 
> ### ** Examples
> 
> ## have a look at the defaults
> new("FLXcontrol")
An object of class "FLXcontrol"
Slot "iter.max":
[1] 200

Slot "minprior":
[1] 0.05

Slot "tolerance":
[1] 1e-06

Slot "verbose":
[1] 0

Slot "classify":
[1] "auto"

Slot "nrep":
[1] 1

> 
> ## corce a list
> mycont = list(iter=200, tol=0.001, class="r")
> as(mycont, "FLXcontrol")
An object of class "FLXcontrol"
Slot "iter.max":
[1] 200

Slot "minprior":
[1] 0.05

Slot "tolerance":
[1] 0.001

Slot "verbose":
[1] 0

Slot "classify":
[1] "random"

Slot "nrep":
[1] 1

> 
> 
> 
> cleanEx(); ..nameEx <- "FLXmclust"
> 
> ### * FLXmclust
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: FLXmclust
> ### Title: FlexMix Clustering Demo Driver
> ### Aliases: FLXmclust plotEll
> ### Keywords: cluster
> 
> ### ** Examples
> 
> data(Nclus)
> 
> require("MASS")
Loading required package: MASS
[1] TRUE
> eqscplot(Nclus)
> 
> ## This model is wrong (one component has a non-diagonal cov matrix)
> ex1 <- flexmix(Nclus~1, k=4, model=FLXmclust())
Loading required package: mvtnorm
> print(ex1)

Call:
flexmix(formula = Nclus ~ 1, k = 4, model = FLXmclust())

Cluster sizes:
  1   2   3   4 
149  92 213  96 

convergence after 63 iterations
> plotEll(ex1, Nclus)
Loading required package: ellipse
> 
> ## True model, wrong number of components
> ex2 <- flexmix(Nclus~1, k=6, model=FLXmclust(diag=FALSE))  
> print(ex2)

Call:
flexmix(formula = Nclus ~ 1, k = 6, model = FLXmclust(diag = FALSE))

Cluster sizes:
  1   2   3   4 
150 100  96 204 

convergence after 29 iterations
> 
> plotEll(ex2, Nclus)
> 
> ## Get paramters of first component
> parameters(ex2, component=1)
$center
[1] -2.008795  6.054081

$cov
            [,1]        [,2]
[1,]  2.17825068 -0.02057670
[2,] -0.02057670  1.08106779

> 
> ## Have a look at the posterior probabilies of 10 random observations
> ok <- sample(1:nrow(Nclus), 10)
> p  <- posterior(ex2)[ok,]
> p
              [,1]         [,2]         [,3]          [,4]
 [1,] 1.541356e-07 6.531646e-04 6.650272e-25  9.993467e-01
 [2,] 4.790748e-16 1.000000e+00 4.541016e-21  3.420157e-55
 [3,] 1.000000e+00 4.346009e-26 4.188792e-09  6.383958e-69
 [4,] 2.773366e-05 1.308274e-16 9.976103e-01  2.361919e-03
 [5,] 3.973690e-11 6.321230e-17 9.999970e-01  3.007601e-06
 [6,] 1.000000e+00 1.322897e-26 5.120674e-09  5.238854e-70
 [7,] 6.242923e-25 1.000000e+00 4.636363e-22 2.751679e-130
 [8,] 9.999999e-01 1.764226e-18 6.145835e-08  3.007147e-37
 [9,] 4.496454e-06 2.987152e-05 1.605489e-10  9.999656e-01
[10,] 1.000000e+00 7.613846e-24 5.704144e-10  7.604974e-67
> 
> ## The following two should be the same
> max.col(p)
 [1] 4 2 1 3 3 1 2 1 4 1
> cluster(ex2)[ok]
 [1] 4 2 1 3 3 1 2 1 4 1
> 
> ## Don't show:
>   stopifnot(all.equal(max.col(p), cluster(ex2)[ok]))
> ## End Don't show
> 
> 
> 
> cleanEx(); ..nameEx <- "KLdiv"
> 
> ### * KLdiv
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: KLdiv
> ### Title: Kullback-Leibler Divergence
> ### Aliases: KLdiv KLdiv-methods KLdiv,matrix-method KLdiv,flexmix-method
> ### Keywords: methods
> 
> ### ** Examples
> 
> x = (1:100)/100
> ## Gaussian and Student t are much closer to each other than
> ## to the uniform:
> KLdiv(cbind(u=dunif(x), n=dnorm(x), t=dt(x, df=10)))
          u           n          t
u  0.000000  1.08811353 1.12460761
n -1.066617  0.00000000 0.03529794
t -1.100701 -0.03522973 0.00000000
> 
> 
> 
> cleanEx(); ..nameEx <- "fitted"
> 
> ### * fitted
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: fitted
> ### Title: Extract Model Fitted Values
> ### Aliases: fitted-methods fitted,flexmix-method fitted,FLXrefit-method
> ###   fitted,FLXrefitglm-method
> ### Keywords: methods
> 
> ### ** Examples
> 
> data(NPreg)
> ex1 <- flexmix(yn~x+I(x^2), data=NPreg, k=2)
> matplot(NPreg$x, fitted(ex1), pch=16, type="p")
> points(NPreg$x, NPreg$yn)
> 
> 
> 
> cleanEx(); ..nameEx <- "flexmix"
> 
> ### * flexmix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: flexmix
> ### Title: Flexible Mixture Modeling
> ### Aliases: flexmix flexmix,formula,ANY,ANY,ANY,missing-method
> ###   flexmix,formula,ANY,ANY,ANY,list-method
> ###   flexmix,formula,ANY,ANY,ANY,FLXmodel-method predict,flexmix-method
> ###   show,flexmix-method summary,flexmix-method
> ###   show,summary.flexmix-method
> ### Keywords: regression cluster
> 
> ### ** Examples
> 
> 
> data(NPreg)
> 
> ## mixture of two linear regression models. Note that control parameters
> ## can be specified as named list and abbreviated if unique.
> ex1 <- flexmix(yn~x+I(x^2), data=NPreg, k=2,
+                    control=list(verb=5, iter=100))
Classification: weighted 
   5 Log-likelihood :    -728.2984 
  10 Log-likelihood :    -727.0263 
  15 Log-likelihood :    -709.5686 
  20 Log-likelihood :    -642.5701 
  23 Log-likelihood :    -642.5453 
converged
> 
> ex1

Call:
flexmix(formula = yn ~ x + I(x^2), data = NPreg, k = 2, control = list(verb = 5, 
    iter = 100))

Cluster sizes:
  1   2 
100 100 

convergence after 23 iterations
> summary(ex1)

Call:
flexmix(formula = yn ~ x + I(x^2), data = NPreg, k = 2, control = list(verb = 5, 
    iter = 100))

       prior size post>0 ratio
Comp.1 0.494  100    145 0.690
Comp.2 0.506  100    141 0.709

'log Lik.' -642.5453 (df=9)
AIC: 1303.091   BIC: 1332.775 

> plot(ex1)
> 
> ## now we fit a model with one Gaussian response and one Poisson
> ## response. Note that the formulas inside the call to FLXglm are
> ## relative to the overall model formula.
> ex2 <- flexmix(yn~x, data=NPreg, k=2,
+                model=list(FLXglm(yn~.+I(x^2)), 
+                           FLXglm(yp~., family="poisson")))
> plot(ex2)
> 
> ex2

Call:
flexmix(formula = yn ~ x, data = NPreg, k = 2, model = list(FLXglm(yn ~ 
    . + I(x^2)), FLXglm(yp ~ ., family = "poisson")))

Cluster sizes:
  1   2 
104  96 

convergence after 11 iterations
> table(ex2@cluster, NPreg$class)
   
     1  2
  1 99  5
  2  1 95
> 
> ## for Gaussian responses we get coefficients and standard deviation
> parameters(ex2, component=1, model=1)
$coef
(Intercept)           x      I(x^2) 
-0.14069452  4.73362649  0.04258788 

$sigma
[1] 3.426301

> 
> ## for Poisson response we get only coefficients
> parameters(ex2, component=1, model=2)
$coef
(Intercept)           x 
  1.9391131  -0.1808908 

> 
> ## fitting a model only to the Poisson response is of course
> ## done like this
> ex3 <- flexmix(yp~x, data=NPreg, k=2, model=FLXglm(family="poisson"))
> 
> ## if observations are grouped, i.e., we have several observations per
> ## individual, fitting is usually much faster:
> ex4 <- flexmix(yp~x|id1, data=NPreg, k=2, model=FLXglm(family="poisson"))
> 
> 
> 
> cleanEx(); ..nameEx <- "refit"
> 
> ### * refit
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: refit
> ### Title: Refit a Fitted Model
> ### Aliases: refit refit-methods refit,flexmix-method
> ###   refit,FLXglmmodel-method show,FLXrefit-method summary,FLXrefit-method
> ###   summary,FLXrefitglm-method
> ### Keywords: methods
> 
> ### ** Examples
> 
> data(NPreg)
> ex1 <- flexmix(yn~x+I(x^2), data=NPreg, k=2)
> ex1r <- refit(ex1)
> 
> ## in one component all coefficients should be highly significant,
> ## in the other component only the linear term
> summary(ex1r)

Call:
refit(ex1)

Component 1 :
             Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.208966   0.673895 -0.3101   0.7568
x            4.816982   0.327445 14.7108   <2e-16
I(x^2)       0.036236   0.032545  1.1134   0.2669

-------------
Component 2 :
             Estimate Std. Error t value  Pr(>|t|)
(Intercept) 14.717553   0.890847  16.521 < 2.2e-16
x            9.846132   0.390386  25.221 < 2.2e-16
I(x^2)      -0.968302   0.036951 -26.205 < 2.2e-16


> 
> 
> 
> cleanEx(); ..nameEx <- "seizure"
> 
> ### * seizure
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: seizure
> ### Title: Epileptic Seizure Data
> ### Aliases: seizure
> ### Keywords: datasets
> 
> ### ** Examples
> 
> data(seizure)
> plot(Seizures/Hours~Day, col=as.integer(Treatment),
+      pch=as.integer(Treatment), data=seizure)
> abline(v=27.5, lty=2, col="grey")
> legend(140, 9, c("Baseline", "Treatment"),
+        pch=1:2, col=1:2, xjust=1, yjust=1)
> 
> set.seed(123)
> 
> ## The model presented in the Wang et al paper: two components for
> ## "good" and "bad" days, respectively, each a Poisson GLM with hours of
> ## parental observation as offset
> 
> seizMix <- flexmix(Seizures~Treatment*log(Day),
+                    data=seizure, k=2,
+                    model=FLXglm(family="poisson", offset=log(seizure$Hours)))
> 
> summary(seizMix)

Call:
flexmix(formula = Seizures ~ Treatment * log(Day), data = seizure, 
    k = 2, model = FLXglm(family = "poisson", offset = log(seizure$Hours)))

       prior size post>0 ratio
Comp.1 0.276   37    101 0.366
Comp.2 0.724  103    115 0.896

'log Lik.' -376.1762 (df=9)
AIC: 770.3525   BIC: 796.8272 

> summary(refit(seizMix))

Call:
refit(seizMix)

Component 1 :
                       Estimate Std. Error z value  Pr(>|z|)
(Intercept)            2.844612   0.234145 12.1489 < 2.2e-16
TreatmentYes           1.302730   0.474049  2.7481  0.005994
log(Day)              -0.406097   0.088312 -4.5984 4.257e-06
TreatmentYes:log(Day) -0.431234   0.133471 -3.2309  0.001234

-------------
Component 2 :
                       Estimate Std. Error  z value  Pr(>|z|)
(Intercept)            2.070332   0.089172  23.2174 < 2.2e-16
TreatmentYes           7.432105   0.517962  14.3487 < 2.2e-16
log(Day)              -0.270568   0.038194  -7.0839 1.401e-12
TreatmentYes:log(Day) -2.276297   0.139674 -16.2972 < 2.2e-16


> 
> matplot(seizure$Day, fitted(seizMix)/seizure$Hours, type="l",
+         add=TRUE, col=3:4)
> 
> 
> 
> cleanEx(); ..nameEx <- "stepFlexmix"
> 
> ### * stepFlexmix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: stepFlexmix
> ### Title: Run FlexMix Repeatedly
> ### Aliases: stepFlexmix
> ### Keywords: cluster regression
> 
> ### ** Examples
> 
> data(NPreg)
> 
> ## try 5 times for k=2
> ex1 <- stepFlexmix(yn~x+I(x^2), data=NPreg, K=2, nrep=5)
2 : * * * * *
> ex1

Call:
stepFlexmix(yn ~ x + I(x^2), data = NPreg, K = 2, nrep = 5)

Cluster sizes:
  1   2 
100 100 

convergence after 15 iterations
> 
> ## now for k=2,3,4,5
> ## low nrep to have reasonable execution time of example
> ## even on slow systems
> ex2 <- stepFlexmix(yn~x+I(x^2), data=NPreg, K=2:5, nrep=2)
2 : * *
3 : * *
4 : * *
5 : * *
> ex2
$"2"

Call:
stepFlexmix(yn ~ x + I(x^2), data = NPreg, K = 2:5, nrep = 2)

Cluster sizes:
  1   2 
100 100 

convergence after 20 iterations

$"3"

Call:
stepFlexmix(yn ~ x + I(x^2), data = NPreg, K = 2:5, nrep = 2)

Cluster sizes:
  1   2   3 
  9 102  89 

convergence after 84 iterations

$"4"

Call:
stepFlexmix(yn ~ x + I(x^2), data = NPreg, K = 2:5, nrep = 2)

Cluster sizes:
 1  2  3  4 
12 93 69 26 

convergence after 46 iterations

$"5"

Call:
stepFlexmix(yn ~ x + I(x^2), data = NPreg, K = 2:5, nrep = 2)

Cluster sizes:
 1  2  3  4 
75 22 12 91 

convergence after 61 iterations

> sapply(ex2, logLik)
        2         3         4         5 
-642.5452 -636.2990 -635.3643 -635.3745 
> 
> ## model selection:
> sapply(ex2, AIC)
       2        3        4        5 
1303.090 1300.598 1308.729 1308.749 
> sapply(ex2, BIC)
       2        3        4        5 
1332.775 1346.774 1371.397 1371.417 
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  21.64 0.3 22.01 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')