Homework 1

2a (5 points)

Write a function called even.or.odd. Its parameter x will be a numeric vector. Return a character vector that says “odd” for odd numbers and “even” for the even numbers. The results should correctly classify every value. To determine if a number is even or odd, you can use the modulus operator %% (e.g.: 5%%3 = 2). Note: Try to find a solution that uses vector logic instead of a for loop. In R, this is a good programming practice that will speed up your programs.

Display the results of this function on the vector 1:5.

x <- c ( 1 , 4 , 5 , 7 , 9)
even.or.odd <- function (x) {
  ifelse ( x %% 2 == 0, 'even','odd')
}
even.or.odd(x)

[1] "odd"  "even" "odd"  "odd"  "odd" 

2b (5 points)

Write a function my.sum that computes the sum of a numeric vector x. The user can also specify the.type as a character. If the.type is “even”, the function will compute the sum of only the even values in x. If the.type is “odd”, the function will compute the sum of only the odd values in x. If the.type is “all”, the function will compute the sum of the entire vector x. Within the function, you may use the built-in sum function. The function should omit missing values (NA) from the sum. This can be done using the na.rm argument within the sum function.

Display the results of this function for odd, even, and all values of the vector 1:5.

my.sum <- function(x, the.type) {
  if(the.type == "even"){
    sum(x[(x %%2 == 0)], na.rm = TRUE)
  }
  else if(the.type == "odd"){
    sum(x[(x %%2 == 1)], na.rm = TRUE)
  }
  else if(the.type == "all"){
    sum(x, na.rm = TRUE)
  }
}

my.sum(1:5, the.type ="even")

[1] 6

3a (5 points)

Divide the data into training and testing sets. To do so, let’s create an assignment vector called training_row. Each row of the data set will be assigned to the training set (with training_row set to TRUE) with probability 0.8 or to the test set (with training_row set to FALSE) with probability 0.2. Use the sample function to create the training_row vector of TRUE and FALSE values. The vector should be as long as the number of rows in the iris data set.

Then, divide the iris data set into separate training and test sets according to the training_row assignments.

In order to obtain consistent results, we’ll need to set the seed of R’s pseudo-random number generator. To do so, use set.seed(41) in the code chunk labeled constants above.

library(datasets)
set.seed(41)
##create normalization function 
##nor <-function(x) { (x - min(x)) / (max(x) - min(x))  }
#run normalization function on first 4 columns of df
##df_norm <- as.data.frame (lapply (df [ , c(1,2,3,4) ], nor))

#extract training and testing set
training_row <- sample(1:nrow (iris) , 0.8*nrow (iris))
train <- iris[training_row, ]
test <- iris[- training_row, ]

3b (5 points)

Use the function knn from the package class with k = 2 to classify the data. What proportion of the values are misclassified on the testing set?

Note: In order to use knn, the train and test objects must only include the columns that are used to make the classification. The Species will need to be separated into the cl vector and removed from the train and test objects.

library(class)
#seperate variable of interest into cl vectors
train_labels <- train[, 5] 
test_labels <- test[, 5]  
#predict classication using train_lables for test set
pred <- knn (train[,-5], test[,-5], cl = train_labels, k = 2)
#calculate misclassification rate
mean(test_labels != pred)

[1] 0.1

#alternativley use: sum(test_labels!=pred)/sum(nrow(test)) or use confusion matrix
conf <- table(test_labels,pred)
conf

            pred
test_labels  setosa versicolor virginica
  setosa          9          0         0
  versicolor      0          7         2
  virginica       0          1        11

1-sum(diag(conf))/sum(conf)

[1] 0.1

5a (5 points)

Plot the misclassification rates on the testing set versus the value of k. Use the plot function. Try different values of the arguments (las, xlim, ylim, xlab, ylab, cex, main) to create a nicer display. Use type = “both” to display both the points and a line.

# Plot the misclassification
plot(mis_df$k, mis_df$mis_rate, main = "Misclassification Rates vs Neighbors",  xlab = "k, number of neighbors", ylab = "Misclassification Rates", cex =0.5, type = "both")

5b (5 points)

Now create the same plot placing k on a *logarithmic** scale. Make sure to change the label of the x axis to distinguish this.

plot(mis_df$k, mis_df$mis_rate, log = "x", main = "Misclassification Rates vs Neighbors",  xlab = "K in log-scale", ylab="Misclassification Rate",cex = 0.5, type = "both")

