part-II-writeup.Rmd

---
title: "Part-II-Writeup"
subtitle: "Group 5"
output: pdf_document
---
```{r libraries, include=FALSE}
library(tidyverse)
library(GGally)
library(ggpubr)
library(knitr)
library(BAS)
library(glmnet)
knitr::opts_chunk$set(echo = FALSE)
```

```{r read-data, echo=FALSE}
load("paintings_train.Rdata")
load("paintings_test.Rdata")
load("paintings_validation.Rdata")
```

```{r Useful Functions,echo=FALSE}
# no exp
rmse = function(y, ypred) {
  rmse = sqrt(mean((y - ypred)^2))
  return(rmse)
}

# exp rmse
ermse = function(y, ypred) {
  ermse = sqrt(mean((exp(y) - exp(ypred))^2))
return(ermse)
}

# coverage
coverage <- function(y, cl, cu){
  mean( cl <= y & cu >= y )
}
```


# Part 1: Introduction

For this project, we are being asked to explain which factors played an important role on a selection of paintings sold in Paris in the late 18th century. We will use these factors to identify undervalued/overvalued paintings in the dataset. To facilitate this, we have been provided with auction price data from 1764-1780 containing information on painting/artist characteristics and the sale itself. First, we will use EDA to get to know the data. This will help us get an initial sense of which variables are associated with price, and which variables to include in our preliminary models.

After EDA is complete, we will move on to the model building portion of the analysis. Using the important variables identified in EDA, we will create a multiple linear regression model attempting to explain price as a function of these covariates. After considering potential interactions between the selectors, we will try to further improve/simplify our model via variable selection with AIC/BIC. Finally, we will validate our model by checking the assumptions of linear regression, and assert that we have not made any major violations. How well our model performs will be a result of several factors:

- Does our model correctly capture what the true connection between variables and price is? I.e., does our model have low bias?
- How well does our model fit the data relative to the null model? A common metric for this is the root mean squared error (RMSE), which is a function of the residuals of our model. We are looking for a low RMSE on "new" data (the test dataset).
- Does our model perform well on new data? I.e., does our model have low variance? If we add too many predictors, we risk overfitting our model to the training data and jeopardizing performance on new data.
- Does our model "cover" the predictions well? From our model, we can create prediction intervals for where we expect the price to fall. We would like the true value to be covered most, but not all of the time. If we're covering the true value all of the time, our prediction intervals are too wide, and our prediction intervals are not as precise as we would like.

After our analysis with linear regression, we will proceed forward with a more complicated model. This will be contained in the subsequent writeup.

For the second part of this project, we are focusing on obtaining predictive accuracy as the primary objective. We are allowed to lose some interpretation if the prediction is better. We will investigate splines, BMA, etc. as these are powerful tools that can still be somewhat interpreted.

# Part 2: EDA

To begin with the EDA, first we got a summary of the variables and discard any unneeded portions. Our response variable is `logprice`. As such, we will not use `price` in our analysis. `count` is a column of all 1's, and provides no value. The data is a mix of continuous, binary and categorical variables. The variables `position`, `nfigures` and the ones that are related to the dimensions of the painting(`Surface`, `Width_in` etc) are all continuous. `year` can be taken either as continuous or as factors, we choose to use it as a factor. The variable `lot` is a number so can be used as numeric. The variable pertaining to features of the paintings(like `lrgfont`, `engraved` etc), `diff_origin`, `artist_living` and `Interm` are binary and we will use them as factors. All the other variables are multilevel categorical variables.

## Investigating Binary Variables

Next, we note that a large amount of the variables are binary. We can investigate the association with price via boxplots:

```{r boxplots, warning = F, echo=FALSE, fig.cap="Boxplots between logprice and some binary variables"}
paintings_train %>%
  gather(diff_origin,engraved,prevcoll,paired,finished,lrgfont,lands_sc,lands_elem,Interm,
         othgenre,portrait,still_life,discauth,key = "var",value = "value") %>%
  ggplot(aes(group = value,y = logprice, x = value)) +
  geom_boxplot(na.rm = T) +
  facet_wrap(~ var, scales = "free") +
  theme_bw() + xlab(label = "")
```

