rm(list = ls())
library(mboost)

# download the data
data <- read.table('http://azzalini.stat.unipd.it/Libro-DM/auto.dat', header = TRUE)
data <- data[data$trazione !='4wd', c(2, 3, 5, 7:18)]
colnames(data) <- c('fuelType', 'aspiration', 'driveWheels', 'wheelBase', 'length', 'width',
                    'height', 'weight', 'displacement', 'compressionRatio', 'horsePower',
                    'maxRevolution', 'cityKmL', 'highwayKmL', 'numCylinders')
data$driveWheels <- factor(data$driveWheels)
data <- data[ , c(1:12, 15, 13)]
head(data)

# split the data in training and test set
set.seed(20190922)
train <- sample(1:nrow(data), 2*nrow(data)/3, replace = FALSE)
test <- (1:nrow(data))[-train]
table(data$aspiration[train], data$fuelType[train])

### linear regression with boosting ###
m1 <- glmboost(cityKmL ~ ., data = data[train, ], center = FALSE, control = boost_control(mstop = 1e3))
m2 <- glmboost(cityKmL ~ ., data = data[train, ], center = TRUE, control = boost_control(mstop = 1e3))

# cross-validation to compute the tuning parameter (number of boosting steps)
cv_1 <- cvrisk(m1, folds = cv(model.weights(m1), type = 'kfold', B = 10), papply = lapply)
cv_2 <- cvrisk(m2, folds = cv(model.weights(m2), type = 'kfold', B = 10), papply = lapply)

# check the model
summary(m1[mstop(cv_1)])
summary(m2[mstop(cv_2)])

# evaluation of the error
mean((data$cityKmL[test] - predict(m1[mstop(cv_1)], newdata = data[test, ]))^2)
mean((data$cityKmL[test] - predict(m2[mstop(cv_2)], newdata = data[test, ]))^2)

# to see the model updating step by step
summary(m2[1])
summary(m2[2])
summary(m2[3])
summary(m2[4])
summary(m2[5])
summary(m2[6])
summary(m2[7])
summary(m2[8])
summary(m2[9])
summary(m2[10])
summary(m2[11])
summary(m2[12])

X <- data.frame(data[train , 1:3], scale(data[train, -c(1:3)], scale = FALSE))
mPlot <- glmboost(cityKmL ~ ., data = X, center = FALSE,
               control = boost_control(mstop = 1e3))
plot(mPlot)


# different loss function
m3 <- glmboost(cityKmL ~ ., data = data[train, -15], center = TRUE, control = boost_control(mstop = 1e3),
               family = Laplace())
cv_3 <- cvrisk(m3, folds = cv(model.weights(m3), type = 'kfold', B = 10), papply = lapply)
mean((data$cityKmL[test] - predict(m3[mstop(cv_3)], newdata = data[test, ]))^2)


### generalized additive model with boosting ###
# paste0(colnames(data), collapse = ') + bols(')

m4 <- gamboost(cityKmL ~ bols(fuelType) + bols(aspiration) + bols(driveWheels) + bols(wheelBase) +
                 bols(length) + bols(width) + bols(height) + bols(weight) + bols(displacement) +
                 bols(compressionRatio) + bols(horsePower) + bols(maxRevolution) + bols(numCylinders),
               data = data[train, ], control = boost_control(mstop = 1e3))
cv_4 <- cvrisk(m4, folds = cv(model.weights(m4), type = 'kfold', B = 10), papply = lapply)
mean((data$cityKmL[test] - predict(m4[mstop(cv_4)], newdata = data[test, ]))^2)
coef(m4[mstop(cv_4)])

m4b <- gamboost(cityKmL ~ bols(fuelType, intercept = FALSE) + bols(aspiration, intercept = FALSE) +
                 bols(driveWheels, intercept = FALSE) + bols(wheelBase, intercept = FALSE) +
                 bols(length, intercept = FALSE) + bols(width, intercept = FALSE) +
                 bols(height, intercept = FALSE) + bols(weight, intercept = FALSE) +
                 bols(displacement, intercept = FALSE) + bols(compressionRatio, intercept = FALSE) +
                 bols(horsePower, intercept = FALSE) + bols(maxRevolution, intercept = FALSE) +
                 bols(numCylinders, intercept = FALSE) + bols(X1, intercept = FALSE),
               data = data.frame(1, X), control = boost_control(mstop = 1e3))
