```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
setwd("~/Github/bigdata/lectures/l13/")
library(tidyverse)
```

## Today's Class {.smaller}

1. Review
    - Ensembles
    - Boosting
2. Data Cleaning Basics
3. Dates
4. Wide vs long data
5. Merges/Joins
6. `NA`s, `NaN`, etc. 
7. Predictions 2
    - CI vs PI
    - Sundays
    - My forecast
    - Ensemble PIs

# Review

```{r,echo=F, message=F}
n = 500
set.seed(101)
x1 = runif(n)
x2 = runif(n)
prob  = ifelse(x1 < 0.5 & x1 > 0.25 & x2 > 0.25 & x2<0.75,0.8,0.3)
y  = as.factor(rbinom(n,1,prob))
df = data.frame(y=y,x1=x1,x2=x2)
```

```{r,echo=F}
sequence = seq(0,1,length.out=201)
grid = expand.grid(sequence,sequence)
colnames(grid) = c("x1","x2")
```

## Weighted -- Example

Suppose I have two models: 
- $M_1$ always predicts $\hat y_1  = M_1(X) = \bar y$
- $M_2$ uses $p$ variables and linear regression to predict $\hat y_2 = M_2(X) = \hat \beta_0 + \hat \beta_1 x_1 +...+ \hat \beta_p x_p$

These are not equivalently good. $M_2$ is _likely_ better (not always -- see overfitting again). 

We could plug them both into our naive average, and get a prediction.


## Inverse Variance Weights

The basic idea is that a single measure of the prediction errors of each model can help us generate the best weights.

The models with the smallest average errors should get the most weight, and the models with the largest average errors should get the least weight. 

If you are minimizing MSE, this will look like weights which are proportional to the ``precision" (aka the inverse of the variance). 

## Procedure

1. Estimate $K$ models $M_1,...,M_K$
2. Use cross-validation (k-fold) to estimate OOS errors for each model.
3. Calculate average OOS error $MSE_k$ for each model.
4. Find weights $\textbf{w} \propto MSE_k$ such that $\sum w_k = 1$
5. Make your predictions using those weights

This will be reasonably quick. For $K$ models and $m$ folds in the cross validation, you'll only need to fit $Km$ models. The weights pop out analytically.

# Boosting

## Boosting -- Ensembles 2

Boosting changes things up slightly. It says "we want to do a better job of predicting when we are wrong". 

Thus far, we have been giving each tree slightly different data, and averaging across each tree. 

But each tree still was optimizing the same thing -- mean prediction error\footnote{This could be MSE, etc depending on the model.}
:

$$ \frac1n \sum_{i=1}^n l(\hat{y},y) =  \frac1n \sum_{i=1}^n 1(\hat{y} \neq y)  $$

## Boosting Insight 1

Some observations are harder to predict. 

We want to build a model that predicts everything well. 

We can give _more weight_ to observations that are hard to predict. 

$$ \frac1n \sum_{i=1}^n 1(\hat{y} \neq y) =  \sum_{i=1}^n \frac1n 1(\hat{y} \neq y)$$

The weights ideally sum to 1, but we need not have each weight be $\frac1n$. Feel familiar? 

## Boosting Algorithm

1. Start with a class of models $F$ (e.g. tree, glm, etc). And weights $w_i = \frac1n$
2. Fit a model $M_k$ in that class using those weights.
3. Find the prediction error for each observation from that model.
4. Increase the weights for observations where predictions were most wrong, decrease the weights for observations where predictions were most correct.
5. Repeat steps 2-4 until you've built $K$ models.

Model $M_K$ will be using weights based on how poorly each previous model did at predicting each observation. 

```{r,message=F,warning=F, echo=F}
#Package
library(mboost)
#Tree version
mod1 = blackboost(y~.,data=df,family=Binomial())
#GLM
mod2 = glmboost(y~.,data=df,family=Binomial())
#GAM (relative of trees)
mod3 = gamboost(y~.,data=df,family=Binomial())

#Make predictions
grid$preds1 = predict(mod1,newdata=grid,type="response")
grid$preds2 = predict(mod2,newdata=grid,type="response")
grid$preds3 = predict(mod3,newdata=grid,type="response")
```

## Boosting Trees
```{r}
ggplot(grid,aes(x=x1,y=x2,col=preds1))+geom_point()
```

## Boosting GAM

