---
title: 'Math 50 Fall 2017, Homework #5'
---
<hr>
__NOTE: For your homework download and use the template__ (https://math.dartmouth.edu/~m50f17/HW5.Rmd)

__Read the green comments in the rmd file to see where your answers should go.__


<br><hr><br>


### Normal probability plots 

You can use qqnorm and qqline functions to plot probability plots of residuals. The functions 
rstandard and rstudent
calculate the standardized residuals and R-student residuals, respectively. 

```{r}

prop = read.table("https://math.dartmouth.edu/~m50f17/propellant.csv", header=T, sep=",")
age <- prop$Age
shearS <- prop$ShearS

fitted = lm(shearS ~ age)

stdRes = rstandard(fitted)
rStuRes = rstudent(fitted)

qqnorm(rStuRes, main="Normal Probability Plot (residuals on vertical axis)") 
qqline(rStuRes)

```

However, note that the book uses residuals on the x-axis instead of y-axis. In order to obtain that use the parameter datax as shown below.  In the below graph x-axis denotes the R-student residuals and the y-axis is the theoretical quantiles ( in the book y-axis is probability instead of quantiles).  



```{r}

qqnorm(rStuRes, datax = TRUE ,
       main="Normal Probability Plot") 
qqline(rStuRes, datax = TRUE )


```

### Residual vs predicted response plot 

```{r}
yHat <- predict(fitted)

plot (yHat, rStuRes)
abline(0,0)
```


After observing that the observation points 5 and 6 look like potential outliers, next we delete those points and compare the fitted model of the deleted data with the full data.    

```{r}

plot(age, shearS, xlim=c(0,30), ylim=c(1600,2700))
abline(fitted$coef, lwd = 2, col = "blue")

ageRem <- age[-6] 
ageRem <- ageRem[-5] 

shearSRem <- shearS[-6]
shearSRem <- shearSRem[-5]

fitted2 = lm(shearSRem ~ ageRem)
abline(fitted2$coef, lwd = 2, col = "red")

```


### A note 
There is a dedicated library : 

   MPV: Data Sets from Montgomery, Peck and Vining's Book
   
in order to provide an easy way to load tables from the book. 
To install the library type : 

install.packages("MPV")

Below is an example how to use this library. Check https://cran.r-project.org/web/packages/MPV/MPV.pdf for table names.

```{r}

library(MPV)
data(table.b1)

y <- table.b1$y
x1 <- table.b1$x1
x3 <- table.b1$x3
x8 <- table.b1$x8

y.lm <- lm(y ~ x1 + x3 + x8)
summary(y.lm)


```


<br><hr><br>
<!-- ############################################################################################ --> 
<!-- ################################# NEW QUESTION ############################################# --> 
<!-- ############################################################################################ --> 


## Question-1 
Solve the parts (a), (b) and (c) of Problem 4.1. 
In addition answer the following. 

(d)  Is it possible to perform lack of fit test using the steps (4.20) to (4.24) ?  


### Answer: 
<!-- If question requires a written answer it goes here. R code goes to the chunk below.) -->

```{r}
library("MPV")

fitted = lm(y~x8, data=table.b1)

cat("Part (a). 
    There are some minor deviations but in general it seems OK and we can say no problem with the normality assumptions.")

rStuRes <- rstudent(fitted)
qqnorm(rStuRes, datax = TRUE, main = "Normal Probability Plot")
qqline(rStuRes, datax = TRUE  )


cat("Part (b). 
    Plot seems in a band and it doesn't suggest any violations of our assumptions.")

yHat <- predict(fitted)
plot(yHat, rStuRes)
abline(0, 0) 


cat("Part (c). 
     Plot shows a structure (trend in different parts) so we conclude that adding x2 might contribute to the model.")

plot(table.b1$x2, rStuRes)
abline(0, 0)


cat("Part (d). 
     The lack of fit test explained in steps (4.20) to (4.24) requires repeated observations, therefore not possible to do in this case.")


```


<br><hr><br>
<!-- ############################################################################################ --> 
<!-- ################################# NEW QUESTION ############################################# --> 
<!-- ############################################################################################ --> 


## Question-2 
Chapter 4, Problem 2 all parts. 



### Answer: 
<!-- If question requires a written answer it goes here. R code goes to the chunk below.) -->

