library("BMA")

# e.g. 1 ------------------------------------------------------------------

rates <- round(seq(.1, .9, length.out=15), 2)
rates

set.seed(1234)
x1 <- rbinom(1000, 1, rates[1])
x2 <- rbinom(1000, 1, rates[2])
x3 <- rbinom(1000, 1, rates[3])
x4 <- rbinom(1000, 1, rates[4])
x5 <- rbinom(1000, 1, rates[5])
x6 <- rbinom(1000, 1, rates[6])
x7 <- rbinom(1000, 1, rates[7])
x8 <- rbinom(1000, 1, rates[8])
x9 <- rbinom(1000, 1, rates[9])
x10 <- rbinom(1000, 1, rates[10])
x11 <- rbinom(1000, 1, rates[11])
x12 <- rbinom(1000, 1, rates[12])
x13 <- rbinom(1000, 1, rates[13])
x14 <- rbinom(1000, 1, rates[14])
x15 <- rbinom(1000, 1, rates[15])

x16 <- rnorm(1000)
x17 <- rnorm(1000)
x18 <- rnorm(1000)
x19 <- rnorm(1000)
x20 <- rnorm(1000)
x21 <- rnorm(1000)
x22 <- rnorm(1000)
x23 <- rnorm(1000)
x24 <- rnorm(1000)
x25 <- rnorm(1000)
x26 <- rnorm(1000)
x27 <- rnorm(1000)
x28 <- rnorm(1000)
x29 <- rnorm(1000)
x30 <- rnorm(1000)

y <- rbinom(1000, 1, 0.5)

example1.dat <- data.frame(y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10,
                           x11, x12, x13, x14, x15, x16, x17, x18, x19, x20,
                           x21, x22, x23, x24, x25, x26, x27, x28, x29, x30)
names(example1.dat)
head(example1.dat)

fml <- as.formula(paste0("y ~",
                         paste(names(example1.dat[, -y]), collapse = ' + ')))
output1 <- bic.glm(fml, glm.family = "binomial",
                   data = example1.dat, maxCol = 31)
summary(output1)
output1$postprob
output1$label
output1$probne0
output1$postmean
output1$postsd

output1$mle
output1$se

imageplot.bma(output1)

output1.aic <- glm(fml, data = example1.dat, family = binomial(link = 'logit'))
summary(output1.aic)
step.aic1 <- step(output1.aic)
summary(step.aic1)


# e.g. 2 -------------------------------------------------------------------

rates <- round(seq(.1, .9, length.out=15), 2)
rates

x1 <- rbinom(500, 1, rates[1])
x2 <- rbinom(500, 1, rates[2])
x3 <- rbinom(500, 1, rates[3])
x4 <- rbinom(500, 1, rates[4])
x5 <- rbinom(500, 1, rates[5])
x6 <- rbinom(500, 1, rates[6])
x7 <- rbinom(500, 1, rates[7])
x8 <- rbinom(500, 1, rates[8])
x9 <- rbinom(500, 1, rates[9])
x10 <- rbinom(500, 1, rates[10])

x11 <- rnorm(500)
x12 <- rnorm(500)
x13 <- rnorm(500)
x14 <- rnorm(500)
x15 <- rnorm(500)
x16 <- rnorm(500)
x17 <- rnorm(500)
x18 <- rnorm(500)
x19 <- rnorm(500)
x20 <- rnorm(500)

inv.logit.rate <- exp(x1 + x2 + x3 + x4 + x5 + x11 + x12 + x13 + x14 +x15) /
    (1 + exp(x1 + x2 + x3 + x4 + x5 + x11 + x12 + x13 + x14 +x15))
y <- rbinom(500, 1, inv.logit.rate)

example2.dat <- data.frame(y, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10,
                           x11, x12, x13, x14, x15, x16, x17, x18, x19, x20)

fml <- as.formula(paste0("y ~",
                         paste(names(example2.dat[, -y]), collapse = ' + ')))
output2 <- bic.glm(fml, glm.family = "binomial",
                   data = example2.dat)
summary(output2)
output2$postprob
output2$label
output2$probne0
output2$postmean

output2$mle
output2$se


output2.aic <- glm(fml, data = example2.dat, family = "binomial")
summary(output2.aic)
step.aic <- step(output2.aic)
summary(output2.aic)

# e.g. 3 --------------------------------------------

data("birthwt", package = 'MASS')

birthwt$smoke <- as.factor(birthwt$smoke)
birthwt$race <- as.factor(birthwt$race)
birthwt$ptl <- as.factor(birthwt$ptl)
birthwt$ht <- as.factor(birthwt$ht)
birthwt$ui <- as.factor(birthwt$ui)
birthwt <- subset(birthwt, select = -c(id, bwt))
summary(birthwt)

output3 <- bic.glm(low ~ age + lwt + smoke + ptl + ht + ui +
                       ftv + race, glm.family = binomial, data = birthwt)
summary(output3)
output3$postprob
output3$names
output3$probne0

output3$postmean

output3$mle

output3$se