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("BMA-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('BMA')
Loading required package: survival
Loading required package: splines
Loading required package: leaps
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "MC3.REG"
> 
> ### * MC3.REG
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: MC3.REG
> ### Title: Bayesian simultaneous variable selection and outlier
> ###   identification
> ### Aliases: MC3.REG
> ### Keywords: regression models
> 
> ### ** Examples
> 
> # Example 1:   Scottish hill racing data.
> 
> data(race)
> b<- out.ltsreg(race[,-1], race[,1], 2)
Scalable Robust Estimators with High Breakdown Point (version 0.2-7)
> races.run1<-MC3.REG(race[,1], race[,-1], num.its=20000, c(FALSE,TRUE), rep(TRUE,length(b)), b, 
+    PI = .1, K = 7, nu = .2, lambda = .1684, phi = 9.2)
> races.run1

Call:
MC3.REG(all.y = race[, 1], all.x = race[, -1], num.its = 20000,     M0.var = c(FALSE, TRUE), M0.out = rep(TRUE, length(b)), outs.list = b,     PI = 0.1, K = 7, nu = 0.2, lambda = 0.1684, phi = 9.2)

Model parameters: PI = 0.1 K = 7 nu = 0.2 lambda = 0.1684 phi = 9.2

Models visited: 
   posterior n.visits variables outliers
1  8.159e-01    16499       1 1  0 1 1 1
2  9.105e-02     1739       1 1  0 1 1 0
3  8.956e-02     1696       1 1  1 1 1 1
4  3.192e-03       57       1 1  1 1 1 0
5  2.149e-04        4       0 1  0 1 1 1
6  3.617e-05        2       0 1  0 1 1 0
7  3.251e-05        3       0 1  1 1 1 1
8  1.559e-05        0       1 1  0 0 1 0
9  1.844e-06        0       0 1  1 1 1 0
10 3.443e-07        0       1 1  0 0 1 1
11 2.683e-07        0       1 1  1 0 1 0
12 6.039e-09        0       1 1  1 0 1 1
13 1.415e-09        0       0 1  0 1 0 1
14 8.655e-10        0       1 1  0 1 0 0
15 6.979e-11        0       1 1  0 1 0 1
16 4.441e-11        0       0 1  0 0 1 1
17 1.872e-11        0       1 1  1 1 0 0
18 1.865e-12        0       1 1  1 1 0 1
19 3.082e-23        0       1 0  0 1 1 0
20 6.143e-25        0       1 0  1 1 1 0
21 6.089e-25        0       1 0  0 1 1 1
22 1.625e-25        0       0 0  0 1 1 0
23 2.446e-26        0       0 0  0 1 1 1
24 1.192e-26        0       1 0  1 1 1 1
25 4.010e-28        0       0 0  1 1 1 1
> summary(races.run1)

Call:
MC3.REG(all.y = race[, 1], all.x = race[, -1], num.its = 20000,     M0.var = c(FALSE, TRUE), M0.out = rep(TRUE, length(b)), outs.list = b,     PI = 0.1, K = 7, nu = 0.2, lambda = 0.1684, phi = 9.2)

Model parameters: PI = 0.1 K = 7 nu = 0.2 lambda = 0.1684 phi = 9.2

  25  models were selected
 Best  5  models (cumulative posterior probability =  0.9999 ): 

           prob     model 1    model 2    model 3    model 4    model 5  
variables                                                                
  Climb    0.9997    x          x          x          x          .       
  Time     1.0000    x          x          x          x          x       
outliers                                                                 
  6        0.09278   .          .          x          x          .       
  7        0.99998   x          x          x          x          x       
  18       1.00000   x          x          x          x          x       
  33       0.90570   x          .          x          .          x       
                                                                         
post prob           0.8158986  0.0910511  0.0895564  0.0031922  0.0002149
> 
> # Example 2: Crime data
> library(MASS)
> data(UScrime)
> 
> y.crime.log<- log(UScrime$y)
> x.crime.log<- UScrime[,-ncol(UScrime)]
> x.crime.log[,-2]<- log(x.crime.log[,-2])
> crime.run1<-MC3.REG(y.crime.log, x.crime.log, num.its=2000, rep(TRUE,15), outliers = FALSE)
> crime.run1[1:25,]
     post.prob visit.count M So Ed Po1 Po2 LF M.F Pop NW U1 U2 GDP Ineq Prob
1  0.031524044          50 1  0  1   1   0  0   0   0  1  0  1   0    1    1
2  0.022454627          52 1  0  1   1   0  0   0   0  0  0  1   0    1    1
3  0.020890613          37 1  0  1   1   0  0   0   0  1  0  0   0    1    1
4  0.019953547          15 1  0  1   0   1  0   0   0  1  0  1   0    1    1
5  0.019314528          16 1  0  1   1   0  0   0   0  1  0  1   0    1    1
6  0.019295407          19 1  0  1   1   0  0   0   0  1  0  0   0    1    1
7  0.018984363          22 0  0  1   1   0  0   0   1  1  0  0   0    1    1
8  0.015805206          22 1  0  1   1   0  0   0   0  0  0  0   0    1    1
9  0.014625384          11 1  0  1   1   0  0   0   0  0  0  1   0    1    0
10 0.014307260          20 1  0  1   0   1  0   0   0  0  0  1   0    1    1
11 0.014023784          17 1  0  1   1   1  0   0   0  1  0  1   0    1    1
12 0.012190440           1 1  0  1   1   0  0   0   0  0  0  0   0    1    0
13 0.012087408          25 1  0  1   1   0  0   0   0  1  0  1   1    1    1
14 0.011954754          26 1  0  1   0   1  0   0   0  1  0  0   0    1    1
15 0.011564142          32 1  0  1   1   0  0   0   1  1  0  0   0    1    1
16 0.011208618          15 1  0  1   1   1  0   0   0  0  0  1   0    1    1
17 0.011126362           1 0  0  1   0   1  0   0   1  1  0  0   0    1    1
18 0.010887126          19 1  0  1   1   0  0   0   0  1  0  0   1    1    1
19 0.010505602          17 1  0  1   0   1  0   0   0  1  0  1   1    1    1
20 0.010231427          33 0  0  1   1   0  0   0   0  1  0  0   0    1    1
21 0.009406954          11 1  0  1   1   0  0   0   0  1  0  1   1    1    1
22 0.009097504          41 1  0  1   1   0  0   0   1  1  0  1   0    1    1
23 0.009096343          14 1  0  1   1   0  0   0   0  1  0  0   1    1    1
24 0.008987097          18 1  0  1   1   1  0   0   0  1  0  0   0    1    1
25 0.008882828           2 0  0  1   0   1  0   0   0  1  0  0   0    1    1
   Time
1     0
2     0
3     0
4     0
5     1
6     1
7     0
8     0
9     0
10    0
11    0
12    0
13    0
14    0
15    0
16    0
17    0
18    1
19    0
20    0
21    1
22    0
23    0
24    0
25    0
> summary(crime.run1)

Call:
MC3.REG(all.y = y.crime.log, all.x = x.crime.log, num.its = 2000,     M0.var = rep(TRUE, 15), outliers = FALSE)

Model parameters: PI = 0.1 K = 7 nu = 2.58 lambda = 0.28 phi = 2.85

  967  models were selected
 Best  5  models (cumulative posterior probability =  0.1141 ): 

           prob     model 1  model 2  model 3  model 4  model 5
variables                                                      
  M        0.81448   x        x        x        x        x     
  So       0.14129   .        .        .        .        .     
  Ed       0.99029   x        x        x        x        x     
  Po1      0.76383   x        x        x        .        x     
  Po2      0.44298   .        .        .        x        .     
  LF       0.04621   .        .        .        .        .     
  M.F      0.04939   .        .        .        .        .     
  Pop      0.22944   .        .        .        .        .     
  NW       0.66053   x        .        x        x        x     
  U1       0.08126   .        .        .        .        .     
  U2       0.46491   x        x        .        x        x     
  GDP      0.28059   .        .        .        .        .     
  Ineq     0.99917   x        x        x        x        x     
  Prob     0.87483   x        x        x        x        x     
  Time     0.22291   .        .        .        .        x     
                                                               
post prob           0.03152  0.02245  0.02089  0.01995  0.01931
> 
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "as.data.frame.mc3"
> 
> ### * as.data.frame.mc3
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: as.data.frame.mc3
> ### Title: Type conversion and indexing for mc3 object
> ### Aliases: as.data.frame.mc3 [.mc3
> ### Keywords: methods
> 
> ### ** Examples
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "bic.glm"
> 
> ### * bic.glm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bic.glm
> ### Title: Bayesian Model Averaging for generalized linear models.
> ### Aliases: bic.glm bic.glm.data.frame bic.glm.matrix bic.glm.formula
> ### Keywords: regression models
> 
> ### ** Examples
> 
> 
> ### logistic regression
> library("MASS")
> data(birthwt)
> y<- birthwt$lo
> x<- data.frame(birthwt[,-1])
> x$race<- as.factor(x$race)
> x$ht<- (x$ht>=1)+0
> x<- x[,-9]
> x$smoke <- as.factor(x$smoke)
> x$ptl<- as.factor(x$ptl)
> x$ht <- as.factor(x$ht)
> x$ui <- as.factor(x$ui)
> 
> glm.out.FT<- bic.glm(x, y, strict = FALSE, OR = 20, glm.family="binomial", factor.type=TRUE)
> summary(glm.out.FT)

Call:
bic.glm.data.frame(x = x, y = y, glm.family = "binomial", strict = FALSE,     OR = 20, factor.type = TRUE)


  49  models were selected
 Best  5  models (cumulative posterior probability =  0.3676 ): 

           p!=0   EV        SD         model 1     model 2     model 3   
Intercept  100    0.471634  1.309e+00     1.45068     1.06795     1.20695
age        12.6  -0.007765  2.433e-02       .           .           .    
lwt        68.7  -0.011558  9.703e-03    -0.01865    -0.01692    -0.01881
race        9.6                                                          
    .2            0.115279  3.888e-01       .           .           .    
    .3            0.092740  3.160e-01       .           .           .    
smoke      33.2                                                          
     .1           0.255357  4.281e-01       .           .           .    
ptl        35.1                                                          
   .1             0.617360  8.898e-01       .           .         1.74321
   .2             0.168588  6.150e-01       .           .         0.50107
   .3            -4.910963  5.231e+02       .           .       -13.98561
ht         65.4                                                          
  .1              1.166893  1.034e+00     1.85551     1.96157     1.92368
ui         33.8                                                          
  .1              0.310537  5.057e-01       .         0.93000       .    
ftv         1.5  -0.001339  2.393e-02       .           .           .    
                                                                         
nVar                                        2           3           3    
BIC                                    -753.82285  -753.11035  -752.99778
post prob                                 0.111       0.078       0.073  
           model 4     model 5   
Intercept     1.08354     0.72186
age             .           .    
lwt          -0.01805    -0.01631
race                             
    .2          .           .    
    .3          .           .    
smoke                            
     .1       0.68391     0.65300
ptl                              
   .1           .           .    
   .2           .           .    
   .3           .           .    
ht                               
  .1          1.82202     1.92213
ui                               
  .1            .         0.89627
ftv             .           .    
                                 
nVar            3           4    
BIC        -752.86548  -751.65571
post prob     0.069       0.037  
> imageplot.bma(glm.out.FT)
> 
> glm.out.FF<- bic.glm(x, y, strict = FALSE, OR = 20, glm.family="binomial", factor.type=FALSE)
> summary(glm.out.FF)

Call:
bic.glm.data.frame(x = x, y = y, glm.family = "binomial", strict = FALSE,     OR = 20, factor.type = FALSE)


  35  models were selected
 Best  5  models (cumulative posterior probability =  0.4045 ): 

           p!=0    EV        SD         model 1     model 2     model 3   
Intercept  100     0.652433    1.30285     1.19274     0.72820    -1.04282
age         25.7  -0.017415    0.03480       .           .           .    
lwt         66.5  -0.011511    0.01008    -0.01864    -0.01391       .    
race2       12.9   0.107422    0.33131       .           .           .    
race3        3.2   0.012579    0.10145       .           .           .    
smoke1      14.1   0.078530    0.23451       .           .           .    
ptl1       100.0   1.749160    0.49542     1.73602     1.71252     1.73597
ptl2         1.4   0.007156    0.12710       .           .           .    
ptl3         2.8  -0.397860  147.54575       .           .           .    
ht1         63.3   1.139602    1.05453     1.91085       .           .    
ui1         22.3   0.186405    0.40835       .           .           .    
ftv          1.5  -0.001610    0.02498       .           .           .    
                                                                          
nVar                                         3           2           1    
BIC                                     -762.31510  -760.33661  -760.33635
post prob                                  0.164       0.061       0.061  
           model 4     model 5   
Intercept     0.84764     0.65569
age             .        -0.07545
lwt          -0.01699       .    
race2           .           .    
race3           .           .    
smoke1          .           .    
ptl1          1.66642     1.93649
ptl2            .           .    
ptl3            .           .    
ht1           1.99607       .    
ui1           0.81562       .    
ftv             .           .    
                                 
nVar            4           2    
BIC        -760.30767  -760.23885
post prob     0.060       0.058  
> imageplot.bma(glm.out.FF)
> 
> glm.out.TT<- bic.glm(x, y, strict = TRUE, OR = 20, glm.family="binomial", factor.type=TRUE)
> summary(glm.out.TT)

Call:
bic.glm.data.frame(x = x, y = y, glm.family = "binomial", strict = TRUE,     OR = 20, factor.type = TRUE)


  4  models were selected
 Best  4  models (cumulative posterior probability =  1 ): 

           p!=0   EV       SD         model 1     model 2     model 3   
Intercept  100    0.76092    1.23000     1.45068     0.99831    -1.05711
age         0.0   0.00000    0.00000       .           .           .    
lwt        74.4  -0.01306    0.00964    -0.01865    -0.01406       .    
race        0.0                                                         
    .2            0.00000    0.00000       .           .           .    
    .3            0.00000    0.00000       .           .           .    
smoke       0.0                                                         
     .1           0.00000    0.00000       .           .           .    
ptl        13.3                                                         
   .1             0.23225    0.61790       .           .         1.75026
   .2             0.08647    0.40473       .           .         0.65165
   .3            -1.79255  321.59042       .           .       -13.50896
ht         56.6                                                         
  .1              1.04938    1.06011     1.85551       .           .    
ui          0.0                                                         
  .1              0.00000    0.00000       .           .           .    
ftv         0.0   0.00000    0.00000       .           .           .    
                                                                        
nVar                                       2           1           1    
BIC                                   -753.82285  -751.51602  -750.92334
post prob                                0.566       0.178       0.133  
           model 4   
Intercept    -0.79000
age             .    
lwt             .    
race                 
    .2          .    
    .3          .    
smoke                
     .1         .    
ptl                  
   .1           .    
   .2           .    
   .3           .    
ht                   
  .1            .    
ui                   
  .1            .    
ftv             .    
                     
nVar            0    
BIC        -750.77644
post prob     0.123  
> imageplot.bma(glm.out.TT)
> 
> glm.out.TF<- bic.glm(x, y, strict = TRUE, OR = 20, glm.family="binomial", factor.type=FALSE)
> summary(glm.out.TF)

Call:
bic.glm.data.frame(x = x, y = y, glm.family = "binomial", strict = TRUE,     OR = 20, factor.type = FALSE)


  3  models were selected
 Best  3  models (cumulative posterior probability =  1 ): 

           p!=0    EV       SD        model 1     model 2     model 3   
Intercept  100     0.61697  1.165455     1.19274     0.72820    -1.04282
age          0.0   0.00000  0.000000       .           .           .    
lwt         78.7  -0.01366  0.009516    -0.01864    -0.01391       .    
race2        0.0   0.00000  0.000000       .           .           .    
race3        0.0   0.00000  0.000000       .           .           .    
smoke1       0.0   0.00000  0.000000       .           .           .    
ptl1       100.0   1.73100  0.479102     1.73602     1.71252     1.73597
ptl2         0.0   0.00000  0.000000       .           .           .    
ptl3         0.0   0.00000  0.000000       .           .           .    
ht1         57.4   1.09588  1.096825     1.91085       .           .    
ui1          0.0   0.00000  0.000000       .           .           .    
ftv          0.0   0.00000  0.000000       .           .           .    
                                                                        
nVar                                       3           2           1    
BIC                                   -762.31510  -760.33661  -760.33635
post prob                                0.574       0.213       0.213  
> imageplot.bma(glm.out.TF)
> 
> 
> ### Gamma family 
> library(survival)
> data(veteran)
> surv.t<- veteran$time
> x<- veteran[,-c(3,4)]
> x$celltype<- factor(as.character(x$celltype))
> sel<- veteran$status == 0
> x<- x[!sel,]
> surv.t<- surv.t[!sel]
> 
> glm.out.va<- bic.glm(x, y=surv.t, glm.family=Gamma(link="inverse") ,factor.type=FALSE)
> summary(glm.out.va)

Call:
bic.glm.data.frame(x = x, y = surv.t, glm.family = Gamma(link = "inverse"),     factor.type = FALSE)


  9  models were selected
 Best  5  models (cumulative posterior probability =  0.7968 ): 

                   p!=0    EV         SD         model 1     model 2   
Intercept          100     2.598e-02  3.915e-03   2.641e-02   2.470e-02
trt                  5.6   4.851e-05  4.071e-04       .           .    
celltypelarge       81.2  -5.428e-03  3.425e-03  -6.631e-03       .    
celltypesmallcell   10.9   1.723e-04  1.681e-03       .           .    
celltypesquamous    89.9  -6.411e-03  3.259e-03  -7.431e-03  -3.580e-03
karno              100.0  -2.020e-04  4.984e-05  -1.939e-04  -2.193e-04
diagtime             4.8   8.172e-07  1.948e-05       .           .    
age                  5.5   2.141e-06  1.883e-05       .           .    
prior                6.2  -5.479e-06  3.625e-05       .           .    
                                                                       
nVar                                                3           2      
BIC                                              -4.969e+02  -4.932e+02
post prob                                         0.537       0.087    
                   model 3     model 4     model 5   
Intercept           2.681e-02   2.550e-02   2.308e-02
trt                     .       8.681e-04       .    
celltypelarge      -6.483e-03  -6.412e-03       .    
celltypesmallcell       .           .       5.052e-03
celltypesquamous   -7.171e-03  -7.679e-03       .    
karno              -1.975e-04  -1.993e-04  -2.323e-04
diagtime                .           .           .    
age                     .           .           .    
prior              -8.846e-05       .           .    
                                                     
nVar                  4           4           2      
BIC                -4.925e+02  -4.923e+02  -4.923e+02
post prob           0.062       0.056       0.055    
> imageplot.bma(glm.out.va)
> plot(glm.out.va)
Hit <Return> to see next plot: 
> 
> ### Poisson family
> ### Yates (teeth) data. 
> 
> x<- rbind(
+     c(0, 0, 0),
+     c(0, 1, 0),
+     c(1, 0, 0),
+     c(1, 1, 1))
> 
> y<-c(4, 16, 1, 21)
> n<-c(1,1,1,1)
> 
> models<- rbind(
+     c(1, 1, 0),
+     c(1, 1, 1))
> 
> glm.out.yates <- bic.glm(x,y,n, glm.family = poisson(), factor.type=FALSE) 
> summary(glm.out.yates)

Call:
bic.glm.matrix(x = x, y = y, glm.family = poisson(), wt = n,     factor.type = FALSE)


  4  models were selected
 Best  4  models (cumulative posterior probability =  1 ): 

           p!=0    EV      SD      model 1   model 2   model 3   model 4 
Intercept  100     1.0456  0.5170   0.91629   1.38629   0.91629   0.86750
X1          45.9  -0.3894  0.8906      .     -1.38629      .      0.09531
X2         100.0   1.7892  0.5747   2.00148   1.38629   1.85630   2.00148
X3          51.5   0.5455  0.9722      .      1.65823   0.27193      .   
                                                                         
nVar                                  1         3         2         2    
BIC                                -0.16739   0.00000   0.54115   1.12363
post prob                           0.318     0.292     0.223     0.167  
> 
> ### Gaussian
> library(MASS)
> data(UScrime)
> f<- formula(log(y) ~  log(M)+So+log(Ed)+log(Po1)+log(Po2)+log(LF)+log(M.F)+ 
+    log(Pop)+log(NW)+log(U1)+log(U2)+log(GDP)+log(Ineq)+log(Prob)+log(Time) )
> glm.out.crime<- bic.glm(f, data = UScrime, glm.family = gaussian()) 
> summary(glm.out.crime)

Call:
bic.glm.formula(f = f, data = UScrime, glm.family = gaussian())


  137  models were selected
 Best  5  models (cumulative posterior probability =  0.1849 ): 

           p!=0    EV        SD       model 1     model 2     model 3   
Intercept  100    -23.33177  5.57356   -24.38362   -22.63715   -24.50477
log.M.      95.1    1.33340  0.56520     1.51437     1.47803     1.46061
So          12.2    0.01549  0.05821       .           .           .    
log.Ed.    100.0    2.08869  0.51718     2.38935     2.22117     2.39875
log.Po1.    68.9    0.64171  0.48337     0.91047     0.85244       .    
log.Po2.    34.6    0.28819  0.46013       .           .         0.90689
log.LF.      6.5    0.01945  0.16701       .           .           .    
log.M.F.     5.9   -0.02250  0.36877       .           .           .    
log.Pop.    28.2   -0.01806  0.03576       .           .           .    
log.NW.     81.1    0.08124  0.05376     0.08456     0.10888     0.08534
log.U1.     13.1   -0.02026  0.13174       .           .           .    
log.U2.     72.5    0.24167  0.20227     0.32169     0.28874     0.32977
log.GDP.    26.0    0.15105  0.32081       .           .           .    
log.Ineq.  100.0    1.39354  0.32754     1.23088     1.23775     1.29370
log.Prob.   98.4   -0.23727  0.09784    -0.19062    -0.31040    -0.20614
log.Time.   33.3   -0.09438  0.16214       .        -0.28659       .    
                                                                        
nVar                                       7           8           7    
BIC                                   -108.92725  -108.75890  -108.07840
post prob                                0.051       0.047       0.033  
           model 4     model 5   
Intercept   -22.80644   -26.29204
log.M.        1.26830     1.68239
So              .           .    
log.Ed.       2.17788     2.08662
log.Po1.      0.98597     1.14936
log.Po2.        .           .    
log.LF.         .           .    
log.M.F.        .           .    
log.Pop.     -0.05685       .    
log.NW.       0.09745       .    
log.U1.         .           .    
log.U2.       0.28054     0.33564
log.GDP.        .           .    
log.Ineq.     1.32157     1.57576
log.Prob.    -0.21636    -0.16699
log.Time.       .           .    
                                 
nVar            8           6    
BIC        -107.70043  -107.56949
post prob     0.028       0.026  
> # note the problems with the estimation of the posterior standard deviation 
> # (compare with bicreg example)
> 
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "bic.surv"
> 
> ### * bic.surv
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bic.surv
> ### Title: Bayesian Model Averaging for Survival models.
> ### Aliases: bic.surv bic.surv.data.frame bic.surv.matrix bic.surv.formula
> ### Keywords: regression survival
> 
> ### ** Examples
> 
> 
> ## veteran data
> library(survival)
> data(veteran)
> 
> test.bic.surv<- bic.surv(Surv(time,status) ~ ., data = veteran, factor.type = TRUE)
> summary(test.bic.surv, conditional=FALSE, digits=2)

Call:
bic.surv.formula(f = Surv(time, status) ~ ., data = veteran,     factor.type = TRUE)


  6  models were selected
 Best  5  models (cumulative posterior probability =  0.95 ): 

                    p!=0    EV       SD      model 1   model 2   model 3 
trt                  11.3   0.02909  0.1058      .         .       0.2573
celltype             84.6                                                
        .smallcell          0.61564  0.3538    0.7121      .       0.8196
        .adeno              0.97619  0.4965    1.1508      .       1.1477
        .large              0.28260  0.2830    0.3251      .       0.3930
karno               100.0  -0.03135  0.0052   -0.0309   -0.0332   -0.0311
diagtime              5.4   0.00018  0.0020      .         .         .   
age                   6.1  -0.00036  0.0026      .         .         .   
prior                 5.6   0.00058  0.0054      .         .         .   
                                                                         
nVar                                             2         1         3   
BIC                                          -39.3626  -36.7742  -36.1558
post prob                                      0.562     0.154     0.113 
                    model 4   model 5 
trt                     .         .   
celltype                              
        .smallcell    0.7208    0.7264
        .adeno        1.1643    1.1765
        .large        0.3215    0.3276
karno                -0.0318   -0.0311
diagtime                .         .   
age                  -0.0059      .   
prior                   .       0.0103
                                      
nVar                    3         3   
BIC                 -34.9292  -34.7553
post prob             0.061     0.056 
> plot(test.bic.surv)
Hit <Return> to see next plot: 
> imageplot.bma(test.bic.surv)
> 
> 
> ## pbc data
> data(pbc)
> x<- pbc[1:312,]
> surv.t<- x$time
> cens<- x$status
> x<- x[,-c(6,7,10,17,19)]
> 
> x$bili<- log(x$bili)
> x$alb<- log(x$alb)
> x$protime<- log(x$protime)
> x$copper<- log(x$copper)
Warning in log(x) : NaNs produced
> x$sgot<- log(x$sgot)
> 
> test.bic.surv<- bic.surv(x, surv.t, cens, factor.type=FALSE, strict=FALSE)
Warning in bic.surv.data.frame(x, surv.t, cens, factor.type = FALSE, strict = FALSE) : 
	 There were  2 records deleted due to NA's
> summary(test.bic.surv)

Call:
bic.surv.data.frame(x = x, surv.t = surv.t, cens = cens, strict = FALSE,     factor.type = FALSE)


  39  models were selected
 Best  5  models (cumulative posterior probability =  0.3738 ): 

           p!=0    EV         SD         model 1     model 2     model 3   
age        100.0   3.234e-02  9.177e-03     0.03220     0.02943     0.03066
alb         99.2  -2.750e+00  8.332e-01    -2.81577    -2.44611    -2.38567
alkphos      8.1  -3.719e-06  1.620e-05       .           .           .    
ascites      1.6   2.634e-03  4.282e-02       .           .           .    
bili       100.0   7.765e-01  1.287e-01     0.76060     0.74335     0.78602
edtrt       82.6   7.201e-01  4.385e-01     0.81318     0.82100     1.02642
hepmeg       2.0   4.108e-03  4.259e-02       .           .           .    
platelet     5.4  -5.452e-05  3.194e-04       .           .           .    
protime     80.2   2.565e+00  1.640e+00     3.10698     2.70572       .    
sex          3.4  -8.029e-03  6.846e-02       .           .           .    
sgot        25.2   1.236e-01  2.515e-01       .           .           .    
spiders      1.3   3.737e-04  2.487e-02       .           .           .    
stage       39.7   1.147e-01  1.680e-01       .         0.24386     0.31754
trt          1.6   1.935e-03  2.818e-02       .           .           .    
copper      74.2   2.598e-01  1.953e-01     0.35765     0.34732     0.34781
                                                                           
nVar                                          6           7           6    
BIC                                      -177.95246  -176.44209  -176.21202
post prob                                   0.147       0.069       0.062  
           model 4     model 5   
age           0.03079     0.03595
alb          -3.26805    -2.78021
alkphos         .           .    
ascites         .           .    
bili          0.80120     0.68450
edtrt           .         0.83680
hepmeg          .           .    
platelet        .           .    
protime       3.76955     3.41652
sex             .           .    
sgot            .         0.40676
spiders         .           .    
stage           .           .    
trt             .           .    
copper        0.35676     0.31067
                                 
nVar            5           7    
BIC        -175.84950  -175.53683
post prob     0.051       0.044  
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "bicreg"
> 
> ### * bicreg
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bicreg
> ### Title: Bayesian Model Averaging for linear regression models.
> ### Aliases: bicreg
> ### Keywords: regression models
> 
> ### ** Examples
> 
> library(MASS)
> data(UScrime)
> x<- UScrime[,-16]
> y<- log(UScrime[,16])
> x[,-2]<- log(x[,-2])
> lma<- bicreg(x, y, strict = FALSE, OR = 20) 
> summary(lma)

Call:
bicreg(x = x, y = y, strict = FALSE, OR = 20)


  51  models were selected
 Best  5  models (cumulative posterior probability =  0.2922 ): 

           p!=0    EV         SD       model 1    model 2    model 3  
Intercept  100.0  -23.477820  5.46307  -22.63715  -24.38362  -25.94554
M           97.5    1.401690  0.53084    1.47803    1.51437    1.60455
So           6.3    0.006928  0.03921      .          .          .    
Ed         100.0    2.128169  0.51346    2.22117    2.38935    1.99973
Po1         75.4    0.670723  0.42310    0.85244    0.91047    0.73577
Po2         24.6    0.218588  0.39634      .          .          .    
LF           2.1    0.003718  0.08012      .          .          .    
M.F          5.2   -0.069138  0.44577      .          .          .    
Pop         28.9   -0.018279  0.03609      .          .          .    
NW          89.6    0.091067  0.04961    0.10888    0.08456    0.11191
U1          11.9   -0.036411  0.13634      .          .          .    
U2          83.9    0.274023  0.19116    0.28874    0.32169    0.27422
GDP         33.8    0.197110  0.35442      .          .        0.54105
Ineq       100.0    1.381005  0.33320    1.23775    1.23088    1.41942
Prob        98.8   -0.248303  0.10045   -0.31040   -0.19062   -0.29989
Time        44.5   -0.128942  0.17758   -0.28659      .       -0.29682
                                                                      
nVar                                       8          7          9    
r2                                       0.842      0.826      0.851  
BIC                                    -55.91243  -55.36499  -54.69225
post prob                                0.089      0.067      0.048  
           model 4    model 5  
Intercept  -22.80644  -24.50477
M            1.26830    1.46061
So             .          .    
Ed           2.17788    2.39875
Po1          0.98597      .    
Po2            .        0.90689
LF             .          .    
M.F            .          .    
Pop         -0.05685      .    
NW           0.09745    0.08534
U1             .          .    
U2           0.28054    0.32977
GDP            .          .    
Ineq         1.32157    1.29370
Prob        -0.21636   -0.20614
Time           .          .    
                               
nVar           8          7    
r2           0.838      0.823  
BIC        -54.60434  -54.40788
post prob    0.046      0.042  
> plot(lma)
Hit <Return> to see next plot: 
> imageplot.bma(lma)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "glib"
> 
> ### * glib
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: glib
> ### Title: Model uncertainty in generalized linear models using Bayes
> ###   factors
> ### Aliases: glib glib.matrix glib.data.frame glib.bic.glm as.bic.glm
> ###   as.bic.glm.glib
> ### Keywords: regression models
> 
> ### ** Examples
> 
> 
> ### Finney data
> library(forward)
Loading required package: lqs
Warning: package 'lqs' has been moved back to package 'MASS'
Warning: package 'MASS' has now been loaded
> data(vaso)
> x<- vaso[,1:2]
> y<- vaso[,3]
> n<- rep(1,times=length(y))
> 
> finney.models<- rbind(
+     c(1, 0),
+     c(0, 1),
+     c(1, 1))
> 
> finney.glib <- glib (x,y,n, error="binomial", link="logit", models=finney.models, 
+    glimvar=TRUE, output.priorvar=TRUE, output.postvar=TRUE)
> summary(finney.glib)

Call:
glib.data.frame(x = x, y = y, n = n, error = "binomial", link = "logit",     models = finney.models, glimvar = TRUE, output.priorvar = TRUE,     output.postvar = TRUE)


  3  models were selected
 Best  3  models (cumulative posterior probability =  1 ): 

           p!=0    EV    SD      model 1   model 2   model 3 
Intercept  100                                               
volume     100.0  3.557  1.3740   1.32340      .      3.55840
rate        99.9  2.440  0.8791      .      0.83401   2.43997
                                                             
nVar                                1         1         2    
BIC                               2.48400  -0.01913  17.26707
post prob                         0.001     0.000     0.999  
> 
> finney.bic.glm<- as.bic.glm(finney.glib)
> plot(finney.bic.glm,mfrow=c(2,1))
> 
> ### Yates (teeth) data. 
> 
> x<- rbind(
+     c(0, 0, 0),
+     c(0, 1, 0),
+     c(1, 0, 0),
+     c(1, 1, 1))
> 
> y<-c(4, 16, 1, 21)
> n<-c(1,1,1,1)
> 
> models<- rbind(
+     c(1, 1, 0),
+     c(1, 1, 1))
> 
> glib.yates <- glib (x,y,n,models=models,glimvar=TRUE,output.priorvar=TRUE,
+    output.postvar=TRUE) 
> summary(glib.yates)

Call:
glib.matrix(x = x, y = y, n = n, models = models, glimvar = TRUE,     output.priorvar = TRUE, output.postvar = TRUE)


  2  models were selected
 Best  2  models (cumulative posterior probability =  1 ): 

           p!=0    EV      SD      model 1   model 2 
Intercept  100                                       
X1         100.0  -0.2593  0.6435   0.08995  -0.50413
X2         100.0   1.6540  0.5263   1.88885   1.48933
X3          58.8   0.7489  0.7306      .      0.74887
                                                     
nVar                                  2         3    
BIC                                20.59868  21.30901
post prob                           0.412     0.588  
> 
> ### logistic regression with no models specified
> library("MASS")
> data(birthwt)
> y<- birthwt$lo
> x<- data.frame(birthwt[,-1])
> x$race<- as.factor(x$race)
> x$ht<- (x$ht>=1)+0
> x<- x[,-9]
> x$smoke <- as.factor(x$smoke)
> x$ptl<- as.factor(x$ptl)
> x$ht <- as.factor(x$ht)
> x$ui <- as.factor(x$ui)
> 
> glib.birthwt<- glib(x,y, error="binomial", link = "logit")
> summary(glib.birthwt)

Call:
glib.data.frame(x = x, y = y, error = "binomial", link = "logit")


  256  models were selected
 Best  5  models (cumulative posterior probability =  0.3059 ): 

           p!=0   EV       SD        model 1   model 2   model 3   model 4 
Intercept  100                                                             
age        26.8  -0.06662  0.036059      .         .         .         .   
lwt        63.2  -0.01688  0.007272  -0.01849  -0.01384      .     -0.01686
race2      19.2   0.82273  0.517009      .         .         .         .   
race3      10.7   0.49674  0.448753      .         .         .         .   
smoke1     19.6   0.60547  0.374913      .         .         .         .   
ptl1       98.7   1.73491  0.495829   1.72747   1.70520   1.72890   1.65840
ptl2        5.4   0.53154  0.954573      .         .         .         .   
ptl3        4.4  -5.00251  9.832010      .         .         .         .   
ht1        59.1   1.75873  0.750614   1.89493      .         .      1.98006
ui1        24.6   0.82924  0.452410      .         .         .      0.81246
ftv         6.0  -0.11523  0.173588      .         .         .         .   
                                                                           
nVar                                    3         2         1         4    
BIC                                  11.38862   9.31084   9.32564   9.17692
post prob                             0.130     0.046     0.046     0.043  
           model 5 
Intercept          
age        -0.07499
lwt            .   
race2          .   
race3          .   
smoke1         .   
ptl1        1.92673
ptl2           .   
ptl3           .   
ht1            .   
ui1            .   
ftv            .   
                   
nVar          2    
BIC         9.03798
post prob   0.040  
> 
> glm.birthwt<- as.bic.glm(glib.birthwt)
> 
> imageplot.bma(glm.birthwt)
> 
> plot(glm.birthwt)
Hit <Return> to see next plot: 
> 
> 
> 
> cleanEx(); ..nameEx <- "imageplot"
> 
> ### * imageplot
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: imageplot.bma
> ### Title: Images of models used in Bayesian model averaging
> ### Aliases: imageplot.bma
> ### Keywords: regression models
> 
> ### ** Examples
> 
> # logistic regression using bic.glm
> library("MASS")
> data(birthwt)
> y<- birthwt$lo
> x<- data.frame(birthwt[,-1])
> x$race<- as.factor(x$race)
> x$ht<- (x$ht>=1)+0
> x<- x[,-9]
> x$smoke <- as.factor(x$smoke)
> x$ptl<- as.factor(x$ptl)
> x$ht <- as.factor(x$ht)
> x$ui <- as.factor(x$ui)
> 
> glm.out1<- bic.glm(x, y, strict = TRUE, OR = 20, glm.family="binomial")
> imageplot.bma(glm.out1)
> 
> 
> # logistic regression using glib
> library("MASS")
> data(birthwt)
> y<- birthwt$lo
> x<- data.frame(birthwt[,-1])
> x$race<- as.factor(x$race)
> x$ht<- (x$ht>=1)+0
> x<- x[,-9]
> x$smoke <- as.factor(x$smoke)
> x$ptl<- as.factor(x$ptl)
> x$ht <- as.factor(x$ht)
> x$ui <- as.factor(x$ui)
> 
> glib.birthwt<- glib(x,y, error="binomial", link = "logit")
> glm.birthwt<- as.bic.glm(glib.birthwt)
> imageplot.bma(glm.birthwt)
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "plot.bic"
> 
> ### * plot.bic
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: plot.bicreg
> ### Title: Plots the posterior distributions of coefficients derived from
> ###   Bayesian model averaging
> ### Aliases: plot.bicreg plot.bic.glm plot.bic.surv plot
> ### Keywords: regression models
> 
> ### ** Examples
> 
> library(MASS)
> data(UScrime)
> x<- UScrime[,-16]
> y<- log(UScrime[,16])
> x[,-2]<- log(x[,-2])
> plot( bicreg(x, y)) 
Hit <Return> to see next plot: 
> 
> 
> cleanEx(); ..nameEx <- "summary.bic"
> 
> ### * summary.bic
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: summary.bic
> ### Title: Summaries of Bayesian model averaging objects
> ### Aliases: summary summary.bicreg summary.bic.glm summary.bic.surv
> ###   summary.glib summary.mc3 print print.bicreg print.bic.glm
> ###   print.bic.surv print.glib print.mc3
> ### Keywords: print
> 
> ### ** Examples
> 
> # logistic regression
> library("MASS")
> data(birthwt)
> y<- birthwt$lo
> x<- data.frame(birthwt[,-1])
> x$race<- as.factor(x$race)
> x$ht<- (x$ht>=1)+0
> x<- x[,-9]
> x$smoke <- as.factor(x$smoke)
> x$ptl<- as.factor(x$ptl)
> x$ht <- as.factor(x$ht)
> x$ui <- as.factor(x$ui)
> 
> glm.out1<- bic.glm(x, y, OR = 20, glm.family="binomial", factor.type=TRUE)
> summary(glm.out1, conditional = TRUE)

Call:
bic.glm.data.frame(x = x, y = y, glm.family = "binomial", OR = 20,     factor.type = TRUE)


  49  models were selected
 Best  5  models (cumulative posterior probability =  0.3676 ): 

           p!=0   EV        SD         cond EV    cond SD    model 1   
Intercept  100    0.471634  1.309e+00    0.47163  1.309e+00     1.45068
age        12.6  -0.007765  2.433e-02   -0.06168  3.707e-02       .    
lwt        68.7  -0.011558  9.703e-03   -0.01683  6.953e-03    -0.01865
race        9.6                                                        
    .2            0.115279  3.888e-01    1.20068  5.473e-01       .    
    .3            0.092740  3.160e-01    0.96593  4.610e-01       .    
smoke      33.2                                                        
     .1           0.255357  4.281e-01    0.76820  3.966e-01       .    
ptl        35.1                                                        
   .1             0.617360  8.898e-01    1.75847  8.211e+00       .    
   .2             0.168588  6.150e-01    0.48020  7.592e+00       .    
   .3            -4.910963  5.231e+02  -13.98826  8.828e+02       .    
ht         65.4                                                        
  .1              1.166893  1.034e+00    1.78544  7.287e-01     1.85551
ui         33.8                                                        
  .1              0.310537  5.057e-01    0.91842  4.449e-01       .    
ftv         1.5  -0.001339  2.393e-02   -0.08735  1.728e-01       .    
                                                                       
nVar                                                              2    
BIC                                                          -753.82285
post prob                                                       0.111  
           model 2     model 3     model 4     model 5   
Intercept     1.06795     1.20695     1.08354     0.72186
age             .           .           .           .    
lwt          -0.01692    -0.01881    -0.01805    -0.01631
race                                                     
    .2          .           .           .           .    
    .3          .           .           .           .    
smoke                                                    
     .1         .           .         0.68391     0.65300
ptl                                                      
   .1           .         1.74321       .           .    
   .2           .         0.50107       .           .    
   .3           .       -13.98561       .           .    
ht                                                       
  .1          1.96157     1.92368     1.82202     1.92213
ui                                                       
  .1          0.93000       .           .         0.89627
ftv             .           .           .           .    
                                                         
nVar            3           3           3           4    
BIC        -753.11035  -752.99778  -752.86548  -751.65571
post prob     0.078       0.073       0.069       0.037  
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  176.25 2.66 180.92 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')