08-Exploring-state-data-set.Rmd

# The state dataset

## Reading in and manipulating data


The state data sets include state information in the early years around the 1970s in the USA. We pick state.abb, state.x77, and state.region to form our data files. The detailed information is listed here.

**state.abb**: A vector with 2-letter abbreviations for the state names.

**state.x77**:
A matrix with 50 rows and 8 columns giving the following statistics in the respective columns.

- Population: population estimate as of July 1, 1975

- Income: per capita income (1974)

- Illiteracy: illiteracy (1970, percent of population)

- Life Exp: life expectancy in years (1969-71)

- Murder: murder and non-negligent manslaughter rate per 100,000 population (1976)

- HS Grad: percent high-school graduates (1970)

- Frost: mean number of days with minimum temperature below freezing (1931-1960) in capital or large city

- Area: land area in square miles

**state.region**:
A factor containing the regions (Northeast, South, North Central, West) that each state belongs to.

For the convenience of analyzing the state data, let's merge the three data sets into a single data set "sta".  It is a data frame with 10 columns and 50 rows. 

```{r}
tem <- data.frame(state.x77)               # Transform matrix into data frame
sta <- cbind(state.abb, tem, state.region) # Combine the three data sets
colnames(sta)[1] <- "State"                # Rename first column
colnames(sta)[10] <- "Region"              # Rename the 10th column
head(sta)
str(sta)
summary(sta)
```

## Visualizing data

Let's start by visualizing the distributions of numeric variables. In many cases, we want to know if our data follows a normal distribution or not so that we can decide whether some methods are suitable or not. Here are some ways that we can check the normality of a variable.

- Histogram: Does the histogram approach a normal density curve? If yes, then the variable more likely follows a normal distribution.

- Q\-Q plot: Do the sample quantiles almost fall into a straight line? If yes, then the variable more likely follows a normal distribution.

- Shapiro-Wilk test: This is a widely used normality test. The null hypothesis is that a variable follows a normal distribution. Small p-value indicates a non-normality of the variable.

```{r message=FALSE, warning=FALSE, fig.width=8, fig.height=8}
a <- colnames(sta)[2:9] # Pick up all numeric columns/variables according to the names 
par(mfrow = c(4, 4))    # Layout outputs in 4 rows and 4 columns
for (i in 1:length(a)){ 
  sub = sta[a[i]][,1]   # Extract corresponding variable a[i] in sta
                        
  hist(sub, main = paste("Hist. of", a[i], sep = " "), xlab = a[i])
  qqnorm(sub, main = paste("Q-Q Plot of", a[i], sep = " ")) #
  qqline(sub)           # Add a QQ plot line. 
 
   if (i == 1) {
     s.t <- shapiro.test(sub) # Normality test for population
   } else {
     s.t <- rbind(s.t, shapiro.test(sub)) # Bind a new test result to previous row
 }
}
s.t <- s.t[, 1:2]       # Take the first two columns of shapiro.test result
s.t <- cbind(a, s.t)    # Add variable name for the result
s.t
```

From the histograms and QQplots we can see that the distributions of Population, Illiteracy, and Area skew to the left. Income and Life.Exp are distributed close to normal. The shapiro tests show that Income, Life.Exp and Frost are normally distributed with p-values greater than 0.05, while Murder and HS.Grad are almost normally distributed with p-values very close to 0.05. There is no evidence that Population, Illiteracy, and Area having normal distribution.

In the state data, there is a categorical variable *region* which contains 4 observations. What is the distribution of the categorical variable? Let's take a look at the number of observations(states) in each region and the corresponding percentage.

(ref:state-region) State count in each region.

```{r state-region, fig.cap='(ref:state-region)', fig.align='center'}
counts <- sort(table(sta$Region), decreasing = TRUE)  # Number of states in each region
percentages <- 100 * counts / length(sta$Region)      
barplot(percentages, ylab = "Percentage", col = "lightblue") 
text(x=seq(0.7, 5, 1.2), 2, paste("n=", counts))      # Add count to each bar
```

The barplot tells us that we have relatively more states in the South(16) and fewer states in the Northeast(9). North Central and West have a similar number of states(12 and 13).

A lollipop plot is a hybrid of a scatter plot and a barplot. Here is the lollipop chart shows the relationship between state and population.

(ref:state-pop) Loppipop plot of the population in each state.

