## Eyeballing your Models

- Plots to assess model fit 
- Plots generated by base R
- Using fortify() and ggplot()
- Customize your graphs 


## Linear regression assumptions and plots to diagnose them

Assumption | Plot
----------|:----------
Independence of residuals | It has to do with study design or data collection. No plot for cross-sectional data
Homoscedasticity (constant variance of the residuals) | Residuals or standardized residuals vs fitted values
Residuals are normally distributed | Q-Q plot of residuals
Linear relationship between dependent and independent variables | Plot residuals versus individual independent variables and/or residuals vs fitted values





## Example 1 - well-behaved data 

Generating the data
```{r}
set.seed(13)
x = runif(n = 100, min = 0, max = 10)
y = 2 + 0.8*x + rnorm(n = 100, mean = 0, sd = 1)
mydata = data.frame(x,y)
```

Fitting a linear regression
```{r}
mylinreg = lm(y ~ x, data = mydata)
```

## Example 1: well-behaved data 

```{r, tidy=TRUE}
summary(mylinreg)
```

## Plots generated by base R
```{r}
par(mfrow = c(2, 2))
plot(mylinreg)
par(mfrow = c(1, 1))
```

## Plots generated by base R
Linearity and homoscedasticity assumptions
```{r, echo=FALSE}
plot(mylinreg, which = 1)
```

## Plots generated by base R
Normality assumption
```{r, echo=FALSE}
plot(mylinreg, which = 2)
```

## Plots generated by base R
Homoscedasticity assumption
```{r, echo=FALSE}
plot(mylinreg, which = 3)
```

## Plots generated by base R
Detect extreme values
```{r, echo=FALSE}
plot(mylinreg, which = 5)
```


## Using fortify() and ggplot()
- fortify(): Adds the variables listed below to the original dataset
  
Assumption | Plot
----------|:----------
.hat | Measure of the leverage of a data point
.sigma | Estimate of residual standard deviation when corresponding observation is dropped from the model
.cooksd | Cook’s distances
.fitted | Fitted values
.resid | Residuals
.stdresid | Standardized residuals


## Using fortify() and ggplot()

With fortified data we can recreate the plots given by plot() and we can even modify them
```{r}
library(ggplot2)
myfortdata = fortify(mylinreg)
head(myfortdata)
```

## Homoscedasticity assumption
```{r, fig.height=4, fig.width=6}
ggplot(data = myfortdata, aes(x = .fitted, y = .resid)) +
  geom_hline(yintercept = 0, colour = "firebrick3") +
  geom_point() 
```

## Linearity assumption
```{r, fig.height=4, fig.width=6, message=FALSE}
ggplot(data = myfortdata, aes(x = .fitted, y = .resid)) +
  geom_hline(yintercept = 0, colour = "firebrick3") +
  geom_point() +
  geom_smooth(se = FALSE)
```

## Normality assumption
```{r, fig.height=4, fig.width=6}
ggplot(data = myfortdata, aes(sample = .stdresid)) +
  stat_qq() +
  geom_abline(colour = "firebrick3")
```

## Standardized residuals vs fitted values
```{r, fig.height=4, fig.width=6}
ggplot(data = myfortdata, aes(x = .fitted, y = .stdresid)) +
  geom_hline(yintercept = 0, colour = "firebrick3") +
  geom_point() 
```

## Residuals vs. leverages
```{r, warning=FALSE, fig.height=4, fig.width=6, message=FALSE}
ggplot(data = myfortdata, aes(x = .hat, y = .stdresid)) +
  geom_point() +
  geom_smooth(se = FALSE) 
```

## Customized plots
```{r, fig.height=4, fig.width=6}
# Start with a scatter plot of your data
ggplot(data = myfortdata, aes(x = x, y = y)) +
  geom_point() 
```

## Customized plots
```{r, fig.height=4, fig.width=6}
# Explore the relationship using "loess"
ggplot(data = myfortdata, aes(x = x, y = y)) +
  geom_point() +
  stat_smooth(method = "loess", se = FALSE) 
```

## Customized plots
```{r, fig.height=4, fig.width=6}
# Add a regression line 
ggplot(data = myfortdata, aes(x = x, y = y)) +
  geom_point() +
  stat_smooth(method = "lm", formula = y ~ x, se = FALSE)
```

## Customized plots
```{r, fig.height=4, fig.width=6}
# Adding observation number
ggplot(data = myfortdata, aes(x = .fitted, y = .resid)) +
  geom_hline(yintercept = 0, colour = "firebrick3") +
  geom_text(label = rownames(myfortdata))
```

## Customized plots
```{r, fig.height=4, fig.width=6, warning=FALSE}
# Points size reflecting Cook's distance
ggplot(data = myfortdata, aes(x = .fitted, y = .resid, size = .cooksd)) +
  geom_hline(yintercept = 0, colour = "firebrick3") +
  geom_point() +
  scale_size_area("Cook’s distance") 
```

