George and McCulloch (1997), “Approaches to Bayesian Variable Selection”, Statistica Sinica, 7, 339 – 373.
Ghosh and Clyde (2011) “Rao-Blackwellization for Bayesian variable selection and model averaging in linear and binary regression: A novel data augmentation approach”, Journal of the American Statistical Association, 106 1041-1052. http://homepage.stat.uiowa.edu/~jghsh/ghosh_clyde_2011_jasa.pdf
BoomSpikeSlab packge.
niter <- 1000
model <- lm.spike(yy ~ mls(yy, 1, 1)+ mls(xx,0:5, 3)+ fmls(zz, 34, 66),
niter=niter)
## =-=-=-=-= Iteration 0 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 100 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 200 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 300 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 400 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 500 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 600 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 700 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 800 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
## =-=-=-=-= Iteration 900 Wed Mar 04 16:12:03 2015
## =-=-=-=-=
# plot.ts(model$beta[1:20])
# opar <- par(ask=TRUE)
# on.exit(par(opar))
# hist(model$sigma) ## should be near ?
plot(model, inclusion.threshold = .001)
#plot(model, "fit")
#plot(model, "residuals")
plot(model, "size")
summary(model)
## coefficients:
## mean sd mean.inc sd.inc inc.prob
## mls(yy, 1, 1) 7.02e-01 0.083900 0.70200 0.0839 1.000
## fmls(zz, 34, 66)X.18/m -3.93e-02 0.091400 -0.23700 0.0594 0.166
## (Intercept) -2.82e-02 0.085800 -0.26600 0.0785 0.106
## fmls(zz, 34, 66)X.21/m -1.46e-02 0.046500 -0.15100 0.0409 0.097
## mls(xx, 0:5, 3)X.5/m 3.33e-03 0.044100 0.55400 0.1480 0.006
## mls(xx, 0:5, 3)X.1/m 2.84e-03 0.037700 0.47300 0.1310 0.006
## fmls(zz, 34, 66)X.32/m -3.72e-04 0.006830 -0.12400 0.0158 0.003
## fmls(zz, 34, 66)X.34/m -9.95e-05 0.002320 -0.04980 0.0204 0.002
## fmls(zz, 34, 66)X.31/m -8.70e-05 0.002360 -0.04350 0.0422 0.002
## fmls(zz, 34, 66)X.30/m -5.89e-05 0.001860 -0.05890 0.0000 0.001
## fmls(zz, 34, 66)X.29/m -1.75e-04 0.005540 -0.17500 0.0000 0.001
## fmls(zz, 34, 66)X.26/m -1.38e-04 0.004360 -0.13800 0.0000 0.001
## fmls(zz, 34, 66)X.24/m -1.77e-05 0.000558 -0.01770 0.0000 0.001
## fmls(zz, 34, 66)X.23/m 4.33e-06 0.000137 0.00433 0.0000 0.001
## fmls(zz, 34, 66)X.13/m -1.34e-04 0.004230 -0.13400 0.0000 0.001
## fmls(zz, 34, 66)X.2/m 5.13e-05 0.001620 0.05130 0.0000 0.001
## mls(xx, 0:5, 3)X.2/m -9.60e-05 0.003040 -0.09600 0.0000 0.001
## fmls(zz, 34, 66)X.33/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.28/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.27/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.25/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.22/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.20/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.19/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.17/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.16/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.15/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.14/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.12/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.11/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.10/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.9/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.8/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.7/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.6/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.5/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.4/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.3/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.1/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## fmls(zz, 34, 66)X.0/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## mls(xx, 0:5, 3)X.4/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## mls(xx, 0:5, 3)X.3/m 0.00e+00 0.000000 0.00000 0.0000 0.000
## mls(xx, 0:5, 3)X.0/m 0.00e+00 0.000000 0.00000 0.0000 0.000
##
## residual.sd =
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.5332 0.6288 0.6605 0.6613 0.6920 0.8265
##
## r-square =
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.1617 0.1856 0.2580 0.2527 0.3275 0.5166
# xxx <- cbind( mls(yy, 1, 1),mls(xx,0:5, 3), fmls(zz, 55, 66))
# #head(xxx)
# ncol(xxx)
# prior <- IndependentSpikeSlabPrior(x = xxx,
# y = yy,
# expected.r2 = .5,
# prior.df = .01, # weaker prior than the default
# expected.model.size = 20, # expect 20 nonzero predictors
# prior.beta.sd = NULL,
# optional.coefficient.estimate = NULL, # shrink
# mean.y = mean(y, na.rm = TRUE),
# sdy = sd(as.numeric(y), na.rm = TRUE),
# sdx = apply(as.matrix(x), 2, sd, na.rm = TRUE),
# prior.inclusion.probabilities = NULL,
# number.of.observations = nrow(x),
# number.of.variables = ncol(x),
# scale.by.residual.variance = FALSE,
# sigma.upper.limit = Inf)
#
#
# model <- lm.spike(yy ~ xxx,
# niter=niter)
# plot.ts(model$beta[1:9])
# opar <- par(ask=TRUE)
# on.exit(par(opar))
# hist(model$sigma) ## should be near ?
# plot(model)
# #plot(model, "fit")
# #plot(model, "residuals")
# plot(model, "size")
# summary(model)
library("spikeslab")
## Loading required package: lars
## Loaded lars 1.2
##
## Loading required package: randomForest
## randomForest 4.6-10
## Type rfNews() to see new features/changes/bug fixes.
## Loading required package: parallel
##
## spikeslab 1.1.5
##
## Type spikeslab.news() to see new features, changes, and bug fixes.
##
obj <- spikeslab(yy~mls(yy, 1, 1)+mls(xx,0:5, 3))
print(obj)
## -------------------------------------------------------------------
## Variable selection method : AIC
## Big p small n : FALSE
## Screen variables : FALSE
## Fast processing : TRUE
## Sample size : 122
## No. predictors : 7
## No. burn-in values : 500
## No. sampled values : 500
## Estimated mse : 0.3438
## Model size : 6
##
##
## ---> Top variables:
## bma gnet bma.scale gnet.scale
## mls(xx, 0:5, 3)X.5/m 0.348 0.520 0.321 0.479
## mls(yy, 1, 1) 0.326 0.308 0.426 0.402
## mls(xx, 0:5, 3)X.0/m 0.288 0.445 0.204 0.315
## mls(xx, 0:5, 3)X.1/m 0.140 0.197 0.111 0.156
## mls(xx, 0:5, 3)X.4/m 0.128 0.221 0.100 0.173
## mls(xx, 0:5, 3)X.2/m -0.032 -0.041 -0.029 -0.038
## -------------------------------------------------------------------
#summary(obj)
plot(obj)
#gnet predictor
#yhat.gnet <- predict(obj)$yhat.gnet
#an equivalent call is...
#yhat.gnet <- predict(obj, x = x)$yhat.gnet
# combine yy xx zz
obj <- spikeslab(yy ~ mls(yy, 1, 1)+
mls(xx,0:5, 3)+ mls(zz, 0:21, 66) )
print(obj)
## -------------------------------------------------------------------
## Variable selection method : AIC
## Big p small n : FALSE
## Screen variables : FALSE
## Fast processing : TRUE
## Sample size : 122
## No. predictors : 29
## No. burn-in values : 500
## No. sampled values : 500
## Estimated mse : 0.301
## Model size : 21
##
##
## ---> Top variables:
## bma gnet bma.scale gnet.scale
## mls(yy, 1, 1) 0.341 0.328 0.446 0.428
## mls(xx, 0:5, 3)X.5/m 0.256 0.358 0.236 0.330
## mls(xx, 0:5, 3)X.0/m 0.244 0.326 0.173 0.231
## mls(zz, 0:21, 66)X.21/m -0.172 -0.179 -0.122 -0.128
## mls(xx, 0:5, 3)X.4/m 0.079 0.115 0.061 0.090
## mls(zz, 0:21, 66)X.18/m -0.077 -0.088 -0.083 -0.095
## mls(xx, 0:5, 3)X.1/m 0.067 0.105 0.053 0.083
## mls(zz, 0:21, 66)X.0/m 0.028 0.048 0.036 0.062
## mls(zz, 0:21, 66)X.10/m -0.024 -0.047 -0.031 -0.061
## mls(zz, 0:21, 66)X.1/m 0.018 0.036 0.015 0.031
## mls(zz, 0:21, 66)X.11/m -0.016 -0.023 -0.016 -0.023
## mls(zz, 0:21, 66)X.2/m 0.015 0.028 0.016 0.029
## mls(xx, 0:5, 3)X.3/m 0.012 0.012 0.008 0.008
## mls(zz, 0:21, 66)X.12/m -0.012 -0.020 -0.016 -0.027
## mls(zz, 0:21, 66)X.16/m -0.010 -0.014 -0.010 -0.015
## mls(zz, 0:21, 66)X.6/m 0.009 0.028 0.008 0.025
## mls(zz, 0:21, 66)X.20/m -0.008 -0.012 -0.010 -0.015
## mls(zz, 0:21, 66)X.5/m -0.007 -0.012 -0.009 -0.016
## mls(zz, 0:21, 66)X.4/m 0.006 0.010 0.008 0.013
## mls(zz, 0:21, 66)X.15/m -0.004 -0.005 -0.005 -0.005
## mls(zz, 0:21, 66)X.3/m 0.004 0.006 0.005 0.008
## -------------------------------------------------------------------
#summary(obj)
#par(mfrow=c(1,1))
plot(obj)
interpolation NAs using noclf in zoo package
# combine yy xx zz
# deal with the NAs. randomForest does not allow NAs.
obj <- spikeslab(na.locf(yy, fromLast=TRUE ) ~
mls(na.locf(yy, fromLast=TRUE), 1, 1)+
mls(na.locf(xx),0:8, 3)+ mls(na.locf(zz), 0:131, 66) )
## Error in mls(na.locf(zz), 0:131, 66): Incomplete high frequency data
obj <- spikeslab(x=na.locf(xxx, fromLast=TRUE),
y= na.locf(yy, fromLast=TRUE ),
bigp.smalln = TRUE,
bigp.small.n.factor = 500, max.var = 500)
## Error in na.locf(xxx, fromLast = TRUE): object 'xxx' not found
print(obj)
## -------------------------------------------------------------------
## Variable selection method : AIC
## Big p small n : FALSE
## Screen variables : FALSE
## Fast processing : TRUE
## Sample size : 122
## No. predictors : 29
## No. burn-in values : 500
## No. sampled values : 500
## Estimated mse : 0.301
## Model size : 21
##
##
## ---> Top variables:
## bma gnet bma.scale gnet.scale
## mls(yy, 1, 1) 0.341 0.328 0.446 0.428
## mls(xx, 0:5, 3)X.5/m 0.256 0.358 0.236 0.330
## mls(xx, 0:5, 3)X.0/m 0.244 0.326 0.173 0.231
## mls(zz, 0:21, 66)X.21/m -0.172 -0.179 -0.122 -0.128
## mls(xx, 0:5, 3)X.4/m 0.079 0.115 0.061 0.090
## mls(zz, 0:21, 66)X.18/m -0.077 -0.088 -0.083 -0.095
## mls(xx, 0:5, 3)X.1/m 0.067 0.105 0.053 0.083
## mls(zz, 0:21, 66)X.0/m 0.028 0.048 0.036 0.062
## mls(zz, 0:21, 66)X.10/m -0.024 -0.047 -0.031 -0.061
## mls(zz, 0:21, 66)X.1/m 0.018 0.036 0.015 0.031
## mls(zz, 0:21, 66)X.11/m -0.016 -0.023 -0.016 -0.023
## mls(zz, 0:21, 66)X.2/m 0.015 0.028 0.016 0.029
## mls(xx, 0:5, 3)X.3/m 0.012 0.012 0.008 0.008
## mls(zz, 0:21, 66)X.12/m -0.012 -0.020 -0.016 -0.027
## mls(zz, 0:21, 66)X.16/m -0.010 -0.014 -0.010 -0.015
## mls(zz, 0:21, 66)X.6/m 0.009 0.028 0.008 0.025
## mls(zz, 0:21, 66)X.20/m -0.008 -0.012 -0.010 -0.015
## mls(zz, 0:21, 66)X.5/m -0.007 -0.012 -0.009 -0.016
## mls(zz, 0:21, 66)X.4/m 0.006 0.010 0.008 0.013
## mls(zz, 0:21, 66)X.15/m -0.004 -0.005 -0.005 -0.005
## mls(zz, 0:21, 66)X.3/m 0.004 0.006 0.005 0.008
## -------------------------------------------------------------------
plot(obj)
xxx <- as.data.frame( cbind( mls(na.locf(yy, fromLast=TRUE), 1, 1),
mls(na.locf(xx),0:5, 3),
mls(na.locf(zz, fromLast=TRUE), 0:131, 66)
)
)
names(xxx)
## [1] "X.1/m" "X.0/m" "X.1/m" "X.2/m" "X.3/m" "X.4/m" "X.5/m"
## [8] "X.0/m" "X.1/m" "X.2/m" "X.3/m" "X.4/m" "X.5/m" "X.6/m"
## [15] "X.7/m" "X.8/m" "X.9/m" "X.10/m" "X.11/m" "X.12/m" "X.13/m"
## [22] "X.14/m" "X.15/m" "X.16/m" "X.17/m" "X.18/m" "X.19/m" "X.20/m"
## [29] "X.21/m" "X.22/m" "X.23/m" "X.24/m" "X.25/m" "X.26/m" "X.27/m"
## [36] "X.28/m" "X.29/m" "X.30/m" "X.31/m" "X.32/m" "X.33/m" "X.34/m"
## [43] "X.35/m" "X.36/m" "X.37/m" "X.38/m" "X.39/m" "X.40/m" "X.41/m"
## [50] "X.42/m" "X.43/m" "X.44/m" "X.45/m" "X.46/m" "X.47/m" "X.48/m"
## [57] "X.49/m" "X.50/m" "X.51/m" "X.52/m" "X.53/m" "X.54/m" "X.55/m"
## [64] "X.56/m" "X.57/m" "X.58/m" "X.59/m" "X.60/m" "X.61/m" "X.62/m"
## [71] "X.63/m" "X.64/m" "X.65/m" "X.66/m" "X.67/m" "X.68/m" "X.69/m"
## [78] "X.70/m" "X.71/m" "X.72/m" "X.73/m" "X.74/m" "X.75/m" "X.76/m"
## [85] "X.77/m" "X.78/m" "X.79/m" "X.80/m" "X.81/m" "X.82/m" "X.83/m"
## [92] "X.84/m" "X.85/m" "X.86/m" "X.87/m" "X.88/m" "X.89/m" "X.90/m"
## [99] "X.91/m" "X.92/m" "X.93/m" "X.94/m" "X.95/m" "X.96/m" "X.97/m"
## [106] "X.98/m" "X.99/m" "X.100/m" "X.101/m" "X.102/m" "X.103/m" "X.104/m"
## [113] "X.105/m" "X.106/m" "X.107/m" "X.108/m" "X.109/m" "X.110/m" "X.111/m"
## [120] "X.112/m" "X.113/m" "X.114/m" "X.115/m" "X.116/m" "X.117/m" "X.118/m"
## [127] "X.119/m" "X.120/m" "X.121/m" "X.122/m" "X.123/m" "X.124/m" "X.125/m"
## [134] "X.126/m" "X.127/m" "X.128/m" "X.129/m" "X.130/m" "X.131/m"
sum(is.na(na.locf(xxx, fromLast=TRUE)))
## [1] 0
cv.obj <- cv.spikeslab(x = na.locf(xxx, fromLast=TRUE), na.locf(yy, fromLast=TRUE ) , K = 20)
## K-fold: 1
## K-fold: 2
## K-fold: 3
## K-fold: 4
## K-fold: 5
## K-fold: 6
## K-fold: 7
## K-fold: 8
## K-fold: 9
## K-fold: 10
## K-fold: 11
## K-fold: 12
## K-fold: 13
## K-fold: 14
## K-fold: 15
## K-fold: 16
## K-fold: 17
## K-fold: 18
## K-fold: 19
## K-fold: 20
## final analysis (full-data)
## Warning in min(ms.range, na.rm = TRUE): no non-missing arguments to min;
## returning Inf
## -------------------------------------------------------------------
## Variable selection method : cross-validation
## Big p small n : FALSE
## Screen variables : FALSE
## Fast processing : TRUE
## Sample size : 124
## No. predictors : 139
## No. burn-in values : 500
## No. sampled values : 500
## K-fold : 20
## CV mean-squared error : 0.422 +/- 0.066
## Model size : 136
##
##
## Top variables in terms of stability:
## bma bma.cv gnet gnet.cv stability
## x.1 0.378 0.369 0.370 0.282 95
## x.134 -0.014 -0.049 -0.070 -0.075 75
## x.40 -0.010 -0.008 -0.051 -0.016 60
## x.29 -0.013 -0.019 -0.042 -0.025 55
## x.96 0.007 0.006 0.030 0.010 55
## x.99 -0.001 -0.002 -0.006 -0.008 55
## x.132 -0.015 -0.015 -0.055 -0.021 50
## x.66 -0.011 -0.008 -0.035 -0.014 50
## x.26 -0.009 -0.010 -0.036 -0.010 50
## x.115 -0.006 -0.003 -0.036 -0.008 50
## x.98 0.002 0.003 0.010 0.006 50
## x.105 -0.001 -0.005 -0.007 -0.009 50
## x.88 0.001 0.002 0.005 0.005 50
## x.107 0.001 0.001 0.005 0.002 50
## x.101 0.000 0.005 0.002 0.008 50
## x.94 0.000 -0.002 0.000 -0.003 50
## x.19 -0.007 -0.003 -0.037 -0.007 45
## x.92 -0.006 -0.003 -0.016 -0.005 45
## x.68 -0.005 -0.002 -0.021 -0.005 45
## x.32 0.005 0.003 0.029 0.007 45
## x.39 -0.003 -0.002 -0.014 -0.003 45
## x.123 -0.003 -0.005 -0.014 -0.011 45
## x.133 -0.001 -0.002 -0.013 -0.005 45
## x.77 -0.001 -0.001 -0.010 -0.004 45
## x.139 0.001 0.001 0.007 0.004 45
## x.103 -0.001 -0.002 -0.003 -0.002 45
## x.16 -0.001 -0.001 -0.006 -0.001 45
## x.110 0.001 0.001 0.007 0.002 45
## x.42 -0.001 -0.002 -0.001 -0.001 45
## x.108 -0.001 -0.002 -0.003 -0.003 45
## x.45 -0.001 -0.002 -0.006 -0.008 45
## x.24 0.000 0.000 -0.007 -0.003 45
## x.7 0.000 0.003 0.002 0.009 45
## x.20 0.000 -0.001 -0.002 -0.003 45
## x.18 0.000 -0.001 -0.004 -0.003 45
## x.8 0.000 0.002 0.000 0.006 45
## x.48 0.000 0.001 0.000 0.003 45
## x.60 0.000 0.000 0.000 -0.002 45
## x.71 0.000 -0.001 0.000 -0.003 45
## x.69 -0.003 -0.001 -0.012 -0.002 40
## x.33 -0.002 0.000 -0.022 -0.002 40
## x.36 -0.002 -0.001 -0.002 -0.001 40
## x.135 -0.002 -0.002 -0.008 -0.001 40
## x.61 -0.002 -0.001 -0.007 -0.002 40
## x.51 0.001 0.001 0.006 0.002 40
## x.85 0.001 0.001 0.018 0.004 40
## x.4 -0.001 -0.001 -0.006 -0.002 40
## x.125 0.001 0.001 0.002 0.002 40
## x.62 -0.001 -0.001 -0.003 -0.002 40
## x.47 0.000 0.001 0.009 0.002 40
## x.76 0.000 0.000 -0.002 0.000 40
## x.113 0.000 -0.001 -0.002 -0.003 40
## x.10 0.000 0.001 0.002 0.004 40
## x.46 0.000 0.000 0.003 0.002 40
## x.43 0.000 0.001 0.002 0.002 40
## x.83 0.000 -0.001 -0.001 -0.001 40
## x.81 0.000 0.001 0.000 0.002 40
## x.38 0.000 0.000 0.000 -0.001 40
## x.21 0.000 0.000 0.000 -0.002 40
## x.82 0.000 -0.001 0.000 -0.002 40
## x.91 0.000 -0.001 0.000 -0.002 40
## x.130 0.000 0.000 0.000 0.000 40
## x.97 0.004 0.002 0.025 0.006 35
## x.2 0.002 0.003 0.014 0.011 35
## x.128 -0.002 -0.001 -0.005 0.000 35
## x.109 -0.001 0.000 -0.007 -0.001 35
## x.137 -0.001 0.000 0.001 0.000 35
## x.138 0.001 0.001 0.005 0.004 35
## x.30 0.001 0.001 0.003 0.001 35
## x.117 0.001 0.001 0.000 -0.001 35
## x.9 0.001 0.001 0.007 0.003 35
## x.100 0.001 0.000 0.003 0.000 35
## x.93 0.001 0.000 0.007 0.002 35
## x.129 0.001 0.000 0.002 0.000 35
## x.106 -0.001 0.000 -0.003 -0.001 35
## x.54 -0.001 -0.001 -0.002 0.000 35
## x.41 0.000 0.000 -0.010 -0.001 35
## x.124 0.000 0.000 -0.004 -0.002 35
## x.27 0.000 0.000 0.003 0.000 35
## x.121 0.000 0.000 0.007 0.002 35
## x.17 0.000 0.000 -0.001 -0.001 35
## x.58 0.000 0.000 0.007 0.002 35
## x.59 0.000 0.000 0.000 0.000 35
## x.57 0.000 -0.001 0.000 -0.001 35
## x.126 0.000 0.000 0.001 0.001 35
## x.3 0.000 0.001 0.001 0.002 35
## x.90 0.000 0.000 -0.008 -0.001 35
## x.72 0.000 0.000 0.000 0.000 35
## x.114 0.000 0.000 -0.001 -0.001 35
## x.67 0.000 0.000 -0.002 -0.003 35
## x.34 0.000 0.000 -0.001 -0.001 35
## x.31 0.000 0.000 0.000 0.000 35
## x.6 0.000 0.000 0.001 0.000 35
## x.11 0.000 0.000 0.000 0.000 35
## x.55 0.000 -0.001 0.000 0.000 35
## x.56 0.000 0.000 0.000 0.000 35
## x.122 0.000 0.000 0.000 0.000 35
## x.95 -0.006 -0.001 -0.011 -0.001 30
## x.22 0.002 0.000 0.003 0.000 30
## x.49 -0.001 0.000 -0.009 -0.001 30
## x.75 -0.001 -0.001 -0.005 -0.002 30
## x.73 0.001 0.001 0.002 0.001 30
## x.79 -0.001 0.000 0.000 0.000 30
## x.116 0.000 -0.001 -0.002 -0.001 30
## x.111 0.000 0.000 0.002 0.000 30
## x.86 0.000 0.000 0.000 0.000 30
## x.112 0.000 0.000 -0.002 0.000 30
## x.84 0.000 0.000 0.001 0.000 30
## x.64 0.000 0.000 -0.002 -0.001 30
## x.50 0.000 0.000 0.000 -0.001 30
## x.65 0.000 0.000 -0.001 0.000 30
## x.131 0.000 0.000 0.000 -0.001 30
## x.74 0.000 0.000 0.000 0.000 30
## x.87 0.000 0.000 0.000 0.000 30
## x.119 0.000 0.000 0.000 0.000 30
## x.78 0.000 0.000 0.000 0.000 30
## x.25 0.000 0.000 -0.001 0.000 30
## x.104 0.000 0.000 0.000 0.000 30
## x.15 0.000 0.000 0.000 0.000 30
## x.118 0.000 0.000 0.000 -0.001 30
## x.37 0.000 -0.002 0.000 0.000 30
## -------------------------------------------------------------------
plot(cv.obj, plot.type = "cv")
# combine yy xx zz
obj <- spikeslab(na.locf(yy, fromLast=TRUE ) ~
mls(na.locf(yy, fromLast=TRUE), 1, 1)+
mls(na.locf(xx),0:8, 3)+ mls(na.locf(zz), 0:65, 66) )
## Error in mls(na.locf(zz), 0:65, 66): Incomplete high frequency data
print(obj)
## -------------------------------------------------------------------
## Variable selection method : AIC
## Big p small n : FALSE
## Screen variables : FALSE
## Fast processing : TRUE
## Sample size : 122
## No. predictors : 29
## No. burn-in values : 500
## No. sampled values : 500
## Estimated mse : 0.301
## Model size : 21
##
##
## ---> Top variables:
## bma gnet bma.scale gnet.scale
## mls(yy, 1, 1) 0.341 0.328 0.446 0.428
## mls(xx, 0:5, 3)X.5/m 0.256 0.358 0.236 0.330
## mls(xx, 0:5, 3)X.0/m 0.244 0.326 0.173 0.231
## mls(zz, 0:21, 66)X.21/m -0.172 -0.179 -0.122 -0.128
## mls(xx, 0:5, 3)X.4/m 0.079 0.115 0.061 0.090
## mls(zz, 0:21, 66)X.18/m -0.077 -0.088 -0.083 -0.095
## mls(xx, 0:5, 3)X.1/m 0.067 0.105 0.053 0.083
## mls(zz, 0:21, 66)X.0/m 0.028 0.048 0.036 0.062
## mls(zz, 0:21, 66)X.10/m -0.024 -0.047 -0.031 -0.061
## mls(zz, 0:21, 66)X.1/m 0.018 0.036 0.015 0.031
## mls(zz, 0:21, 66)X.11/m -0.016 -0.023 -0.016 -0.023
## mls(zz, 0:21, 66)X.2/m 0.015 0.028 0.016 0.029
## mls(xx, 0:5, 3)X.3/m 0.012 0.012 0.008 0.008
## mls(zz, 0:21, 66)X.12/m -0.012 -0.020 -0.016 -0.027
## mls(zz, 0:21, 66)X.16/m -0.010 -0.014 -0.010 -0.015
## mls(zz, 0:21, 66)X.6/m 0.009 0.028 0.008 0.025
## mls(zz, 0:21, 66)X.20/m -0.008 -0.012 -0.010 -0.015
## mls(zz, 0:21, 66)X.5/m -0.007 -0.012 -0.009 -0.016
## mls(zz, 0:21, 66)X.4/m 0.006 0.010 0.008 0.013
## mls(zz, 0:21, 66)X.15/m -0.004 -0.005 -0.005 -0.005
## mls(zz, 0:21, 66)X.3/m 0.004 0.006 0.005 0.008
## -------------------------------------------------------------------
plot(obj)
Bayes.poisson.out <- MCMCregress(yy~mls(yy, 1, 1)+mls(xx,0:5,3))
summary((Bayes.poisson.out))
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) -0.20659 0.16688 0.0016688 0.0016688
## mls(yy, 1, 1) 0.37336 0.08143 0.0008143 0.0008375
## mls(xx, 0:5, 3)X.0/m 0.40263 0.13601 0.0013601 0.0014351
## mls(xx, 0:5, 3)X.1/m 0.20032 0.15675 0.0015675 0.0015675
## mls(xx, 0:5, 3)X.2/m -0.06689 0.09840 0.0009840 0.0009840
## mls(xx, 0:5, 3)X.3/m -0.03787 0.11042 0.0011042 0.0010864
## mls(xx, 0:5, 3)X.4/m 0.27495 0.14929 0.0014929 0.0014975
## mls(xx, 0:5, 3)X.5/m 0.59358 0.14505 0.0014505 0.0014505
## sigma2 0.35251 0.04755 0.0004755 0.0005040
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -0.52758 -0.31919 -0.20802 -0.0974584 0.1256
## mls(yy, 1, 1) 0.21511 0.31900 0.37252 0.4274306 0.5345
## mls(xx, 0:5, 3)X.0/m 0.13370 0.31334 0.40183 0.4927874 0.6715
## mls(xx, 0:5, 3)X.1/m -0.10497 0.09541 0.19916 0.3040547 0.5156
## mls(xx, 0:5, 3)X.2/m -0.26306 -0.13237 -0.06584 -0.0007112 0.1262
## mls(xx, 0:5, 3)X.3/m -0.25704 -0.11169 -0.03801 0.0368691 0.1746
## mls(xx, 0:5, 3)X.4/m -0.01652 0.17466 0.27494 0.3741174 0.5720
## mls(xx, 0:5, 3)X.5/m 0.30899 0.49678 0.59448 0.6903994 0.8792
## sigma2 0.27267 0.31851 0.34789 0.3812831 0.4560
#plot(Bayes.poisson.out)
model1 <- MCMCregress(yy~ mls(yy, 1, 1)+mls(xx,0:5, 3))
summary(model1)
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) -0.20659 0.16688 0.0016688 0.0016688
## mls(yy, 1, 1) 0.37336 0.08143 0.0008143 0.0008375
## mls(xx, 0:5, 3)X.0/m 0.40263 0.13601 0.0013601 0.0014351
## mls(xx, 0:5, 3)X.1/m 0.20032 0.15675 0.0015675 0.0015675
## mls(xx, 0:5, 3)X.2/m -0.06689 0.09840 0.0009840 0.0009840
## mls(xx, 0:5, 3)X.3/m -0.03787 0.11042 0.0011042 0.0010864
## mls(xx, 0:5, 3)X.4/m 0.27495 0.14929 0.0014929 0.0014975
## mls(xx, 0:5, 3)X.5/m 0.59358 0.14505 0.0014505 0.0014505
## sigma2 0.35251 0.04755 0.0004755 0.0005040
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -0.52758 -0.31919 -0.20802 -0.0974584 0.1256
## mls(yy, 1, 1) 0.21511 0.31900 0.37252 0.4274306 0.5345
## mls(xx, 0:5, 3)X.0/m 0.13370 0.31334 0.40183 0.4927874 0.6715
## mls(xx, 0:5, 3)X.1/m -0.10497 0.09541 0.19916 0.3040547 0.5156
## mls(xx, 0:5, 3)X.2/m -0.26306 -0.13237 -0.06584 -0.0007112 0.1262
## mls(xx, 0:5, 3)X.3/m -0.25704 -0.11169 -0.03801 0.0368691 0.1746
## mls(xx, 0:5, 3)X.4/m -0.01652 0.17466 0.27494 0.3741174 0.5720
## mls(xx, 0:5, 3)X.5/m 0.30899 0.49678 0.59448 0.6903994 0.8792
## sigma2 0.27267 0.31851 0.34789 0.3812831 0.4560
#plot(model1)
model2 <- MCMCregress(yy~ mls(yy, 1, 1)+mls(xx,0:8, 3))
summary(model2)
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) -0.08837 0.22376 0.0022376 0.0022142
## mls(yy, 1, 1) 0.32661 0.10336 0.0010336 0.0010336
## mls(xx, 0:8, 3)X.0/m 0.40867 0.13945 0.0013945 0.0013945
## mls(xx, 0:8, 3)X.1/m 0.22706 0.17381 0.0017381 0.0016879
## mls(xx, 0:8, 3)X.2/m -0.03451 0.13572 0.0013572 0.0013572
## mls(xx, 0:8, 3)X.3/m -0.02870 0.14713 0.0014713 0.0014713
## mls(xx, 0:8, 3)X.4/m 0.36038 0.17571 0.0017571 0.0017571
## mls(xx, 0:8, 3)X.5/m 0.57955 0.14859 0.0014859 0.0014859
## mls(xx, 0:8, 3)X.6/m -0.11465 0.12651 0.0012651 0.0012651
## mls(xx, 0:8, 3)X.7/m 0.13981 0.17562 0.0017562 0.0017556
## mls(xx, 0:8, 3)X.8/m 0.06541 0.15713 0.0015713 0.0015713
## sigma2 0.35787 0.04942 0.0004942 0.0005487
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -0.52183 -0.23857 -0.08991 0.06201 0.3540
## mls(yy, 1, 1) 0.11993 0.25793 0.32650 0.39384 0.5330
## mls(xx, 0:8, 3)X.0/m 0.13104 0.31768 0.40853 0.49929 0.6866
## mls(xx, 0:8, 3)X.1/m -0.11624 0.11219 0.22466 0.34242 0.5698
## mls(xx, 0:8, 3)X.2/m -0.30009 -0.12424 -0.03565 0.05831 0.2309
## mls(xx, 0:8, 3)X.3/m -0.32100 -0.12769 -0.02814 0.07262 0.2542
## mls(xx, 0:8, 3)X.4/m 0.01149 0.24108 0.36248 0.47809 0.6977
## mls(xx, 0:8, 3)X.5/m 0.28532 0.48121 0.57944 0.67807 0.8760
## mls(xx, 0:8, 3)X.6/m -0.36380 -0.20035 -0.11416 -0.02883 0.1340
## mls(xx, 0:8, 3)X.7/m -0.20317 0.02196 0.14013 0.26009 0.4826
## mls(xx, 0:8, 3)X.8/m -0.24280 -0.04103 0.06742 0.17034 0.3709
## sigma2 0.27372 0.32315 0.35327 0.38813 0.4651
model3 <- MCMCregress(yy~ mls(yy, 1, 1)+mls(xx,0:11, 3))
summary(model3)
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) -0.14576 0.25464 0.0025464 0.0025464
## mls(yy, 1, 1) 0.30373 0.09831 0.0009831 0.0010507
## mls(xx, 0:11, 3)X.0/m 0.45689 0.14431 0.0014431 0.0014431
## mls(xx, 0:11, 3)X.1/m 0.38830 0.16788 0.0016788 0.0016788
## mls(xx, 0:11, 3)X.2/m 0.01854 0.15394 0.0015394 0.0015394
## mls(xx, 0:11, 3)X.3/m 0.05380 0.13779 0.0013779 0.0013779
## mls(xx, 0:11, 3)X.4/m 0.44459 0.16973 0.0016973 0.0016997
## mls(xx, 0:11, 3)X.5/m 0.50938 0.14534 0.0014534 0.0014534
## mls(xx, 0:11, 3)X.6/m -0.07929 0.13078 0.0013078 0.0013078
## mls(xx, 0:11, 3)X.7/m -0.04997 0.17414 0.0017414 0.0017414
## mls(xx, 0:11, 3)X.8/m 0.07245 0.15390 0.0015390 0.0015390
## mls(xx, 0:11, 3)X.9/m 0.34158 0.13104 0.0013104 0.0013104
## mls(xx, 0:11, 3)X.10/m -0.10533 0.16330 0.0016330 0.0015973
## mls(xx, 0:11, 3)X.11/m -0.36078 0.14117 0.0014117 0.0014117
## sigma2 0.30560 0.04312 0.0004312 0.0004808
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -0.64204 -0.31783 -0.14745 0.023926 0.36203
## mls(yy, 1, 1) 0.11032 0.23836 0.30327 0.371170 0.49385
## mls(xx, 0:11, 3)X.0/m 0.17466 0.36255 0.45866 0.552361 0.73796
## mls(xx, 0:11, 3)X.1/m 0.05272 0.27646 0.38817 0.501666 0.71823
## mls(xx, 0:11, 3)X.2/m -0.28268 -0.08359 0.01760 0.123762 0.31787
## mls(xx, 0:11, 3)X.3/m -0.21308 -0.03938 0.05313 0.145660 0.32474
## mls(xx, 0:11, 3)X.4/m 0.11857 0.32886 0.44269 0.556880 0.77371
## mls(xx, 0:11, 3)X.5/m 0.22404 0.41353 0.51135 0.607280 0.79142
## mls(xx, 0:11, 3)X.6/m -0.33827 -0.16604 -0.07862 0.007870 0.17746
## mls(xx, 0:11, 3)X.7/m -0.39895 -0.16515 -0.05107 0.067865 0.28982
## mls(xx, 0:11, 3)X.8/m -0.22989 -0.03212 0.07059 0.176392 0.37474
## mls(xx, 0:11, 3)X.9/m 0.08574 0.25329 0.34180 0.429835 0.60223
## mls(xx, 0:11, 3)X.10/m -0.42488 -0.21584 -0.10540 0.005763 0.21174
## mls(xx, 0:11, 3)X.11/m -0.63535 -0.45560 -0.36123 -0.265503 -0.08521
## sigma2 0.23225 0.27476 0.30199 0.332431 0.40091
Error message:
BF <- BayesFactor(model1, model2, model3)
mod.probs <- PostProbMod(BF)
Error in BF.logmarglike[i] <- attr(model.list[[i]], "logmarglike") :
replacement has length zero
Bayes.poisson.out <- MCMCregress(yy~mls(yy, 1, 1)+
mls(xx,0:5, 3)+ fmls(zz, 21, 66))
summary((Bayes.poisson.out))
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## (Intercept) -0.237395 0.17721 0.0017721 0.0017721
## mls(yy, 1, 1) 0.345379 0.08900 0.0008900 0.0009100
## mls(xx, 0:5, 3)X.0/m 0.405199 0.14448 0.0014448 0.0014448
## mls(xx, 0:5, 3)X.1/m 0.236058 0.16906 0.0016906 0.0016906
## mls(xx, 0:5, 3)X.2/m -0.045597 0.10685 0.0010685 0.0010859
## mls(xx, 0:5, 3)X.3/m -0.018291 0.11186 0.0011186 0.0011679
## mls(xx, 0:5, 3)X.4/m 0.224002 0.15970 0.0015970 0.0015970
## mls(xx, 0:5, 3)X.5/m 0.665406 0.15903 0.0015903 0.0015903
## fmls(zz, 21, 66)X.0/m 0.089753 0.07998 0.0007998 0.0007998
## fmls(zz, 21, 66)X.1/m 0.052828 0.06279 0.0006279 0.0006279
## fmls(zz, 21, 66)X.2/m 0.061710 0.06881 0.0006881 0.0006881
## fmls(zz, 21, 66)X.3/m 0.064267 0.08695 0.0008695 0.0008695
## fmls(zz, 21, 66)X.4/m 0.084036 0.08382 0.0008382 0.0008382
## fmls(zz, 21, 66)X.5/m -0.086786 0.07929 0.0007929 0.0007929
## fmls(zz, 21, 66)X.6/m 0.139246 0.06462 0.0006462 0.0006479
## fmls(zz, 21, 66)X.7/m -0.006835 0.06726 0.0006726 0.0006726
## fmls(zz, 21, 66)X.8/m 0.046615 0.07977 0.0007977 0.0007844
## fmls(zz, 21, 66)X.9/m 0.023698 0.07017 0.0007017 0.0007017
## fmls(zz, 21, 66)X.10/m -0.157532 0.08684 0.0008684 0.0008684
## fmls(zz, 21, 66)X.11/m -0.009440 0.06218 0.0006218 0.0006537
## fmls(zz, 21, 66)X.12/m -0.023983 0.07957 0.0007957 0.0008054
## fmls(zz, 21, 66)X.13/m 0.024205 0.07475 0.0007475 0.0007475
## fmls(zz, 21, 66)X.14/m -0.006538 0.07080 0.0007080 0.0007080
## fmls(zz, 21, 66)X.15/m -0.079532 0.06854 0.0006854 0.0006854
## fmls(zz, 21, 66)X.16/m -0.017451 0.06831 0.0006831 0.0006995
## fmls(zz, 21, 66)X.17/m 0.014686 0.07756 0.0007756 0.0007756
## fmls(zz, 21, 66)X.18/m -0.084527 0.07420 0.0007420 0.0007420
## fmls(zz, 21, 66)X.19/m 0.004848 0.07878 0.0007878 0.0007878
## fmls(zz, 21, 66)X.20/m -0.044359 0.08587 0.0008587 0.0008587
## fmls(zz, 21, 66)X.21/m -0.157944 0.05367 0.0005367 0.0005271
## sigma2 0.310751 0.04727 0.0004727 0.0006024
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## (Intercept) -0.58146 -0.357808 -0.237987 -0.11829 0.11361
## mls(yy, 1, 1) 0.16698 0.285385 0.346610 0.40592 0.51704
## mls(xx, 0:5, 3)X.0/m 0.12395 0.307136 0.404801 0.50105 0.69152
## mls(xx, 0:5, 3)X.1/m -0.09870 0.122976 0.235516 0.35054 0.56355
## mls(xx, 0:5, 3)X.2/m -0.25409 -0.116220 -0.045917 0.02663 0.16561
## mls(xx, 0:5, 3)X.3/m -0.23750 -0.094939 -0.017832 0.05920 0.19610
## mls(xx, 0:5, 3)X.4/m -0.08822 0.117578 0.223305 0.33040 0.54042
## mls(xx, 0:5, 3)X.5/m 0.35327 0.560729 0.664809 0.77172 0.97356
## fmls(zz, 21, 66)X.0/m -0.06886 0.036352 0.091068 0.14316 0.24692
## fmls(zz, 21, 66)X.1/m -0.06936 0.010961 0.053264 0.09378 0.17543
## fmls(zz, 21, 66)X.2/m -0.07344 0.016525 0.062092 0.10818 0.19683
## fmls(zz, 21, 66)X.3/m -0.10493 0.006640 0.063899 0.12132 0.23625
## fmls(zz, 21, 66)X.4/m -0.08114 0.027874 0.083516 0.13991 0.24999
## fmls(zz, 21, 66)X.5/m -0.24224 -0.140241 -0.086860 -0.03323 0.06945
## fmls(zz, 21, 66)X.6/m 0.01090 0.096926 0.139834 0.18219 0.26666
## fmls(zz, 21, 66)X.7/m -0.13838 -0.051709 -0.007277 0.03764 0.12604
## fmls(zz, 21, 66)X.8/m -0.11115 -0.006065 0.046108 0.10039 0.20301
## fmls(zz, 21, 66)X.9/m -0.11176 -0.023422 0.024097 0.07061 0.15987
## fmls(zz, 21, 66)X.10/m -0.32767 -0.215501 -0.157873 -0.09985 0.01401
## fmls(zz, 21, 66)X.11/m -0.13124 -0.051206 -0.010213 0.03220 0.11233
## fmls(zz, 21, 66)X.12/m -0.18123 -0.076524 -0.023984 0.02966 0.13414
## fmls(zz, 21, 66)X.13/m -0.12258 -0.025553 0.023860 0.07401 0.17114
## fmls(zz, 21, 66)X.14/m -0.14706 -0.053223 -0.006911 0.04081 0.13062
## fmls(zz, 21, 66)X.15/m -0.21267 -0.124249 -0.080754 -0.03328 0.05380
## fmls(zz, 21, 66)X.16/m -0.15080 -0.063138 -0.016692 0.02782 0.11662
## fmls(zz, 21, 66)X.17/m -0.13607 -0.037147 0.014160 0.06657 0.16896
## fmls(zz, 21, 66)X.18/m -0.22879 -0.134557 -0.084968 -0.03395 0.06233
## fmls(zz, 21, 66)X.19/m -0.14782 -0.047395 0.003921 0.05674 0.16160
## fmls(zz, 21, 66)X.20/m -0.21204 -0.102015 -0.044172 0.01321 0.12384
## fmls(zz, 21, 66)X.21/m -0.26379 -0.193874 -0.157642 -0.12174 -0.05318
## sigma2 0.23294 0.277030 0.306348 0.33891 0.41552
#plot(Bayes.poisson.out)
Next step: user defined likelihood function, prior.
## log of inverse gamma density
logdinvgamma <- function(sigma2, a, b){
logf <- a * log(b) - lgamma(a) -(a+1) * log(sigma2) -b/sigma2
return(logf)
}
## log of multivariate normal density
logdmvnorm <- function(theta, mu, Sigma){
d <- length(theta)
logf <- -0.5*d * log(2*pi) - 0.5*log(det(Sigma)) -
0.5 * t(theta - mu) %*% solve(Sigma) %*% (theta -mu)
return(logf)
}
k <- length(theta)
beta <- theta[1:(k-1)]
alpha <- exp(theta[k])
mu <- exp(X %*% beta)
## evaluate log-likelihood at (alpha, beta)
log.like <- sum(
lgamma(y+alpha) - lfactorial(y) - lgamma(alpha) +
alpha * log(alpha/(alpha+mu)) +
y * log(mu/(alpha+mu))
)
## evaluate log prior at (alpha, beta)
## note Jacobian term necessary b/c of transformation
log.prior <- logdinvgamma(alpha, c0, d0) + theta[k] +
logdmvnorm(beta, b0, B0)