#################################################################################
### BMTRY 790: Machine Learning								###
### SPRING 2023										###
###													###						
### Lecture 16: Multivariate Adaptive Regression Splines in R using the mda 	###
###		    and earth packaged								###
#################################################################################

### MARS Example
# LOAD LIBRARIES AND DATA
library(mda)
library(earth)

immresp<-read.csv("H:\\public_html\\BMTRY790_Spring2023\\Datasets\\EnvironContamImmuneResp2.csv")
head(immresp)

### Fitting an ADDITIVE MARS model using the mda package
mars.fit1<-mars(immresp[,-10], immresp[,10], degree=1, prune=TRUE, forward.step = TRUE)
mars.fit1$gcv
mars.fit1$coef
mars.fit1$all.terms
mars.fit1$selected.terms
mars.fit1$cuts
mars.fit1$factor

	### PLOTTING ALL INDIVIDUAL PREDICTORS
par(mfrow = c(3, 3), mar=c(4,3,3,2), pty="s")
for (i in 1:9)
  {
    xp <- matrix(sapply(immresp[,-10], mean), nrow(immresp), ncol(immresp) - 1, byrow = TRUE)
    xr <- sapply(immresp, range)
    xp[, i] <- seq(xr[1, i], xr[2, i], len=nrow(immresp))
    xf <- predict(mars.fit1, xp)
    plot(xp[, i], xf, xlab = colnames(immresp)[i], ylab = "", type = "l")
  }


### Fitting an INTERACTION MARS model using the mda package (2-way intx only)
mars.fit2 <- mars(immresp[,c(2:6,9)], immresp[,10], degree = 2, prune = TRUE, forward.step = TRUE)
mars.fit2$gcv
mars.fit2$coef
mars.fit2$all.terms
mars.fit2$selected.terms
mars.fit2$cuts
mars.fit2$factor

	### PLOTTING ALL INDIVIDUAL PREDICTORS
par(mfrow = c(2, 3), mar=c(4,3,3,2), pty="s")
for (i in 1:6)
  {
    xp <- matrix(sapply(immresp[,c(2:6,9)], mean), nrow(immresp), 6, byrow = TRUE)
    xr <- sapply(immresp[,c(2:6,9)], range)
    xp[, i] <- seq(xr[1, i], xr[2, i], len=nrow(immresp))
    xf <- predict(mars.fit2, xp)
    plot(xp[, i], xf, xlab = colnames(immresp)[c(2:6,9)][i], ylab = "", type = "l")
  }



### Fitting an ADDITIVE MARS model using the earth package
mars.fit3<-earth(immresp[,-10], immresp[,10], degree=1)
names(mars.fit3)
mars.fit3$cuts
plot(mars.fit3)


### We can convert a MARS model developed in mda into the format observed in earth using the earth package
mars.fit4<-mars.to.earth(mars.fit2)
mars.fit4
update(mars.fit4)