***Figure 1*** provide a great visualization on how a given binary variable is associated with changes in logprice. It's hard to quantify how much of a difference in price should warrant variable inclusion, but as we will be performing variable selection later, it shouldn't hurt to include more than necessary. We see the following potential predictors from the graphs:

- `prevcoll`
- `lrgfont`
- `diff_origin`
- `finished`
- `engraved`
- `lands_sc`
- `portrait`
- `paired`
- `still_life`
- `lands_elem`
- `Interm`
- `figures`
- `discauth`
- `othgenre`

We note that many of these binary variables seem to be strongly associated with `logprice`, so we will consider them further when we are building the initial model.

The variable `Interm` has many missing values corresponding to the winning bidder being "Unknown", and we found that including that "missingness" as a level, the logprice was lower than the when we know whether there is an intermediary or not. Surprisingly, for the same observations, `type_intermed` is just blank, same as when `Interm` is 0, not missing. So we can just use `type_intermed` instead of `Interm` and assume that the blanks in place of NAs are imputations, or maybe that's how it was supposed to be coded.

We also note that the variable `peasant` and `othgenre` are like two levels of a factor(they cannot be 1 together as described in the codebook, and both describe the scene of the painting). So, while we don't think `peasant` has any significant impact on `logprice`, we may include it anyway.

## Investigating categorical variables

Several variables have multiple categories (country origin, type of endbuyer, etc.). We can visualize these with box plots, there will just be more boxes than a binary variable.

```{r boxplots for categorical vars, warning = F, echo=FALSE, fig.cap="Boxplots between logprice and some categorical variables"}
paintings_train %>%
  gather(dealer,origin_author,origin_cat,school_pntg,winningbiddertype,mat,type_intermed, key = "var",value = "value") %>%
  ggplot(aes(group = value,y = logprice, x = value)) +
  geom_boxplot() +
  facet_wrap(~ var, scales = "free") +
  theme_bw() + xlab("")
```

```{r boxplots for year, warning = F, echo=FALSE, fig.cap="Boxplots between logprice and year", fig.height=4, fig.width=6}
paintings_train %>%
  gather(year, key = "var",value = "value") %>%
  ggplot(aes(group = value,y = logprice, x = value)) +
  geom_boxplot() +
  facet_wrap(~ var, scales = "free") +
  theme_bw() + xlab("year")
```

Several of the categories for variables have very small sample sizes, and as such cannot be used for our analysis (`mat`, `school_pntg`, `winningbiddertype`, `type_intermed`). The other variables initially seem feasible. It would make sense that nationality could be an important factor; perhaps certain countries are known more for their art. Similarly, it makes sense that `year` could be associated with price. Perhaps based on the economic circumstances, people would be more willing to spend different amounts of money on art. There doesn't seem to be a linear trend, however, so we will need to treat `year` as a categorical variable rather than numeric.  

We will also choose only `origin_author` from `origin_author`, `origin_cat` and `school_pntg`, because they are likely to be the same in many observations and out of the three, `origin_author` has all the categories and has reasonable amount of observations for each category compared to `school_pntg`.

## Continuous Variables

Finally, certain variables are continuous in nature. We can investigate the relationship with price via scatterplots.

* Note: For the `position` variable, there were a few observations whose values were not between 0 and 1. We assumed this was an error caused by percentage being entered instead of fractions. These values were corrected by dividing by 100.

```{r data preprocess, echo=FALSE}
# clean position(Only R1771 has position more than 1)
ind <- which(paintings_train$sale == "R1771")
paintings_train$position[ind] <- paintings_train$position[ind]/100
# test
ind <- which(paintings_test$sale == "R1771")
paintings_test$position[ind] <- paintings_test$position[ind]/100
```

```{r scatterplots of logprice with those continuous varaibles, message=FALSE, warning = F, echo=FALSE, fig.height=4, fig.width=8, fig.cap="Scatterplots between logprice and position as well as log(Surface)"}
ggs <- ggpairs(data = paintings_train %>% mutate(logSurface = log(Surface)), columns = c(3,29,10), progress = F)
pts <- lapply(1:2, function(x) getPlot(ggs,3,x))
ggarrange(plotlist = pts)
```

