###################################################
# Predicting fat content from spectrometric curves: 
# R/S-plus commandlines
###################################################
#
# Entering spectrometric data 
#############################
SPECDAT <- as.matrix(read.table("npfda-spectrometric.dat"))
attributes(SPECDAT)$dimnames <- NULL
SPECURVES <- SPECDAT[,1:100]           # sample of curves
Response <- SPECDAT[,101]              # sample of responses
learning <- 1:160
testing <- 161:215
SPECURVES1 <- SPECURVES[learning,]     # learning sample of curves
SPECURVES2 <- SPECURVES[testing,]      # testing sample of curves
Specresp1 <- Response[learning]        # learning sample of responses
Specresp2 <- Response[testing]         # testing sample of responses
##
# Computing predicted fat content
#################################
result.reg <- funopare.knn.lcv(Specresp1, SPECURVES1, SPECURVES2, 2, nknot=20, c(0,1))
result.mode <- funopare.mode.lcv(Specresp1, SPECURVES1, SPECURVES2, 2, nknot=20, c(0,1))
result.quantile <- funopare.quantile.lcv(Specresp1, SPECURVES1, SPECURVES2, 2, nknot=20, c(0,1), alpha=0.5)
##
# Computing square errors
#########################
Se.reg <- (Specresp2 - result.reg$Predicted.values)^2
Se.mode <- (Specresp2 - result.mode$Predicted.values)^2
Se.quantile <- (Specresp2 - result.quantile$Predicted.values)^2
mse.reg <- round(sum(Se.reg)/length(Specresp2),2)
mse.mode <- round(sum(Se.mode)/length(Specresp2),2)
mse.quantile <- round(sum(Se.quantile)/length(Specresp2),2)
##
# Plotting predicted responses (don't forget to active your graphics device!!)
##############################################################################
par(mfrow=c(2,2), pty="s") # figure will be drawn row-by-row in an 1 by 3
                           # matrix on the current graphical device,
                           # pty = "s"  generates a square plotting region.
xlimits <- range(Specresp2)
ylimits <- range(c(result.mode$Predicted.values, result.quantile$Predicted.values, result.reg$Predicted.values))
plot(Specresp2, result.reg$Predicted.values, xlim=xlimits, ylim=ylimits, xlab="Responses of testing sample", ylab="Predicted responses")
mtext(paste("Cond. Expect.: MSE=",mse.reg,""), cex=1.15, line=1)
abline(0,1) 
plot(Specresp2, result.mode$Predicted.values, xlim=xlimits, ylim=ylimits, xlab="Responses of testing sample", ylab="Predicted responses")
mtext(paste("Cond. Mode: MSE=",mse.mode, sep=""), cex=1.15, line=1)
abline(0,1)
plot(Specresp2, result.quantile$Predicted.values, xlim=xlimits, ylim=ylimits, xlab="Responses of testing sample", ylab="Predicted responses")
mtext(paste("Cond. Median: MSE=",mse.quantile,""), cex=1.15, line=1)
abline(0,1)
##
# Improving the predictions by averaging the 
# three predicted values (three methods)
############################################
result.multimethods <- (result.mode$Predicted.values + result.quantile$Predicted.values + result.reg$Predicted.values)/3
Se.multimethods <- (Specresp2 - result.multimethods)^2
mse.multimethods <- round(sum(Se.multimethods)/length(Specresp2),2)
error.names <- c(paste("Cond. Expect.\n\n MSE=",mse.reg,sep=""), paste("Cond. Mode\n\n MSE=",mse.mode,sep=""), paste("Cond. Median\n\n MSE=",mse.quantile,sep=""), paste("multimethods\n\n MSE=",mse.multimethods,sep=""))
par(mfrow=c(1,1))
boxplot(Se.reg, Se.mode, Se.quantile, Se.multimethods, names = error.names)
title("Square Error")