cv_4b <- cvrisk(m4b, folds = cv(model.weights(m4b), type = 'kfold', B = 10), papply = lapply)
coef(m4b[mstop(cv_4b)])


m5 <- gamboost(cityKmL ~ bols(fuelType) + bols(aspiration) + bols(driveWheels) + bols(wheelBase) +
                 bols(length) + bols(width) + bols(height) + bols(weight) + bols(displacement) +
                 bols(compressionRatio) +
                 bols(horsePower) + bbs(horsePower, center = FALSE) +
                 bols(maxRevolution) + bols(numCylinders),
               data = data[train, ], control = boost_control(mstop = 1e3))
cv_5 <- cvrisk(m5, folds = cv(model.weights(m5), type = 'kfold', B = 10), papply = lapply)
mean((data$cityKmL[test] - predict(m5[mstop(cv_5)], newdata = data[test, ]))^2)


### different outcome ###

# logistic regression
dataBin <- read.csv('https://archive.ics.uci.edu/ml/machine-learning-databases/00451/dataR2.csv',
                    header = TRUE)
head(dataBin)
# split in training and test data 
train <- c(sample(which(dataBin$Classification == '1'), 2*sum(dataBin$Classification == '1')/3),
           sample(which(dataBin$Classification == '2'), 2*sum(dataBin$Classification == '2')/3))

# model with exponential loss (AdaBoost)
dataBin$Classification <- as.factor(dataBin$Classification)
m6 <- glmboost(Classification ~ ., data = dataBin[train, ], family = AdaExp(),
               control = boost_control(mstop = 1e3))
cv_6 <- cvrisk(m6, folds = cv(model.weights(m6), type = 'kfold', B = 10), papply = lapply)
summary(m6[mstop(cv_6)])

mean(as.factor(dataBin$Classification)[-train] != predict(m6[mstop(cv_6)],
                                                          newdata = dataBin[-train, ], type = 'class'))

# logistic
m7 <- glmboost(Classification ~ ., data = dataBin[train, ], family = Binomial(),
               control = boost_control(mstop = 1e3))
cv_7 <- cvrisk(m7, folds = cv(model.weights(m7), type = 'kfold', B = 10), papply = lapply)
summary(m7[mstop(cv_7)])

mean(as.factor(dataBin$Classification)[-train] != predict(m7[mstop(cv_7)],
                                                          newdata = dataBin[-train, ], type = 'class'))

# survival analysis
library('survival')
download.file('https://www.imbi.uni-freiburg.de/imbi/Royston-Sauerbrei-book/Multivariable_Model-building/downloads/datasets/giloma.zip',
              'dataGlioma.zip')
glioma <- read.csv(unz('dataGlioma.zip', 'giloma/glioma.csv'), header = TRUE)
head(glioma)
# transform the response variable in a survival object
dataSurv <- glioma[ , 1:18]
dataSurv[ , 2] <- dataSurv[ , 2] - 1

# Cox model with linear predictor
m8 <- glmboost(Surv(survtime, cens) ~ ., data = dataSurv, family = CoxPH())
cv_8 <- cvrisk(m8, papply = lapply)
summary(m8[mstop(cv_8)])

# additive Cox model
m9 <- gamboost(Surv(survtime, cens) ~ btree(sex) + btree(time) + btree(gradd1) + btree(gradd2) +
                 bbs(age) + btree(kard1) + btree(kard2) + btree(surgd1) + btree(surgd2) + btree(convul) +
                 btree(cort) + btree(epi) + btree(amnesia) + btree(ops) + btree(aph) + btree(trt),
               data = dataSurv, family = CoxPH())
cv_9 <- cvrisk(m9, papply = lapply)
summary(m9[mstop(cv_9)])


### high-dimensional case ###
library(golubEsets)
data("Golub_Merge")
dim(Golub_Merge)
Golub_Merge
X <- t(exprs(Golub_Merge))
# extract clinical data
y <- pData(Golub_Merge)[ , 2]

X <- scale(X, scale = FALSE)
mHD <- glmboost(y ~ ., data = data.frame(y, X), family = Binomial(type = 'glm', link = 'logit'))
cv_HD <- cvrisk(mHD, folds = cv(model.weights(mHD), type = 'kfold', B = 10), papply = lapply)
summary(mHD[mstop(cv_HD)])