In no way does `position` seem to be associated with price, so we will not consider it further. We expected `Surface` to play a much larger role in price. Bigger paintings have more expensive materials and take artists longer to finish. We see from the graph, however, that there is only a very slight increase in `logprice` as `Surface` increases. As such, we will not include for now as it doesn't seem to be as related as the other variables. Additionally, the test data contains several NA values for `Surface` compared to the other variables, and improper imputation could be misleading.

## EDA between Predictors

Apart from the EDA for the predictors vs. `logprice`, on looking at the relationships between the predictors, we found a few things which seem important:

- Only dealer "R" can have `lrgfont` = 1. It also looks like the dealer is most likely "R" when `engraved` or `prevcoll` are 1.
- The variable `lands_sc` and `lands_elem` are never 1 together in the training data, but there very few cases in the test data. We observe the same in the variables `portrait` and `still_life`.
- There is only one observation with `type_intermed` "EB" along with `winningbiddertype` being "EBC", same is the case in test data. Since, there is only one observation of "EB", it will be difficult to use `type_intermed` unless we change something.
- Whenever there is an intermediary involved, the end buyer is most likely a collector.
- The dealers are not dealing each year, every dealer has certain years where they sell paintings, in other years, they don't. This is observed in both the training and the test data.

## Variable Exclusion

Further explanation of why certain variables were included/excluded are included in the formal modal analysis below. For brevity, they will not be included here.


# Part 3: Discussion of preliminary model Part I

We were fairly happy with our initial model. As we kept removing predictors, every criteria improved (RMSE, coverage, maxDeviation, and MeanAbsDeviation). However, removing predictors has the effect of increasing the bias of our model. We found the bias-variance tradeoff to be worth it, but our model started performing poorly on the training set.

From this, we believe that the training and test data sets are different. We believe this for 2 reasons:

1. The null model outperformed every other model, suggesting that adding predictors made things worse.
2. We used cross-validation to split our training data set into a sub-training and sub-testing data set for several different seeds. The final model fit using the sub-training data performed well on the sub-testing data, suggesting that it should perform well on the test data if the training and testing data were similar. However, we noted much poorer performance on the test data than the sub-test data, again suggesting that the 2 data sets are different.

From these results, we would like to fit a model that performs well on the training data set, then use shrinkage to improve test performance. We need to choose a shrinkage parameter $\lambda$ that balances fitting the training data well with good test performance.

However, we were later given the correct test set and found our previous model performed much better, so we thought we might not need to constrain the model anymore.

# Part 4: Development of the final model

## Moving from Initial to Final

```{r modify dataset, include=FALSE}
inds <- unlist(map(paintings_train, function(pred){
  if(is.numeric(pred)){
    if(all(pred == 1 | pred == 0 | is.na(pred)))
      return(TRUE)
    else
      return(FALSE)
  }
  else
    return(FALSE)
}))
# numeric binary to factors
paintings_train[inds] <- map_df(paintings_train[inds], as.factor)
paintings_test[inds] <- map_df(paintings_test[inds], as.factor)
paintings_validation[inds] <- map_df(paintings_validation[inds], as.factor)

# year to factors
paintings_train$year <- as.factor(paintings_train$year)
paintings_test$year <- as.factor(paintings_test$year)
paintings_validation$year <- as.factor(paintings_validation$year)

# Converting type_intermed "EB" to "E"
ind_tmp <- which(paintings_train$type_intermed == "EB")
paintings_train$type_intermed[ind_tmp] <- "E"
ind_tmp <- which(paintings_test$type_intermed == "EB")
paintings_test$type_intermed[ind_tmp] <- "E"
ind_tmp <- which(paintings_validation$type_intermed == "EB")
paintings_validation$type_intermed[ind_tmp] <- "E"
```

```{r}
model1 = lm(logprice ~ (dealer + endbuyer + type_intermed)*(diff_origin + discauth + lrgfont + engraved + prevcoll + paired) + (endbuyer + type_intermed)*(figures + finished + portrait + still_life + peasant + othgenre + lands_sc + lands_elem) + year + origin_author, data = paintings_train)

model2 <- step(model1, trace = F)
```

