house5a.Rmd

---
title: "House Price Prediciton with Linear Ensemble and SVM"
output:
  html_document:
    keep_md: yes
---

Set this variable appropriately before running the rest:
```{r}
running_as_kaggle_kernel = FALSE

```

This approach uses a simple 50/50 logarithmic average of my linear predictor
and mtyxwp's simple support vector machine (copied from "svm_simple" kernel).
My linear predictor does quite a bit of pre-processing using ordinary least squares
as a tool and then feeds the resulting model into an ensemble of fitting methods,
with weights for combining the methods determined from the fit on a validation set.

Start by loading libraries.
```{r}
library(plyr)
library(caret)
library(Metrics)
library(parallel)
library(doParallel)
library(ggplot2) # Data visualization
library(readr) # CSV file I/O, e.g. the read_csv function
library(Amelia)
library(mice)
library(lattice)
library(rpart)
library(xgboost)
library(e1071)

```

Here are the linear fitting methods.  Many of these will be excluded from the final
predictions, if they don't seem to help with performance on the validation set.
And one will be excluded for performance reasons if running on Kaggle.

```{r}
modelnames = c("lars2",     # Least Angle Regression
               "cubist",    # Cubist Regression Tree
               "glmboost",  # Boosted Generalized Linear Model
               "glmnet",    # Generalized linear model via penalized maximum likelihood
               "lasso",     # Least absolute shrinkage & selection operator (L1 penalty)
               "bayesglm",  # Bayesian Generalized Linear Model
               "ridge",     # Ridge Regression (L2 penalty)
               "xgbLinear", # eXtreme Gradient Boosting, linear method
               "nnls",      # Non-Negative Least Squares
               "icr",       # Independent Component Regression
               "gbm")       # Stochastic Gradient Boosting

too_slow <- c("xgbLinear")

```

Read the data files, and make any environment-specific adjustments.
```{r}
if (running_as_kaggle_kernel) {
  trainfile <- "../input/train.csv"
  testfile <- "../input/test.csv"
  keep = ! (modelnames %in% too_slow)
  modelnames <- modelnames[ keep ]
} else {
  trainfile <- "train.csv"
  testfile <- "test.csv"
}

```


## Functions

For factors that have an explicit or implicit order (e.g. Poor, Fair, Good, etc.), I am going to recode them as integers at first, so the order doesn't get lost.  Also,
when I call this function, I will be careful to choose the ordering that is expected
(based on intuition, but usually obvious) to be positively associated with sale price.
You'll see why later.
```{r}
# Function to recode levels to numeric in specified order and add ".n" to name
recode <- function( df, var, lev ) { 
  to <- as.character( 0:(length(lev)-1) )
  newvar <- as.numeric( as.character( mapvalues(df[[var]], from=lev, to=to) ) )
  newname <- paste0(var,".n")
  df <- cbind( df, newvar )
  names(df)[ncol(df)] <- newname
  df[var] <- NULL
  df
}

```

