STAT270 - Final Project

library(ggplot2)
library(tidyverse)
library(graphics)
library(mdsr) # install package if not installed
library(discrim) # install package if not installed
library(klaR) # install package if not installed
library(kknn) # install package if not installed
library(utils) # install package if not installed
library(sp) # install package if not installed
library(fs)

Note: If you Rmd file submission knits you will receive total of (5 points Extra Credit)

Directions: Complete Task 1 and one of the Task 2 or 3!

Task 1 (Total 60 pts):

High-earners in the 1994 United States Census

A marketing analyst might be interested in finding factors that can be used to predict whether a potential customer is a high-earner. The 1994 United States Census provides information that can inform such a model, with records from 32,561 adults that include a binary variable indicating whether each person makes greater or less than $50,000 (more than $80,000 today after accounting for inflation). This is our response variable.

A bit of data preparation

library(tidyverse)
library(mdsr)
url <-
"http://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data"

census <- read_csv(
  url,
  col_names = c(
    "age", "workclass", "fnlwgt", "education", 
    "education_1", "marital_status", "occupation", "relationship", 
    "race", "sex", "capital_gain", "capital_loss", "hours_per_week", 
    "native_country", "income"
  )
) %>%
  mutate(income = factor(income), income_ind = as.numeric(income == ">50K")) # create indicator variable income_ind (0 - low, 1 - high earner)

Rows: 32561 Columns: 15── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (9): workclass, education, marital_status, occupation, relationship, race, sex, native_country, income
dbl (6): age, fnlwgt, education_1, capital_gain, capital_loss, hours_per_week
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# look at the structure of the data
glimpse(census)

Rows: 32,561
Columns: 16
$ age            <dbl> 39, 50, 38, 53, 28, 37, 49, 52, 31, 42, 37, 30, 23, 32, 40, 34, 25, 32, 38, 43, 40, 54, 35, 43, 59, 56, 19, 54, 39…
$ workclass      <chr> "State-gov", "Self-emp-not-inc", "Private", "Private", "Private", "Private", "Private", "Self-emp-not-inc", "Priva…
$ fnlwgt         <dbl> 77516, 83311, 215646, 234721, 338409, 284582, 160187, 209642, 45781, 159449, 280464, 141297, 122272, 205019, 12177…
$ education      <chr> "Bachelors", "Bachelors", "HS-grad", "11th", "Bachelors", "Masters", "9th", "HS-grad", "Masters", "Bachelors", "So…
$ education_1    <dbl> 13, 13, 9, 7, 13, 14, 5, 9, 14, 13, 10, 13, 13, 12, 11, 4, 9, 9, 7, 14, 16, 9, 5, 7, 9, 13, 9, 10, 9, 9, 12, 10, 1…
$ marital_status <chr> "Never-married", "Married-civ-spouse", "Divorced", "Married-civ-spouse", "Married-civ-spouse", "Married-civ-spouse…
$ occupation     <chr> "Adm-clerical", "Exec-managerial", "Handlers-cleaners", "Handlers-cleaners", "Prof-specialty", "Exec-managerial", …
$ relationship   <chr> "Not-in-family", "Husband", "Not-in-family", "Husband", "Wife", "Wife", "Not-in-family", "Husband", "Not-in-family…
$ race           <chr> "White", "White", "White", "Black", "Black", "White", "Black", "White", "White", "White", "Black", "Asian-Pac-Isla…
$ sex            <chr> "Male", "Male", "Male", "Male", "Female", "Female", "Female", "Male", "Female", "Male", "Male", "Male", "Female", …
$ capital_gain   <dbl> 2174, 0, 0, 0, 0, 0, 0, 0, 14084, 5178, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ capital_loss   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2042, 0, 0, 0, 0, 0, 0, 0, 0, 1408, 0, 0, 0, …
$ hours_per_week <dbl> 40, 13, 40, 40, 40, 40, 16, 45, 50, 40, 80, 40, 30, 50, 40, 45, 35, 40, 50, 45, 60, 20, 40, 40, 40, 40, 40, 60, 80…
$ native_country <chr> "United-States", "United-States", "United-States", "United-States", "Cuba", "United-States", "Jamaica", "United-St…
$ income         <fct> <=50K, <=50K, <=50K, <=50K, <=50K, <=50K, <=50K, >50K, >50K, >50K, >50K, >50K, <=50K, <=50K, >50K, <=50K, <=50K, <…
$ income_ind     <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …

