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("mix-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('mix')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "da.mix"
> 
> ### * da.mix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: da.mix
> ### Title: Data Augmentation for Unrestricted General Location Model
> ### Aliases: da.mix
> ### Keywords: models
> 
> ### ** Examples
> 
> data(stlouis)
> s <- prelim.mix(stlouis,3)  # preliminary manipulations
> thetahat <- em.mix(s) # find ML estimate
Steps of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...101...102...103...104...105...106...107...108...109...110...111...112...113...114...115...116...117...118...119...120...121...122...123...124...125...126...127...128...129...130...131...132...133...134...135...136...137...138...139...140...141...142...143...144...145...146...147...148...149...150...151...152...153...154...155...156...157...158...159...160...161...162...163...164...165...166...167...168...169...170...171...172...173...174...175...176...177...178...179...180...181...
> rngseed(1234567)   # set random number generator seed
> newtheta <- da.mix(s, thetahat, steps=100, showits=TRUE)  # take 100 steps
Steps of Data Augmentation: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...
> ximp1 <- imp.mix(s, newtheta) # impute under newtheta
> 
> 
> 
> cleanEx(); ..nameEx <- "dabipf.mix"
> 
> ### * dabipf.mix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: dabipf.mix
> ### Title: Data Augmentation/Bayesian IPF Algorithm for Restricted General
> ###   Location Models
> ### Aliases: dabipf.mix
> ### Keywords: models
> 
> ### ** Examples
> 
> data(stlouis)
> s <- prelim.mix(stlouis,3)      # do preliminary manipulations
> margins <- c(1,2,3)       # saturated contingency table model
> design <- diag(rep(1,12))  # identity matrix  D=no of cells
> thetahat <- ecm.mix(s,margins,design) # find ML estimate
Steps of ECM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...101...102...103...104...105...106...107...108...109...110...111...112...113...114...115...116...117...118...119...120...121...122...123...124...125...126...127...128...129...130...131...132...133...134...135...136...137...138...139...140...141...142...143...144...145...146...147...148...149...150...151...152...153...154...155...156...157...158...159...160...161...162...163...164...165...166...167...168...169...170...171...172...173...174...175...176...177...178...179...180...181...
> rngseed(1234567)       # random generator seed
> newtheta <- dabipf.mix(s,margins,design,thetahat,steps=200)
> ximp <- imp.mix(s,newtheta,stlouis)   # impute under newtheta
> 
> 
> 
> cleanEx(); ..nameEx <- "ecm.mix"
> 
> ### * ecm.mix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ecm.mix
> ### Title: ECM Algorithm for Restricted General Location Model
> ### Aliases: ecm.mix
> ### Keywords: models
> 
> ### ** Examples
> 
> data(stlouis)
> s <- prelim.mix(stlouis,3)         # preliminary manipulations
> margins <- c(1,2,3)          # saturated loglinear model
> design <- diag(rep(1,12))     # identity matrix, D=no of cells
> thetahat <- ecm.mix(s,margins,design)  # should be same as em.mix(s)
Steps of ECM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...101...102...103...104...105...106...107...108...109...110...111...112...113...114...115...116...117...118...119...120...121...122...123...124...125...126...127...128...129...130...131...132...133...134...135...136...137...138...139...140...141...142...143...144...145...146...147...148...149...150...151...152...153...154...155...156...157...158...159...160...161...162...163...164...165...166...167...168...169...170...171...172...173...174...175...176...177...178...179...180...181...
> loglik.mix(s,thetahat)             # loglikelihood at thetahat
[1] -110.3436
> 
> 
> 
> cleanEx(); ..nameEx <- "em.mix"
> 
> ### * em.mix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: em.mix
> ### Title: EM Algorithm for Unrestricted General Location Model
> ### Aliases: em.mix
> ### Keywords: models
> 
> ### ** Examples
> 
> data(stlouis)
> s <- prelim.mix(stlouis,3)    # do preliminary manipulations
> thetahat <- em.mix(s)   # compute ML estimate
Steps of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...101...102...103...104...105...106...107...108...109...110...111...112...113...114...115...116...117...118...119...120...121...122...123...124...125...126...127...128...129...130...131...132...133...134...135...136...137...138...139...140...141...142...143...144...145...146...147...148...149...150...151...152...153...154...155...156...157...158...159...160...161...162...163...164...165...166...167...168...169...170...171...172...173...174...175...176...177...178...179...180...181...
> getparam.mix(s,thetahat, corr=TRUE) # look at estimated parameters
$pi
, , D2=1

          D1=1       D1=2
G=1 0.14040543 0.09654434
G=2 0.01627643 0.04617068
G=3 0.03103847 0.01452370

, , D2=2

         D1=1       D1=2
G=1 0.1358792 0.01847542
G=2 0.1076896 0.17768936
G=3 0.1119740 0.10333338


$mu
        [,1]     [,2]      [,3]     [,4]      [,5]      [,6]     [,7]     [,8]
R1 109.77629 77.22294 115.80413 110.7962  93.52564  56.09284 123.7264 104.1379
V1 130.76154 65.92582 137.70299 134.1615 113.81725  58.14789 161.0442 135.9854
R2  99.71291 89.54078  82.91913 103.9516 126.02480  88.08945 118.0722 109.4170
V2 116.75389 76.68685  96.22231 132.3583 142.64726 105.19174 132.6509 109.1568
       [,9]    [,10]    [,11]    [,12]
R1 105.6808 119.6342 106.9637 107.4552
V1 127.0207 141.4109 103.5737 107.6212
R2 100.7178 136.9676  97.0306 106.9518
V2 128.0930 181.0136 102.4299 104.9008

$sdv
      R1       V1       R2       V2 
13.02237 21.05411 10.67396 22.98537 

$r
          R1        V1        R2        V2
R1 1.0000000 0.8024177 0.6985010 0.8217321
V1 0.8024177 1.0000000 0.6803922 0.8186402
R2 0.6985010 0.6803922 1.0000000 0.7818385
V2 0.8217321 0.8186402 0.7818385 1.0000000

> 
> 
> 
> cleanEx(); ..nameEx <- "getparam.mix"
> 
> ### * getparam.mix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: getparam.mix
> ### Title: Present Parameters of General Location Model in an
> ###   Understandable Format
> ### Aliases: getparam.mix
> ### Keywords: models
> 
> ### ** Examples
> 
> data(stlouis)
> s <- prelim.mix(stlouis,3)    # do preliminary manipulations
> thetahat <- em.mix(s)   # compute ML estimate
Steps of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...101...102...103...104...105...106...107...108...109...110...111...112...113...114...115...116...117...118...119...120...121...122...123...124...125...126...127...128...129...130...131...132...133...134...135...136...137...138...139...140...141...142...143...144...145...146...147...148...149...150...151...152...153...154...155...156...157...158...159...160...161...162...163...164...165...166...167...168...169...170...171...172...173...174...175...176...177...178...179...180...181...
> getparam.mix(s, thetahat, corr=TRUE)$r   # look at estimated correlations 
          R1        V1        R2        V2
R1 1.0000000 0.8024177 0.6985010 0.8217321
V1 0.8024177 1.0000000 0.6803922 0.8186402
R2 0.6985010 0.6803922 1.0000000 0.7818385
V2 0.8217321 0.8186402 0.7818385 1.0000000
> 
> 
> 
> cleanEx(); ..nameEx <- "imp.mix"
> 
> ### * imp.mix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: imp.mix
> ### Title: Impute Missing Data Under General Location Model
> ### Aliases: imp.mix
> ### Keywords: models
> 
> ### ** Examples
> 
> data(stlouis)
> s <- prelim.mix(stlouis,3)    # do preliminary manipulations
> thetahat <- em.mix(s)   # ML estimate for unrestricted model
Steps of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...101...102...103...104...105...106...107...108...109...110...111...112...113...114...115...116...117...118...119...120...121...122...123...124...125...126...127...128...129...130...131...132...133...134...135...136...137...138...139...140...141...142...143...144...145...146...147...148...149...150...151...152...153...154...155...156...157...158...159...160...161...162...163...164...165...166...167...168...169...170...171...172...173...174...175...176...177...178...179...180...181...
> rngseed(1234567)     # set random number generator seed
> newtheta <- da.mix(s,thetahat,steps=100) # data augmentation
> ximp <- imp.mix(s, newtheta, stlouis)  # impute under newtheta
> 
> 
> 
> cleanEx(); ..nameEx <- "loglik.mix"
> 
> ### * loglik.mix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: loglik.mix
> ### Title: Loglikelihood for Incomplete Data under the General Location
> ###   Model
> ### Aliases: loglik.mix
> ### Keywords: models
> 
> ### ** Examples
> 
> data(stlouis)
> s <- prelim.mix(stlouis,3)       # preliminary manipulations
> thetahat <- em.mix(s)      # MLE under unrestricted general location model
Steps of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...22...23...24...25...26...27...28...29...30...31...32...33...34...35...36...37...38...39...40...41...42...43...44...45...46...47...48...49...50...51...52...53...54...55...56...57...58...59...60...61...62...63...64...65...66...67...68...69...70...71...72...73...74...75...76...77...78...79...80...81...82...83...84...85...86...87...88...89...90...91...92...93...94...95...96...97...98...99...100...101...102...103...104...105...106...107...108...109...110...111...112...113...114...115...116...117...118...119...120...121...122...123...124...125...126...127...128...129...130...131...132...133...134...135...136...137...138...139...140...141...142...143...144...145...146...147...148...149...150...151...152...153...154...155...156...157...158...159...160...161...162...163...164...165...166...167...168...169...170...171...172...173...174...175...176...177...178...179...180...181...
> loglik.mix(s, thetahat)  # loglikelihood at thetahat
[1] -110.3436
> 
> 
> 
> cleanEx(); ..nameEx <- "prelim.mix"
> 
> ### * prelim.mix
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: prelim.mix
> ### Title: Preliminary Manipulations on Matrix of Incomplete Mixed Data
> ### Aliases: prelim.mix
> ### Keywords: models
> 
> ### ** Examples
> 
> data(stlouis)
> s <- prelim.mix(stlouis, 3)  # do preliminary manipulations 
> s$nmis              # look at nmis
 G D1 D2 R1 V1 R2 V2 
 0 28 28 21 30 16 17 
> s$r                 # look at missing data patterns
   G D1 D2 R1 V1 R2 V2
12 1  1  1  1  1  1  1
1  1  0  1  1  1  1  1
8  1  1  0  1  1  1  1
6  1  0  0  1  1  1  1
1  1  1  1  0  1  1  1
1  1  0  1  0  1  1  1
1  1  1  0  0  1  1  1
3  1  0  0  0  1  1  1
3  1  1  1  1  0  1  1
1  1  0  0  1  0  1  1
1  1  1  1  0  0  1  1
6  1  0  1  0  0  1  1
2  1  0  0  0  0  1  1
4  1  1  1  1  1  0  1
1  1  0  1  1  0  0  1
1  1  0  0  0  0  0  1
3  1  1  1  1  0  1  0
1  1  1  1  0  0  1  0
2  1  0  1  0  0  1  0
1  1  0  0  0  0  1  0
1  1  1  0  1  1  0  0
1  1  0  0  1  1  0  0
4  1  1  1  1  0  0  0
2  1  1  0  1  0  0  0
1  1  0  0  1  0  0  0
1  1  0  1  0  0  0  0
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  1.79 0.08 1.86 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')