```{r}
ggplot(grid,aes(x=x1,y=x2,col=preds3))+geom_point()
```

# Data Cleaning

## Data Cleaning 101

Rules:

> 1. Keep looking at the data
> 2. Make small changes. 
> 3. Test your changes before you overwrite variables.
> 4. Don't overwrite actual files unless you're certain.
> 5. Don't throw away potentially useful data.

## Basics

What is data cleaning?


> (A) We get raw data -- which consists of all kinds of nonsense and structure  

> (B) We want to have a matrix $x$ to model an outcome $y$ with.

> To get from $A\rightarrow$B$, we may need to remove some structure, build our own structure, deal with missing values, get rid of irrelevant data, add in relevant variables, deal with sampling design, and more. This is cleaning. 

# Formats

## CSVs aren't everything

Data you want comes in many forms.[E.g. See SCF](https://www.federalreserve.gov/econres/scfindex.htm) 

- flat: CSV, arrow, parquet, etc
- compressed: .zip, .RData, etc
- proprietary: .xlsx, STATA, SAS, etc
- Internet-standards (scraped): html, json, etc

Converting between these is non-trivial and critical. 

> - But there are mostly packages for this -- so a google search will typically solve your problems, and I'll assume you can get to a data.frame or two.

## Other Formats

Most other non-trivial data formats are an entire courses worth of material. E.g.How to deal with:

    - paragraphs of text
    - Images

At a high level. My advice for those situations _at the moment_ is to train a model (or borrow a pre-trained one) that can deal with those things. 

## Images

For example, Images are really just a matrix with a color at each point. Detecting a face or something in that matrix is well explored at this point, and well beyond what we can cover here.

- We can use other people's plug-and-play models for face detection as inputs.
    - E.g. run that model to determine "is there a face", then use "face-presence" as a variable.
  
## Text

Text models are similar. Well explored, mostly beyond the time we have in this course. But you can take someone else's sentiment analysis, or subject analysis as an input to your models very easily.

## Models are Legos

I want you to see models as building blocks. We have some goal, we can combine 10 models to build inputs that we throw into an 11th model to help us with our goals.

1. This is a multi-step thing.
2. Most of the time, good use of data involves doing a good job constructing other pieces from it.

# Dates

## Seconds since 1970

Most dates are stored as an integer value representing the number of seconds since 12:00:00 AM on Jan 1, 1970.

```{r}
time = Sys.time()
time
as.integer(time)
```

## Days since 1970

```{r}
date = Sys.Date()
date
as.integer(date)
```

This is the number of days since Jan 1, 1970. 

## Advantages

The great thing about setting up variables like this, is that we can easily determine the gap between two dates.

```{r}
newtime = Sys.time()
newtime-time
```

That is how long it took to run the code in the middle there.

Beyond this -- it is also a specific, definitely identified time. 

## Problems

I don't care about the number of seconds or days since 1970. 

For basically this reason, dates are notorious to work with. They are usally stored in a format which is focused on being good at things that we don't care about. 

> Subsecond accuracy also going to require other things...

## Other Formats

You may be used to this format:

```{r}
date = "2021-05-11 8:05:32 AM"
is.character(date)
```

This is a character. I can understand what it is saying, but it can't be added and subtracted.

Moreover, it is underidentified. 

- Is this date in November or May? Need more info.
- Is this time during Business hours in NYC?
    - Need timezone info. 
    - Including Daylight Savings

## Dates in R

Tidyverse has the only good package for working with dates I've encountered. 

```{r,message=F}
library(lubridate)
```

Package comes with a host of useful functions for dealing with dates. E.g. I can convert characters to "Dates".

```{r}
date = "2021-05-11"
ymd(date) #This doesn't look different, but it is.
```
## Example

```{r}
as.integer(ymd(date))
as.integer(date)
```


## Example: CDC Data

[Posted Data from my prediction](https://codowd.com/bigdata/predictions/cdc_cases.csv)

```{r}
cases = read_csv("https://codowd.com/bigdata/predictions/cdc_cases.csv",skip=2)
head(cases$Date)
```

## Example

```{r}
cases = cases %>%
  mutate(Date=mdy(Date))
head(cases$Date)
head(as.integer(cases$Date))
```

## Example

I care about the day of week -- going to focus on Sundays
```{r}
cases = cases %>%
  mutate(dow = wday(Date))
head(cases$dow)
```

These are "Thurs","Weds","Tues","Mond","Sun","Sat" respectively. 

## Example Etc

Other functions for year, month, etc. 

```{r}
cases = cases %>% mutate(months = month(Date))
head(cases$months)
```

## Fixed Effects 

We haven't discussed Fixed Effects. 

But if you wanted Month and Year fixed effects, using the month and year functions to create columns, converting those columns to factors and using those factors as inputs would do it.

```{r}
cases = cases %>% mutate(years = year(Date))
head(cases$years)
```

## Fixed Effects
But if you wanted Month by Year (AKA Month x Year) fixed effects, you need to combine both.
```{r}
cases = cases %>% mutate(m_y = paste0(months,"-",years))
sample(cases$m_y,5)
```

That will give us a character again. Which can become a factor, etc.
    

# Wide vs Long

## Panel Data

A reasonably common data structure is known as 'panel data'.

In this setting, we have a number of units (e.g. States) which we observe. We also have a number of times at which we observe those units. There are some number of variables we observe. 

There are many ways you could store this information.

## Wide

Wide makes each unit a column and each time a row. Then it stores our variable of interest in the cells.

```{r,message=F}
cigs = read_csv("cigarettes.csv")
cigs
```

## Long

Long has many rows. In particular, if there are $n$ units, $T$ times, and 1 variable we observe, Long has $nT$ rows and 3 columns.

```{r,echo=F}
cigs_long = cigs %>% pivot_longer(!year,names_to="State",values_to="purchases")
```

```{r}
cigs_long
```


## Converting Wide-to-Long and back

There are _many_ different ways to pivot. Tidvyverse has a number of builtin tools. You could write something to do it by hand. Etc. 

```{r, eval=F}
cigs_long =cigs %>%
      pivot_longer(cols = !year,
                  names_to="State",
                  values_to="purchases")
```

```{r, eval=F}
cigs_wide = cigs_long %>%
    pivot_wider(names_from = State,
                values_from = purchases)
```

## Complications

That looked nice and smooth.

As soon as you have complications, like missing values, or other nonsense, it gets worse.

# Merges/Joins

## Multiple Sources

Sometimes you have two different sources of information. 

E.g. Zillow home data -- we have a pile of characteristics for each home. We also have sale prices and sale dates for each home. 

Why are these separate? 

## Looking at Zillow Data

```{r,echo=F,message=F,cache=T}
load("~/Github/bigdata/lectures/l10/zillow/prop_16_A.RData")
load("~/Github/bigdata/lectures/l10/zillow/prop_16_B.RData")
prop = rbind(propA,propB)
rm(propA,propB)
sales = read_csv("~/Github/bigdata/lectures/l10/zillow/train_2016_v2.csv")
```
```{r}
sales
```
Each place could be sold multiple times. 

## Looking at Zillow Data

```{r}
prop
```

But they have a lot of properties already, with stable characteristics. Much easier to store these things separately.

## Merges

But we want to use the characteristics to predict sales qualities. So we need all the data combined.

This is a merge. 

## Merge Problems: Non-perfect Match

```{r}
nrow(prop)
nrow(sales)
```

We do not have a sale for every property. In many situations, we may not have property characteristics for every sale (not a problem here). 

## Merge Problems: Non-perfect Match

In this case, we don't care about non-matched properties. We want a training data set. We will subset to the data in the sales (our outcome).

Sales data without characteristics would pose a larger concern. But only accepting "matched" values helps. 

## Merge Problems: Duplicates

```{r}
sum(duplicated(sales$parcelid))
```

And we have some duplicates. 

The easiest thing to do with duplicates is to throw them away. 

# STOP

## Merge Problems: Duplicates

Throwing away duplicates could cause problems. What kind of problems?

If homes that sell frequently are fundamentally different from those that don't, and are important to our target questions, then throwing them out may bias our whole procedure. 

What then? Keep all sales _if it won't break your model_. 

## Joins

There are a few basic Joins between Data A and B:

- Left Join: Keep all data in A, add in data from B when available.
- Right Join: Keep all data in B, add in data from A when available
- Inner Join: Keep all observations that are in both A and B
- Outer Join: Keep all observations. Set `NA` for missing columns for data that isn't in both datasets.

We can't build a model using characteristics to predict sales errors without data in A and B, so we want an inner join. 

## Joining

```{r}
data = inner_join(sales,prop,by=c(parcelid="parcelid"))
data
```

## IDs

How did we do this? 

We had a unique ID for each property. The `parcelid`. This ID was present in both datasets. So if we saw the same ID in each dataset, we knew the property was the same and we could match it.

Sometimes this happens. And when it does, it becomes very easy to combine information across different sources. 

- IDFA: unique ID for each apple device
- SSN: unique ID for most US taxpayers
- URL: unique ID for most websites
- IP Address: unique ID for most web-capable devices

## IDs pt 2

But often we don't have this information for some of the data we are using.

Methods:

1. Try to reconstruct an Identifier. 
     - SSNs historically were relatively easy to identify if you knew birthdate and place of birth. 
2. Match on other dimensions.
    - Maybe you don't observe SSN in any of your data. Or place of birth. 
    - But you see a billing address, a name, and a phone number.
    - We can match on those dimensions too. 

## IDs

Matching on other characteristics becomes tricky. The ideal is to say you have a match when multiple characteristics are all identical. 

In practice, even this is unlikely. 

- Is "Connor Dowd" the same as "Connor J Dowd", or "Connor M Dowd" or "connor dowd"?
- Billing addresses may leave off apartment numbers, etc. 
- One of the datasets may be missing one or more details. Does "1(800)333-2283" = NA?

## IDs

But these things add up. Even if the address is an apartment building with 5k people in it, you've narrowed things down from 7 billion to 5k. Now how many are named "Connor" or "Dowd"? 

Matching on these things rapidly becomes more art than science. There are many resources I can point you to if you have more Qs.

## Privacy In two minutes

Privacy in data is mostly out of bounds for this class.

- If you are publishing data, you need to be careful about what you include that people could use to match with. 
  - "Anonymising" by dropping e.g. names, is not enough
- The recommended procedure is to add some noise to your data before publishing if you're worried
  - E.g. Randomly change 2 digits in each phone number, randomly adjust numbers by small amounts, randomly drop some data and add in a few fake data points

# NAs

## What is NA? 

Usually `NA` denotes missingness. We don't know whether or not home 1235532 has a fireplace. So we don't put `TRUE` and we don't put `FALSE`, we put `NA`. 

Once again, much like duplicates, the standard advice is to drop observations with some `NA` values. 

## Dropping NAs

Dropping `NA` observations tends to be justified by the following assumption:

- The data is missing _at random_

That is -- the observations where we don't know are more like a clerical error than a selection issue. It is not that missing observations are fundamentally different, it is that someone forgot to check a box. 

This is frequently implausible. 

## Dropping NAs Nevertheless

Even under that assumption, dropping those observations can be devastating to your sample size.

```{r}
data %>% drop_na()
```

## Dropping NAs

```{r}
colMeans(apply(data,2,is.na))[1:24]
```

## NAs

Some of those columns _were entirely missing observations_. We can probably ditch those columns.

But some of them just had a lot of missing observations. Like `fireplacecnt`, which is missing in 89% of homes. 

It seems likely that it is not missing at random.

```{r}
table(data$fireplacecnt,useNA="ifany")
```

It seems likely that people without fireplaces didn't enter a number for fireplace count. Dropping those observations would be a mistake. Most of them are just 0s. 

## NAs 

Instead, what you should do is try to also model the `NA`s. 

For some model types and variables this is difficult. We could probably just make `NA` a 0 for `fireplacecnt`. But what about for `airconditioningtype` which is just a factor anyhow? Or for `Square footage`, where it definitely exists?

## NAs

For factor variables, this is straightforward. Just make one of your levels "NA". But what about Square footage?

A good rule of thumb _can_ be to replace NA with a value that is impossible. E.g. replace `NA` with `-1` for square footage. 

1. This will be easily spotted by others using the data, so you don't screw them up. 
2. For flexible model types, the model can now just try to make predictions when square footage is negative.
    - E.g. for tree/forest model or KNN, if the model wants it can easily partition those observations away from everything else. 
    
## NAs

```{r,eval=F}
na_helper = function(x) {
  if (is.factor(x)) {
    levels(x) = c("-1",levels(x))
    x[is.na(x)] = -1
  } 
  if (is.logical(x)) {
    x = x*1
    x[is.na(x)] = -1
  }
  if (is.character(x)) {
    x[is.na(x)] = "-1"
  }
  if (is.numeric(x)) {
    x[is.na(x)] = -1
  }
  x
}
```


## NAs 
```{r,eval=F}
prop2na = prop2 %>% mutate(fips = as.numeric(fips))
prop2na = prop2na %>% mutate(across(where(is.character),as.factor))
prop2na = prop2na %>% mutate(across(everything(),na_helper))
```

This will not work as well for a linear model -- where its predictions for square footage `-1` affect its predictions for square footage of 2000.

# Predictions 2

## My forecast

My forecast was very straightforward. 
```{r}
small = cases %>% filter(dow == 1 | dow == 5) #Pull out Thursdays and sundays
small %>% filter(Date > "2021-04-25") %>% select(`New Cases`, Date)
```

```{r}
37885*(44766/60196)
```

## Log Diff-in-Diff

Diff in Diff says that the change between April 29th and May 6th is going to be the same as the change between May 2 and May 9. This is a presumption that the trends hold.

Log Diff-in-Diff says that the percent change from April 29th to May 6th will be the same as the percent change from May 2 to May 9. 

```{r}
logdiff = log(44766)-log(60196)
logsun = log(37885)
logpred =  logsun + logdiff
pred = exp(logpred)
pred
```

## Picking That
I also checked to see if the log diff-in-diff model using Thursday's 7-day averages was better. It wasn't. 

This was _extremely_ easy to do out of sample _because_ there was nothing estimated. Each week, I purely take 3 numbers from the prior 8 days to predict sunday. 

Additionally, this was automatically adaptive to sundays. Its a week-on-week change times the last sunday. Any "sunday fixed effect" gets incorporated automatically. 

## Sundays are Important

A lot of models looked tuned to predict an average next day (or 3rd day away), rather than a Sunday. 
```{r}
ggplot(cases,aes(x=Date,y=`New Cases`,fill=dow==1))+
         geom_col(size=0.1)
```


## Price-Is-Right

There was no real wrong answer here. I decided that I would use the 20th percentile for my prediction. A lot of people seemed to have similarly ad hoc procedures. 

- _predicting is different in competitions_

But a big difference was confidence interval percentiles vs predictive interval percentiles. 

## CI vs PI
Let me illustrate the difference with some made up data.
```{r}
n = 100
x = rnorm(n)
y = rnorm(n)+2*x
df = data.frame(x=x,y=y)
```
```{r}
ggplot(df,aes(x=x,y=y))+
  geom_point()+
  geom_smooth(method="lm",se=T)
```

## Probabilities

When you have an unbiased distribution of errors (i.e. out-of-sample distribution) for a model, getting event probabilities becomes somewhat straightforwards.

Take your prediction. Add those errors. See how often thing happens.

Almost all of you did something like this. Well done.

## Ensemble Probabilities
What do you do if you have 10 models though? Find probability _under_ each model -- then combine.

1. Build each model.
2. Get honest error distributions for each model.
3. Predict P[event] given error distributions for each model
4. Combine $P[event|M_1],...,P[event|M_{10}]$ in a "sensible way"

## Sensible? 

If you're averaging these models for your predictions, you can average these probabilities. 

What I'm trying to avoid is the following procedure:

1. Build each model.
2. Average to get ensemble prediction 
3. Look at distribution of predictions around average to determine uncertainty

And other similar mistakes. We are uncertain about the best model, and each model has beliefs about $P[event]$. We need to combine them well. 

## Model Uncertainty

What if we only built one model? 

- How much do you trust it?

If you had a model telling you that tomorrow the stock price for Apple was guaranteed to go up 50%, does that tell you about Apple or about the model? 

If a guy at a casino tells you he has a system, and you should bet it all on red for _certain_ winnings, are you now certain of winnings if you bet it all on red? 

## Model Uncertainty

More broadly, we want to know about the ways in which our models might fail that are relevant. 

So an important question is "how often does this model fail in ways relevant to this questions?"

For models that are reasonably certain about an event -- this will dominate our uncertainty. E.g. we don't know the guy at the casino, so our difficulty trusting him dominates the chance that we don't win on red.

## Model Uncertainty 

Why does this matter?

Models in big data settings can become reasonably certain about events not happening, or happening. They have a lot of data, they have a model, the two combine for a lot of certainty.

But the model being certain does not mean _you_ need to be certain. _The data wants to trick you_

And in competitive settings, like a casino, or the stock market, "the data wants to trick you" is less a metaphor than you might think. 

# Wrap up

## Things to do

Homework 5 is due tomorrow.

See you Thursday.  


# Bye!