Function for initial data cleaning.  (Much of this could be done more elegantly
and consistently, but I confess I didn't start out this project by making sure 
the data were clean, as I should have.)
```{r}
cleanData <- function( df ) {
  # Convert pseudo-numeric "type of dwelling" identifier to a factor
  df$MSSubClass <- as.factor( df$MSSubClass )

  # Deal with numeric variables that have missing values
  df$LotFrontage = as.character( df$LotFrontage )
  df$HasLotFrontage = ifelse( df$LotFrontage=="NA", 0, 1 )
  df$LotFrontage = ifelse( df$LotFrontage=="NA", "0", df$LotFrontage ) 
  df$LotFrontage = as.numeric( df$LotFrontage )

  df$MasVnrArea = as.character( df$MasVnrArea )
  df$HasMasVnr = ifelse( df$MasVnrArea=="NA", 0, 1 )
  df$MasVnrArea = ifelse( df$MasVnrArea=="NA", "0", df$MasVnrArea ) 
  df$MasVnrArea = as.numeric( df$MasVnrArea )

  df$GarageYrBlt = as.character( df$GarageYrBlt )
  df$HasGarageYr = ifelse( df$GarageYrBlt=="NA", 0, 1 )
  df$GarageYrBlt = ifelse( df$GarageYrBlt=="NA", "0", df$GarageYrBlt ) 
  df$GarageYrBlt = as.numeric( df$GarageYrBlt )

  # Dummy for "has basement"
  df$HasBasement = ifelse( df$BsmtQual=="NA", 0, 1 )

  # Recode ordered factors as pseudo-continuous numerical variables
  qualcats  = c( "Po",  "Fa",  "TA",   "Gd",   "Ex" )
  qualcats2 = c( "NA",  qualcats )
  funcats   = c( "Sal", "Sev", "Maj2", "Maj1", "Mod", "Min2", "Min1", "Typ" )
  basecats  = c( "NA",  "Unf", "LwQ",  "Rec",  "BLQ", "ALQ",  "GLQ"         )
  df <- recode( df, "ExterCond",    qualcats  )
  df <- recode( df, "ExterQual",    qualcats  )
  df <- recode( df, "HeatingQC",    qualcats  )
  df <- recode( df, "KitchenQual",  qualcats  )
  df <- recode( df, "BsmtCond",     qualcats2 )
  df <- recode( df, "FireplaceQu",  qualcats2 )
  df <- recode( df, "GarageQual",   qualcats2 )
  df <- recode( df, "GarageCond",   qualcats2 )
  df <- recode( df, "Functional",   funcats   )
  df <- recode( df, "BsmtFinType1", basecats  )
  df <- recode( df, "BsmtFinType2", basecats  )
  df <- recode( df, "PavedDrive",   c("N",   "P",      "Y"                     ) )                                         
  df <- recode( df, "Utilities",    c("ELO", "NoSeWa", "NoSewr", "AllPub"      ) )
  df <- recode( df, "LotShape",     c("IR3", "IR2",    "IR1",    "Reg"         ) )                                         
  df <- recode( df, "BsmtExposure", c("NA",  "No",     "Mn",     "Av",    "Gd" ) )
  df <- recode( df, "PoolQC",       c("NA",  "Fa",     "TA",     "Gd",    "Ex" ) )
  df <- recode( df, "GarageFinish", c("NA",  "Unf",    "RFn",    "Fin"         ) )

  # BsmtHeight needs special treatment, since it's really a categorized continuous variable
  from <- c("NA", "Po", "Fa", "TA", "Gd", "Ex"  ) 
  to   <- c("0",  "50", "75", "85", "95", "120" )                                          
  df$BsmtHeight <- as.numeric( mapvalues(df$BsmtQual, from=from, to=to) )
  df$BsmtQual <- NULL

  # Fix numeric variables that will get read as factors in test set
  df$BsmtFinSF1 = as.numeric( as.character( df$BsmtFinSF1 ) )
  df$BsmtFinSF1[is.na(df$BsmtFinSF1)] = mean(df$BsmtFinSF1, na.rm=TRUE)
  df$BsmtFinSF2 = as.numeric( as.character( df$BsmtFinSF2 ) )
  df$BsmtFinSF2[is.na(df$BsmtFinSF2)] = mean(df$BsmtFinSF2, na.rm=TRUE)
  df$BsmtUnfSF = as.numeric( as.character( df$BsmtUnfSF ) )
  df$BsmtUnfSF[is.na(df$BsmtUnfSF)] = mean(df$BsmtUnfSF, na.rm=TRUE)
  df$BsmtFullBath = as.numeric( as.character( df$BsmtFullBath ) )
  df$BsmtFullBath[is.na(df$BsmtFullBath)] = mean(df$BsmtFullBath, na.rm=TRUE)
  df$BsmtHalfBath = as.numeric( as.character( df$BsmtHalfBath ) )
  df$BsmtHalfBath[is.na(df$BsmtHalfBath)] = mean(df$BsmtHalfBath, na.rm=TRUE)
  df$GarageCars = as.numeric( as.character( df$GarageCars ) )
  df$GarageCars[is.na(df$GarageCars)] = mean(df$GarageCars, na.rm=TRUE)
  df$GarageArea = as.numeric( as.character( df$GarageArea ) )
  df$GarageArea[is.na(df$GarageArea)] = mean(df$GarageArea, na.rm=TRUE)
  
  # Fix missing values
  df$Utilities.n[is.na(df$Utilities.n)] = 3       # Modal value
  df$Functional.n[is.na(df$Functional.n)] = 7     # Modal value
  df$KitchenQual.n[is.na(df$KitchenQual.n)] = 3   # Modal value
  df$Electrical[df$Electrical=="NA"] = as.factor("SBrkr")

  # Take logarithms where appropriate
  df$X1stFlrSF = log( df$X1stFlrSF )
  names(df)[names(df)=="X1stFlrSF"] <- "Ln1stFlrSF"
  df$GrLivArea = log( df$GrLivArea )
  names(df)[names(df)=="GrLivArea"] <- "LnLivArea"
  df$OFHEO = log( df$OFHEO )
  names(df)[names(df)=="OFHEO"] <- "LnOFHEO"

  # Normalize dependant variable (if this is a training set)  
  if (!is.null(df$SalePrice)) {
    df$SalePrice = log( df$SalePrice ) - df$LnOFHEO 
    names(df)[names(df)=="SalePrice"] <- "RelPrice"
  }

  df
}

```
OK, Wait!  We interrupt this function to bring you a special announcement!
That last section does something non-standard.  It normalizes the target
variable by dividing it by a house price index.  The new target variable
is called "RelPrice" and will have to be un-normalized at the end.  More on this issue
when I get to "Analysis" section of this notebook.

Now on to the fancy stuff.  The following function takes a baseline OLS model,
adds a factor to it, and returns the factor's OLS coefficients, to be used later
to recode the factor as a continuous variable.  (Note that the factor enters into
the OLS as an exhaustive set of dummies -- one for every level present -- because
the constant term is removed.)
```{r}
# Function to get coefficients to be used to make factor continuous given baseline model
getCoeffs <- function( df, basemodel, factor ) {
  mod <- paste0( basemodel, "+", factor, "-1" )
  lm <- lm(formula=mod, data=df)
  fnames <- grep( factor, names(lm$coefficients), fixed=TRUE )
  lm$coefficients[fnames]
}

```

...and the next function takes those dummy coefficients and uses them to make
the factor into a continuous variable:
```{r}
# Function to make factor continuous (given dummy coefficients) and add "_r" to name
makeContinuous <- function( df, factor, coeffs ) {
  outvar <- 0*(1:nrow(df))
  fact <- df[[factor]]
  for ( n in levels(fact) ) {
     outvar[fact==n] <- coeffs[paste0(factor,n)]
  }   
  df <- cbind( df, outvar )  
  names(df)[ncol(df)] <- paste0( factor, "_r" )
  df[factor] <- NULL
  df
}

```

That's how I process many-valued unordered factors.  As for the ordered factors,
I started out by making them into integers, but the following set of functions
makes them more genuinely continuous, using a procedure somewhat similar to
what I did with the unordered factors.  These next functions also apply to variables
that were originally presented as integers (e.g. 10-point scale) but don't
represent phenomena where the ratios of different values are meaningful.