```{r}

cat("Part (a). 
    There are some deviations from the normality in the first part. First half above and second half seems to be below the normal distribution line (qq-line).")

fitted <- lm(y ~ x2 + x7 + x8, data=table.b1)
rStuRes <- rstudent(fitted)
qqnorm(rStuRes, datax = TRUE, main = "Normal Probability Plot")
qqline(rStuRes, datax = TRUE)


cat("Part (b). 
 This also seems in a band with no apparent violations in our assumptions from this plot.")

yhat <- predict(fitted)
plot(yhat, rStuRes)
abline(0, 0)

cat("Part (c).
    The plots of residuals of x7 shows non-constant variance (variance increasing in the first half). Plot against x2 has also some narrowing in the middle and suggests non-constant variance. On the other hand plot against x8 is a little better, although it shows some narrowing towards ends there are too few points (only two points in each end). ")

plot(table.b1$x2, rStuRes)
abline(0, 0)

plot(table.b1$x7, rStuRes)
abline(0, 0)

plot(table.b1$x8, rStuRes)
abline(0, 0)

cat("Part (d).
    We can see the linear relation between y and the regressors x2 and x8. For the regressor x7 this is less apparent. ")

fit.y_x2deleted <- lm(y ~ x7 + x8, data=table.b1)
res.y_x2deleted <- y - predict(fit.y_x2deleted)
fit.x2 <- lm(x2 ~ x7 + x8, data=table.b1)
res.x2 <- table.b1$x2 - predict(fit.x2)
plot(res.y_x2deleted, res.x2)

fit.y_x7deleted <- lm(y ~ x2 + x8, data=table.b1)
res.y_x7deleted <- y - predict(fit.y_x7deleted)
fit.x7 <- lm(x7 ~ x2 + x8, data=table.b1)
res.x7 <- table.b1$x7 - predict(fit.x7)
plot(res.y_x7deleted, res.x7)

fit.y_x8deleted <- lm(y ~ x2 + x7, data=table.b1)
res.y_x8deleted <- y - predict(fit.y_x8deleted)
fit.x8 <- lm(x8 ~ x2 + x7, data=table.b1)
res.x8 <- table.b1$x8 - predict(fit.x8)
plot(res.y_x8deleted, res.x8 )


library(MASS)
stuRes <- studres(fitted)

cat("Part (e).
    No R-Student is above our threshold 3, the observation point with 2.54 R-student residual can be considered as an outlier.  The difference between the studentized and R-student residuals is small (see below). This shows Si^2 and MSres are not significantly different than from each other and this suggest that there is no influential point with a strong  influence on variance estimation.")
summary(rStuRes)
summary(rStuRes-stuRes)

 
```


<br><hr><br>
<!-- ############################################################################################ --> 
<!-- ################################# NEW QUESTION ############################################# --> 
<!-- ############################################################################################ --> 


## Question-3 
Chapter 4, Problem 19 all parts. In addition answer the following. 


(c) Find the point with largest (in absolute value) r-student residual as a potential outlier. Repeat the regression analysis after deleting that point from the observation data. Construct the probability plot and residual vs predicted response plot. Calculate the differences (deleted vs full data) in fitted coefficients,  $MS_{res}$ and $R^2$. Comment on the differences in the plots and the values.   Do you think it is an influential point?  Do they imply any improvement? 


### Answer: 
<!-- If question requires a written answer it goes here. R code goes to the chunk below.) -->

```{r}

data(p4.19)
fitted <- lm( y ~ x1+x2+x3, data=p4.19)
stdRes <- rstandard(fitted)
rStuRes <- rstudent(fitted)

cat("Part a.
    The normality assumtions seems to be valid : ")

qqnorm(rStuRes, datax = TRUE , main="Normal Probability Plot") 
qqline(rStuRes, datax = TRUE )

cat("However in the below plot we see a violation of the linearity assumption, (as the behaviour in the middle suggests that). Similar to above one of the observations even though not exactly above our cut-off value 3, it is separated from the rest (absolute value of its R-Student) and therefore can be considered as an outlier. ")
 
yHat <- predict(fitted)
plot(yHat, rStuRes)
abline(0, 0)

# Only 6 data repeated at:  x=0
# So we can calculate directly 
repeatedY_atX0 <- c(133, 133, 140, 142, 145, 142)
SS_PE <- sum((repeatedY_atX0 - mean(repeatedY_atX0))^2)
SSres <- sum( (p4.19$y - yHat)^2 );
SS_LOF <- SSres - SS_PE

m <- 9
p <- 4
n <- length(yHat)
F_lof <- ((SS_LOF/(m-p)))/(SS_PE/(n-m))
fQuant <- qf(0.95,  m-p, n-m) 

cat ("Part b.
      Fo value is ",F_lof ,"  and this is slighlt greater than the critical quantile value (at significance level 0.05)",fQuant, "
       Thus we conclude that lack-of-fit test suggest that linear model is not adequate.") 


deletedData <- p4.19[-(which(rStuRes==max(rStuRes))), ]
delFit <- lm(y~x1+x2+x3, data=deletedData)
rStuResD = rstudent(delFit)

qqnorm(rStuResD,  datax = TRUE , main="Normal Probability Plot") 
qqline(rStuResD,  datax = TRUE )

yHatD <- predict(delFit)
plot(yHatD, rStuResD)
abline(0,0)

MSres <- sum( (p4.19$y - yHat)^2)/(length(yHat)-4)
MSresD <- sum( (deletedData$y-yHatD)^2)/(length(yHatD)-4)

summary(fitted)
summary(delFit)

cat("Part c.
     As we see from summary above R2 values dropped. MSres before was ", MSres, " 
     and after deletion it is  ", MSresD, "
    , There is a considerable drop in MSres. 
     We can say removing this observation point have some improvement in MSres, 
     but this is expected since it was a potential outlier. Howver in
     x-space removed point is not well separated from the rest (see table), thus not an influential point.")



```

