#Neural network on spam data

library(MASS)
library(nnet)
library(class)

fil = "http://www.uio.no/studier/emner/matnat/math/STK2100/data/spam_data.txt"
spam = read.table(fil,header=T)
summary(spam)
#Transform input data due to long-tailed distribution
spam[,1:57] = log(spam[,1:57]+0.1)

#Only first three covariates

fit.glm3 = glm(y~x1+x2+x3,subset=spam$train,data=spam,family=binomial)
pred.glm3 = predict(fit.glm3,spam[!spam$train,],type="response")>0.5
err.glm3 = mean(spam$y[!spam$train]!=pred.glm3)
library(gam)
fit.gam3 = gam(y~s(x1)+s(x2)+s(x3),subset=spam$train,data=spam,family=binomial)
pred.gam3 = predict(fit.gam3,spam[!spam$train,],type="response")>0.5
err.gam3 = mean(spam$y[!spam$train]!=pred.gam3)

show(c(err.glm3,err.gam3))

#All covariates

fit = glm(y~.-train,data=spam,family=binomial,subset=spam$train)
pred = predict(fit,spam[!spam$train,],type="response")>0.5
err.logist=mean(spam$y[!spam$train]!=pred)

nam = names(spam)[1:57]
k = 57
formula = as.formula(paste("y",paste(paste("s(",nam[1:k],")",sep=""), 
                                     collapse="+"), sep="~"))
fit.gam = gam(formula,data=spam,subset=spam$train,family=binomial)
pred.gam = predict(fit.gam,spam[!spam$train,],type="response")>0.5
err.gam=mean(spam$y[!spam$train]!=pred.gam)

show(c(err.logist,err.gam))

#Using variables selected in textbook
ind = c(5:8,16:17,25:27,37,45:46,52:53,56:57)

formula = as.formula(paste("y",paste(paste(nam[ind],sep=""), 
                                     collapse="+"), sep="~"))
fit.glm.var = glm(formula,data=spam,subset=spam$train,family=binomial)
pred.glm.var = predict(fit.glm.var,spam[!spam$train,],type="response")>0.5
err.glm.var=mean(spam$y[!spam$train]!=pred.glm.var)

formula = as.formula(paste("y",paste(paste("s(",nam[ind],")",sep=""), 
                                     collapse="+"), sep="~"))
fit.gam.var = gam(formula,data=spam,subset=spam$train,family=binomial)
pred.gam.var = predict(fit.gam.var,spam[!spam$train,],type="response")>0.5
err.gam.var=mean(spam$y[!spam$train]!=pred.gam.var)

show(c(err.glm.var,err.gam.var))

#Changing loss
table(spam$y[!spam$train],pred.gam.var)
L01=10;L10=1
crit = L10/(L01+L10)
pred.gam.var2 = predict(fit.gam.var,spam[!spam$train,],type="response")>crit
err.gam.var2=mean(spam$y[!spam$train]!=pred.gam.var2)
table(spam$y[!spam$train],pred.gam.var2)

#Using principal components
x.prcomp = prcomp(spam[,1:57],retx=TRUE,scale=TRUE)
d = data.frame(x.prcomp$x,y=spam$y,train=spam$train)
k = 15
fit.pc = glm(y~.-train,data=d[,c(1:k,58,59)],family=binomial,subset=d$train)
pred.pc = predict(fit.pc,d[!d$train,],type="response")>0.5
err.pca = mean(d$y[!d$train]!=pred.pc)
show(err.pca)

nam = names(d)[1:57]
formula = as.formula(paste("y",paste(paste("s(",nam[1:k],")",sep=""), 
                                     collapse="+"), sep="~"))
fit.gam.pca = gam(formula,data=d,family=binomial,subset=d$train)
pred.gam.pca = predict(fit.gam.pca,d[!d$train,],type="response")>0.5
err.gam.pca = mean(d$y[!d$train]!=pred.gam.pca)
plot(fit.gam.pca)
plot(fit.gam.pca,se=T)
show(c(err.pca,err.gam.pca))


#c)
plot(fit.gam3,se=T)
#d)
fit.pc.glm3 = glm(y~PC1+PC2+PC3,subset=d$train,data=d,family=binomial)
pred.pc.glm3 = predict(fit.pc.glm3,d[!d$train,],type="response")>0.5
err.pc.glm3 = mean(d$y[!d$train]!=pred.pc.glm3)
library(gam)
fit.pc.gam3 = gam(y~s(PC1)+s(PC2)+s(PC3),subset=d$train,data=d,family=binomial)
pred.pc.gam3 = predict(fit.pc.gam3,d[!d$train,],type="response")>0.5
err.pc.gam3 = mean(d$y[!d$train]!=pred.pc.gam3)

plot(fit.pc.gam3,se=T)

plot.ts(x.prcomp$rotation[,1])
plot.ts(x.prcomp$rotation[,2])
plot.ts(x.prcomp$rotation[,3])
plot.ts(x.prcomp$rotation[,4])
plot.ts(x.prcomp$rotation[,5])
plot.ts(x.prcomp$rotation[,6])
plot.ts(x.prcomp$rotation[,7])
plot.ts(x.prcomp$rotation[,8])
#e)
k=20
nam = names(d)[1:57]
formula = as.formula(paste("y",paste(paste("s(",nam[1:k],")",sep=""), 
                                     collapse="+"), sep="~"))
fit.pc.gam = gam(formula,subset=d$train,data=d,family=binomial)
pred.pc.gam = predict(fit.pc.gam,d[!d$train,],type="response")>0.5
err.pc.gam = mean(d$y[!d$train]!=pred.pc.gam)



#f)
err.pca.ns = rep(NA,57)
for(k in 1:57)
{
  print(k)
  formula = as.formula(paste("y",paste(paste("s(",nam[1:k],")",sep=""), 
                                          collapse="+"), sep="~"))
  fit2 = gam(formula,data=d,subset=d$train,family=binomial)
  pred2 = predict(fit2,d[!d$train,],type="response")>0.5
  err.pca.ns[k] = mean(d$y[!d$train]!=pred2)
  rm(fit2,pred2)
  gc()
}
plot(1:57,err.pca.ns)
kmin = which.min(err.pca.ns)

show(err.pca.ns[kmin])


plot(fit2,se=T)