5c (10 points)

Let’s examine how the results would change if we were to run the knn classifier multiple times. Perform the following steps:

1. Re-perform the previous work 3 more times. Each time, you should create a new training and test set, apply knn on each value of k from 1 to 50, and compute the misclassification rates on the testing set.

2. Plot the results of the earlier work along with the 3 new iterations on a single plot. Use the lines function to add additional lines to the earlier plot from 5a (using the linear scale). Use different colors, line types (lty), and point characters (pch) to distinguish the lines.

3. Use the legend command to place a legend in the top left corner (x = “topleft”) of the plot. Use the same colors and point characters to display which line is which. Label the iterations 1 through 4.

#create function knn_output(k_values)
knn_output <- function(k_values = 1:50){
  training_row = sample(1:nrow (iris) , 0.8*nrow (iris))
  train = iris[ training_row, ]
  test = iris[- training_row, ]
  train_labels <- train[, 5] 
  test_labels <- test[, 5]  
  k_values <- 1 : 50
  mis_df <- data.frame( k = rep( 0 , length(k_values)), mis_rate = rep (0 , length(k_values)))
  for (i in k_values){
    k <- k_values[i]
    pred.i <- knn(train[,-5], test[,-5], train_labels, k = i)
    mis.i <- mean(test_labels != pred.i)
    mis_df[i,'k']<- k
    mis_df[i,'mis_rate']<- mis.i
  }
  mis_df
}

#run another 3 iterations
knn2 <-knn_output(k_values = 1:50)
knn3 <-knn_output(k_values = 1:50)
knn4 <-knn_output(k_values = 1:50)

#create plot
plot(mis_df$k, mis_df$mis_rate, main = "Misclassification Rates vs Neighbors", xlab = "k, number of neighbors", ylab = "Misclassification Rates", col = "black",lty = 1, pch = 1,type = "b",ylim = c(0,.3))
#add lines of knn2, knn3, knn4
lines(mis_rate ~ k, knn2, col="blue" , lty=2,pch = 2,type = "b",ylim=c(0,.3))
lines(mis_rate ~ k, knn3, col ="red", lty=3,pch = 3,type = "b",ylim=c(0,.3))
lines(mis_rate ~ k, knn4, col ="green", lty=4,pch = 4,type = "b",ylim=c(0,.3))
#add legend
legend("topleft", 
  legend = c('iteration 1','iteration 2','iteration 3','iteration 4'), 
  col = c("black","blue","red","green"), 
  bty = "n", 
  cex = 0.8, 
  text.col = c("black","blue","red","green"), 
  horiz = F,
  lty=c(1,2,3,4),
  pch = c(1,2,3,4))

#My first try:
#knn_mis_rate <- function(k){
#  training_row = sample(1:nrow (iris) , 0.8*nrow (iris))
#  train = iris[ training_row, ]
#  test = iris[- training_row, ]
#  train_labels <- train[, 5] 
#  test_labels <- test[, 5]  
# pred <- knn (train[,-5], test[,-5], train_labels, k)
#  #calculate mis_rate
#  mis_rate <- mean(pred != test_labels)
#  return(mis_rate)
#}

#knn_output <- function(k_values = 1:50){
#  k_values = 1:50
#  mis_rate = c()
#  for (k in k_values){ 
#    mis_rate[k] = knn_mis_rate(k)
#}
#  cbind.data.frame(k_values, mis_rate)
#}

6a (2 points)

What are the dimensions of the data set?

library(ISLR)
dim(Hitters)

[1] 322  20

6b (2 points)

How many salaries are missing (NA)?

sum(is.na(Hitters$Salary))

[1] 59

6c (2 points)

What is the maximum number of career home runs?

max(Hitters$CHmRun)

[1] 548

6d (2 points)

Compute the min, median, mean, and max of Hits, Home Runs, and Runs for a season (not career totals). Remove any missing values from the calculations. Round your results to 1 decimal place.

cal <- function(x) {
  c(min = round(min(x, na.rm = TRUE),1), 
    median = round(median(x,na.rm = TRUE),1),
    mean = round(mean(x,na.rm = TRUE),1),
    max = round(max(x,na.rm = TRUE),1)
  )
}

sapply(Hitters[,c("Hits","HmRun","Runs")], cal)

       Hits HmRun  Runs
min       1   0.0   0.0
median   96   8.0  48.0
mean    101  10.8  50.9
max     238  40.0 130.0