```{r include=FALSE}
set.seed(5)
model3 = bas.lm(formula(model2), data=paintings_train,
                   prior = 'g-prior', modelprior=uniform(),
                   method='MCMC',
                   n.models=5000, MCMC.iterations=1000000, thin=20, 
                   force.heredity = TRUE)
diagnostics(model3)
```

```{r include=FALSE}
image(model3,rotate = F)
plot(model3, cex.lab = 0.7, which=4)
```

```{r include=FALSE}
ytr_bma <- predict(model3, estimator = "BMA", se.fit = T)
ytr_bma <- data.frame(fit = ytr_bma$fit, lwr = ytr_bma$fit - 1.96*ytr_bma$se.bma.pred,
                      upr = ytr_bma$fit + 1.96*ytr_bma$se.bma.pred)
ermse(ytr_bma$fit, paintings_train$logprice)
coverage(paintings_train$logprice, ytr_bma$lwr, ytr_bma$upr)
save(ytr_bma, file="predict-train.Rdata")
```

```{r}
yte_bma <- predict(model3, newdata = paintings_test, estimator = "BMA", se.fit = T)
yte_bma <- data.frame(fit = exp(yte_bma$fit), lwr = exp(yte_bma$fit - 1.96*yte_bma$se.bma.pred),
                      upr = exp(yte_bma$fit + 1.96*yte_bma$se.bma.pred))
predictions <- yte_bma
save(predictions, file="predict-test.Rdata")
```

### Using the old test data

Using the "incorrect" test data set, we followed the same steps as fitting the initial model up to the first step function. Here we noted that the model was performing exceptionally well on the training data compared to the test data. In part I we started to manually remove predictors, but it was difficult to check which order of removal would be optimal. So instead we now opt for Bayesian Model Averaging (BMA) to see if the model can shrink and which variables are unimportant. This also has the added benefit of checking for correlation between the variables.

BMA removed a couple variables, but it showed practically every other predictor was associated with the response with high inclusion probability. Additionally, we didn't find much evidence of correlation between the predictors. The testing performance for this model is still bad, so we will further shrink the model using constraints. Lasso is the technique we used for this.

We applied Lasso regression on the best predictive model (BPM), giving an optimal value of $\lambda = 0$ (no shrinkage). As stated previously, this did not result in good performance, so we wanted to manually choose a lambda. After testing, we found $\lambda = 0.3$ to offer similar performances on the training and testing data sets.

### Using the updated test data

With the updated test data set, we found that our model started to perform better than the null model and it's no longer overfitting. So we changed our approach from removing predictors and adding constrains to including more predictors to reduce the bias. We started with initial model again then used a step function but with AIC instead of BIC to gain better prediction performance. Then we used BMA again on the output from the step function and found the model was better than other models we have tried.


## Residual Discussion

```{r}
par(mfrow= c(1,2))
plot(model3, which = 1)
resd <- -fitted(model3)+paintings_train$logprice
qqnorm(resd)
qqline(resd)

# checking the linear model with bpm predictors
# bpm_formula = as.formula("logprice~ dealer + endbuyer + type_intermed + diff_origin +
#    lrgfont + engraved + prevcoll + paired + finished +
#    portrait + still_life + lands_sc + year + origin_author +
#     dealer:diff_origin")
# plot(lm(bpm_formula,data = paintings_train),which = 4)
# plot(lm(bpm_formula,data = paintings_train))[4]
```

From these two plots, we found that the residuals of our model were normally distributed around zero and the variance of residuals across fitted values looks constant. We also didn't find any outliers from above plots. It seems that the assumption of normality is satisfied here.


## Predictive Intervals

We used the `predict` function to get the standard errors of the predicted values, `se.bma.pred`. Then we manually calculated the 95% prediction interval.

## Model Summary

Below we show the coefficient estimates(converted to exponential scale) with their 95% credible intervals and their inclusion probabilities in BMA.

```{r}
df <- as.data.frame( exp(confint(coef(model3))[, 1:3]) )
df$incl_prob <- summary(model3)[1:90,1]
kable(df[order(df$incl_prob, decreasing = T),c(3,1,2,4)], 
      col.names = c("Fit", "2.5%", "97.5%", "Inclusion Probability"), digits =3) 
```

