---
title: Master of Supply Chain, Transport and Mobility - Data Analysis on Transport
  and Logistics - Course 16-17 - Partial Exam
author: 'Lecturer: Lidia Montero'
date: "November, 10th 2016"
output: word_document
---

## Problem 1: All questions account for 1 point

### Dataset: usair

*Data from a dataset of air pollution in US cities. Seven variables were recorded for 41 cities:*

  - SO2:		 Sulphur dioxide content of air in micrograms per cubic meter
  - NegTemp: 	Average annual temperature less than -1 Fahrenheit degrees
  - Manuf:	Number of manufacturing enterprises employing 20 or more workers
  - Pop: 		Population size (1970 census) in thousands
  - Wind: 		Average annual wind speed in miles per hour
  - Precip: 	Average annual precipitation in inches
  - Days: 		Average number of days with precipitation per year.

*Source*

Everitt, B.S. (2005), An R and S-PLUS Companion to Multivariate Analysis, Springer


**Load usair.RData file  in your current R or RStudio session **

**1.	Pop contains the description of thousands of inhabitants for the cities included in the data set. Create a new factor variable consisting on an indicator for small, medium and large cities (named it f.size). Small cities are those with less than half million inhabitants, medium cities are those in the range from half medium to one millium and a half and large cities have a number of inhabitants greater than one million and a half.**

```{r}
library(car)
library(FactoMineR)
load("usair.RData")
summary(usair)
usair$f.size<-factor(cut(usair$Pop,breaks=c(0,500,1500,3500)),labels=c("Small","Medium","Large"))
table(usair$f.size)
100*table(usair$f.size)/nrow(usair)
```


**2.	Our target is defined as SO2. Summarize numerically and graphically the response variable. Make an interpretation of the results. Do you think that SO2 may be considered normally distributed? **

*By graphical inspection throu an histogram with normal curve overlayed, clearly SO2 target is not normally distributed. Shapiro-Wilk test rejects the normality hypothesis.*
```{r}
mm<-mean(usair$SO2)
ssdd<-sd(usair$SO2)
hist(usair$SO2, freq=F)
curve(dnorm(x,mm,ssdd), col="red", lwd=2,add=T)
shapiro.test((usair$SO2))
```

**3.	Calculate the upper threshold to identify severe outliers for SO2. Is there any city satisfying this criteria? Are global outliers retained once city size factor is considered?**

*The upper threshold for severe outliers is at 3 times IQR from Q3, thus according to the summary of S02, 101 mg/m3. Chicago is the only severe outlier and the biggest city. Once the target is examined for each group defined by the city factor size, Chicago is not an outlier. Providence is a small city with a very high sulphur dioxide content in air, clearly a severe outlier in its group. In the medium size group, three group outliers in the high SO2 contents appear: Cleveland, Pittsburgh and St Louis*

```{r}
par(mfrow=c(1,2))
ss<-summary(usair$SO2);ss
# Upper threshold
utso2<-ss[5]+3*(ss[5]-ss[2]);utso2
Boxplot(usair$SO2,labels=row.names(usair))
abline(h=101,col="red",lwd=3)
Boxplot(usair$SO2~usair$f.size,labels=row.names(usair),col=heat.colors(3))
```

**4.	Which are the numerical variables statistically associated with the response (SO2)? Indicate the suitable measure of association and/or tests that support your answer. Assess linearity association to SO2 for available variables.**

*Spearman correlations between target and the rest of the variables have to be calculated. The most associated variables are NegTemp (average annual temperature less than -1F) and average number of days with rain/snow. Null hypothesis are postulated for each of the six possible pairs and tested with cor.test(), pvalues under 0.05 are found for Neg.Temp, Manuf, Pop and Days, thus significant statistical association is found.*

```{r}
round(cor(usair[,1:7],method="spearman"),dig=3)
cortest<-rep(0,6)
for (j in 2:7 ) {cortest[j-1]<-cor.test(usair[,1],usair[,j],type="spearman")$p.value}
names(cortest)<-names(usair)[2:7]
sort(cortest)
plot(usair[,1:7])
```

**5.	Describe the profile for the SO2 numeric target using the tools available in FactoMineR package.**

*An alternative answer could be given using condes() procedure in FactoMineR. Target is the first variable in dataset. Global association with the target is significant for Manuf, Pop, Neg.Temp and Days. Global association between SO2 target and the city size factor is significant. In the category output, the mean of SO2 for the group of Large cities (3 obs) is significantly 29.7 units over the grand mean*
```{r}
library(FactoMineR)
condes(usair,1)
```

**6.	The average SO2 in the cities can be argued to be the same for all city size levels (f.size)? Which are the groups that show a significant greater average SO2 than the others?**