6e (2 points)

What percentage of these players’s seasons had at least 100 hits and 20 home runs? Use the percent function in the scales package to convert a decimal proportion to a percentage.

library('scales')
n <- nrow(subset(Hitters,Hits >= 100 & HmRun >= 20))
percent(n/nrow(Hitters))  

[1] "16.8%"

6f (2 points)

What is the relationship between different pairs of variables? Let’s look at Salary, Hits, Runs, HmRun, Errors, and Assists. Use the pairs function to display scatterplots of each pair of these variables.

pairs(Hitters[,c("Salary","Hits","Runs","HmRun","Errors","Assists")], lower.panel = NULL )

6g (2 points)

Based on these scatterplots, which variables appear to be correlated with Salary, and which ones appear to have little or no correlation with Salary? Provide a short explanation for your assessment.

Based on these scatterplots, Hits, Runs, HmRun appear to be positively correlated with Salary because the there’re clear positive associate for Salary-Hits, Salary-Runs, Salary-HmRun, which represent moderate linear relationships. Errors and Assists appear to have little or no correlation with Salary as the points are dispersed without clear cirection and form. To quantify the correlations, a correlation matrix can be applied using the function: cor(df[,c(“Salary”,“Hits”,“Runs”,“HmRun”,“Errors”,“Assists”)], use = “complete.obs”)

6h (2 points)

Create a new variable called HighRBI for those players with at least 75 RBI (TRUE). Players with less than 75 RBI should have the value FALSE.

library(dplyr)
H <- mutate(Hitters , HighRBI = ifelse ( RBI >= 75, 'TRUE', "FALSE"))
H$HighRBI <- as.logical(H$HighRBI)

#H <- Hitters %>% 
#  mutate(HighRBI = as.factor(case_when(
#  RBI >=75 ~ 'TRUE',
#  RBI <75 ~ 'FALSE'
#  )
#  )
#)

6i (2 points)

What percentage of hitters qualified as HighRBI during these seasons?

percent(mean(H$HighRBI == TRUE))

[1] "18.6%"

6j (2 points)

What is the correlation of HighRBI, Home Runs, Hits, Runs, Assists, and Errors with Salary? Use only the cases in which both variables are measured. Round the answer to two decimal places.

##calculate correlation coefficient matrix for numeric variables, round to 2 decimal
round(cor(H[,c("Salary","Hits","Runs","HmRun","Assists","Errors","HighRBI")], use = "complete.obs") , 2)

        Salary Hits Runs HmRun Assists Errors HighRBI
Salary    1.00 0.44 0.42  0.34    0.03  -0.01    0.37
Hits      0.44 1.00 0.91  0.53    0.30   0.28    0.55
Runs      0.42 0.91 1.00  0.63    0.18   0.19    0.54
HmRun     0.34 0.53 0.63  1.00   -0.16  -0.01    0.70
Assists   0.03 0.30 0.18 -0.16    1.00   0.70    0.00
Errors   -0.01 0.28 0.19 -0.01    0.70   1.00    0.06
HighRBI   0.37 0.55 0.54  0.70    0.00   0.06    1.00

6k (2 points)

How did the salaries differ for players with and without HighRBI? Use the boxplot function and split the salary data by HighRBI status. Do HighRBI players have a higher median salary?

#careate a boxplot
boxplot(Salary~HighRBI,data= H, xlab = "RBI greater than 75")

##Based on the Boxplot, HighRBI players have a higher median salary than non-highRBI players.

6l (2 points)

Show a histogram of home runs using the hist function with breaks = 20 and freq = FALSE.

hist(H$HmRun,breaks = 20, freq = FALSE, main = "Home Runs", xlab = "Home Runs" )

7a (2 points)

What is the mean and standard deviation of Hits, Runs, Home Runs, RBI, Assists, Errors, and Salaries? Remove any missing values from the calculations. Round the answers to 1 decimal place.

# calculate mean for Hits, Runs, Home Runs, RBI, Assists, Errors, and Salaries

v <- c("Hits", "Runs", "HmRun", "RBI", "Assists","Errors","Salary")

multi.fun <- function(x) {
      c(mean = round(mean(x, na.rm = TRUE),1), sd = round(sd(x,na.rm = TRUE),1))
}

sapply(H[, v], multi.fun)

      Hits Runs HmRun  RBI Assists Errors Salary