As the first step, I have a function that takes an integer variable and creates
an above/below dummy for each of its values.
```{r}
# Function to convert integer variable to above/below dummies for each of its values
dummify <- function( df, var ) { 
  v <- df[[var]]
  vals <- sort( unique(v) )
  n <- length(vals) - 1
  if (n==0) { return (NULL) }
  for (i in vals[1:n]) {
    newname <- paste0(var,".gt",as.character(i))
    newvar <- as.numeric(v > i)
    df <- cbind( df, newvar )
    names(df)[ncol(df)] <- newname
  }
  df
}

```

As the second step, I have a function that adds those dummies to the model.
But it then removes any dummies whose coefficients fit with the wrong sign.
```{r}
# Function to add to model any of the created dummies that have positive coefficients
addDummiesToModel <- function( df, var, mod ) {
  vars <- c()
  for (name in names(df)) {
    if ( grepl(paste0(var,".gt"), name, fixed=TRUE) ) {
      vars <- c(vars,name)
    }
  }
  newmod <- paste0( mod, " + ", paste(vars, collapse=" + ") )
  lmfull <- lm( formula=newmod, data=df )
  count = length(vars)
  for (var in vars) {
    if ( is.na( lmfull$coefficients[var] ) ) {
      print( paste0( "Variable ", var, " removed due to rank deficiency"))
      newmod <- gsub( paste0(" + ",var), "", newmod, fixed=TRUE )
      count = count - 1
    }
    else if ( lmfull$coefficients[var] < 0 ) {
      newmod <- gsub( paste0(" + ",var), "", newmod, fixed=TRUE )
      count = count - 1
    }
  }
  if (count==0) { return(NULL) }
  newmod
}

```

The next function runs the resulting model and extracts the dummy coefficients.
```{r}
# Function to get coefficients to be used to make ordered variable continuous
#   (given a model that includes any dummies found to have positive coefficients)
getOrderedCoeffs <- function( df, newmod, var ) {
  lm1 <- lm(formula=newmod, data=df)
  fnames <- grep( var, names(lm1$coefficients), fixed=TRUE )
  coeffs <- lm1$coefficients[fnames]
  names(coeffs) <- gsub( paste0(var,".gt"), "", names(coeffs), fixed=TRUE )
  coeffs
}

```

But it would get verbose having to run the last three functions separately, so I have
another function that runs them all and checks the results.
```{r}
# Function to make coefficients to be used to make ordered variable continuous
makeOrderedCoeffs <- function( df, mod, var ) {
  df <- dummify( df, var )
  if (is.null(df)) { 
    print( paste0("dummify returned NULL for ", var) )
    return(NULL) 
  }
  mod <- addDummiesToModel( df, var, mod )
  if (is.null(mod)) {
    print("addDummiesToModel returned NULL")
    return(NULL) 
  }
  coeffs <- getOrderedCoeffs( df, mod, var )
  if (is.null(coeffs)) {
    print("getOrderedCoeffs returned NULL")
    return(NULL) 
  }
  coeffs
}

```

And finally, the next function uses the resulting coefficients to make
the integer variable into a continuous one.  (This often has to be done separately
so that coefficients generated from training data can be applied to test data.)
```{r}
# Function to make ordered variable continuous (given coefficients) and add "_r" to name
makeOrderedContinuous <- function( df, var, coeffs ) {
  outvar <- 0*(1:nrow(df))
  v <- df[[var]]
  outvar <- 0
  for (n in names(coeffs)) {
    outvar <- outvar + ifelse( v>as.numeric(n), coeffs[n], 0 )
  }
  df <- cbind( df, outvar )  
  names(df)[ncol(df)] <- paste0( var, "_r" )
  df[var] <- NULL
  df
}

```

But even with the functions I've already written, the whole process can get verbose,
and it may have to be done several times, so I have the following function
that goes through a bunch of ordered variables and makes them all continuous.
Note that, for each variable, it uses as a base model one that includes all of
the previously processed variables but not the ones yet to be processed (with one
exception -- well, two, but the second one would have been included either way).
I could have included the yet-to-be-processed variables as pseudocontinuous
integers, but I chose not to.  Or I could have just used the original baseline
model for each variable.  (This issue also came up with the many-valued factors.)
It's not clear that what I chose was the best choice.  But it's also not clear
that it wasn't the best choice.
```{r}
# Function to make a bunch of ordered variables continuous using functions above
makeOrderedVariablesContinuous <- function(df, mod, varlist) {

  # Note: Ordered variables should be defined so that positive coefficients are expected

  # Start with a baseline model.  For each ordered variable (in sequence, so it matters
  #   how the variables are arranged in the variable list), create above/below dummies
  #   for each level.  Delete any dummies with wrong sign.  Use OLS coefficients on
  #   remaining dummies to define a continuous version of the ordered variable. 
  
  # "OverallQual" and "OverallCond" need special treatment because they are already
  #    in the baseline model so must be removed before adding corresponding dummies.
  already_in_baseline_model <- c("OverallQual", "OverallCond")
  
  orderedCoeffs <- list() # List that will contain coefficients for ordered variables
  varsToDrop <- c() # List that will contain variable names dropped because all wrong sign
  i <- 0 
  for ( var in varlist ) {
    if ( var %in% already_in_baseline_model ) {
      mod <- gsub( paste0("+ ", var), "", mod, fixed=TRUE ) # Delete from model
    } 
    co <- makeOrderedCoeffs( df, mod, var )
    if ( is.null(co) ) {
      varsToDrop <- c(varsToDrop, var)
      df[var] <- NULL
    }
    else {
      df <- makeOrderedContinuous( df, var, co )
      mod <- paste0( mod, " + ", var, "_r")
      i <- i + 1
      orderedCoeffs[[i]] <- co
      names(orderedCoeffs)[[i]] <- var
    }
  }
  output <- list( df, mod, orderedCoeffs, varsToDrop )
  names(output) <- c("df", "mod", "coeffs", "drop")
  output
}

```