*Non-parametric test on means for SO2 according to f.size have to be used. The null hypotesis of equal SO2 in all groups has a p.value of 0.054 according to Kruskal-Wallis test. It is on the limit, but since sample size is small we have to be cauteous and reject the null hypothesis, thus there is some average mean group different from the others. With the pairwise Wilcoxon tests pvalues for Mean S02(Large)=Mean S02(Small) and Mean S02(Large)=Mean S02(Medium) are rejected, but means for Small and Medium cities are clearly equal, so the mean sulphur dioxide content in air for large cities is clearly different (and higher) than those for small and medium cities.*

```{r}
Boxplot(SO2~f.size,data=usair)
kruskal.test(SO2~f.size,data=usair)
oneway.test(SO2~f.size,data=usair)  # Not suitable

with(usair,pairwise.wilcox.test(SO2,f.size))
with(usair,pairwise.t.test(SO2,f.size))

```

**7.	The variance of SO2 in the cities can be argued to be the same for all city size levels (f.size)? Which are the groups that are likely to have a greater dispersion of SO2 than the others?**

*Non-parametric test on variances for SO2 according to f.size has to be used. The null hypotesis of equal variance for SO2 in all groups has a p.value of 0.367 according to Fligner-Killeen test. Thus we do not have evidence to reject the null hypothesis, and thus, it is accepted. Dispersion of the target does not depend on city size group*

```{r}
Boxplot(SO2~f.size,data=usair)
fligner.test(SO2~f.size,data=usair)
bartlett.test(SO2~f.size,data=usair)  # Not suitable
```

**8.	Consider a multiple regression model (m1) for target SO2 on all numeric variables in the dataset. Assess the quality of the model.**

*The coefficient of determination is 67%, indicating that 2/3 of the target S02 variability is explained by the model (all numeric variables included in the dataset). Some of the variables have a p.value for the null *
```{r}
names(usair)
m1<-lm(SO2~.,data=usair[,1:7])
summary(m1)
```

**9.	Consider model (m1), check significance for all variables.**

*Some of the variables have a p.value for the null hypothesis test of coefficient equal 0 greater than 0.05, thus apparently, we have few observations - power of the test is low, Precip and Days are not significant. If we check for net-effects of the variables with Anova() method, perfectly suited for net-effect testing, the same conclussions appear: Precip and Days net-effects when the rest of the variable have already been included in the linear predictor are not significant*
```{r}
summary(m1)
Anova(m1)
```

**10.	Assess default residual plots in R for model (m1): are there any atypical residuals? Which one/s?**

*The first plot depicts the raw residuals vs fitted values according to the model, a noise pattern has to shown for valid models. In this case no patter is present, but some large residuals on the positive Y axis (Providence, Pittsburgh). Normality of the residuals is checked with a QQPlot showing 2-3 cities that are not on the QQline, again these cities are Providence, Pittsburgh and a new city Phoenix; at least for the first two cities, residuals are too large to follow a normal distribution. On the scale-location plot, the smoother line is not flat, indicating non constant variance, but, but with 41 observations one has to be cautelous. The last plot on the right-down part shows an atypical city according to its leverage, so far away from the multidimensional cloud of points included in the design matrix, that does not seem relevant because the residual is close to 0.*
```{r}
par(mfrow=c(2,2))
plot(m1,id.n=5)

```

**11.	For your model, determine the presence of outliers in residuals. Specify city names, selected criteria and behavioral discrepancy.**

*Atypical residuals appear for Providence and Pittsburgh, so lack of fit for these 2 cities are remarkable: the observed SO2 is much, much greater than the predicted value according to m1 model. Pittsburgh is a medium city with observed SO2 of 61 micrograms per cubic meter  while the model predicts 30.92  and Providence is a small city showing observed SO2 of 94 micrograms per cubic meter while the model predicts 45.24, so a large lack of fit is found for these two cities*
```{r}
par(mfrow=c(1,1))
llist<-Boxplot(resid(m1),labels=row.names(usair))  # For assessing atypical residuals
llist
predict(m1)[llist]
usair[llist,]
```

**12.	For your final model, determine the presence of actual influent data. Specify city names, selected criteria and behavior.**

*Boxplot for Cook's distance identifies two cities with large distances that lay far away from the distances for the rest of the observations: Providence (too high SO2 air contents) and Phoenix (low SO2 air contents, with a negative prediction according to m1. Both observations are influent data that show a lack of fit with the current model m1*
```{r}
sort(cooks.distance(m1),decreasing=T)
llist<-Boxplot(cooks.distance(m1),labels=row.names(usair))  # For assessing actual influent observations
llist <-llist[c(1,5)];llist
predict(m1)[llist]
usair[llist,]
```