mean 101.0 50.9  10.8 48.0   106.9    8.0  535.9
sd    46.5 26.0   8.7 26.2   136.9    6.4  451.1

7b (3 points)

Some players only get to play part-time. Show the mean and standard deviations for the same variables as in the previous question only for players with at least 300 AtBat.

sapply(H[H$AtBat >= 300, v], multi.fun)

      Hits Runs HmRun  RBI Assists Errors Salary
mean 126.3 64.3  13.8 60.3   133.4    9.4  617.1
sd    35.2 21.3   9.0 23.2   155.0    6.8  467.7

7c (3 points)

Show a scatter plot of Salary versus Home Runs for players with at least 300 AtBat.

plot(Salary ~ HmRun, H[H$AtBat>=300,])

7d (2 points)

There is a player with zero home runs and a salary over 2,000 (more than 2 million dollars). Who is this player? What does it look like happened during the season? Are these numbers accurate? Use the internet to search for this player’s results in 1986 and 1987.

#find the player with HmRun === 0 and Salary > 2000
subset(Hitters, HmRun == 0 & Salary>2000 )

              AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun
-Mike Schmidt    20    1     0    0   0     0     2     41     9      2
              CRuns CRBI CWalks League Division PutOuts Assists Errors
-Mike Schmidt     6    7      4      N        E      78     220      6
                Salary NewLeague
-Mike Schmidt 2127.333         N

#the player is Mike Schmidt. These numbers are inaccurate. He had 552 AtBat, 160 Hits, 37 HmRuns and 119 RBI in 1986 season.

8a (7 points)

Build a linear regression model and explain how (or why) you choose certain predictors in your model. Use 70% of the valid data for training and the remaining 30% of the valid data for testing. Please report the Room Mean Squared Error of the model on both the training and testing sets. Note that, what data are considered as “valid” is up to you based on your data exploration. For example, you can exclude certain data because of either missing data or outliers. But please explain how you determine your validate dataset.

round(cor(H[,c(1:13,16:19)], use = "complete.obs") , 2)

        AtBat Hits HmRun  Runs  RBI Walks Years CAtBat CHits CHmRun CRuns
AtBat    1.00 0.96  0.56  0.90 0.80  0.62  0.01   0.21  0.23   0.21  0.24
Hits     0.96 1.00  0.53  0.91 0.79  0.59  0.02   0.21  0.24   0.19  0.24
HmRun    0.56 0.53  1.00  0.63 0.85  0.44  0.11   0.22  0.22   0.49  0.26
Runs     0.90 0.91  0.63  1.00 0.78  0.70 -0.01   0.17  0.19   0.23  0.24
RBI      0.80 0.79  0.85  0.78 1.00  0.57  0.13   0.28  0.29   0.44  0.31
Walks    0.62 0.59  0.44  0.70 0.57  1.00  0.13   0.27  0.27   0.35  0.33
Years    0.01 0.02  0.11 -0.01 0.13  0.13  1.00   0.92  0.90   0.72  0.88
CAtBat   0.21 0.21  0.22  0.17 0.28  0.27  0.92   1.00  1.00   0.80  0.98
CHits    0.23 0.24  0.22  0.19 0.29  0.27  0.90   1.00  1.00   0.79  0.98
CHmRun   0.21 0.19  0.49  0.23 0.44  0.35  0.72   0.80  0.79   1.00  0.83
CRuns    0.24 0.24  0.26  0.24 0.31  0.33  0.88   0.98  0.98   0.83  1.00
CRBI     0.22 0.22  0.35  0.20 0.39  0.31  0.86   0.95  0.95   0.93  0.95
CWalks   0.13 0.12  0.23  0.16 0.23  0.43  0.84   0.91  0.89   0.81  0.93
PutOuts  0.31 0.30  0.25  0.27 0.31  0.28 -0.02   0.05  0.07   0.09  0.06
Assists  0.34 0.30 -0.16  0.18 0.06  0.10 -0.09  -0.01 -0.01  -0.19 -0.04
Errors   0.33 0.28 -0.01  0.19 0.15  0.08 -0.16  -0.07 -0.07  -0.17 -0.09
Salary   0.39 0.44  0.34  0.42 0.45  0.44  0.40   0.53  0.55   0.52  0.56
         CRBI CWalks PutOuts Assists Errors Salary