OK, I'm finished with the "feature generation" functions, but once they've run
there will be more data cleaning to do (which, again, one could probably do
more elegantly and robustly than I did).
```{r}
# Function to do final data cleaning after variables have been processed into features
finalCleaning <- function( df ) {

  # Fix numeric variables with missing values in test and/or validation set
  df$MSSubClass_r[is.na(df$MSSubClass_r)] = mean(df$MSSubClass_r, na.rm=TRUE)
  df$Exterior1st_r[is.na(df$Exterior1st_r)] = mean(df$Exterior1st_r, na.rm=TRUE)
  df$Exterior2nd_r[is.na(df$Exterior2nd_r)] = mean(df$Exterior2nd_r, na.rm=TRUE)
  df$Condition2_r[is.na(df$Condition2_r)] = mean(df$Condition2_r, na.rm=TRUE)
  
  # Collapse sale condition categories
  salecon <- as.character(df$SaleCondition)
  df$SaleMisc <- ifelse( salecon=="Family" | salecon=="Partial", 1, 0 )
  df$SaleAbnormal <- ifelse( salecon=="Abnorml", 1, 0 )
  df$SaleCondition <- NULL

  # Collapse sale type categories
  st <- as.character(df$SaleType)
  con <- c("Con", "ConLw", "ConLI", "ConLD")
  wd <- c("WD", "CWD", "VWD")
  df$Contract <- ifelse( st %in% con, 1, 0 )
  df$WrntyDeed <- ifelse( st %in% wd, 1, 0 )
  df$NewSale <- ifelse( st=="New", 1, 0 )
  df$SaleType <- NULL

  # Only one kind of building type seems to be different
  df$SingleFam <- ifelse( as.character(df$BldgType)=="1Fam", 1, 0 )
  df$BldgType <- NULL

  # It matters if you have a garage, but this is captured by "HasGarageYear"
  # It also matters if it's a real garage or just a car port, so:
  df$CarPort <- ifelse( as.character(df$GarageType)=="CarPort", 1, 0 )
  df$GarageType <- NULL

  # Residential vs. nonresidential seems to be only relevant aspect of zoning
  zo <- as.character(df$MSZoning)
  res_zone <- c( "FV", "RH", "RL", "RP", "RM" )
  df$Residential <- ifelse( zo %in% res_zone, 1, 0 )
  df$MSZoning <- NULL

  # Get rid of RoofMatl. It is an overfit dummy for one case.
  # Earlier analysis showed all levels got OLS coefficients that were
  # very significantly different from zero but not different from one another.
  # "ClyTile" was the omitted category and was only one case.
  df$RoofMatl <- NULL

  # Get rid of MiscFeature. Per earlier analysis, it's a mess. Don't want to deal with it.
  df$MiscFeature <- NULL

  # Factors that earlier analyses didn't like and too much of a pain in the neck to keep
  df$Fence <- NULL
  df$RoofStyle <- NULL
  df$Heating <- NULL

  # I didn't see any residual seasonal pattern, so:
  df$MoSold <- NULL

  # These nonlinearitiesn seem to matter
  df$LotFrontage2 <- df$LotFrontage^2
  df$SinceRemod <- df$YrSold - df$YearRemodAdd
  df$SinceRemod2 <- df$SinceRemod^2
  df$YrSold <- NULL
  df$YearRemodAdd <- NULL
  df$BsmtFinSF1sq <- df$BsmtFinSF1^2

  # The following turn out to be redundant. But may want to bring them back later.
  df$TotalBsmtSF <- NULL
  df$HasMasVnr <- NULL
  df$KitchenAbvGr_r <- NULL
  df$GarageCond.n_r <- NULL
  
  df
}

```

Phew!  So much for data processing functions.  Now on to the model fitting.

The following function fits the models using the "train" function from the "caret" package.
```{r}
# Function to fit formula "fo" using several fitting methods "models"
fitModels <- function( df, fo, models, runParallel, seed ) {

  if (runParallel) {
    # Set up for multiple cores
    cluster <- makeCluster(detectCores() - 1) # convention to leave 1 core for OS
    registerDoParallel(cluster)
    }

  # Fit the models
  modelfits = list()
  for (m in models) {
    print ( paste("Training model:", m) )
    set.seed(seed)
    fit <- train( as.formula(fo), data=df, method=m )
    modelfits = c(modelfits, list(fit))
  }

  if (runParallel) {
    # Go back to sequential processing
    stopCluster(cluster)
    registerDoSEQ()
  }
  
  names(modelfits) <- models

modelfits 
}

```

And the following function uses those fits to make predictions.
```{r}
# Function to make predictions given several model fits
makePredictions <- function( df, modelfits, basepred, is_test ) {

  modelnames <- names(modelfits)
  rmses <- list()
  predicted <- list()
  for (fi in modelfits) {
    writeLines ( paste("\n\n\nPredicting for model:", fi[[1]]) )
    p <- predict(fi, newdata=df)
    if (!is.null(basepred)) {
      p[is.na(p)] <- basepred[is.na(p)]
      }
    predicted <- c(predicted, list(p))
    if (!is_test) {
      rmses <- c(rmses, rmse(df$RelPrice, p))
    }
  }
  names(predicted) <- modelnames
  if (!is_test) {
    names(rmses) <- modelnames
    print( rmses )
  }
  
  predicted
}

```