```{r state-pop, fig.cap='(ref:state-pop)', fig.align='center'}
library(ggplot2)
ggplot(sta, aes(x = State, y = Population)) +
  geom_point(size = 3, color = "red") + 
  geom_segment(aes(x = State, xend = State, y = 0, yend = Population)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) # Rotate axis label
```

From the plot we can see that even in the early days, California and New York are the top two states in population. South Dakota had little population even in the 1970s.

Other questions we may ask are: how about the murder rate distribution in the early days? Is it the same for different states and different regions? What are the main effect factors for the murder rate? Can we use a model of other factors to explain their contribution to the murder rate? 

A choropleth map may give us an overall view.

(ref:state-map) Map of murder rate distribution.

```{r state-map, warning=FALSE, fig.cap='(ref:state-map)', fig.align='center'}
library(maps)
library(ggplot2)
sta$region <- tolower(state.name)  # Lowercase states' names
states <- map_data("state")        # Extract state data
map <- merge(states, sta, by = "region", all.x = T)  # Merge states and state.x77 data
map <- map[order(map$order), ]     # Must order first
ggplot(map, aes(x = long, y = lat, group = group)) +  
  geom_polygon(aes(fill = Murder)) +   
  geom_path() + 
  scale_fill_gradientn(colours = rev(heat.colors(10))) +
  coord_map() +
  labs(x = "Longitude", y = "Latitude") +
  guides(fill = guide_legend(title = "Murder Rate"))   
```

We can see from the map that the bottom and right of the map are close to red while the top middle and left are yellow. There is an area on top-right are yellow too, which means that murder rate are higher in south and east states but less in north central, northwest and northeast states.

(ref:state-murder) Ridgeline plot for murder rate in each region.

```{r state-murder, warning=FALSE, message=FALSE, fig.cap='(ref:state-murder)', fig.align='center'}
library(ggridges)
ggplot(sta, aes(x = Murder, y = Region, fill = Region)) +
  geom_density_ridges() +
  theme_ridges() +                                 # No color on backgroud
  theme(legend.position = "none",                  # No show legend
        axis.title.x = element_text(hjust = 0.5),  # x axis title in the center
        axis.title.y = element_text(hjust = 0.5))  # y axis title in the center
```

The ridgeline plot tells us that the murder rate skewed to the left for region west, northeast, and north central, but skewed to the right for the region south, which confirms with the map above that south has a bigger murder rate than the other regions. 

>
```{exercise}

>
Similar to Figure \@ref(fig:state-map), use choropleth map to obtain the Illiteracy map in state.x77 data set and give a brief interpretation. Hint: You can combine state.abb and state.x77 or use the row names of state.x77 data set directly. You can start from importing the data:
>
tem <- data.frame(state.x77) 
>
sta <- cbind(state.abb, tem, state.region)
>
colnames(sta)[10] <- "Region" 
>
```

>
```{exercise}

>
Similar to Figure \@ref(fig:state-murder), use ridgeline plot to explore the regional distribution of Illiteracy for state.x77 and state.region data sets and interpret your figure.
```

## Analyzing the relationship among variables

A scatter matrix is a pair-wise scatter plot of multiple variables presented in a matrix format. To visualize the linear relationship among variables in a plot, a scatter matrix is the best choice. The range of the correlation coefficient is [-1, 1]. The coefficient -1 implies two variables are strictly negatively related, such as $y=-x$. And coefficient 1 implies positive related, such as $y=2x+1$. Here is the scatter matrix for our state data.

(ref:state-corrplot) Corrplot for numeric variables.

```{r state-corrplot, warning=FALSE, message=FALSE, fig.width=6, fig.height=6, fig.cap='(ref:state-corrplot)', fig.align='center'}
st <- sta[, 2:9] # Take numeric variables as goal matrix
library(ellipse) 
library(corrplot)
corMatrix <- cor(as.matrix(st)) # Calculate correlation matrix
col <- colorRampPalette(c("red", "yellow", "blue"))  # 3 colors to represent coefficients -1 to 1.
corrplot.mixed(corMatrix, order = "AOE", lower = "number", lower.col = "black", 
               number.cex = .8, upper = "ellipse",  upper.col = col(10), 
               diag = "u", tl.pos = "lt", tl.col = "black") # Mix plots of "number" and "ellipse"
```

