# NOT RUN {
#-------------------linear model----------------------#
# Generate simulated data
n <- 200
p <- 20
k <- 5
rho <- 0.4
seed <- 10
Tbeta <- rep(0, p)
Tbeta[1:k*floor(p/k):floor(p/k)] <- rep(1, k)
Data <- gen.data(n, p, k, rho, family = "gaussian", beta = Tbeta, seed = seed)
x <- Data$x[1:140, ]
y <- Data$y[1:140]
x_new <- Data$x[141:200, ]
y_new <- Data$y[141:200]
lm.bss <- bess(x, y)
lambda.list <- exp(seq(log(5), log(0.1), length.out = 10))
lm.bsrr <- bess(x, y, type = "bsrr", method = "pgsection")
coef(lm.bss)
coef(lm.bsrr)
print(lm.bss)
print(lm.bsrr)
summary(lm.bss)
summary(lm.bsrr)
pred.bss <- predict(lm.bss, newx = x_new)
pred.bsrr <- predict(lm.bsrr, newx = x_new)
# generate plots
plot(lm.bss, type = "both", breaks = TRUE)
plot(lm.bsrr)
#-------------------logistic model----------------------#
#Generate simulated data
Data <- gen.data(n, p, k, rho, family = "binomial", beta = Tbeta, seed = seed)
x <- Data$x[1:140, ]
y <- Data$y[1:140]
x_new <- Data$x[141:200, ]
y_new <- Data$y[141:200]
logi.bss <- bess(x, y, family = "binomial")
lambda.list <- exp(seq(log(5), log(0.1), length.out = 10))
logi.bsrr <- bess(x, y, type = "bsrr", family = "binomial", lambda.list = lambda.list)
coef(logi.bss)
coef(logi.bsrr)
print(logi.bss)
print(logi.bsrr)
summary(logi.bss)
summary(logi.bsrr)
pred.bss <- predict(logi.bss, newx = x_new)
pred.bsrr <- predict(logi.bsrr, newx = x_new)
# generate plots
plot(logi.bss, type = "both", breaks = TRUE)
plot(logi.bsrr)
#-------------------poisson model----------------------#
Data <- gen.data(n, p, k, rho=0.3, family = "poisson", beta = Tbeta, seed = seed)
x <- Data$x[1:140, ]
y <- Data$y[1:140]
x_new <- Data$x[141:200, ]
y_new <- Data$y[141:200]
poi.bss <- bess(x, y, family = "poisson")
lambda.list <- exp(seq(log(5), log(0.1), length.out = 10))
poi.bsrr <- bess(x, y, type = "bsrr",
family = "poisson", lambda.list = lambda.list)
coef(poi.bss)
coef(poi.bsrr)
print(poi.bss)
print(poi.bsrr)
summary(poi.bss)
summary(poi.bsrr)
pred.bss <- predict(poi.bss, newx = x_new)
pred.bsrr <- predict(poi.bsrr, newx = x_new)
# generate plots
plot(poi.bss, type = "both", breaks = TRUE)
plot(poi.bsrr)
#-------------------coxph model----------------------#
#Generate simulated data
Data <- gen.data(n, p, k, rho, family = "cox", scal = 10, beta = Tbeta)
x <- Data$x[1:140, ]
y <- Data$y[1:140, ]
x_new <- Data$x[141:200, ]
y_new <- Data$y[141:200, ]
cox.bss <- bess(x, y, family = "cox")
lambda.list <- exp(seq(log(5), log(0.1), length.out = 10))
cox.bsrr <- bess(x, y, type = "bsrr", family = "cox", lambda.list = lambda.list)
coef(cox.bss)
coef(cox.bsrr)
print(cox.bss)
print(cox.bsrr)
summary(cox.bss)
summary(cox.bsrr)
pred.bss <- predict(cox.bss, newx = x_new)
pred.bsrr <- predict(cox.bsrr, newx = x_new)
# generate plots
plot(cox.bss, type = "both", breaks = TRUE)
plot(cox.bsrr)
#----------------------High dimensional linear models--------------------#
# }
# NOT RUN {
data <- gen.data(n, p = 1000, k, family = "gaussian", seed = seed)
# Best subset selection with SIS screening
lm.high <- bess(data$x, data$y, screening.num = 100)
# }
# NOT RUN {
#-------------------group selection----------------------#
beta <- rep(c(rep(1,2),rep(0,3)), 4)
Data <- gen.data(200, 20, 5, rho=0.4, beta = beta, seed =10)
x <- Data$x
y <- Data$y
group.index <- c(rep(1, 2), rep(2, 3), rep(3, 2), rep(4, 3),
rep(5, 2), rep(6, 3), rep(7, 2), rep(8, 3))
lm.group <- bess(x, y, s.min=1, s.max = 8, type = "bss", group.index = group.index)
lm.groupbsrr <- bess(x, y, type = "bsrr", s.min = 1, s.max = 8, group.index = group.index)
coef(lm.group)
coef(lm.groupbsrr)
print(lm.group)
print(lm.groupbsrr)
#'summary(lm.group)
summary(lm.groupbsrr)
pred.group <- predict(lm.group, newx = x_new)
pred.groupl0l2 <- predict(lm.groupbsrr, newx = x_new)
#-------------------include specified variables----------------------#
Data <- gen.data(n, p, k, rho, family = "gaussian", beta = Tbeta, seed = seed)
lm.bss <- bess(Data$x, Data$y, always.include = 2)
# }
# NOT RUN {
#-------------------trim32 data analysis in doi: 10.18637/jss.v094.i04----------------------#
# import trim32 data by:
load(url('https://github.com/Mamba413/bess/tree/master/data/trim32.RData'))
# or manually downloading trim32.RData in the github page:
# "https://github.com/Mamba413/bess/tree/master/data/" and read it by:
load('trim32.RData')
X <- trim32$x
Y <- trim32$y
dim(X)
# running bess with argument method = "sequential".
fit.seq <- bess(X, Y, method = "sequential")
summary(fit.seq)
# the bess function outputs an 'lm' type of object bestmodel associated
# with the selected best model
bm.seq <- fit.seq$bestmodel
summary(bm.seq)
pred.seq <- predict(fit.seq, newdata = data$x)
plot(fit.seq, type = "both", breaks = TRUE)
# We now call the function bess with argument method = "gsection"
fit.gs <- bess(X, Y, family = "gaussian", method = "gsection")
bm.gs <- fit.gs$bestmodel
summary(bm.gs)
beta <- coef(fit.gs, sparse = TRUE)
class(beta)
pred.gs <- predict(fit.gs, newdata = X)
# }
Run the code above in your browser using DataLab