Finally, the following function uses the predictions (presumably from a cross-validation
set predicted with a model fitted on a training set) to create an ensemble of fitting
procedures.  It starts out by including all the procedures that have been used
to fit the models.  To generate the first candidate for a new set of predictions,
it runs a no-constant OLS with coefficients constrained to sum to one.
Then, one by one, it eliminates procedures with negative coefficients, until
it ends up with an OLS where all coefficients are positive.  That OLS then generates
the weights for the ensemble.  (I guess this is what machine learning people call
"stacking" an OLS on top of the other models.)
```{r}
# Function to choose an ensemble of models and weight them
chooseEnsemble <- function( df, predicted ) {

  preddf <- cbind( as.data.frame(predicted), df$RelPrice )
  colnames(preddf) <- c(modelnames, "actual")
  bestcoef = 0
  i = 1

  while ( bestcoef <= 0 ) {

    # Run full equation
    predeq <- lm(actual~.-1, data=preddf)
    summary( predeq )
    cof <- predeq$coefficients
    bestmod <- names( cof[order(cof,decreasing=TRUE)] )[i]
    i = i + 1

    # Force coefficients to sum to 1
    preddf2 <- preddf
    predmod <- "actual~-1"
    for (n in names(preddf2)) {
      preddf2[n] <- preddf2[n] - preddf[bestmod]
      if (!(n=="actual"||n==bestmod)) {
        predmod = paste0(predmod, "+", n)
      }
    }

    # Keep dropping variables until all coefficients are positive
    eq <- lm( predmod, data=preddf2 )
    while( min(eq$coefficients) < 0 ) {
      dropv <- names( which.min(eq$coefficients) )
      predmod <- gsub( paste0("+",dropv), "", predmod, fixed=TRUE)
      eq <- lm(predmod, data=preddf2)
    }

    # Calculate missing coefficient
    bestcoef = 1 - sum(eq$coefficients)
    names(bestcoef) <- bestmod
    weights <- c(eq$coefficients, bestcoef)
  }
  
  weights
}

```

We're almost done with functions.  The following is a function that I copied from
mtyxwp's "svm_simple" kernel.  It deals with missing values much more elegantly
than I do in my own cleaning functions.  (Since I'd already written and tested them,
I didn't want to risk introducing new bugs by replacing parts of them with this,
but maybe in a future version....)
This function is called in the last section of this notebook, where it starts
over from scratch to fit an SVM and then averages its result with the linear one.
```{r}
# mtyxwp's function to deal with missing data (copied from "svm_simple" kernel
deal_missing <- function(simpledf){
  for(i in 1:ncol(simpledf))
  {
    u <- simpledf[,i]
    if (is.numeric(u))
    {
      simpledf[is.na(u),i] <- median(simpledf[!is.na(u),i])
    } else
    {
      u[is.na(u)] <- "Not Available"
      simpledf[,i] <- as.factor(u)
    }
  }
  return(simpledf)
}

```

OK, so much for functions.  Now on to data.

## Data