On the top-right of the correlation figure, we can see that the red and narrow shape between Murder and Life.Exp shows a high negative correlation, the blue narrow shape between Murder and Illiteracy shows a high positive correlation, the red-orange narrow shapes between Murder and Frost, HS.Grad show a mild negative correlation, also the orange shape between Murder and Income shows small negative correlation and light-blue shapes between Murder and both Area and Population show a small positive correlation.

The pearson and spearman correlation matrix on the bottom-left gives us the r-values between each pair of the variables, which confirm the correlation shape on the top-right.

Positive correlation between Murder and Illiteracy with an r-value of 0.70 means that the lower education level the state has, the higher chance of murder rate the state will have; Negative correlations between Murder and Life.Exp, Frost, with r-values of -0.78, and -0.54 illustrate that the more occurrence of murder, the shorter life the state will expect; And the colder the weather, the lower chance the murder will occur. It's interesting! Is it not easy for a murderer living in a cold area?

>
```{exercise}

>
Similar to Figure \@ref(fig:state-corrplot), plot a scatter matrix among 7 variables: *mpg*, *cyl*, *disp*, *hp*, *drat*, *wt* and *qsec* in the data set *mtcars*. Give a brief interpretation of the scatter matrix plot.
```

Now let's see the cluster situation of these variables.

(ref:state-dendrogram) Cluster dendrogram for state numeric variables.

```{r state-dendrogram, fig.cap='(ref:state-dendrogram)', fig.align='center'}
plot(hclust(as.dist(1 - cor(as.matrix(st)))))  # Hierarchical clustering
```

The cluster Dendrogram tells us that there are two clusters for these variables. Murder is mostly close to Illiteracy, and then to Population and Area. Similar situation, HS.Grad is mostly close to Income, and then to Life.Exp and Frost. Though illiteracy and HS.Grad are in the different clusters, we know for the same state, illiteracy is highly correlated with high school graduation rate; the lower the illiteracy, the higher the high school graduation rate. An r-value of -0.66 between Illiteracy and HS.Grad in the corrplot tells the same story.

>
```{exercise}

>
Similar to Figure \@ref(fig:state-dendrogram), plot a cluster dendrogram of the 7 variables: *mpg*, *cyl*, *disp*, *hp*, *drat*, *wt* and *qsec* in the data set *mtcars*. Give a brief interpretation of your output.
```


we can use density plot to see the distribution of Illiteracy by region.

(ref:state-illiteracy) Illiteracy distribution by region.

```{r state-illiteracy, fig.cap='(ref:state-illiteracy)', fig.align='center'}
ggplot(sta, aes(x = Illiteracy, fill = Region)) + geom_density(alpha = 0.3)
```

We can see that the north-central region has narrow density distribution with most Illiteracy less than 1 percent of the population. While the south region has an open and bell-shaped distribution with illiteracy covered from 0.5 to 3. Though region west has a spread out distribution too, it's left-skewed; There are lots of west states with illiteracy of less than 1% of its population. Most northeast region states have illiteracy less than 1.5% of their population. 

>
```{exercise}

>
Similar to Figure \@ref(fig:state-illiteracy), use density plot to see the distribution of *mpg* by *cyl* in the data set *mtcars*. 
```

Because of the relationship of Murder with both Population and Area, We add one more column of Pop.Density for the population per square miles of area to see the correlation between Murder and this density.

(ref:state-population) Box plot of population density by region.

```{r state-population, fig.cap='(ref:state-population)', fig.align='center'}
sta$Pop.Density <- sta$Population/sta$Area
boxplot(sta$Pop.Density ~ sta$Region, xlab = "Region", ylab = "Population Density")
model <- aov(sta$Pop.Density ~ sta$Region, sta)
summary(model)
```

The box plot shows that the mean Pop.Density of Northeast is much more than the other regions, while West has the lowest mean Pop.Density. ANOVA test with a p-value of 6.3e-06 also give us the evidence to reject the null hypothesis that the mean Pop.Densities are the same for different regions, which means at least one of the regional mean population densities is different from the others.

Here is the scatterplot for Illiteracy and Murder with Population per area.

(ref:state-IlliteracyMurder) Scatterplot for illiteracy and murder sized by population density and colored by region.

```{r state-IlliteracyMurder, fig.cap='(ref:state-IlliteracyMurder)', fig.align='center'}
ggplot(sta, aes(x = Illiteracy, y = Murder)) + 
  geom_point(aes(size = Pop.Density, color = Region)) + 
  geom_smooth(method = 'lm',formula = y ~ x)  # Add regression line
```

