```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
download.file('https://wkristan.github.io/template.docx', 'template.docx', mode = 'wb')
```

## Classic model selection

We will work with data from the USDA National Nutrient Database. Import the data set:

```{r import.data}



```

### Graphing the data

To get an idea of the distribution of each of the variables, we can plot a histogram for each. To save space, we'll plot them in a 4x4 array.

```{r plot.untransformed.histograms, results='hide'}



```


Now plot the transformed histograms:

```{r plot.transformed.histograms, results='hide'}



```

Run this chunk to add a function, panel.cor(), that will put the correlation coefficient and p-value into the lower panel of the correlation matrix plot.

```{r panel.cor.function}

panel.cor <- panel.cor <- function(x, y, cex.cor = 0.8, method = "pearson", ...) {
options(warn = -1)                              # Turn off warnings (e.g. tied ranks)
usr <- par("usr"); on.exit(par(usr))            # Saves current "usr" and resets on exit
par(usr = c(0, 1, 0, 1))                        # Set plot size to 1 x 1
r <- cor(x, y, method = method, use = "pair")   # correlation coef
p <- cor.test(x, y, method = method)$p.val      # p-value
txt <- format(r, digits = 3)                    # Format r-value
txt1 <- format(p, digits = 3)                   # Format p-value
txt2 <- paste0("r= ", txt, '\n', "p= ", txt1)   # Make panel text
text(0.5, 0.5, txt2, cex = cex.cor, ...)        # Place panel text
options(warn = 0)                               # Reset warning
}

```

Next we should see how the variables relate to one another. First, we will make a scatterplot matrix using kcal, along with the un-transformed versions of the predictors:

```{r kcal.untrans.pairs.plots}



```

Now make one for kcal against the log-transformed versions of the predictors:

```{r kcal.trans.pairs.plots}



```

Using the log of kcal might help - plot it against the un-transformed predictors:

```{r log.kcal.untrans.pairs.plots}



```

Finally, plot the transformed log.kcal against the transformed predictors:

```{r log.kcal.trans.pairs.plot}



```

#### Question: why is it better to use the transformed versions of the variables instead of the raw variables, based on the scatterplots?

> 

### Fitting models

Make an empty list to hold the models:

```{r make.model.list}



```

Add all of the single predictor models to the list:

```{r add.single.predictor.models.to.list}



```

Get the r.squared values for each model:

```{r r.squared}



```

Get the adj.r.squared values:

```{r adj.r.squared}



```

Combine the r.squared and adj.r.squared into a data frame to make them easier to see:

```{r p.values}



```


#### Question: which single predictor variable model has the highest adjusted r-squared?

> 

The two predictors log.ash and log.protein are strongly correlated. Fit a model with both included and get a summary here:

```{r log.ash.log.protein.model}



```

#### Question: did the model with both log.ash and log.protein have a higher r.squared than than either of them by themselves? Did it have a higher adj.r.sqaured than either by themselved?

> 

Fit a model with all of the predictors included:

```{r all.variables}



```

#### Question: did the model with all predictors included have a higher adj.r.squared than all the others so far?

> 

Use a Type II SS ANOVA table to decide whether there are variables that are not contirbuting statistically significant amounts of explained variation to the model:

```{r type.II.ss.anova}



```

#### Question: which variables are statistically significant?

> 

Fit a model with only the statistically significant predictors from the ANOVA:

```{r signif.predictors.model}



```

#### Question: did the model with only the significant predictors have a better adj.r.squared than all the others?

> 

#### Question: log.protein did not make it into the best model. Does this mean that protein has no calories? Explain.

> 

Double check that the three significant predictors that you dropped are not just statistically redundant with one another by testing the model without them against the model that included them:

```{r test.reduced.against.full.model}



```

#### Question: how did your test of the reduced model against the full model show you that log.ash, log.protein, and log.fiber are not just statistically redundant with one another?

> 

#### Question: look at the p-values for each of the models. Are any of the models non-significant? If you had just looked at a single model, would you have considered it to be a poor model based on the p-value alone for any of these? What does that tell you about the importance of fitting multiple models and comparing their adjusted R^2?

> 

To make sure we aren't selecting a "best" model that is in reality terrible, we should make sure that the model is statistically significant. Get the summary() of signif.pred:

```{r summary.of.signif.pred}



```

#### Question: is the model with only the statistically significant predictors significant overall? How do you know (what part of the output tells you this)?

>


## Models as hypotheses

Once you have your own hypotheses created in the Excel spreadsheet provided, import it as a dataset called food.hypotheses:

``` {r import.food.hypotheses}



```

To start, make an empty list for the models you will use (my H.example, and your two hypotheses), called food.hypotheses.list:

```{r make.food.hypotheses.list}



```

Fit a model with my hypothesis (H.example) as a predictor, and assign it to the food.hypotheses.list:

```{r make.H.example.model}



```

Code chunks for your own models should follow.