## Customized plots
```{r, warning=FALSE, fig.height=4, fig.width=6, message=FALSE}
# Residuals vs. leverages with observation number
ggplot(data = myfortdata, aes(x = .hat, y = .stdresid)) +
  geom_point() +
  geom_smooth(se = FALSE) +
  geom_text(label = rownames(myfortdata), size = 4) +
  geom_hline(yintercept = 2, lty = "dotted", colour = "slateblue1") +
  geom_hline(yintercept = -2, lty = "dotted", colour = "slateblue1") 
```


## Example 2: midwest data 
```{r}
mylinreg2 = lm(poptotal ~ popwhite + popblack, data = midwest)
summary(mylinreg2)
```

## Explore your data
```{r, fig.height=4, fig.width=6}
myfortdata2 = fortify(mylinreg2)
ggplot(data = myfortdata2, aes(x = popwhite, y = poptotal)) +
  geom_point() +
  stat_smooth(method = "loess", se = FALSE)
```

## Explore your data
```{r, fig.height=4, fig.width=6}
ggplot(data = myfortdata2, aes(x = popwhite, y = poptotal)) +
  geom_point() +
  stat_smooth(method = "lm", formula = y ~ x + I(x^2), colour = "red", se = FALSE) +
  stat_smooth(method = "lm", formula = y ~ x + I(x^2) + I(x^3), colour = "blue", se = FALSE)
```

## Plots generated by base R
```{r}
par(mfrow = c(2, 2))
plot(mylinreg2)
par(mfrow = c(1, 1))
```

## Homoscedasticity assumption
```{r, fig.height=4, fig.width=6}
ggplot(data = myfortdata2, aes(x = .fitted, y = .resid)) +
  geom_hline(yintercept = 0, colour = "firebrick3") +
  geom_point() 
```

## Normality assumption
```{r, fig.height=4, fig.width=6}
ggplot(data = myfortdata2, aes(sample = .stdresid)) +
  stat_qq() +
  geom_abline(colour = "firebrick3")
```

## Standardized residuals vs fitted values
```{r, fig.height=4, fig.width=6}
ggplot(data = myfortdata2, aes(x = .fitted, y = .stdresid)) +
  geom_hline(yintercept = 0, colour = "firebrick3") +
  geom_point() 
```

## Residuals vs. leverages
```{r, warning=FALSE, fig.height=4, fig.width=6, message=FALSE}
ggplot(data = myfortdata2, aes(x = .hat, y = .stdresid)) +
  geom_point() +
  geom_smooth(se = FALSE)
```



## Next examples
In the previous examples we focused in plots to assess model assumptions once we have fitted the model. In the following two examples we will explore the data and assess linear relationships before fitting any model, specially for correlated data.


## Example 3: iris data
```{r}
data(iris)
summary(iris)
```

## Explore your data
```{r, fig.height=4, fig.width=6}
# Scatter plot and smoother
ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length)) +
  geom_point() +
  stat_smooth(method = "loess", se = FALSE)
```

## Explore your data
```{r, fig.height=3.85, fig.width=6}
# Faceting by Species
ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length)) +
  geom_point() +
  stat_smooth(method = "loess", se = FALSE) +
  facet_wrap( ~ Species)
```

## Explore your data
```{r, fig.height=4, fig.width=6}
# Pattern within each group in the same graph
ggplot(data = iris, aes(x = Petal.Length, y = Sepal.Length, colour = factor(Species))) +
  geom_point() +
  stat_smooth(method = "loess", se = FALSE)
```



## Example 4: BodyWeight data
```{r, fig.height=4, fig.width=6}
library(nlme)
data(BodyWeight)
summary(BodyWeight)
```

## Explore your data
```{r, fig.height=4, fig.width=6}
# Scatter plot and smoother
ggplot(data = BodyWeight, aes(x = Time, y = weight)) +
  geom_point() +
  stat_smooth(method = "loess", se = FALSE)
```

## Explore your data
```{r, fig.height=4, fig.width=6}
# Individual trajectories over time
ggplot(data = BodyWeight, aes(x = Time, y = weight, group = Rat)) +
  geom_point() +
  geom_line()
```

## Explore your data
```{r, fig.height=3.85, fig.width=6}
# Faceting by Rat
ggplot(data = BodyWeight, aes(x = Time, y = weight, group = Rat)) +
  geom_point() +
  geom_line() +
  facet_wrap( ~ Rat, ncol = 4)
```

## Explore your data
```{r, fig.height=4, fig.width=6}
# Observed mean over time for all the data
ggplot(data = BodyWeight, aes(x = Time, y = weight)) +
  geom_point() +
  stat_summary(fun.y = mean, geom="line", colour = "red")
```


## Explore your data
```{r, fig.height=3.85, fig.width=6}
# Observed mean over time by Diet
ggplot(data = BodyWeight, aes(x = Time, y = weight)) +
  geom_point() +
  stat_summary(fun.y = mean, geom="line", colour = "red") +
  facet_wrap( ~ Diet)
```


## Additional Resources
- [Color options in ggplot](http://sape.inf.usi.ch/quick-reference/ggplot2/colour)
- [Smooth options in ggplot](http://www.ats.ucla.edu/stat/r/faq/smooths.htm)
- [Practical Regression and Anova using R by Julian J. Faraway](http://www.ats.ucla.edu/stat/r/sk/books_pra.htm)