The plot shows that murder and illiteracy are positively correlated. All states in the other three regions have murder rates of less than 12 per 100,000 population except some of the south states(green). All north-central states(blue) have illiteracy less than 1, all northeast states(red) have less than 1.5 of illiteracy. The illiteracy of west(purple) and south states have much bigger variance, some western states have illiteracy less than 0.5, while some south states have illiteracy more than 2.75. Many Northeast states have big population densities but middle illiteracy rates compared with the states in the other three regions.

Because of the high correlation of murder and Life.Exp, we will take a look at the distribution of Life.Exp.

(ref:state-LifeExp) Regional life expectancy.

```{r state-LifeExp, fig.cap='(ref:state-LifeExp)', fig.align='center'}
ggplot(sta, aes(x = Region, y = Life.Exp, fill = Region)) + 
  geom_violin(trim = FALSE) + 
  geom_boxplot(width = 0.1)
```

On average, the south has a lower life expectancy than the other three regions. North Central has the highest Life.Exp. West has extending distribution with two long tails on each end, which means some west states have very long life expectancy, while some states expect short life expectancy though they are in the same region.

According to the corrplot, there are other variables in the other cluster, such as life expectancy, high school graduation and income. We are interested in whether they affect the murder rate. Here is the plot for relationship between murder and these variables. Firstly we group the income into 5 IncomeTypes.

(ref:state-LifeMurder) Relationship between murder rate and life expectancy, high school graduation, and income.

```{r state-LifeMurder, warning=FALSE, message=FALSE, fig.cap='(ref:state-LifeMurder)', fig.align='center'}
library(dplyr)
# group income into IncomeType first
sta.income <- sta %>% mutate(IncomeType = factor(ifelse(Income < 3500, "Under3500",  
                             ifelse(Income < 4000 & Income >= 3500, "3500-4000",
                               ifelse(Income < 4500 & Income >= 4000, "4000-4500",
                                 ifelse(Income < 5000 & Income >= 4500, "4500-5000", 
                                   "Above5000"))))))
ggplot(sta.income, aes(x = Murder, y = Life.Exp)) + 
  geom_point(aes(shape = IncomeType, color = Region, size = HS.Grad)) + 
  geom_smooth(method = 'lm', formula = y ~ x)
```

Murder is negatively correlated with Life.Exp. Some states with higher murder rates over 12 have relatively small symbols, which means their high school graduation rates are close to 40%; And these small symbols with murder rates bigger than 12 are all colored as green, which means they all belong to the south region. 

It looks like the income type does not affect the murder rate a lot because all different symbols scatter around in different murder rates, especially between murder rates 8 and 10.

Most southern states have lower HS.Grad, lower Life.Exp but higher murder frequency, while states in the other three regions have relatively higher high school graduation rate and income but lower murder rate.

>
```{exercise}

>
Use a scatter plot to analyze the correlation between Illiteracy and those variables in the other cluster shown in Figure \@ref(fig:state-dendrogram). Interpret your plot.
```

## The whole picture of the data set 

We analyzed the relationship of murder rate with the variables in both clusters. It looks like all variables are correlated, more or less. Now let's see the whole picture of these variables. Heat map and Segment diagram are the popular methods.

(ref:state-heatmap) Heat map for whole state data set.

```{r state-heatmap, message=FALSE, warning=FALSE, fig.width=8, fig.height=10, fig.cap='(ref:state-heatmap)', fig.align='center'}
library(gplots)
st.matrix <- as.matrix(st)    # Transfer the data frame to matrix
s <- apply(st.matrix, 2, function(y)(y - mean(y)) / sd(y))  # Standardize data
a <- heatmap.2(s, 
          col = greenred(75), # Color green red
          density.info = "none", 
          trace = "none", 
          scale = "none", 
          RowSideColors = rainbow(4)[sta$Region],  
          srtCol = 45,        # Column labels at 45 degree
          margins = c(5, 8),  # Bottom and right margins
          lhei = c(5, 15)     # Relative heights of the rows
) 
legend("topright", levels(sta$Region), fill = rainbow(4), cex = 0.8)  # Add legend 
```

Same as cluster Dendrogram plot, Heatmap shows that Life.Exp, Income, HS.Grad, together with Frost build one cluster, while Illiteracy, Murder, Population, and area build another cluster.