From the table above, we see that lots of varaibles have inclusion probability that is close to 1, these predictors are `dealer`, `endbuyer`, `type_intermed`, `origin_author`, `year`, `prevcoll`, `paired`, `finished`, `lrgfont`, `engraved`, `diff_origin`, `lands_sc` and the interaction term `dealer:diff_origin`. Out of these, `paired` and `lands_sc` reduce the price while all the other painting features(the binary predictors) increase it. `portrait` and `still_life` have high probabilities but are not too close to 1 and both seem to reduce the price.  
Some predictors, especially all the other interaction terms, have inclusion probability close to 0, and very few predictors have intermediate inclusion probability, suggesting presence of correlation among them.

Since we have a lot of predictors with high inclusion probabilities, we discuss a few of the predictors below which seem to have a large impact on the price of the paintings. We discuss the effects while ignoring the interaction effects, even though most of them do not seem to have any impact(coefficient of 1) anyway.

- For paintings sold in year 1777, we expect the median price in livres increase by 1000% (with a credible interval of (723%, 1420%)), compared to the baseline year 1764. For paintings sold in year 1769, we expect the median price in livres increase by 680% (with a credible interval of (364%, 1237%)), compared to the baseline year 1764. Some of other years also saw an increased price but these two years of sale seemed to be the largest.
- We expect the median price of paintings sold by DealerL in livres increase by 440% (with a credible interval of (192%, 966%)), compared to the baseline dealerJ, when `paired`,`diff_origin`,`discauth`,`prevcoll` are zero. We expect the median price of paintings sold by DealerR in livres increase by 279% (with a credible interval of (163%, 466%)), compared to the baseline dealerJ, when `paired`,`diff_origin`,`discauth`,`prevcoll` are zero. Dealer L supposedly sold the highest priced paintings followed by Dealer R then Dealer P then Dealer J.
- We expect the median price of paintings bought by a buyer in Livres to be 276% (with a credible interval of (99%, 611%)) more than the baseline when the endbuyer is completely unknown(no information, not even the name) when `finished`,`diff_origin`,`land_sc` are zero.
- We expect the median price of paintings with the origin of the author unknown to be 66%(with a credible interval of (19%,86%)) lower than when the author is Austrian. From the results, it seems that Austrian, Dutch/Flemish and Spanish authors' paintings are more expensive.
- We expect the median price of paintings with an additional paragraph included in livres increase by 155% (with a credible interval of (99%, 226%)).
- We expect the median price of paintings in livres to increase by 130% (with a credible interval of (74%, 206%)) if the dealer mentions engravings done after the painting.

```{r include=F}
plot(coef(model3), ask = F)
```

We also looked at the posterior densities of the model coefficients. Most of the densities either had a large point mass at 0, or were normally distributed and had a single mode. The few coefficients that were multimodal were `diff_origin`(very small bump), `discauth` and `paired`. Looking into correlations between these and the other predictors may help us improve the model or at least give us some insight into their relationship with price. 

```{r}
yva_bma <- predict(model3, newdata = paintings_validation, estimator = "BMA", se.fit = T)
yva_bma <- data.frame(fit = exp(yva_bma$fit), lwr = exp(yva_bma$fit - 1.96*yva_bma$se.bma.pred),
                      upr = exp(yva_bma$fit + 1.96*yva_bma$se.bma.pred))
yva_predictions <- yva_bma
save(yva_predictions, file="predict-validation.Rdata")
```


```{r include=FALSE}
BPM_predict=predict(model3,estimator="BPM")
variable.names(BPM_predict)
```

# Part 5: Assessment of Final Model

## Model Evaluation

First, let's get a sense of how much of the variation in price our model is explaining

```{r}
cor(paintings_train$logprice, ytr_bma[, 1])^2
```

Our model is explaining approximately 66% of the variation in the log price of paintings. It is important to note that an $R^2$ cannot be used in a vacuum to tell if the model is good or bad, but it's useful to note. We observe an RMSE of 1538.9, and a coverage of 0.956(using a 95% prediction interval). As stated previously from the residual plot, it looks like all of the assumptions for MLR have been satisfied, so our standard errors + confidence intervals will be appropriate.