AtBat    0.22   0.13    0.31    0.34   0.33   0.39
Hits     0.22   0.12    0.30    0.30   0.28   0.44
HmRun    0.35   0.23    0.25   -0.16  -0.01   0.34
Runs     0.20   0.16    0.27    0.18   0.19   0.42
RBI      0.39   0.23    0.31    0.06   0.15   0.45
Walks    0.31   0.43    0.28    0.10   0.08   0.44
Years    0.86   0.84   -0.02   -0.09  -0.16   0.40
CAtBat   0.95   0.91    0.05   -0.01  -0.07   0.53
CHits    0.95   0.89    0.07   -0.01  -0.07   0.55
CHmRun   0.93   0.81    0.09   -0.19  -0.17   0.52
CRuns    0.95   0.93    0.06   -0.04  -0.09   0.56
CRBI     1.00   0.89    0.10   -0.10  -0.12   0.57
CWalks   0.89   1.00    0.06   -0.07  -0.13   0.49
PutOuts  0.10   0.06    1.00   -0.04   0.08   0.30
Assists -0.10  -0.07   -0.04    1.00   0.70   0.03
Errors  -0.12  -0.13    0.08    0.70   1.00  -0.01
Salary   0.57   0.49    0.30    0.03  -0.01   1.00

#From the correlation matrix, CAtBat,CHits,CHmRun,CRuns,CRBI,CWalks are highly correalted with #Salary, with correlation coefficient all above 0.50. Besides, a player's salary is based on his #career records, not very related to stats for the current season. Therefore, I will use thoese #variables as predictors for my model.

# there are missing values in set, which will hamper the linear models. Remove NA values:
H <- na.omit(Hitters[,c("Salary","CAtBat","CHits","CHmRun","CRuns","CRBI","CWalks")])
sum(is.na(H))

[1] 0

#find outliers using boxplot outcome
findOutliners <- function(x){
  outliers <- boxplot(x, plot=FALSE)$out
  which(x %in% outliers)
}

#remove all outliers 
dim(H)

[1] 263   7

H <- H[ -findOutliners(H$CHits), ]
H <- H[ -findOutliners(H$CRBI), ]
H <- H[ -findOutliners(H$CAtBat), ]
H <- H[ -findOutliners(H$CHmRun), ]
H <- H[ -findOutliners(H$CRuns), ]
H <- H[ -findOutliners(H$CWalks), ]
H <- H[ -findOutliners(H$Salary), ]
dim(H) #212 records left, 263-212 =51 observations removed

[1] 212   7

#subset train and test set
set.seed(55)
train_index = sample(1:nrow (H) , 0.7*nrow (H))
H_train = H[train_index, ]
H_test = H[- train_index, ]

#build linear regression model
lm <- lm(Salary ~ CRBI + CAtBat + CHmRun + CRuns + CHits + CWalks,data = H_train)
#prediction on training set
pred_train <- predict(lm)
#calculate rmse
lm_train_rmse <- sqrt(mean((H_train$Salary - pred_train)^2, na.rm = TRUE))
#prediction on test set
pred_test <- predict(lm,H_test)
lm_test_rmse <- sqrt(mean((H_test$Salary - pred_test)^2, na.rm = TRUE))

cbind(Model = "LM", `Training RMSE` = lm_train_rmse, `Test RMSE` = lm_test_rmse)

     Model Training RMSE     Test RMSE        
[1,] "LM"  "210.02047935331" "196.27687440876"

8b (7 points)

Repeat question 8a using KNN with 5 neighbors.

#subset data, assign cl labels
knn_train = H_train[ , -1]
knn_test = H_test[, -1]
knn_train_labels <- H_train[, 1] 
knn_test_labels <- H_test[, 1] 

#build knn model on train and test sets
knn_training_pred <- knn(knn_train, knn_train, knn_train_labels,  k = 5)
knn_test_pred <- knn(knn_train, knn_test,knn_train_labels,  k = 5)

#calculate rmse
knn_training_rmse <- sqrt(mean((as.numeric(knn_training_pred) - knn_train_labels)^2, na.rm = TRUE))
knn_test_rmse <- sqrt(mean((as.numeric(knn_test_pred) - knn_test_labels)^2, na.rm = TRUE))

cbind(Model = "KNN",`Training RMSE` = knn_training_rmse, `Test RMSE` = knn_test_rmse)

     Model Training RMSE      Test RMSE         
[1,] "KNN" "486.129796527746" "493.416545941693"