Compared with other states, lots of southern states with lower Life.Exp, Income, and HS.Grad have higher Murder rates and Illiteracy, like Mississippi and Alabama. On the contrary, some northern and western states which have higher Life.Exp, Income, and HS.Grad show lower Area, Population, Murder, and Illiteracy, like Nebraska and South Dakota. Though the income of South Dakota show a little bit green.

(ref:state-segment) Segment diagram for all states.

```{r state-segment, fig.width=8, fig.height=8, fig.cap='(ref:state-segment)', fig.align='center'}
row.names(st) <- sta$State
stars(st, key.loc = c(13, 1.5), draw.segments = T)   
```

The segment Diagram shows us different aspects of each state. For example, South Dakota has big Frost(yellow), big Life Expectancy(blue), a relatively high percentage of high school graduation rate(pink) and good income(red), but has a small area, and very tiny population, illiteracy, and murder rate compared with the other states.

We use principal components analysis to explore the data a little bit more!
```{r}
pca = prcomp(st, scale = T)  # scale = T to normalize the data
pca
plot(pca)     # Plot the amount of variance each principal components captures
summary(pca)  # Shows the importance of the components
percentVar <- round(100 * summary(pca)$importance[2, 1:7], 0) # Compute % variances
percentVar
```

The first two components contribute 45% and 20%, together 65% of the variance. The third component is a little bit less but still over 10% of the variance. The bar plot of each component's variance shows how the each component dominate. 

(ref:state-bioplot) Biplot for PCA. 

```{r state-bioplot, warning=FALSE, fig.cap='(ref:state-bioplot)', fig.align='center'}
library(ggfortify)
row.names(sta) <- sta$State
autoplot(prcomp(st,  scale = T), data = sta, 
         colour = 'Region', shape = FALSE, label = TRUE, label.size = 3.5,
         loadings = TRUE, loadings.colour = 'blue', loadings.label = TRUE, 
         loadings.label.size = 4, loadings.label.colour = 'blue')
```

The Biplot illustrates the special role of these variables to the first and second components of the variance. Illiteracy positively contributes to the component of the variance PC1, while Life.Exp and Frost negatively contribute to the component of the variance PC1. Area positively contributes to the component of the variance PC2. The other four variables contribute to both components of the variance PC1 and PC2 positively or negatively. 
From the figure, we also find that many states in the south region, such as Louisiana(LA) and Mississippi(MS) are well-known for their big Illiteracy and high murder rate, while some north-central states like Minnesota(MN) and North Dakota(ND) are noticed for their high life expectancy and long frost time. Big area is the special feature for two states in the west region, Alaska(AK) and California(CA). 

## Linear model analysis
According to the analysis above, we try to find a model to explain the murder rate. Because of the high correlation of HS.Grad with Illiteracy, Life.Exp, and Income, we will not put HS.Grad in the model. For a similar reason, we leave Frost out too.
```{r}
lm.data <- sta[, c(2:6, 9:10)]
lm.data <- within(lm.data, Region <- relevel(Region, ref = "South"))  # Set region South as reference
model <- lm(Murder ~ .,  data = lm.data)
summary(model)
```

Murder rate is most related to Life Expectancy and Population of the state, also it is affected by Illiteracy of the state. The region is another factor contributing to murder rate. The estimates illustrate that every unit of increase in Life.Exp expects 1.445 units lower of murder rate, while every unit of increase in population and illiteracy will increase 0.000259 and 1.861 units of the murder rate. At the same time, if the state belongs to the northeast region, the murder rate will be 2.673 units less. The model will explain 82% of the variance of the murder rate. If we know the population, Life.Exp, Illiteracy of the certain state in those years, we can estimate murder rate as follow:
$Murder = 105.9 - 1.445 * Life.Exp + 0.000259 * Population + 1.861 * Illiteracy - 2.673 * RegionNortheast$ 

>
```{exercise}

>
Do a linear model analysis for Illiteracy and interpret your result. Hint: Check the corrplot figure \@ref(fig:state-corrplot) and pay attention to the high correlation between murder rate and life expectancy.
```

## Conclusion

-Southern region shows a higher murder rate with lower life expectancy, income, and high school graduation rate but higher illiteracy, while northern region shows a lower murder rate with higher population density, life expectancy, income, and high school graduation rate but lower illiteracy. 

-The information of life expectancy, population, illiteracy of the state in the 1970s, and whether the state belongs to the northeast region will help to estimate the murder rate of the state at that time.