################
#
# Swiss fertility data
#   47 French speaking provinces of Switzerland around 1888
#
#           Fertility     - standardized fertility measure for each province 
#           Agriculture   - population involved in agriculture 
#           Examiniation  - draftees receiving highest mark in army exam 
#           Education     - population educated beyond primary school 
#           Catholic      - population who are catholic 
#           Infant.Mort   - live births living less than 1 year 
    
# The interest is to relate Fertility (the response) to the other variables
# (the predictors).  
# The hypothesis is that the predictors are good proxies for the true causes
# of high and low fertility.
# Note our cautious formulation: these variables are probably not direct
# causes of fertility but part of a cultural mix that overall determines
# fertility rates.
    par(mfrow=c(1,1))
	par(mar=c(4,4,2,2))
    data(swiss)               # The dataset comes with R.
    options(width=100)        # I like a longer line for a wide window.
    swiss                     # It's a small dataset, so we can look at it.
    names(swiss)              # The variables.
    dim(swiss)                # The size of the data.
# For my own interest, I wanted to know the names of some of these provinces,
# especially those that are high in Education and Fertility.
# To this end, we obtain an index vector "ord" with the "order" function
# to sort the data, first according to Education, then Fertility:
# - High Education:
    ord <- rev(order(swiss[,"Education"]))
    rownames(swiss)[ord][1:7]                 # 7 highest on education
# ...includes Geneva, Neuchatel, Lausanne, Vevey, which are "larger" cities
# in French speaking Switzerland.  I believe Rive Gauche and Rive Droite
# are suburbs of Geneva (left and right of the river Rhone), but I'm not sure.
# - High Fertility:
    ord <- rev(order(swiss[,"Fertility"]))
    rownames(swiss)[ord][1:10]               # 10 highest on fertility
# ...names that seem rural or alpine to me (I forgot a lot about Switzerland).
# Look at a map of Switzerland.  "   
# - Low Fertility:
    ord <- order(swiss[,"Fertility"])
    rownames(swiss)[ord][1:10]
# ...includes the larger cities.   
# Now some basic plots: response (Fertility versus predictors):
   par(pch=16, mfrow=c(1,5), mar=c(2,2,1.5,0), mgp=c(1.5,.5,0),cex=1, cex.axis=.5, cex.main=.6)
    for(i in 2:6)
      plot(swiss[,i], swiss[,1], main=colnames(swiss)[i], xlab="", ylab="")
# If these marginal plots are any indication, we can expect
# Education and its proxy, Examination, to be negatively correlated
# with Fertility.
# On the other hand, Agriculture and Infant.Mort seem to be
# positively correlated; which may make sense because these are
# generally indicators of a low developmental state ("backwardness").
# An interesting connection exists with Catholic: highly catholic
# provinces have higher fertililty than those that are highly protestant,
# with the exception of three provinces that are highly mixed, and probably
# more cosmopolitan: they are Geneva, Rive Gauche and Rive Droite:    
# The easiest approach is with a scatterplot matrix:
    pairs(swiss, pch=16, cex.labels=1.5, gap=.1)
#   (What a nuisance!  The extreme points fall virtually on the borders.
#    Here is a version that allocates some space between 
#    the points and the borders: download it from the class website.)
	source("pairsplus")
    pairs.plus(swiss, pch=16, cex.labels=1.5, gap=.1)
# Right away we see that Examination and Education are indeed
# positively correlated, that Agriculture negatively correlates
# with both, although more so with Examination, and that
# Infant.Mort does not correlate with anything very much, except for a
# positive correlation with Catholic.
# Finally, here is an application of the "tree()":
library(tree)
source("Trees541")
swiss.Scart <- tree(Fertility ~ ., data=swiss,control=tree.control(47, mincut=3, minsize=6))
show.tree(swiss.Scart, type="u")
# Observations:
# - The top split breaks out a group of 7 provinces (=14.9% of 47) with
#   high education levels and low fertility.
# - In the low Education bucket, the tree splits on Catholic with the
#   expected direction: the high Catholic bucket has higher fertility.
# - The next level shows splits on Infant.Mort on both sides, with
#   the expected direction: higher Infant.Mort, higher Fertility.
#   (This latter relation may not be a causal left to right,
#   equally likely the opposite: right to left, or both ways.)
# - The last split on Agriculture is not very believable as it has the
#   unexpected direction: higher Agriculture, lower Fertility;
#   But then the difference is the smallest of any split in the tree,
#   which adds evidence that the split is probably not "real".
    
# Overall, though, the standard tree function gave us in this case a
# very interpretable tree.  This may be due to the small size of the dataset.
    
# Our own CART routine does not exactly reproduce the tree generated
# by the standard library (there are differences in the code):
    swiss.CART <- regression.tree(swiss[,c(2:5,1)], method="CART",
                                  minsize=5, title="Swiss Fertility Data")
    show.tree(swiss.CART, type="u")
# The tree is interesting because it seems to suggest that
# Agriculture does in fact have an unexpected correlation with
# fertility after adjustment for Education and Catholic!
    
# Now we use another option in our routine to look for high Fertility buckets:
    swiss.XM <- regression.tree(swiss[,c(2:5,1)], method="XM",
                                  minsize=5, title="Swiss Fertility Data")
    show.tree(swiss.XM, type="u") 
# We see that for finding high Fertility, one should look for low Examination 
# above all and high Catholic next.
# Recall that Examination is an education proxy,
# so it is surprising that deeper down Education also appears.
# Overall we learn that for high Fertility we should look for the
# lowest levels of the education proxies Examination and Education.
    
# Next we look for low fertility:
    swiss.NM <- regression.tree(swiss[,c(2:5,1)], method="NM",
                                  minsize=5, title="Swiss Fertility Data")
    show.tree(swiss.NM, type="u") 
    row.names(swiss)[swiss[,"Education"] > 23.2]
# Again, at the top a high Education bucket is broken off,
# followed by a high Examination bucket,
# followed by two buckets with low Catholic.
#
# Interesting:
# There are two splits on Agriculture, once again in the unexpected
# direction: high Agriculture - low Fertility.
# Do we have to conclude that all things being equal (education, religion),
# agriculture in 1888 was such a taxing environment that it
# depressed fertility?
# Note that we say "all things being equal", which means in formal terms:
#   "after adjusting for other variables",
# or equivalently for trees,
#   "after forming buckets with small variability in the other variables". 
    
# Looking ahead, let's use a linear model fit to see whether the
# overall qualitative behaviors are confirmed, in particular that
# of Agriculture:
    summary(lm(Fertility ~ ., data=swiss))
# (We'll say more about the linear model function "lm" later.)
#------------------------ output begin
     ...
    Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
    (Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
    Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
    Examination      -0.25801    0.25388  -1.016  0.31546    
    Education        -0.87094    0.18303  -4.758 2.43e-05 ***
    Catholic          0.10412    0.03526   2.953  0.00519 ** 
    Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
     ...
#------------------------ output end
# Yes, Agriculture has a significant negative slope!
# Note another feature of the linear model, though:
# Examination is insignificant!
# After looking at the trees and the plots we know better than that.
# Let us remind ourselves of a couple of plots:
    windows(width=7, height=14)
    par(mfrow=c(2,1), mgp=c(1.5,.5,0), mar=c(2.5,2.5,1,1))
    plot(swiss[,c("Examination","Fertility")], pch=16)
    plot(swiss[,c("Education","Fertility")], pch=16)
# Education comes in so strong because of the 7 high Education provinces
# By comparison, Examination has a weaker but more robust correlation
# with Fertility: no 7 provinces would account for the correlation.
# More on this below.