Regression Modelling Using R Part-2

title: "Regression Modelling Part-2"
author: "Dr Nisha Arora"
date: "Monday, June 06, 2016"
output: html_document
---

# Regression Diagnostics
## Playing with different Models

To import csv data file & Understand the data
```{r}
adver=read.csv("Advertising.csv"); names(adver);
adver=adver[, -1]
dim(adver); summary(adver); fix(adver)
par(mfrow=c(2,2))
plot(adver$TV, adver$Sales, xlab="TV", Ylab = "Sales", main="Sales Vs TV")
plot(adver$Radio, adver$Sales, xlab="Radio", Ylab = "Sales", main="Sales Vs Radio")
plot(adver$Newspaper, adver$Sales, xlab="Newspaper", Ylab = "Sales", main="Sales Vs Newspaper")
hist(adver$Sales)
```

To fit different multiple regression models
```{r}
summary(lm(Sales~Radio, data=adver)); 
attach(adver); summary(lm(Sales~Newspaper))
m1=lm(Sales~TV, data=adver); summary (m1)

m2=lm(Sales~TV+Radio, data=adver); summary(m2)

m3=lm(Sales~., data=adver); summary(m3)
```
We can observe a strange behaviour of predictor 'Newspaper'. While the newspaper regression **coefficient estimate in simple regression is significantly non-zero, the coefficient estimate for newspaper in the multiple regression model is close to zero, and the corresponding p-value is no longer significant**, with a value around 0.86.
**The simple and multiple regression coefficients can be quite different.**
This difference stems from the fact that in the simple regression case, the slope term represents the average effect of a $1,000 increase in newspaper advertising, ignoring other predictors such as TV and radio. In contrast, in the multiple regression setting, the coefficient for newspaper represents the average effect of increasing newspaper spending by $1,000 while holding TV and radio fixed.

To compare these models (Nested models)
```{r}
anova(m1, m2, m3)
```
Results shows that adding third predictor 'Newspaper' does not significantly reduces RSS.
NOte: See more about coxtest, jtest, encomptest & waldtest for comparing non-nested models.

To look at diagnostics of model 'm2'

```{r}
plot(m2)
std_resid <- (resid(m2) - mean(resid(m2))) / sd(resid(m2))
plot(adver$TV, std_resid, xlab="TV", Ylab = "Std_resid", main="Std_resid Vs TV")
plot(adver$Radio, std_resid, xlab="TV", Ylab = "Std_resid", main="Std_resid Vs Radio")
hist(std_resid); curve(dnorm, add = TRUE)
probDist <- pnorm(std_resid)
plot(ppoints(length(std_resid)), sort(probDist), main = "PP Plot", xlab = "Observed Probability", ylab = "Expected Probability")
abline(0,1)

```
Plots shows presence of extreme points, deviation of residual from normality and presence of herteroskedasticity

```{r}
Pairs(adver)
install.packages("PerformanceAnalytics"); library(PerformanceAnalytics)
chart.Correlation(adver, histogram=TRUE, pch=19)
```
It shows statistically significant correlation between Radio & Newspaper

## Applying statistical tests for model diagnostics
### To test if residuals have zero mean [using one sample t-test]
```{r}
t.test(resid(m2), alternative="two.sided", mu = 0, conf.level = 0.95)
```
It shows that expected value of error term is zero. [Assumption satisfied]

### To test if residuals have constant variation [Homoskedasticity]

```{r}
library(lmtest); gqtest(m2)
```
Result of  Goldfeld-Quandt test also shows presence of heteroskedasticity at 10 % level of significance

```{r}
library(lmtest); hmctest(m2)
```
Result of  Harrison-McCabe test also shows presence of heteroskedasticity at 10 % level of significance

```{r}
library(lmtest); bptest(m2)
```
Result of  Breusch-Pagan test shows presence of heteroskedasticity at 10 % level of significance

```{r}
library(car); ncvTest(m2)
```
Result of  NCV test shows presence of heteroskedasticity at 5 % level of significance

### To test if residuals have no auto-correlation [Durbin Whatson test]
```{r}
dwtest(m2)
```
Results of Durbin Watson test shows presence of auto-correlation in the data

### To test if residuals are normally distribued
```{r}
hist(resid(m2)); hist(std_resid)
shapiro.test(resid(m2)); shapiro.test(std_resid)
shapiro.test(adver$Sales)
```
Results shows deviation from normality of residuals

### To test if error terms & predictors are truly independent
```{r}
library(car); vif(m2); vif(m3)
```
Results shows there is no multi-collinearity in the data

### To test if there is model specification error
```{r}
resettest(m2)
```
Results of Ramsey’s RESET test shows that we can reject the null hypothesis "that the functional form of the model is correctly specified" It shows that we have model specification error

### 
```{r}
harvtest(m2)
```
Result of Harvey-Collier test suggests that the true relationship is linear

## Let's understand the relatioship
Notice that the **correlation between radio and newspaper is 0.35.** This reveals a tendency to spend more on newspaper advertising in markets where more is spent on radio advertising.
Now suppose that the multiple regression is correct and newspaper advertising has no direct impact on sales, but radio advertising does increase sales. Then in markets where we spend more on radio our sales will tend to be higher, and as our correlation matrix shows, we also tend to spend more on newspaper advertising in those same markets. Hence, in a simple linear regression which only examines sales versus newspaper, we will observe that higher values of newspaper tend to be associated with higher values of sales, even though newspaper advertising does not actually affect sales. **So newspaper sales are a surrogate for radio advertising; newspaper gets “credit” for the effect of radio on sales.**

Also, **Suppose that spending money on radio advertising actually increases the effectiveness of TV advertising, so that the slope term for TV should increase as radio increases.**

## Let's include interaction term in the model
```{r}
m4=lm(Sales~TV*Radio, data=adver); summary(m4);plot(m4)
```
## To remove extreme points
```{r}
adver1=adver[-c(6,131, 179),]; 
m5=lm(Sales~(TV+Radio+TV*Radio), data=adver1); summary(m5); plot(m5)
```
## To use heteroskedasticity correction
See 'sandwich package' or 'car package' for various Heteroskedasticity & Autocorrelation-consistent estimation of the covariance matrix of the coefficient estimates in regression models

```{r}
library(sandwich); coeftest(m5, vcov(m5))
coeftest(m5, vcovHC(m5, "HC3"))
library(car); coeftest(m5,vcov=hccm(m5))
library(car); summaryR(m5, type="hc0") # White SE (1980)
summaryR(m5, type="hc1") # STATA default
summaryR(m5, type="hc2") # Unbiased under heteroskedasticity
summaryR(m5, type="hc3") # Default (Conservative)
```

**"Essentially, all models are wrong, but some are useful."**

--- Box, George E. P.; Norman R. Draper (1987). Empirical Model-Building and Response Surfaces, p. 424, Wiley. ISBN 0471810339.