ggplot(census, aes(x = education, fill = income)) +
  geom_bar(position = "fill") +
  coord_flip() +
  labs(y = "Proportion", title = "Proportion of High vs Low Earners by Education")

a) (10 pts) Split the data set into two pieces by separating the rows at random. A sample of 70% of the rows will become the training data set, with the remaining 30% set aside as the testing (or “hold-out”) data set. Use set.seed(364) in the beginning of your code. How many records are in the testing set? 9768 Hint: One possible way to do this is to use the function initial_split(prop = 0.7) and call the functions training() and testing().

library(tidymodels)
set.seed(364)

# Perform the split
census_split <- initial_split(census, prop = 0.7)
train <- training(census_split)
test <- testing(census_split)

# Check how many rows are in the test set
nrow(test)

[1] 9769

Note: You should get around 24% of those in the sample make more than $50k. Thus, the accuracy of the null model is about 76%, since we can get that many right by just predicting that everyone makes less than $50k.

# if your training set is called `train` this code will produce the correct percentage
pi_bar <- train %>%
  count(income) %>%
  mutate(pct = n / sum(n)) %>%
  filter(income == ">50K") %>%
  pull(pct)

print(c("Percent >50K", pi_bar)) 

[1] "Percent >50K"      "0.241751491751492"

Pro Tip: Always benchmark your predictive models against a reasonable null model.

b) (10 pts) Use KNN algorithm to classify the earners (<=50K, >50K, low/high) for High-earners in the 1994 United States Census data above. Select only the quantitative variables age,education_1, capital_gain, capital_loss, hours_per_week. Use mode = "classification" and k=1(use the closest neighbor) in the nearest_neighbor function arguments. Print the confusion matrix. State the accuracy.

Hint: See Programming exercises Week 13 for details about KNN implementation.

library(kknn)

# Select only the numeric columns (excluding 'fnlwgt') and 'income' as the target variable
train_q <- train %>% dplyr::select(income, where(is.numeric), -fnlwgt)

# Define the KNN classifier (k=1, classification task)
knn_model <- kknn(income ~ age + education_1 + capital_gain + capital_loss + hours_per_week, 
                  train = train_q, 
                  test = train_q, 
                  k = 1, 
                  distance = 2, # Euclidean distance by default
                  kernel = "rectangular") # Rectangular kernel (default)

Warning: variable 'income' is absent, its contrast will be ignored

# Predict the income using the KNN classifier
train_q$income_knn <- predict(knn_model)

# Print the confusion matrix
confusion_matrix <- table(train_q$income, train_q$income_knn)
print(confusion_matrix)

       
        <=50K  >50K
  <=50K 15426  1856
  >50K   1627  3883

# Calculate accuracy
accuracy <- sum(diag(confusion_matrix)) / sum(confusion_matrix)
print(paste("Accuracy is:", round(accuracy, 4)))

[1] "Accuracy is: 0.8472"

Accuracy is:0.8472

# Find the Accuracy = (true positive and true negative)/total or use the `accuracy()` function.

Load necessary libraries

library(tidyverse)

Define the formula for the logistic regression model

form <- as.formula( “income_ind ~ age + education_num + sex + marital_status” )

Fit logistic regression model using glm function

logit_model <- glm(form, data = train, family = binomial)

Summary of the logistic regression model to see coefficients

summary(logit_model)

Predict probabilities for the test set