For my analysis I use a price index taken from the Federal Housing Finance Agency
(which used to be called the Office of Federal Housing Enterprise Oversight,
and I use the old initials, mostly because I didn't know the name had changed
when I first wrote this code, but I also just like the old acronym better).
The data come from [FHFA's website](https://www.fhfa.gov/DataTools/Downloads/pages/house-price-index.aspx),
but I just cut and pasted the relevant portion into my program.  It refers
to the West North Central House Price Index, monthly from 2006 through 2010.
(Ames, Iowa, where the competition data are from, is in that region, and the
sales span that range of time.)
```{r}
OFHEO = c(209.32, 210.23, 211.68, 212.71, 214.37, 215.37,
          216.37, 216.22, 215.45, 214.48, 214.73, 211.92,
          212.23, 214.99, 215.82, 216.99, 217.89, 218.28,
          218.69, 216.78, 217.27, 212.78, 212.72, 211.6, 
          208.58, 208.62, 209.68, 210.28, 209.78, 210.87,
          209.68, 208.77, 206.08, 206.07, 200.51, 201.47,
          201.78, 204.24, 201.05, 203.8,  205.1,  206.55,
          205.27, 204.63, 203.47, 204.22, 202.74, 199.78,
          196.35, 197.64, 198.89, 202.13, 204.25, 204.61,
          200.13, 201.76, 198.03, 197.87, 195.11, 193.46 )

Year = c( 2006, 2006, 2006, 2006, 2006, 2006,
          2006, 2006, 2006, 2006, 2006, 2006,
          2007, 2007, 2007, 2007, 2007, 2007,
          2007, 2007, 2007, 2007, 2007, 2007,
          2008, 2008, 2008, 2008, 2008, 2008,
          2008, 2008, 2008, 2008, 2008, 2008,
          2009, 2009, 2009, 2009, 2009, 2009,
          2009, 2009, 2009, 2009, 2009, 2009,
          2010, 2010, 2010, 2010, 2010, 2010,
          2010, 2010, 2010, 2010, 2010, 2010 )

Month = c( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
           1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
           1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
           1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
           1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 )
ofheo <- data.frame( Month, Year, OFHEO )

```

Then I merge that data with the competition data
```{r}
rawdata <- read.csv(trainfile, na.strings="")
data1 <- merge(rawdata, ofheo, by.x=c("YrSold","MoSold"), by.y=c("Year","Month"))
data1 <- cleanData( data1 )

```
...and I got me a training set

My approach was to divide these data into 3 subsets:
- "train1"   (60%)  for  primary training
- "validate" (20%)  for  cross-validation
- "testing"  (20%)  for  initial testing
```{r}
set.seed(999)
inTrain <- createDataPartition(y=data1$RelPrice, p=0.8, list=FALSE)
training <- data1[inTrain,]
testing <- data1[-inTrain,]
inTrain1 <- createDataPartition(y=training$RelPrice, p=0.75, list=FALSE)
train1 <- training[inTrain1,]
validate <- training[-inTrain1,]

```

## Initial Analysis and Feature Processing of Primary Training Data

I start by making a working copy of my initial training set.  (In retrospect,
I was probably being too careful, but for now, this is the code.)
```{r}
da <- train1

```

I'm going to come clean here.  I'm convinced that house prices (even in a place
like Ames, Iowa that's far from the craziness in New York and Tokyo) depend heavily
on broad market conditions.  For Heaven's sake, did we not all live through 2008?
I don't care what the training data say:  if you tried to sell a house in 2009,
even in Iowa, you didn't sell it for the same price you could have gotten in 2006.
In some ad hoc exploratory analysis I ran before writing this code, the data didn't
seem to like the OFHEO (FHFA, whatever) price index variable very much.
But I'm damned if I'm going to let some officious regularization procedure 
tell me it doesn't belong in the model.  Until you can show me test data where
the variable makes the predictions less accurate, it stays in the damn model!

Nonetheless, I'm going to employ regularization procedures, because that's how
you get good predictions.  So I'm forcing them to accept my OFHEO variable.  How?
But putting it on the left-hand side, in the denominator.  Of course it will also
be on the right-hand side, but it will need a coefficient on -1 to be effectively
removed from the model.  (I'm using logarithms, so it effectively has a coefficient
of -1 on the left-hand side, which would wash out if it also has one on the right-hand
side.)  The regularization procedures won't be able to take it out of the model,
because they only affect the right-hand side.  If they take it out of the right-hand
side, it will be reflected one-for-one in the predictions.

So here's my baseline model
```{r}
basemod <- "RelPrice ~ LnOFHEO + Ln1stFlrSF + LnLivArea + OverallQual + OverallCond"
```

Very simple, but it captures the most important stuff.  If I were in a big hurry,
I'd just use this model and be done with it.  No regularization or crazy-ass support
vector machines.  Just good, old-fashioned ordinary least squares with a small and
obvious set of variables.  It won't win any competitions, but it will generate pretty
good predictions with little effort.

The following wasn't in my original code, but just to make a point, I'm going to run this with my normalization taken out:
```{r}
simple <- lm( RelPrice+LnOFHEO ~ Ln1stFlrSF + LnLivArea + OverallQual + OverallCond, da)
summary( simple )

```
It explains well over 3/4 of the variance.  With 4 variables!

Just for the hell of it:
```{r}
rmse( validate$RelPrice+validate$LnOFHEO, predict(simple, validate) )
```
Behold the power of OLS!

We can fight like feral dogs over the next 4 percentage points of RMSE,
but FWIW, these 4 variables do the heavy lifting.  The difference between simple OLS
and some of the best models in this competition might not be enough to pay
the realtor's commission on one of the houses we're predicting.

That said, some realtors make a damn good living on those few percentage points, 
so let's continue.

Process the many-valued factors.  (Don't worry about warnings that there are
missing levels.)
```{r}
# Make factors continuous and add continuous versions to baseline model one by one
factors <- c( "Neighborhood", "MSSubClass", "Condition1", "Exterior1st", "Condition2", 
              "Exterior2nd",  "LotConfig",  "Foundation")
mod <- basemod
coeffs <- list()
i <- 0 
for (f in factors) {    
  co <- getCoeffs( da, mod, f )
  i <- i + 1
  coeffs[[i]] <- co
  names(coeffs)[i] <- f
  da <- makeContinuous( da, f, co )
  mod <- paste0( mod, "+", f, "_r" ) 
}

```

We're up to 87% of the variance now (but the adjusted R-squared is not meaningful,
since my "make factors continuous" procedure stole a bunch of degrees of freedom
and hid them under the bed)
```{r}
# Show output of augmented model
print( summary( lm( formula=mod, data=da ) ) )
```

Now we make the ordered variables continuous
```{r}
ordered <- c("OverallQual",    "OverallCond",  # 1st 2=special cases bcuz in baseline mod  
             "Functional.n", "Fireplaces",     "KitchenQual.n", 
             "BsmtExposure.n", "HeatingQC.n",  "Utilities.n",    "FullBath",  
             "HalfBath",       "GarageCars",   "BsmtFullBath",   "GarageQual.n", 
             "BsmtFinType1.n", "PavedDrive.n", "BsmtCond.n",     "GarageCond.n", 
             "FireplaceQu.n",  "ExterQual.n",  "TotRmsAbvGrd",   "LotShape.n", 
             "BsmtHalfBath",   "PoolQC.n",     "BsmtFinType2.n", "ExterCond.n", 
             "BedroomAbvGr",   "BsmtHeight",   "KitchenAbvGr",   "GarageFinish.n") 

out <- makeOrderedVariablesContinuous(da, mod, ordered)
da <- out$df
mod <- out$mod
orderedCoeffs <- out$coeffs
varsToDrop <- out$drop

```
Clean the data
```{r}
da = finalCleaning( da )

```
Explore the result by looking at an OLS with everything but the kitchen sink.  (Oh,
wait, maybe there is a variable about the kitchen sink.  I'll have to check the
data description again.)
```{r}
regmodel = lm( formula="RelPrice ~ .", data=da )
print( summary( regmodel ) )

```

## Quick and Dirty Look at Validation Set

Make a working copy.
```{r}
da2 <- validate

```

Apply the transformations that we got from the training set.
```{r}
for (f in factors) {    
  co <- coeffs[[f]]
  da2 <- makeContinuous( da2, f, co )
}
for ( var in ordered ) {
  if ( var %in% varsToDrop ) {
    da2[var] <- NULL
  }
  else {
    da2 <- makeOrderedContinuous( da2, var, orderedCoeffs[[var]] )
  }
}
da2 <- finalCleaning( da2 )

```

Predict using the "Kitchen Sink" OLS Model
```{r}
# Make predictions
prediction <- predict(regmodel, da2, type="response")

# Fill in missing values
baselm <- lm(formula=basemod, data=train1)
basepred <- predict( baselm, validate, type="response")
prediction[is.na(prediction)] <- basepred[is.na(prediction)]

# RMSE
rmse(da2$RelPrice,prediction)

```
That's implausibly good.  Which means CARET's random assignment has given us
an unusually easy-to-predict
validation set.  (It also casts doubt on my earlier ranting about realtors and such,
but I think my basic point still stands.)

Also why I should use k-fold cross-validation in the future, but I'm going to continue
with this one for now, since the code is already written and tested.

Inspecting the earlier OLS output, I threw out the variables that looked like
they weren't helping, and this is what I came up with:
```{r}
fo = "RelPrice ~ LotFrontage + LotArea + Alley + LandContour + LandSlope "
fo = paste0(fo, "+ YearBuilt + MasVnrArea + BsmtFinSF1 + BsmtFinSF2 ")
fo = paste0(fo, "+ BsmtUnfSF + CentralAir + Ln1stFlrSF + X2ndFlrSF ")
fo = paste0(fo, "+ LowQualFinSF + LnLivArea + GarageArea + WoodDeckSF ")
fo = paste0(fo, "+ OpenPorchSF + EnclosedPorch + X3SsnPorch + ScreenPorch ")
fo = paste0(fo, "+ LnOFHEO + HasLotFrontage + Neighborhood_r + Condition1_r ")
fo = paste0(fo, "+ Exterior1st_r + Condition2_r + Exterior2nd_r + LotConfig_r ")
fo = paste0(fo, "+ Foundation_r + OverallQual_r + OverallCond_r + Functional.n_r ")
fo = paste0(fo, "+ Fireplaces_r + KitchenQual.n_r + BsmtExposure.n_r + HeatingQC.n_r ")
fo = paste0(fo, "+ FullBath_r + HalfBath_r + GarageCars_r + BsmtFullBath_r ")
fo = paste0(fo, "+ GarageQual.n_r + BsmtFinType1.n_r + ExterQual.n_r + TotRmsAbvGrd_r ")
fo = paste0(fo, "+ PoolQC.n_r + BsmtFinType2.n_r + ExterCond.n_r + SaleAbnormal ")
fo = paste0(fo, "+ Contract + NewSale + SingleFam + Residential + LotFrontage2 ")
fo = paste0(fo, "+ SinceRemod + SinceRemod2 + BsmtFinSF1sq ")

mymodel = lm( formula=fo, data=da )
prediction <- predict(mymodel, da2, type="response")
prediction[is.na(prediction)] <- basepred[is.na(prediction)]
rmse(da2$RelPrice,prediction)

```
Not a big improvement, but, whatev.


# Choosing an Ensemble of Models

Fit the models.  (Go have lunch while this section runs, if you're using Kaggle or
have a slow computer.)
```{r}
modelfits <- fitModels( da, fo, modelnames, !running_as_kaggle_kernel, 998 )

```

Create the ensemble.
```{r}
predicted <- makePredictions ( da2, modelfits, basepred, FALSE )
weights <- chooseEnsemble( da2, predicted )
bestmodels <- names(weights)

```

Take a look at the ensemble results for the validation set.  (Note that these are
in-sample predictions from the point of view of the ensemble-maker, so take with
a grain of salt.)
```{r}
preddf <- cbind( as.data.frame(predicted), da2$RelPrice )
colnames(preddf) <- c(modelnames, "actual")
p_ensemble <- as.data.frame( as.matrix(preddf[bestmodels]) %*% weights )
names( p_ensemble ) <- "ensemble"
rmse( da2$RelPrice, p_ensemble )

```

# Initial Testing

Redo the feature creation process on data that include the validation set
```{r}
# Make working copy of data
da3 <- training

# Make factors continuous
mod <- basemod
coeffs <- list()
i <- 0 
for (f in factors) {    
  co <- getCoeffs( da3, mod, f )
  i <- i + 1
  coeffs[[i]] <- co
  names(coeffs)[i] <- f
  da3 <- makeContinuous( da3, f, co )
  mod <- paste0( mod, "+", f, "_r" ) 
}

# Make ordered variables continuous
out <- makeOrderedVariablesContinuous(da3, mod, ordered)
da3 <- out$df
mod <- out$mod
orderedCoeffs <- out$coeffs
varsToDrop <- out$drop

# Make data nice
da3 = finalCleaning( da3 )

```

Process the initial test data accordingly
```{r}
# Make working copy of data
da4 <- testing

# Make factors continuous and add continuous versions to baseline model one by one
for (f in factors) {    
  co <- coeffs[[f]]
  da4 <- makeContinuous( da4, f, co )
}

# Make ordered variables continuous
for ( var in ordered ) {
  if ( var %in% varsToDrop ) {
    da4[var] <- NULL
  }
  else {
    da4 <- makeOrderedContinuous( da4, var, orderedCoeffs[[var]] )
  }
}

# Make data nice
da4 <- finalCleaning( da4 )

```

Fit models.  (Go have a snack.)
```{r}
modelfits <- fitModels( da3, fo, bestmodels, !running_as_kaggle_kernel, 997 )
```

Make individual model predictions.
```{r}

predicted <- makePredictions ( da4, modelfits, NULL, FALSE )

```

Combine results into an ensemble prediction, and check performance.
```{r}

p_ensemble <- 0*da4$RelPrice
preddf <- as.data.frame( predicted )
p_ensemble <- as.data.frame( as.matrix(preddf[bestmodels]) %*% weights )
names(p_ensemble) <- "ensemble"

# Estimated forecast error
print( rmse( da4$RelPrice, p_ensemble ) )

```

Not as good as I would have hoped, but maybe CARET gave us a hard-to-predict
initial test set, just as it gave us an easy-to-predict validation set.
Again, the future is k-fold.

# Training the Final Model

Feature creation again
```{r}
# Make working copy of data
da5 <- data1

# Make factors continuous
mod <- basemod
coeffs <- list()
i <- 0 
for (f in factors) {    
  co <- getCoeffs( da5, mod, f )
  i <- i + 1
  coeffs[[i]] <- co
  names(coeffs)[i] <- f
  da5 <- makeContinuous( da5, f, co )
  mod <- paste0( mod, "+", f, "_r" ) 
}

# Make ordered varibles continuous
out <- makeOrderedVariablesContinuous(da5, mod, ordered)
da5 <- out$df
mod <- out$mod
orderedCoeffs <- out$coeffs
varsToDrop <- out$drop

# Make data nice
da5 = finalCleaning( da5 )

```

Let's take a look to the kitchen sink regression for the full training set
```{r}
regmodel = lm( formula="RelPrice ~ .", data=da5 )
print( summary( regmodel ) )

```

Based on these results, I'm going to make some changes to the feature list
```{r}
fo <- gsub( "+ ExterQual.n_r ", "", fo, fixed=TRUE  )
fo <- gsub( "+ ExterCond.n_r ", "", fo, fixed=TRUE )
fo <- gsub( "+ PoolQC.n_r ", "", fo, fixed=TRUE )
fo <- paste0(fo, "+ FireplaceQu.n_r + BsmtHeight_r + SaleMisc ")

regmodel = lm( formula=fo, data=da5 )
print( summary( regmodel ) )

```

Here comes the final fit.  (Time for dinner yet?)
```{r}
modelfits <- fitModels( da5, fo, bestmodels, !running_as_kaggle_kernel, 996 )

```

# Predicting the Test Data with Linear Model Ensemble

Read in and process the test data.  (This looks easy, but it was hard
the first time I did it, because there was considerable untidiness in the test data
that I foolishly failed to anticipate.)
```{r}
# Read data
testdat <- read.csv(testfile, na.strings="")
data2 <- merge(testdat, ofheo, by.x=c("YrSold","MoSold"), by.y=c("Year","Month"))
data2 <- cleanData( data2 )

# Make factors continuous and add continuous versions to baseline model one by one
for (f in factors) {    
  co <- coeffs[[f]]
  data2 <- makeContinuous( data2, f, co )
}

# Make ordered variables continuous
for ( var in ordered ) {
  if ( var %in% varsToDrop ) {
    data2[var] <- NULL
  }
  else {
    data2 <- makeOrderedContinuous( data2, var, orderedCoeffs[[var]] )
  }
}

data2 <- finalCleaning( data2 )

```

Generate predictions
```{r}
predicted <- makePredictions ( data2, modelfits, NULL, TRUE )

preddf <- as.data.frame( predicted )
p_ensemble <- as.data.frame( as.matrix(preddf[bestmodels]) %*% weights )
names(p_ensemble) <- "ensemble"
```

Transform predictions back into the form that Kaggle wants them
```{r}
prediction <- p_ensemble + data2$LnOFHEO

result <- data.frame( cbind( data2$Id, exp(prediction) ) )
names(result) <- c("Id", "SalePrice")
sorted_result <- result[order(result$Id),]

```

And one could just store the results, but I'm going to buy another 40 basis points
or so of RMSE improvement by averaging these predictions with those of a completely
different, but comparably performing, model.

# Training Simple SVM and Predicting Test Data

This code is from mtyxwp's "svm_simple" kernel.  mtyxwp's approach to the data is
completely different from mine, and the model is completely different, and
diversity is valuable when it comes to making predictions.  Since I didn't write
this code, I'm not sure I can throw much light on it, beyond what is self-explanatory.
Here's the data processing, which makes no attempt to create new features beyond
what is necessary for tidiness.
```{r}
train <- read.csv(trainfile, stringsAsFactors = F)
test <- read.csv(testfile, stringsAsFactors = F)
price = train[,c("SalePrice")]

simpledf <- rbind(train[,-81], test)

simpletrain <- deal_missing(simpledf)
sdf = simpletrain

sdftrain = sdf[1:1460,]
sdftrain = cbind(sdftrain,price)
sdftest = sdf[1461:2919,]
sdftrain = sdftrain[,-1]
id = sdftest[,1]
sdftest = sdftest[,-1]

straincpy = sdftrain
straincpy[sapply(straincpy, is.factor)] <- lapply(straincpy[sapply(straincpy, is.factor)], as.numeric)
stestcpy = sdftest
stestcpy[sapply(stestcpy, is.factor)] <- lapply(stestcpy[sapply(stestcpy, is.factor)], as.numeric)

```

And then we fit probably the simplest possible support vector machine, from
the user's point of view.  (A lot is going on under the hood.)
```{r}
model.svm <- svm(price ~ ., data = sdftrain, cost = 1)
```

And we make predictions
```{r}
price.svm = predict(model.svm, sdftest)

svmResult = cbind(Id = id, SalePrice = price.svm)
colnames(svmResult) = c("Id","SalePrice")

```


# Avearging Linear Results with SVM Results

One could try to do some sophisticated stacking, and maybe I will in the future,
but for now, just take the (logarithmic) average.
```{r}
final_answer <- exp (   ( log(sorted_result$SalePrice) + log(svmResult[,"SalePrice"]) ) / 2   )
final_result <- data.frame( cbind( sorted_result$Id, final_answer ) )
names(final_result) <- c("Id", "SalePrice")
write.csv(final_result, file="kaggleSubmission5a.csv", row.names=FALSE)
```

Ready to submit.