rm( list = ls( ) ) # in this linear regression example, there's one row in # the data set for each day in 1976 (n = 366) # the outcome variable is # log( daily maximum of 24 hourly average ozone values ) # measured at a location in the Los Angeles basin, # standardized (mean 0, sd 1 ) # the original predictor variables were # day_of_year (1 to 366), to pick up time trends # temp_sandburg at sandburg air force base # pressure_500 atmospheric air pressure # humidity relative humidity # inversion_height height of inversion layer # wind wind speed # pressure_gradient strength of air pressure gradient # inversion_temp inversion layer temperature difference # visibility km of visibility at sandberg air force base # we standardized all of the original x variables # loess (nonparametric local regression) indicated cubic # relationships between the outcome and both temp_sandburg # and inversion_temp, and lots of quadratic relationships # between the outcome and the other predictors # with 9 main effects there are 9 * 8 / 2 two-way interaction # terms formed by multiplying, e.g., humidity * wind # our final model had 9 main effects, 9 quadratic terms, # 2 cubic terms, and 36 two-way interactions, for a total # of k = 56 predictors ozone <- matrix( scan( "ozone-data.txt" ), ncol = 57, byrow = T ) dim( ozone ) # 366 57 no.missing <- ( apply( is.na( ozone ), 1, sum ) == 0 ) ozone <- ozone[ no.missing, ] dim( ozone ) # 330 57 # after omitting rows with any missing data, we have n = 330 # and k = 56 y <- ozone[ , 1 ] x <- ozone[ , -1 ] summary( lm( y ~ x ) ) # Call: # lm(formula = y ~ x) # Residuals: # Min 1Q Median 3Q Max # -1.0216 -0.2087 0.0066 0.2322 0.9819 # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 0.91897 0.09399 9.777 < 2e-16 *** # x1 -0.22681 0.04458 -5.088 6.75e-07 *** # x2 -0.06742 0.03021 -2.232 0.026422 * # x3 0.10640 0.07341 1.449 0.148386 # x4 0.01912 0.04416 0.433 0.665389 # x5 0.17349 0.10701 1.621 0.106113 # x6 -0.07188 0.05932 -1.212 0.226650 # x7 -0.21929 0.05833 -3.759 0.000208 *** # x8 0.27355 0.14262 1.918 0.056150 . # x9 -0.17680 0.04053 -4.362 1.83e-05 *** # x10 -0.46300 0.05409 -8.559 8.36e-16 *** # x11 -0.03924 0.02578 -1.522 0.129145 # x12 -0.13051 0.10503 -1.243 0.215115 # x13 -0.15101 0.05168 -2.922 0.003769 ** # x14 -0.04322 0.12975 -0.333 0.739326 # x15 -0.09093 0.10105 -0.900 0.368953 # x16 -0.19070 0.05515 -3.458 0.000631 *** # x17 0.26892 0.36730 0.732 0.464716 # x18 0.05603 0.02937 1.908 0.057442 . # x19 0.01560 0.02396 0.651 0.515647 # x20 -0.06217 0.02635 -2.359 0.019014 * # x21 -0.01098 0.03547 -0.310 0.757120 # x22 -0.02493 0.07386 -0.338 0.735967 # x23 -0.05955 0.04933 -1.207 0.228425 # x24 0.23065 0.08668 2.661 0.008258 ** # x25 -0.12236 0.06329 -1.933 0.054236 . # x26 -0.09895 0.06337 -1.561 0.119614 # x27 -0.32071 0.12318 -2.604 0.009730 ** # x28 -0.05382 0.03322 -1.620 0.106389 # x29 -0.07011 0.07221 -0.971 0.332397 # x30 -0.08908 0.04562 -1.952 0.051902 . # x31 -0.04572 0.09086 -0.503 0.615264 # x32 0.03426 0.06770 0.506 0.613266 # x33 0.08159 0.04623 1.765 0.078719 . # x34 0.17544 0.12640 1.388 0.166254 # x35 0.02575 0.03579 0.719 0.472566 # x36 0.02479 0.08707 0.285 0.776075 # x37 0.20944 0.18908 1.108 0.268972 # x38 -0.05278 0.16494 -0.320 0.749244 # x39 -0.14402 0.12287 -1.172 0.242160 # x40 -0.02437 0.31608 -0.077 0.938607 # x41 -0.03208 0.06822 -0.470 0.638582 # x42 0.19448 0.12151 1.601 0.110627 # x43 -0.07026 0.09153 -0.768 0.443380 # x44 0.01430 0.08323 0.172 0.863685 # x45 -0.10521 0.17319 -0.607 0.544032 # x46 -0.02132 0.04411 -0.483 0.629265 # x47 -0.08792 0.17750 -0.495 0.620762 # x48 -0.00652 0.14125 -0.046 0.963218 # x49 -0.29814 0.37437 -0.796 0.426512 # x50 0.18773 0.09608 1.954 0.051745 . # x51 -0.02976 0.08638 -0.345 0.730690 # x52 0.20866 0.34358 0.607 0.544164 # x53 -0.03527 0.06724 -0.525 0.600311 # x54 0.12947 0.20030 0.646 0.518584 # x55 -0.03962 0.05515 -0.718 0.473104 # x56 -0.15294 0.13456 -1.137 0.256702 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # Residual standard error: 0.3995 on 273 degrees of freedom # Multiple R-squared: 0.866, Adjusted R-squared: 0.8385 # F-statistic: 31.5 on 56 and 273 DF, p-value: < 2.2e-16 library( BMA ) bma.1 <- bicreg( x, y ) summary( bma.1 ) # Call: # bicreg(x = x, y = y) # 71 models were selected # Best 5 models (cumulative posterior probability = 0.1994 ): # p!=0 EV SD model 1 model 2 model 3 # Intercept 100.0 7.201e-01 0.073979 0.75756 0.67720 0.65138 # X1 100.0 -2.137e-01 0.041125 -0.17844 -0.18245 -0.17364 # X2 97.6 -8.711e-02 0.031015 -0.08389 -0.09561 -0.08738 # X3 24.5 3.136e-02 0.061417 . . . # X5 19.7 3.027e-02 0.068327 . . . # X7 61.4 -7.559e-02 0.071298 -0.10513 . . # X8 100.0 5.148e-01 0.084637 0.57599 0.50468 0.60810 # X9 100.0 -2.081e-01 0.034004 -0.21563 -0.20835 -0.20797 # X10 100.0 -4.365e-01 0.053268 -0.47598 -0.41316 -0.40365 # X11 0.4 -9.768e-05 0.001928 . . . # X12 44.7 -2.718e-02 0.032129 -0.06077 . -0.05715 # X13 100.0 -1.271e-01 0.032759 -0.12968 -0.11496 -0.10245 # X14 57.0 -4.437e-02 0.042084 . -0.07398 . # X16 100.0 -1.786e-01 0.025776 -0.18352 -0.17111 -0.18484 # X17 0.4 -1.222e-04 0.002622 . . . # X18 100.0 7.391e-02 0.021187 0.07476 0.07545 0.07835 # X20 64.4 -3.355e-02 0.029665 -0.05866 . -0.05569 # X23 2.1 -9.201e-04 0.007252 . . . # X24 0.0 0.000e+00 0.000000 . . . # X25 0.4 -9.527e-05 0.002095 . . . # X26 61.9 -4.563e-02 0.042477 . . . # X27 0.0 0.000e+00 0.000000 . . . # X28 0.5 -1.524e-04 0.002700 . . . # X29 0.0 0.000e+00 0.000000 . . . # X30 57.8 -3.939e-02 0.041816 . . . # X33 18.0 1.267e-02 0.030583 . . . # X34 2.8 1.501e-03 0.010048 . . . # X38 0.3 -3.277e-04 0.006327 . . . # X42 100.0 1.421e-01 0.036862 0.12303 0.16318 0.12866 # X51 100.0 -1.298e-01 0.027232 -0.13854 -0.13273 -0.14224 # X52 2.7 1.878e-03 0.013508 . . . # nVar 13 11 12 # r2 0.832 0.826 0.829 # BIC -513.99051 -513.86943 -513.58803 # post prob 0.047 0.044 0.039 # model 4 model 5 # Intercept 0.77381 0.72259 # X1 -0.24957 -0.23843 # X2 -0.08792 -0.07988 # X3 0.14697 . # X5 . 0.15642 # X7 -0.11579 -0.15490 # X8 0.47680 0.47375 # X9 -0.20527 -0.21509 # X10 -0.45601 -0.44289 # X11 . . # X12 . -0.06377 # X13 -0.14885 -0.14011 # X14 -0.07858 . # X16 -0.18174 -0.18275 # X17 . . # X18 0.07068 0.07821 # X20 -0.04126 -0.06119 # X23 . . # X24 . . # X25 . . # X26 -0.08467 -0.07389 # X27 . . # X28 . . # X29 . . # X30 -0.05796 . # X33 . . # X34 . . # X38 . . # X42 0.14876 0.12205 # X51 -0.12411 -0.11753 # X52 . . # nVar 16 15 # r2 0.841 0.838 # BIC -513.41332 -513.31995 # post prob 0.035 0.034 # the next two commands produce posterior means and SDs for the regression # coefficients averaged over the models with high posterior probability bma.1$postmean # (Intercept) X1 X2 X3 X5 # 7.200944e-01 -2.136823e-01 -8.711448e-02 3.136008e-02 3.027459e-02 # X7 X8 X9 X10 X11 # -7.558820e-02 5.148150e-01 -2.080596e-01 -4.365438e-01 -9.768129e-05 # X12 X13 X14 X16 X17 # -2.717735e-02 -1.271128e-01 -4.437123e-02 -1.785687e-01 -1.222472e-04 # X18 X20 X23 X24 X25 # 7.390978e-02 -3.355270e-02 -9.201258e-04 0.000000e+00 -9.527104e-05 # X26 X27 X28 X29 X30 # -4.563153e-02 0.000000e+00 -1.523711e-04 0.000000e+00 -3.939289e-02 # X33 X34 X38 X42 X51 # 1.266717e-02 1.500571e-03 -3.277389e-04 1.420571e-01 -1.297645e-01 # X52 # 1.877559e-03 bma.1$postsd # [1] 0.073979278 0.041125476 0.031015122 0.061417007 0.068326999 0.071297621 # [7] 0.084636693 0.034004268 0.053267819 0.001927955 0.032128644 0.032759296 # [13] 0.042084181 0.025776226 0.002622097 0.021187159 0.029664507 0.007251732 # [19] 0.000000000 0.002094577 0.042477121 0.000000000 0.002699860 0.000000000 # [25] 0.041815973 0.030582957 0.010047634 0.006327234 0.036862084 0.027231981 # [31] 0.013508010 plot( bma.1, mfrow = c( 3, 3 ) ) library( BMS ) ozone.bms <- bms( ozone, mprior = "uniform", g = "UIP", user.int = F, iter = 1000000 ) coef( ozone.bms ) # PIP Post Mean Post SD Cond.Pos.Sign Idx # beta 1 1.000000 -0.1917502107 0.036923352 0.00000000 1 # beta 10 1.000000 -0.4636109780 0.059200361 0.00000000 10 # beta 16 1.000000 -0.1832016709 0.031045534 0.00000000 16 # beta 9 0.999638 -0.1890370279 0.040348140 0.00000000 9 # beta 8 0.964896 0.4457047341 0.137676418 1.00000000 8 # beta 13 0.960852 -0.1160656618 0.043425677 0.00000000 13 # beta 2 0.948431 -0.0912197092 0.034655971 0.00000000 2 # beta 15 0.909843 -0.1457807915 0.065757171 0.00000000 15 # beta 18 0.904739 0.0670822286 0.030777727 1.00000000 18 # beta 7 0.834935 -0.1290581989 0.075418183 0.00000000 7 # beta 43 0.761267 -0.1066395093 0.071895885 0.00013005 43 # beta 20 0.623837 -0.0319245697 0.029428835 0.00000000 20 # beta 42 0.502652 0.0586289608 0.072364375 0.99921815 42 # beta 12 0.486047 -0.0289511407 0.036732317 0.00204713 12 # beta 5 0.388395 0.0660952034 0.098141219 0.99993563 5 # beta 51 0.372360 -0.0403003472 0.061435172 0.00092652 51 # beta 37 0.280159 -0.0200372175 0.044623221 0.05299491 37 # beta 39 0.269753 -0.0272328187 0.055434754 0.00514174 39 # beta 26 0.268726 -0.0177986183 0.035659523 0.00116103 26 # beta 45 0.221403 0.0069587549 0.067387672 0.69219478 45 # beta 38 0.190548 -0.0143177273 0.037292220 0.00945693 38 # beta 30 0.179247 -0.0088017717 0.023970764 0.00000000 30 # beta 28 0.159528 -0.0062295089 0.017772130 0.00000000 28 # beta 3 0.157536 0.0156311322 0.044965688 0.99861619 3 # beta 40 0.153334 0.0043349827 0.045563031 0.60527346 40 # beta 36 0.144694 0.0106702154 0.037472850 0.87779037 36 # beta 44 0.140328 0.0093350074 0.033632561 0.91618209 44 # beta 14 0.131673 -0.0060199935 0.023209209 0.10231407 14 # beta 22 0.119782 -0.0061905806 0.023497997 0.04118315 22 # beta 54 0.117474 0.0040991838 0.039845863 0.65697090 54 # beta 27 0.115953 -0.0075217398 0.035422445 0.02346640 27 # beta 6 0.109729 -0.0090399240 0.036484740 0.06927066 6 # beta 19 0.109211 0.0013680988 0.009319821 0.75231433 19 # beta 48 0.100468 -0.0009614776 0.027962408 0.43542222 48 # beta 23 0.096300 -0.0014460572 0.012743566 0.22501558 23 # beta 33 0.095918 0.0033092228 0.016424696 0.91524010 33 # beta 52 0.093240 0.0024579032 0.034958508 0.50128700 52 # beta 49 0.093058 -0.0013568244 0.023659500 0.39254014 49 # beta 35 0.092969 0.0022322686 0.010177038 0.99362153 35 # beta 50 0.082451 0.0027436910 0.015775076 0.92281476 50 # beta 47 0.081719 -0.0018859296 0.018440761 0.26871352 47 # beta 17 0.079831 -0.0006258293 0.021901876 0.49555937 17 # beta 24 0.079100 0.0029400344 0.025426775 0.60108723 24 # beta 11 0.076583 -0.0012133683 0.006474708 0.00257237 11 # beta 4 0.071125 -0.0018047884 0.013541730 0.12808436 4 # beta 25 0.070811 -0.0006558970 0.014147125 0.61481973 25 # beta 41 0.070330 -0.0009189372 0.011234257 0.33426703 41 # beta 34 0.064014 0.0009513981 0.010052167 0.83847283 34 # beta 56 0.063533 0.0003261082 0.012437870 0.67012419 56 # beta 46 0.063316 -0.0006849353 0.009087031 0.24834165 46 # beta 31 0.063297 0.0005548212 0.008808805 0.73932414 31 # beta 29 0.062114 0.0003185423 0.007992786 0.63607238 29 # beta 53 0.061746 0.0001330240 0.008527331 0.45390795 53 # beta 32 0.058892 -0.0006051994 0.008397498 0.30474428 32 # beta 55 0.058851 0.0001427340 0.008174164 0.55876706 55 # beta 21 0.054707 -0.0003266687 0.006181832 0.24879814 21 summary( ozone.bms ) # Mean no. regressors Draws Burnins # "17.2613" "1e+06" "1000" # Time No. models visited Modelspace 2^K # "1.980647 mins" "195470" "7.2e+16" # % visited % Topmodels Corr PMP # "2.7e-10" "4.2" "0.8846" # No. Obs. Model Prior g-Prior # "330" "uniform / 28" "UIP" # Shrinkage-Stats # "Av=0.997" topmodels.bma( ozone.bms )[, 1:4 ] # 0c3db5000002000 0cbdb5000002000 0c3db5000006000 0cbdb5000006000 # beta 1 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 2 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 3 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 4 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 5 0.0000000000 1.0000000000 0.0000000000 1.0000000000 # beta 6 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 7 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 8 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 9 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 10 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 11 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 12 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 13 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 14 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 15 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 16 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 17 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 18 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 19 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 20 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 21 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 22 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 23 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 24 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 25 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 26 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 27 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 28 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 29 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 30 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 31 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 32 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 33 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 34 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 35 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 36 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 37 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 38 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 39 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 40 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 41 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 42 0.0000000000 0.0000000000 1.0000000000 1.0000000000 # beta 43 1.0000000000 1.0000000000 1.0000000000 1.0000000000 # beta 44 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 45 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 46 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 47 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 48 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 49 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 50 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 51 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 52 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 53 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 54 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 55 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # beta 56 0.0000000000 0.0000000000 0.0000000000 0.0000000000 # PMP (Exact) 0.0008655389 0.0008396589 0.0007421381 0.0006594533 # PMP (MCMC) 0.0008040000 0.0007420000 0.0007340000 0.0005420000 image( ozone.bms ) plotModelsize( ozone.bms )