logit_pred_prob <- predict(logit_model, newdata = test, type = “response”)

Assign 1/0 based on a threshold of 0.5 for high-earner prediction

pred_y_logit <- as.numeric(logit_pred_prob > 0.5)

Create confusion matrix comparing predicted values and actual values

confusion_logit <- table(pred_y_logit, test$income_ind)

Print confusion matrix

print(confusion_logit)

Calculate accuracy

accuracy_logit <- (confusion_logit[1, 1] + confusion_logit[2, 2]) / sum(confusion_logit) print(paste(“Accuracy is:”, accuracy_logit))

`# Fit the logistic regression model with all chosen predictors form <- as.formula( “income_ind ~ age + education_1 + sex + marital_status” )

logit_model <- glm(form, data = train, family = binomial)

Summarize the logistic regression model to see which predictors are significant

summary(logit_model)

Fit the logistic regression model

logreg <- glm(income_ind ~ age + sex + education_1, family = “binomial”, data = train)

Load broom package to use tidy function

library(broom)

Summarize the model using tidy()

tidy_logreg <- tidy(logreg)

Print the tidy output to see the coefficients, standard errors, and p-values

print(tidy_logreg)

**FYI:** The predicted probabilities and predicted values for `income_ind` can be found using the code, uncomment the code lines to use (highlight chink and press **CTRL + SHIFT + C**):

``# Predict probabilities using the logreg model on the test set
logit_pred_prob_test <- predict(logreg, newdata = test, type = "response")

# Convert predicted probabilities to class predictions (1 for high earner, 0 for low earner)
pred_y_test <- as.numeric(logit_pred_prob_test > 0.5)

# Create confusion matrix comparing predicted and actual values from the test set
confusion_test <- table(pred_y_test, test$income_ind)

# Print confusion matrix
print(confusion_test)

# Calculate accuracy by comparing predicted values with actual values from the test set
accuracy_test <- mean(pred_y_test == test$income_ind, na.rm = TRUE)

# Print accuracy
print(paste("Accuracy on the test set is:", accuracy_test))

e) (10 pts) Assessing the Logit model from part d) using the test set saved in test R-object.

What is the accuracy of the model?

# Fit logistic regression model using the training set
logreg <- glm(income_ind ~ age + sex + education_1, family = "binomial", data = train)

# Display model summary
summary(logreg)

Call:
glm(formula = income_ind ~ age + sex + education_1, family = "binomial", 
    data = train)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)    
(Intercept) -7.674689   0.116165  -66.07   <2e-16 ***
age          0.042087   0.001344   31.31   <2e-16 ***
sexMale      1.304009   0.043565   29.93   <2e-16 ***
education_1  0.365419   0.007782   46.96   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 25212  on 22791  degrees of freedom
Residual deviance: 20234  on 22788  degrees of freedom
AIC: 20242

Number of Fisher Scoring iterations: 5

f) (10 pts) Which one of the classification models achieved the highest accuracy? The KNN model achieved the highest accuracy of 84.72% ## TASK 2 (Total 40 pts, 20 pts each part a & b below):

Let us consider the unsupervised learning process of identifying different types of cars. The United States Department of Energy maintains automobile characteristics for thousands of cars: miles per gallon, engine size, number of cylinders, number of gears, etc.

Here, we use the 2016 FEGuide.xlsx file that contains fuel economy rating for the 2016 model year.

Next, we will use the readxl package to read this file into R, clean up some of the resulting variable names, select a small subset of the variables, and filter for distinct models of Toyota vehicles. The resulting data set contains information about 75 different models that Toyota produces.

Store the data file “2016 FEGuide.xlsx” in a subfolder by the name data in your working directory.

Note: You may need to adjust the code below to specify the “2016 FEGuide.xlsx” file location if you opt out to store the data file in different location.

# Load necessary libraries
library(readxl)
library(dplyr)
library(janitor)
library(ape)

Warning: package ‘ape’ was built under R version 4.4.3
Attaching package: ‘ape’

The following object is masked from ‘package:rsample’:

    complement

The following object is masked from ‘package:dials’:

    degree

The following object is masked from ‘package:dplyr’:

    where

# Load and clean data for Honda
filename <- "data/2016 FEGuide.xlsx"

# Read the data and clean column names
cars <- read_excel(filename) %>% 
  clean_names()  # Clean column names

Error: `path` does not exist: ‘data/2016 FEGuide.xlsx’

As a large automaker, Toyota has a diverse lineup of cars, trucks, SUVs, and hybrid vehicles. Can we use unsupervised learning to categorize these vehicles in a sensible way with only the data we have been given?

For an individual quantitative variable, it is easy to measure how far apart any two cars are: Take the difference between the numerical values. The different variables are, however, on different scales and in different units. For example, gears ranges only from 1 to 8, while city_mpg goes from 13 to 58. This means that some decision needs to be made about rescaling the variables so that the differences along each variable reasonably reflect how different the respective cars are. There is more than one way to do this, and in fact, there is no universally “best” solution—the best solution will always depend on the data and your domain expertise. The dist() function takes a simple and pragmatic point of view: Each variable is equally important.

The output of dist() gives the distance from each individual car to every other car.

Create distance matrix object from the car_diffs

Choose one of the car makers: General Motors, Nissan, Ford Motor Company, Honda, Mercedes-Benz, BMW, Kia - preferably maker that you are familiar with the models but not necessarily.

a) (Total 20 pts) Create a tree constructed by hierarchical clustering that relates carmaker car models to one another.

Hint: You can use the code above and make necessary modification.

YOUR CODE HERE:

# Load necessary libraries
library(readxl)
library(dplyr)
library(janitor)
library(ape)

# Load and clean data for Honda
filename <- "data/2016 FEGuide.xlsx"

cars <- read_excel(filename) %>% 
  clean_names() %>%
  select(
    make = mfr_name,
    model = carline,
    displacement = eng_displ,
    number_cyl = cylinders,
    number_gears = gears,
    city_mpg = city_fe_guide_conventional_fuel,
    hwy_mpg = hwy_fe_guide_conventional_fuel
  ) %>%
  distinct(model, .keep_all = TRUE) %>%
  filter(make == "Honda")  # Filtering for Honda only

# View the cleaned data
glimpse(cars)

# Scaling the data to normalize the features
cars_scaled <- cars %>%
  select(displacement, number_cyl, number_gears, city_mpg, hwy_mpg) %>%
  scale()

# Compute the distance matrix
car_diffs <- dist(cars_scaled)

# Create hierarchical clustering
hc <- hclust(car_diffs)

# Plot the dendrogram (tree)
plot(hc, main = "Hierarchical Clustering of Honda Car Models", xlab = "Car Models", sub = "", cex = 0.6, labels = cars$model)

b) (Total 20 pts) Attempt to interpret the tree, how the models in same cluster are similar and how clusters differ.

YOUR COMMENTS:

TASK 3 (Total 40 pts, 20 pts each part a & b below) K-means clustering (Geospatial data example)

Another way to group similar cases is to assign each case to one of several distinct groups, but without constructing a hierarchy. The output is not a tree but a choice of group to which each case belongs. (There can be more detail than this; for instance, a probability for each group that a specific case belongs to the group.) This is like classification except that here there is no response variable.

Geospatial data example:

Consider the cities of the world (in WorldCities, in the mdsr package). Cities can be different and similar in many ways: population, age structure, public transportation and roads, building space per person, etc. The choice of features (or variables) depends on the purpose you have for making the grouping.

Our purpose is to show you that clustering via machine learning can actually identify genuine patterns in the data.

We will choose features that are utterly familiar: the latitude and longitude of each city.

You already know about the location of cities. They are on land. And you know about the organization of land on earth: most land falls in one of the large clusters called continents.

But the WorldCities data doesn’t have any notion of continents. Perhaps it is possible that this feature, which you long ago internalized, can be learned by a computer that has never even taken grade-school geography.

Consider the 4,000 biggest cities in the world and their longitudes and latitudes.

BigCities <- world_cities %>% arrange(desc(population)) %>% 
  head(4000) %>% 
  dplyr::select(longitude, latitude)
glimpse(BigCities)

Note that in these data, there is no ancillary information—not even the name of the city. However, the k-means clustering algorithm will separate these 4,000 points—each of which is located in a two-dimensional plane—into k clusters based on their locations alone.

set.seed(15)
# install the package first if not installed
#install.packages("mclust")
library(mclust) 

# form 6 cluster iteratively
city_clusts <- BigCities %>%
kmeans(centers = 6) %>% fitted("classes") %>% as.character()

# form 6 cluster iteratively, by forming initially 10 random sets
km <- kmeans(BigCities, centers = 6, nstart = 10)
# inspect the structure of the kmeans output cluster object
str(km)

# access two important features of cluster, their size and centers
km$size
km$centers

BigCities <- BigCities %>% mutate(cluster = city_clusts) 

# graph the clusters, using the cluster variable to pick the color in standard cartesian coordinate system
BigCities %>% ggplot(aes(x = longitude, y = latitude)) +
geom_point(aes(color = cluster), alpha = 0.5)

a) Total (10 pts) What did the above clustering algorithm seems to have identified?

Answer:

b) Total (30 pts)

Projections: The Earth happens to be an oblate spheroid—a three-dimensional flattened sphere. Yet we would like to create two-dimensional representations of the Earth that fit on pages or computer screens. The process of converting locations in a three-dimensional geographic coordinate system to a two-dimensional representation is called projection.

A coordinate reference system (CRS) is needed to keep track of geographic locations. Every spatially-aware object in R can have a projection string, encoded using the PROJ.4 map projection library. These can be retrieved (or set) using the proj4string() command.

There are many CRSs, but a few are most common. A set of EPSG (European Petroleum Survey Group) codes provides a shorthand for the full PROJ.4 strings (like the one shown above). The most commonly-used are:

EPSG:4326 Also known as WGS84, this is the standard for GPS systems and Google Earth.

EPSG:3857 A Mercator projection used in maps tiles4 by Google Maps, Open Street Maps, etc.

EPSG:4269 NAD83, most commonly used by U.S. federal agencies.

R-Code below uses EPSG:4326. Use the other two standards EPSG:3857, EPSG:4269

Use K-means algorithm for each of the three (3) projections above and compare the three projections to the standard cartesian coordinates used in the example. Which one is best in identifying the continents?

Note: Graphing the clusters in each projection is worth 10 pts

# assign the BigCities data.frame to a working data.frame object d 
d <- BigCities #or BigCities[,c('longitude', 'latitude')]

# create spatial object from d
coordinates(d) <- 1:2

# Set WGS 84 (EPSG:4326) standard for projecting longitude latitude coordinates
proj4string(d) <- CRS("+init=epsg:4326")

# coordinate reference system using the EPSG:4326 standard
CRS.new <- CRS("+init=epsg:4326")

# the d object in the new CRS, you may print out few records to see how it looks in the new CRS
d.new <- spTransform(d, CRS.new)

# just for information review the 
proj4string(d.new) %>% strwrap()


# form 6 cluster iteratively
city_clusts <- as.data.frame(d.new) %>%
kmeans(centers = 6) %>% fitted("classes") %>% as.character()

# add a variable for the newly formed clusters
df.new <- as.data.frame(d.new) %>% mutate(cluster = city_clusts, longitude = coords.x1, latitude = coords.x2)

# graph the clusters, using the cluster variable to pick the color
df.new %>% ggplot(aes(x = coords.x1, y = coords.x2)) +
geom_point(aes(color = cluster), alpha = 0.5) +
scale_color_brewer(palette = "Set3")

YOUR CODE HERE: