################################################################
################################################################
#
#    CART -- Classification and Regression Trees
#
# The idea:
#
#   When data are heterogeneous, clustered, messy,
#     it might not be the best idea to fit smooth functions;
#   in particular, it might not be a good idea to fit linear functions,
#     as is done in linear regression.
#   Instead, one might consider breaking the data up into
#     subsets that are homogeneous.
#   Once we have homogeneous subsets, we might do the simplest possible thing:
#     fit a constant, that is, a mean.
#   Therefore, tree-based regression does local averaging, like smoothers,
#     but this time on non-overlapping subsets.
#
#   Q: How do we arrive at homogeneous subsets?
#   A: By repeatedly breaking the data up along axis splits.
#
#   But how would we choose such axis splits?
#   Example: In the following plot, where would we reasonably place a split?
#
     library(MASS)
	  data(Boston)
     attach(Boston)  #rende disponibile il data set nella worksppace
	  par(mfrow=c(1,1),mar=c(4,4,2,2))
	  nvar <- ncol(Boston) 
      plot(Boston[,"lstat"], Boston[,"medv"],pch=21, xlab="lstat", ylab="medv",col="blue")
#   Here is a handchosen split at lstat=10:
      spl <- 10
      sel.left  <- Boston[,"lstat"] < spl
      sel.right <- Boston[,"lstat"] >= spl
      plot(Boston[,"lstat"], Boston[,"medv"], xlab="lstat", ylab="medv",col=c("red","green")[sel.left+1], pch=16)
#   We said we will fit only the simplest possible model, a  mean in each subset:
      lines(c(spl,spl), c(0,50),lwd=2,lty=2,col="blue")
      lines(c(0,spl),  rep(mean(Boston[sel.left,"medv"]), 2),  lwd=2,lty=2,col="seagreen4")
      lines(c(spl,40), rep(mean(Boston[sel.right,"medv"]), 2), lwd=2,lty=2,col="cyan4")
# 
# How do we formalize a rule for selecting the axis splits? That's the big question.
# Here is a recipe that is used most often:
#
#   For a given split on a variable "xj" (="lstat" in the plot) and threshold "th" (=10 in the plot),
#   evaluate a criterion, such as the residual sum of squares,
#   after fitting a separate mean left and right of the threshold:
      sum((Boston[sel.left,"medv"]  - mean(Boston[sel.left,"medv"]))^2) +
      sum((Boston[sel.right,"medv"] - mean(Boston[sel.right,"medv"]))^2)
#   which should be the same as
      var(Boston[sel.left,"medv"])  * (sum(sel.left)-1) +
      var(Boston[sel.right,"medv"]) * (sum(sel.right)-1)
#   Indeed!  Then:
#   For all possible splits (505 for "boston")  of all x-variables (13 for "Boston"),
#   evaluate the criterion, and pick the split that has the lowest value of the criterion:
#   For example
      lstat.vals <- sort(Boston[,"lstat"])
      crit.vals <- rep(NA,length(lstat.vals))
      for (i in 1:length(lstat.vals)){
        spl <- lstat.vals[i]
        sel.left  <- Boston[,"lstat"] < spl; sel.right <- !sel.left
        crit.vals[i] <-
        sum((Boston[sel.left,"medv"]  - mean(Boston[sel.left,"medv"]))^2) +
        sum((Boston[sel.right,"medv"] - mean(Boston[sel.right,"medv"]))^2)
      }
    plot(lstat.vals, crit.vals, pch=16, cex=.4)
#   Look at that: spl=10 was near the minimum of the criterion values!
#   Then:
#   Repeat on both sides of the split.
#   repeat on and on and on..., recursively.
#   Stop when the subsets get too small.
# 	It acts like a divisive clustering
#
#    Cavalli-Sforza split
	  Spli.lev<-Inf;m1<-nvar-1
	  for (j in 1:m1){
	  	    lstat.vals <- sort(Boston[,j])
     		crit.vals <- Inf
     		 for (i in 1:length(lstat.vals)){
        			spl <- lstat.vals[i]
        			sel.left  <- Boston[,j] < spl; sel.right <- !sel.left
       				Temp <-
         		sum((Boston[sel.left,"medv"]  - mean(Boston[sel.left,"medv"]))^2) +
            	sum((Boston[sel.right,"medv"] - mean(Boston[sel.right,"medv"]))^2)
            	if (Temp<crit.vals) {crit.vals<-Temp;ipos<-i}
      			}
      		if (crit.vals<Spli.lev) {Spli.lev<-crit.vals;ivar<-j}
      		}
      	cat("Best split at ",Spli.lev,"  for variable ",ivar,"\n")
#      Nuovo grafico
      a<-names(Boston)
      spl <- Spli.lev
      sel.left  <- Boston[,ivar] < spl
      sel.right <- Boston[,ivar] >= spl
      plot(Boston[,14],Boston[,ivar],pch=21,cex=.8,xlab=a[ivar],ylab="medv",col=c("red","green")[sel.left+1])
#
# There are many complications to the idea, but let's leave things simple.
# We are going to use trees not for prediction but for interpretation.
#   
# To be teached after a while
library(tree)
source("Trees541.r")
#
# Here is the conventional "CART" method:
boston.CART <- regression.tree(Boston, method="CART",minsize=23, title="Boston Housing Data")
show.tree(boston.CART, type="u", cex=.8) 
# [Stretch this window as much as possible to avoid much of the overplotting.]
#
# The tree looks interesting and we should stare at it for a while
# to get a sense of what it might tell us:
#   The first split is on Rm, the average number of rooms:
#     85% with Rm < 6.94,      15% with Rm > 6.94
#     [small/normal homes]     [large homes]
#     The 85% of smaller homes is split on Lstat, % of lower status people,
#       51% with Lstat < 14.55,  15% with Lstat > 14.55
#       [^^few poor people^^^]   [^^more poor people^^]
#     The 15% of larger homes is split on Rm again,
#       9.1% with Rm < 7.42,    5.9% with Rm > 7.42
#       [medium large homes]    [largest homes] 
# This makes some sense:
#   among large homes, size drives the price;
#   among small/normal homes, social environment drives the price.
# We had found something like this when studying scatterplots of Medv, Rm, Lstat.
# The variable of interest was really Nox, but it is only of  low importance in the tree.
#  Let's move on to one of our new methods:
# It tries to repeatedly peel off subsets that have HIGH response values:
	boston.XM <- regression.tree(Boston, method="XM", minsize=23, title="Boston Housing Data")
	show.tree(boston.XM, type="u", cex=.8) 
# The story becomes clearer: Lstat is the most important variable, but
#   to predict where the highest prices are, one has to look at Rm.  
#
# 
# And here is another one of our new methods:
# It does the opposite by trying to repeatedly peel off subsets with LOW response values:
    boston.NM <- regression.tree(Boston, method="NM",minsize=23, title="Boston Housing Data")
    show.tree(boston.NM, type="u", cex=.8) 
# Another story emerges, complementary to the first story, and more complete, more intuitive!!
#   In order to ascend housing prices, from lowest to highest, peel off high crime,
#     then racially segregated areas,
#     then descend down the % of lower status people,
#     then climb up the number of rooms.
# Among the three trees we constructed, this one makes the most sense.
#    