Most of the assignment includes tasks that you have done in previous assignments.
I’ve written most of the challenging code for you. When I’ve asked you to write code, I’ve given you an example to follow.
Later code chunks depend on computations I’d like you to do, so those have been commented out by moving the {R} inside the first line of the code chunk. Simply move the {R} back after the first triple-tick to let those compute.
Overview
Goal: Develop a logistic regression model to predict the source of fish based on PCA components of the elemental isotopes and compare that with a QDA model.
I clean and transform the data for you.
I break the data into a “training” and “test” set so that we can develop (train) models and then see how they perform (test).
We do PCA for each element (Ba, Ca, Mg, and Sr) by combining the isotopes into a single principal component, each explaining more than 99% of the variability of the set of isotopes.
I do one, you have to do the other three. Copy/paste, then rename the variables.
Logistic regression on the training data
I do an example main-effects model that doesn’t fit well.
Put in two-way interactions and quadratic terms and it will fit well, even after backward selection.
The classification metrics on the training set are pretty good (Sensitivity and Specificity).
Classification on the test data
I do the hard part of transforming the test data into the same PCA space as the training data. (This is handy code if you find you ever need to do this.)
I have all the code to do the classification using both Logistic regression and QDA. You only need to interpret the results of the classification.
Finally, consider how well the models worked and what could be done to improve them.
You don’t need to do additional improvements, just brainstorm a little.
San Juan River Razorback Suckers
Peer mentor (Spring 2016) Adam Barkalow’s wife uses stable isotope ratios to analyze fish.
Razorback Suckers were collected in 2014 on the San Juan River. Elemental isotopic ratios from finrays were analyzed for Ba (Barium 56), Ca (Calcium 20), Mg (Magnesium 12), and Sr (Strontium 38). Finrays are non-lethally obtained and are used to detect natal origin since the material in the reys are partially developed early in life.
One issue is that hatchery fish can get into the river and lose their tags. It is important for environmental resource managers to know whether untagged fish are wild or hatchery fish. However, this is actually quite easy to determine with the isotope data that we have; I’d like you to solve a harder problem.
Another issue is distinguishing between the two Wild fish, NAP from Ponds and SJR from the River.
There are five fish sources in the dataset.
4 known sources, 1 mix of unknown sources
Hatchery
DEX = Dexter National Fish Hatchery
GJH = Ouray National Fish Hatchery, Grand Valley Unit
Wild
NAP = NAPI ponds
SJR = San Juan River
Unknown
UNK = untagged Razorback Suckers captured in the San Juan River
these could be from any of the above sources
Our goal is to classify the Wild fish into their location sources of Pond or River. First, we will use PCA to reduce the sets of highly correlated isotopes of the same element to a single PC feature for each element. Then, using the binary response of “Pond” versus “River”, and we’ll fit a logistic regression and do model selection. Finally, we’ll assess the expected classification accuracy using logistic regression for classification, and predict some observations I excluded from the model building.
Clean and transform data
Looking at the scatterplot matrix below, clean and/or transform the data if you think it will be helpful. Note that measurement error can be an issue in complicated biological measurements. Furthermore, a transformation might help separate observations that are tightly grouped in space.
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# First, download the data to your computer,# save in the same folder as this Rmd file.dat_sjrs <-read_csv("ADA2_CL_21_Clustering_SanJuanRazorbackSuckers_data2014.csv" ) %>%# the last set of observations only have a few isotopes, so excludena.omit() %>%# select the columns to use for analysisselect( Source, Ba137:Sr88 ) %>%# include only the Wild groupsfilter( Source %in%c("NAP", "SJR")## There are a few unual observations, remove those assuming measurement errors# remove two small Ca43 values , Ca43 >0.5 ) %>%mutate(# change to character so we can easily change the labelsSource = Source %>%as.character()# Simplify Source to be "Pond" and "River" , Source =case_when( Source =="NAP"~"Pond" , Source =="SJR"~"River" )# refactor with new labels , Source =factor(Source)# transforming the Ba values separates the tight clustering on the boundary , Ba137 =log10(Ba137) , Ba138 =log10(Ba138) )
Rows: 1512 Columns: 13
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): Station, Source, Type
dbl (10): Sort Key, Ba137, Ba138, Ca43, Mg24, Mg25, Mg26, Sr86, Sr87, Sr88
ℹ 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.
## NOTE HERE## Subset for classification later# Start random number generator in same place for everyone# and so that random partitions are the same each time code is runset.seed(3)# sample a subset of observation indices to predictind_pred <-sample.int(nrow(dat_sjrs) , size =100 ) %>%sort()ind_pred
# prediction subsetdat_sjrs_pred <- dat_sjrs %>%slice( ind_pred )# remove observations to predict from data to develop the modeldat_sjrs <- dat_sjrs %>%slice(-ind_pred )# data sizesdat_sjrs %>%dim()
[1] 275 10
dat_sjrs_pred %>%dim()
[1] 100 10
Known fish scatterplot
Note that this plot can take a while to generate. You’re welcome to subset the data further for this plot if some of the variables are redundant (highly correlated). You could probably get away with 5 columns of data without any loss of interpretation. If you want to do this, replace the dat_sjrs in the ggpairs() function with dat_sjrs %>% select(col1, col2, ...) and specify the columns to plot.
In this section, we’ll reduce the number of variables in the dataset by using PCA to generate new features which are linear combinations of selected variables. In this way, we can greatly reduce the dimension of the problem while retaining most of the information in the data.
(2 p) PCA of selected sets of variables
I expect that you’ll have four features at the end of this part.
I’ll do the first one as an example and you can do the rest.
Importance of components:
Comp.1 Comp.2
Standard deviation 0.2064719 0.009226185
Proportion of Variance 0.9980072 0.001992764
Cumulative Proportion 0.9980072 1.000000000
pca_Ba %>%loadings() %>%print(cutoff =0)
Loadings:
Comp.1 Comp.2
Ba137 0.714 0.700
Ba138 0.700 -0.714
Comp.1 Comp.2
SS loadings 1.0 1.0
Proportion Var 0.5 0.5
Cumulative Var 0.5 1.0
# If the loadings for Comp.1 are negative,# then switch the signs of the scores (observations on the PCA scale)# so that positive still indicates larger values.# For Ba, we need to use a positive sign in front of the scores to do this.dat_sjrs <- dat_sjrs %>%mutate(PC1_Ba = pca_Ba$scores[, "Comp.1"] %>%as.numeric() )
Note that Comp.1 explains 99.801% of the variability of both the Ba variables.
Calculate the remaining features to use with PCA below, and report the proportion of variance explained by the first component, Comp.1.
Solution
[answer]
Ca variables:
Mg variables:
Sr variables:
(1 p) Plot a scatterplot matrix of the new PCA variables
This plot should have five variables with a title indicating what is being plotted.
Hint: Use the dat_sjrs %>% select(Source, PC1_Ba, ...) command in the first argument of the ggpairs() function.
Solution
[answer]
Logistic Regression
In this section, we’ll use logistic regression to develop a model using the PCA features to calculate the probability that a given fish is from the Pond (versus River).
(2 p) Fit a logistic regression model
We will model the probability that a fish came from a pond. First we need a variable indicating whether it is from a pond or not.
Fit the logistic regression model below. If it does not fit, consider a more complex model (interactions and quadratic terms) until you find that the model fits. Perform backward selection and make sure reduced model also fits.
{R}
glm_pond <-
glm(
cbind(Pond, 1 - Pond) ~ PC1_Ba + PC1_Ca + PC1_Mg + PC1_Sr
, family = binomial
, data = dat_sjrs
)
summary(glm_pond)
# Test residual deviance for lack-of-fit (if > 0.10, little-to-no lack-of-fit)
dev_p_val <- 1 - pchisq(glm_pond$deviance, glm_pond$df.residual)
dev_p_val
# option: trace = 0 doesn't show each step of the automated selection
glm_pond_red_AIC <-
step(
glm_pond
, direction = "both"
, trace = 0
)
# the anova object provides a summary of the selection steps in order
glm_pond_red_AIC$anova
summary(glm_pond_red_AIC)
# Test residual deviance for lack-of-fit (if > 0.10, little-to-no lack-of-fit)
dev_p_val <- 1 - pchisq(glm_pond_red_AIC$deviance, glm_pond_red_AIC$df.residual)
dev_p_val
Note that the model doesn’t fit well since the lack-of-fit p-value < 0.10. Adding higher-order terms, such as two-way interactions and squared terms, may help.
In logistic regression, we have a prediction probability of success. We can find a threshold of that prediction probability (e.g., \(\hat{p}=0.3\)) as a boundary to classify two groups. Below, we summarize all possible threshold with an ROC curve. We also extract the optimal threshold (giving the jointly best Sensitivity and Specificity).
In this section, we’ll will compare the predictions using the logistic regression model to the discriminant analysis from last class. We’ll predict the observations that were held out (after first projecting them into the PCA feature space). Then we’ll create the confusion matrix (table of which observations were classified from which populations into which other populations) then compare the error rates between the two methods.
Projecting the “test” set of observations into the PCA space
In order to project the “test” subset of the data into the PCA space, we need to perform the same centering and PCA rotation (the loadings) as was done on the “training” subset of data. All this information is in the PCA objects calculated above. Below, we first subtract the centering values (the means of the training data) then use matrix multiplication with the loadings to calculate the linear combinations of the isotope data for the rotation. As before, we choose the same \(+\) or \(-\) sign as above for each isotopic element. Now the “test” subset is mapped onto the same PC1 axis as the “training” subset.
{R}
## The equation first subtracts the mean of the variables (in $center),
## then calculates the linear combination for PC1 via matrix multiplication (%*%).
# A function to perform the transformation
f_pca_pred <-
function(
dat
, var_list
, pca_obj
) {
## TESTING ## dat = dat_sjrs_pred; var_list = c("Ba137", "Ba138"); pca_obj = pca_Ba
out <-
(as.matrix(dat %>% select(var_list)) -
matrix(rep(pca_obj$center, nrow(dat)), byrow = TRUE, nrow = nrow(dat), ncol = length(pca_obj$center))) %*%
as.matrix(pca_obj$loadings[,1], ncol = 1) %>%
as.numeric()
return(out)
}
# Do the transformation for each element
dat_sjrs_pred$PC1_Ba <- f_pca_pred(dat_sjrs_pred, c("Ba137", "Ba138") , pca_Ba)
dat_sjrs_pred$PC1_Ca <- f_pca_pred(dat_sjrs_pred, c("Ca43") , pca_Ca)
dat_sjrs_pred$PC1_Mg <- f_pca_pred(dat_sjrs_pred, c("Mg24", "Mg25", "Mg26"), pca_Mg)
dat_sjrs_pred$PC1_Sr <- f_pca_pred(dat_sjrs_pred, c("Sr86", "Sr87", "Sr88"), pca_Sr)
The test data should look similar to the training data because it was randomly sampled from the whole.
Use these values for the Sensitivity and Specificity for the interpretation below.
{R}
qda_roc_pred$confusion_matrix
Summarize the error rates and note and differences you observe between the logistic regression and QDA methods.
Decide which method is preferred.
If you’d like, change the set.seed() value to draw a different random sample and see how the results change from sample to sample, but please return the seed to the original value before turning in your solutions.
Solution
[answer]
(1 p) Accuracy assessment
Comment on whether you are surprised by the prediction accuracy on the “test” subset. Did the classifier work much better, worse, or roughly what you expected?
Do you have any ideas for how classification could be improved between these two difficult-to-distinguish sources?