## Model Testing

For the testing data, we obtained a bias of 220.59. Relative to other groups, this seems like a reasonable bias. We could attempt to include more variables, but we think we already captured the relevant variables associated with log price. Adding more variables here could be beneficial, but likely wouldn't result in too much of a bias improvment. The other statistics here seem reasonable as well:

- Our coverage is 0.956, which is very good
- Our maxDeviation is 13429.63, almost half of some of the other groups. Therefore there are no egregiously different values that our model is predicting.
- Our MeanAbsDeviation is 455.68.
- Finally, our RMSE is 1262.07. For some reason, the RMSE is lower on the test data than the training data. It's possible that this is due to the smaller sample size of the test data. Or it could have "easier" cases to predict.

We believe this model is reasonable due to its large improvement over the null model (no predictors).

## Model Result

```{r}
top_ind <- order(yva_predictions[,1], decreasing = T)[1:10]
top_paintings <- paintings_validation[top_ind,]
top_paintings[, c("dealer", "authorstandard", "year")]
```

Here are some observations from this:

- The dealer for ALL of the top 10 predicted paintings is R.
- Only one painting was not sold in 1769 or 1777. This observation was only different by one year (1776). These years must have had a good market, or R might have been dealing more than the others.
- One author has 2 paintings in the top 10: Rembrandt Harmenszoon van Rijn. This is the famous Rembrandt, so it's expected. One of the paintings, Arquebusiers / Ronde de nuit, is "The Night Watch", a very famous painting. The other is difficult to find due to the French subject with a incorrectly translated character. It's likely one of the "Holy Family" paintings, a series of 5 famous paintings.
- Most of the endbuyers were collectors that all used dealers as intermediaries.
- All of the top predicted paintings had artists that were either of Dutch/Flemish or French origin. For these, the dealer's catalog also correctly labeled their origin.
- All of these paintings had "lrgfont": the dealer devotes an additional paragraph to each painting

# Part 6: Conclusion

## Summary of the results

`Dealer`, `Year`, `Endbuyer` and `Origin_author`(when the author is unknown) seem to impact the price of the paintings most(in terms of the coefficient values). We found that the people involved are more important than the painting features, when the people are unknown, the features may become more important though. For example, when the artist in unknown, painting price drops a lot. We also found that certain dealers were more active in certain years. We were also surprised that whether the author was living or not was not important to the price.

If we just want to focus on the features of the paintings and not the people involved, the best features for a painting to have are highly polished finishing, an additional paragraph and engravings done after the painting. The worst features a painting could have are if it is described as a plain landscape, if it is described as a portrait and if its description indicates still life elements. It would be advisable to look for the presence(or absence) of these features to get an idea about the price, however, our model also includes the effects of the people involved in the sale, so we need to be careful if we look at these features on their own. A model to just predict the price of the painting using the features may be better if we only want to look at their effects, but this model would not consider facts such as the reputation of the artist or the dealer etc.

On the old test dataset, we tried various different models including tree models and found that using BMA followed by lasso on the best predictive model gave us a model that could give us reasonable insight into the training data while also not poorly predicting on the test data. 

On the new dataset, we used just the BMA model without adding any constraints using lasso. If we had more time, we could have tested tree models like random forests to see if the results are any better. Since, we are just using categorical variables in our model anyway, a tree model might be better suited to this problem in terms of predictions since it can capture many more interactions than a linear model can. The BMA model is still useful as it can give easily interpretable results and we can compare the effects of various predictors.

## Things learned

- We should do EDA before we try to build a model and don't forget to consider collinearity bewteen variables.
- Data is not always as clean as the ones we see in homework, we can't always remove observations with missing values. There are some methods to do the missing value imputation, for example, using the MICE package. As far as this project concerned, we found the `Surface` variable had random missing values, but we didn't think we need to impuate the missing values of this variable because it's not very important.
- It's easy to overcomplicated with models, since we have learned lots of techniques. But sometimes keeping it simple may be the best stratagy.