Model-based visualizations: Machine learning
Fitting a straight line to a curvilinear relationship
set.seed(4)
library(survival)
library(rpart)
dat = data.frame(X1 = sample(x = c(1,2,3,4,5), size = 1000, replace=T))
dat$t = rexp(1000, rate=dat$X1)
dat$t = dat$t / max(dat$t)
dat$e = rbinom(n = 1000, size = 1, prob = 1-dat$t )
# Tree Fit:
tfit = rpart(formula = Surv(t, event = e) ~ X1 , data = dat, control=rpart.control(minsplit=30, cp=0.01))
library(partykit)
tfit2 <- as.party(tfit)
plot(tfit2)
Using Pima diabetes data
::::
Using Pima diabetes data
More details here
import pandas as pd
from matplotlib import pyplot as plt
from scipy.cluster import hierarchy
import numpy as np
import seaborn as sns
# Data set
iris = sns.load_dataset('iris')
df = iris.iloc[:,:4]
# Calculate the distance between each sample
Z = hierarchy.linkage(df, 'average')
ind = hierarchy.cut_tree(Z,3)
# Plot with Custom leaves
hierarchy.set_link_color_palette(['m','b','g'])
hierarchy.dendrogram(Z, leaf_rotation=90, leaf_font_size=8, color_threshold=1.8, above_threshold_color='grey');
plt.show()
d = iris %>% janitor::clean_names() %>% select(-species) %>% dist()
cl1 <- hclust(d, method='single')
cl2 = hclust(d, method='average')
cl3 = hclust(d, method='complete')
iris1 <- iris %>% janitor::clean_names() %>% mutate(ind1=factor(cutree(cl1, 3)), ind2=factor(cutree(cl2,3)),
ind3=factor(cutree(cl3,3)))
plt1 <- ggplot(iris1, aes(x=sepal_length,y=sepal_width, color=ind1))+geom_point()+
stat_ellipse() + labs(x = 'sepal length',y='sepal width', title='Linkage: single')+
theme(legend.position='none')
plt2 <- ggplot(iris1, aes(x=sepal_length,y=sepal_width, color=ind2))+geom_point()+
stat_ellipse()+ labs(x = 'sepal length',y='sepal width', title='Linkage: average')+
theme(legend.position='none')
plt3 <- ggplot(iris1, aes(x=sepal_length,y=sepal_width, color=ind3))+geom_point()+
stat_ellipse()+ labs(x = 'sepal length',y='sepal width', title='Linkage: complete')+
theme(legend.position='none')
library(cowplot)
plot_grid(plt1, plt2, plt3, ncol=2)
Typically k-means clustering looks at spherical clusters
iris %>% janitor::clean_names() %>%
select(-species) %>%
as.matrix() %>%
kmeans(centers=3)-> cl
iris %>% janitor::clean_names() %>%
select(-species) %>%
mutate(cl = as.factor(cl$cluster)) %>%
ggplot(aes(sepal_length, sepal_width, color=cl))+
geom_point() +
stat_ellipse() +
labs(x = 'Sepal length', y = 'Sepal Width')+
theme(legend.position='none')
Silhouette statistic
The silhouette value is a measure of how similar an object is to its own cluster (cohesion) compared to other clusters (separation). The silhouette ranges from −1 to +1, where a high value indicates that the object is well matched to its own cluster and poorly matched to neighboring clusters. If most objects have a high value, then the clustering configuration is appropriate. If many points have a low or negative value, then the clustering configuration may have too many or too few clusters. (Wikipedia)
from sklearn.cluster import KMeans
from yellowbrick.cluster import silhouette_visualizer
from yellowbrick.datasets import load_credit
# Load a clustering dataset
X, y = load_credit()
# Specify rows to cluster: under 40 y/o and have either graduate or university education
X = X[(X['age'] <= 40) & (X['edu'].isin([1,2]))]
# Use the quick method and immediately show the figure
silhouette_visualizer(KMeans(5, random_state=42), X, colors='yellowbrick')
SilhouetteVisualizer(ax=<Axes: title={'center': 'Silhouette Plot of KMeans Clustering for 18941 Samples in 5 Centers'}, xlabel='silhouette coefficient values', ylabel='cluster label'>,
colors='yellowbrick',
estimator=KMeans(n_clusters=5, random_state=42))In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. SilhouetteVisualizer(ax=<Axes: title={'center': 'Silhouette Plot of KMeans Clustering for 18941 Samples in 5 Centers'}, xlabel='silhouette coefficient values', ylabel='cluster label'>,
colors='yellowbrick',
estimator=KMeans(n_clusters=5, random_state=42))KMeans(n_clusters=5, random_state=42)
KMeans(n_clusters=5, random_state=42)
from sklearn.cluster import KMeans
from yellowbrick.cluster import SilhouetteVisualizer
from yellowbrick.datasets import load_nfl
# Load a clustering dataset
X, y = load_nfl()
# Specify the features to use for clustering
features = ['Rec', 'Yds', 'TD', 'Fmb', 'Ctch_Rate']
X = X.query('Tgt >= 20')[features]
# Instantiate the clustering model and visualizer
model = KMeans(5, random_state=42)
visualizer = SilhouetteVisualizer(model, colors='yellowbrick')
visualizer.fit(X) # Fit the data to the visualizer
visualizer.show() # Finalize and render the figure
KElbowVisualizer(ax=<Axes: title={'center': 'Distortion Score Elbow for KMeans Clustering'}, xlabel='k', ylabel='distortion score'>,
estimator=KMeans(n_clusters=9, random_state=4), k=(2, 10))In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook. KElbowVisualizer(ax=<Axes: title={'center': 'Distortion Score Elbow for KMeans Clustering'}, xlabel='k', ylabel='distortion score'>,
estimator=KMeans(n_clusters=9, random_state=4), k=(2, 10))KMeans(n_clusters=9, random_state=4)
KMeans(n_clusters=9, random_state=4)
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
from yellowbrick.cluster import KElbowVisualizer
# Generate synthetic dataset with 8 random clusters
X, y = load_nfl()
# Instantiate the clustering model and visualizer
model = KMeans()
visualizer = KElbowVisualizer(
model, k=(2,10), metric='calinski_harabasz', timings=False
)
visualizer.fit(X) # Fit the data to the visualizer
visualizer.show() # Finalize and render the figure
# Libraries
import seaborn as sns
import pandas as pd
from matplotlib import pyplot as plt
# Data set
url = 'https://raw.githubusercontent.com/holtzy/The-Python-Graph-Gallery/master/static/data/mtcars.csv'
df = pd.read_csv(url)
df = df.set_index('model')
# Prepare a vector of color mapped to the 'cyl' column
my_palette = dict(zip(df.cyl.unique(), ["orange","yellow","brown"]))
row_colors = df.cyl.map(my_palette)
# plot
sns.clustermap(df, metric="correlation", method="single", cmap="Blues", standard_scale=1, row_colors=row_colors);
PCA (principal components analysis) is a method to “rotate” the data based on variability of the data cloud, resulting in fewer variables that capture the information in the data cloud. It’s a popular method of dimension reduction
PCA first finds the direction in which there is maximum variance
Then finds the directional orthogonal (uncorrelated) to the first direction with the most variance
And so on
In this case, these two new variables capture 99.7% of the variability in the iris measurements
It’s not uncommon that PCA separates groups along one of the axes
tSNE (t-distributed stochastic neighbor embedding) is a machine learning algorithm to visualize high-dimensional data.
It is a non-linear dimension reduction technique, as opposed to PCA, which is a linear technique
tSNE is also probabilistic, so unless you fix the seed of the random number generator, you will get different results on different runs
Using iris
UMAP (uniform manifold approximation and projection) is another non-linear dimension-reduction technique that is increasingly used in bioinformatics to cluster high-dimensional data. It is believed to preserve both local and global structure and give a “fair” represenation of the data.
Using iris
from sklearn.decomposition import PCA
pca = PCA()
prcomp = pca.fit_transform(df)
out = pd.DataFrame(prcomp, columns = ['PCA1','PCA2','PCA3','PCA4'])
out = pd.concat((out, iris.species), axis=1)
fig,ax = plt.subplots(1,3, sharex=False, sharey=False)
sns.scatterplot(data=out, x='PCA1',y='PCA2', hue='species', ax=ax[0])
ax[0].set_ylabel('')
ax[0].set_xlabel('')
ax[0].set_title('PCA')
from sklearn.manifold import TSNE
tsne = TSNE(n_components=2, perplexity=5, n_iter=300)
out = pd.DataFrame(tsne.fit_transform(df), columns = ['TSNE1','TSNE2'])
out = pd.concat((out, iris.species), axis=1)
sns.scatterplot(data = out, x = 'TSNE1', y = 'TSNE2', hue='species', ax=ax[1])
ax[1].set_ylabel('')
ax[1].set_xlabel('')
ax[1].set_title('t-SNE')
from umap import UMAP
u = UMAP()
embedding = pd.DataFrame(u.fit_transform(df), columns = ['UMAP1','UMAP2'])
embedding = pd.concat((embedding, iris.species), axis=1)
sns.scatterplot(data=embedding, x='UMAP1', y='UMAP2', hue='species', ax=ax[2])
ax[2].set_ylabel('')
ax[2].set_xlabel('')
ax[2].set_title('UMAP')
plt.show()
from yellowbrick.datasets import load_credit
from yellowbrick.features import PCA
# Specify the features of interest and the target
X, y = load_credit();
classes = ['account in default', 'current with bills']
visualizer = PCA(scale=True, classes=classes);
visualizer.fit_transform(X, y);
visualizer.show()
from yellowbrick.datasets import load_credit
from yellowbrick.features import Manifold
# Specify the features of interest and the target
X, y = load_credit();
classes = ['account in default', 'current with bills']
visualizer = Manifold(manifold='tsne', classes = classes);
visualizer.fit_transform(X, y);
visualizer.show()
A Voronoi diagram splits up a plane based on a set of original points. Each polygon, or Voronoi cell, contains an original point and all areas that are closer to that point than any other.
This can be used in different contexts, but especially with geospatial data
Here, we’re plotting the locations of airports in the USA
airports <- read.csv(here::here("slides/ml/data/airport-locations.tsv"), sep="\t", stringsAsFactors=FALSE)
source(here::here('slides',"ml/latlong2state.R"))
airports$state <- latlong2state(airports[,c("longitude", "latitude")])
airports_contig <- na.omit(airports)
# Projection
library(mapproj)
airports_projected <- mapproject(airports_contig$longitude, airports_contig$latitude, "albers", param=c(39,45))
par(mar=c(0,0,0,0))
plot(airports_projected, asp=1, type="n", bty="n", xlab="", ylab="", axes=FALSE)
points(airports_projected, pch=20, cex=0.1, col="red")
A Voronoi diagram splits up a plane based on a set of original points. Each polygon, or Voronoi cell, contains an original point and all areas that are closer to that point than any other.
This can be used in different contexts, but especially with geospatial data
library(deldir)
par(mar=c(0,0,0,0))
plot(airports_projected, asp=1, type="n", bty="n", xlab="", ylab="", axes=FALSE)
points(airports_projected, pch=20, cex=0.1, col="red")
vtess <- deldir(airports_projected$x, airports_projected$y)
plot(vtess, wlines="tess", wpoints="none", labelPts=FALSE, add=TRUE, lty=1)
Credit: Nathan Yau
Also, see here for a D3 version)
One of the things we need to understand (beyond patterns we may see in PCA and other dimension reduction methods) is whether we have sets of correlated features that are either explanatory or serendipitious. One way of visualizing this is using a two-dimensional ranking of features.
from yellowbrick.datasets import load_credit
from yellowbrick.features import Rank2D
# Load the credit dataset
X, y = load_credit()
# Instantiate the visualizer with the Pearson ranking algorithm
visualizer = Rank2D()
visualizer.fit(X, y) # Fit the data to the visualizer
visualizer.transform(X) # Transform the data
visualizer.show() # Finalize and render the figure
Something quite similar is seen in genetics in a linkage disequilibrium or LD plot.
coefplot produces a ggplot2 graphic, so we could add to this using what we know about ggplot
There is actually a confidence interval being drawn, but it’s too small to see compared to the dot representing the estimate
Here, I’m taking a 1% sample of the diamonds data for modeling. This increases the magnitude of the standard errors of the coefficients
set.seed(2405)
diamonds2 <- slice_sample(diamonds, prop=0.01) #<<
model2 <- lm(log(price) ~ log(carat) + cut + color + clarity + depth,
data=diamonds2)
model2_tidy <- broom::tidy(model2)
ggplot(model2_tidy %>% filter(term!= '(Intercept)'),
aes(x =term, y = estimate,
ymin = estimate - 2*std.error,
ymax = estimate + 2*std.error))+
geom_pointrange(color='blue')+
labs(x = 'Predictor', y = 'Coefficient',
title = 'Coefficient Plot')+
geom_hline(yintercept = 0, linetype=2)+
coord_flip()
The other variables are held at the following values:
You can also draw marginal effects for categorical predictors
You can also look at marginal effects for one variable conditional on one or more other variables
We expect, in well-fitting models, that the residuals are on average 0, and that they have no pattern with the predicted/fitted values
An assumption of linear (OLS) regression is that the errors are distributed as a Gaussian (normal) distribution.
This can be visually checked using a quantile-quantile (Q-Q) plot.
If the residuals follow a Gaussian distribution, this plot should be close to the diagonal line
We can see both the pattern and histogram of residuals
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from yellowbrick.datasets import load_concrete
from yellowbrick.regressor import ResidualsPlot
# Load a regression dataset
X, y = load_concrete()
# Create the train and test data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
# Instantiate the linear model and visualizer
model = LinearRegression()
visualizer = ResidualsPlot(model, hist=False, qqplot=True)
visualizer.fit(X_train, y_train) # Fit the training data to the visualizer
visualizer.score(X_test, y_test) # Evaluate the model on the test data
visualizer.show() # Finalize and render the figure
Cook’s distance or Cook’s D is a commonly used estimate of the influence of a data point when performing a least-squares regression analysis
Cook’s distance measures the effect of deleting a given observation. Points with a large Cook’s distance are considered to merit closer examination in the analysis.
Leverage is a measure of how far away the independent variable values of an observation are from those of the other observations.
High-leverage points are those observations, if any, made at extreme or outlying values of the independent variables such that the lack of neighboring observations means that the fitted regression model will pass close to that particular observation.
Uses partykit package
Feature importances are typically a measure of how much a predictive performance metric is reduced when the values of the a variable is scrambled (permuted). They represent the predictive power of a variable
Feature importances and p-values aren’t necessarily related; a variable can be predictive and not statistically significant and vice versa. The former depends on the size of the effect, and the latter depends on the signal-to-noise ratio
Feature importances are typically a measure of how much a predictive performance metric is reduced when the values of the a variable is scrambled (permuted). They represent the predictive power of a variable
Feature importances and p-values aren’t necessarily related; a variable can be predictive and not statistically significant and vice versa. The former depends on the size of the effect, and the latter depends on the signal-to-noise ratio
Lasso regression fits a penalized regression with \(L_1\) penalties on the coefficients. This model minimizes
\[ ||y-X\beta||^2_2 + \lambda ||\beta||_1 \] where \(||x||_1 = \sum |x_i|/n\) and \(||x||_2 = \sqrt{\sum x_i^2/n}\)
This forces slope coefficients to 0 as the value of \(\lambda\) changes.
The plot shows the values of the slope coefficients for different values of \(\log(\lambda)\)
library(mlbench)
library(glmnet)
data("BreastCancer")
X = BreastCancer %>%
filter(complete.cases(.)) %>%
select(-Id, -Class) %>%
mutate(across(everything(), ~as.numeric(as.character(.)))) %>%
as.matrix()
y <- BreastCancer %>%
filter(complete.cases(.)) %>%
pull(Class) %>% as.numeric()
fit <- glmnet(X, y, family='binomial')
plot(fit, 'lambda', main = 'Lasso regression', label = T)
The elastic net regression (ELR) minimizes the penalized objective function
\[ ||y-X\beta||^2_2 + \alpha\lambda ||\beta||_1 + \frac{(1-\alpha)\lambda}{2} ||\beta||_2^2 \] Notice that in ELR, the coefficients move to 0 at a slower rate than for lasso regression
Regularized regression is typically used for variable selection in regression models, where \(\lambda\) is estimated via cross-validation
Note that the actual coefficients are provably biased, so they are not interpreted.
library(mlbench)
data("BreastCancer")
X = BreastCancer %>%
filter(complete.cases(.)) %>%
select(-Id, -Class) %>%
mutate(across(everything(), ~as.numeric(as.character(.)))) %>%
as.matrix()
y <- BreastCancer %>%
filter(complete.cases(.)) %>%
pull(Class) %>% as.numeric()
fit <- glmnet(X, y, family='binomial', alpha = 0.5)
plot(fit, 'lambda', main = 'ELR (alpha=0.5)', label = T)
A receiver operating characteristic (ROC) curve is a graphical plot that illustrates the diagnostic ability of a binary classifier system as its discrimination threshold is varied.
The ROC curve is created by plotting the true positive rate (TPR) against the false positive rate (FPR) at various threshold settings. The true-positive rate is also known as sensitivity, recall or probability of detection in machine learning. The false-positive rate is also known as probability of false alarm and can be calculated as (1 − specificity).
A classifier no better than chance will have a ROC curve along the unit diagonal
The precision-recall curve is used for evaluating the performance of binary classification algorithms. It is often used in situations where classes are heavily imbalanced.
A “good” classifier will have high precision and high recall throughout the graph, to be close to the upper right corner
We can predict the probability that an individual is a member of a class, rather than the predicted class. One property of the learning methods we can then assess is calibration.
Calibration curves are used to evaluate how calibrated a classifier is i.e., how the probabilities of predicting each class label differ. We bin the data, and for each bin we compute the average predicted probability (x-axis) and the proportion of positives (y-axis).
We expect bins with low proportion of positives to also have low average probability of being in the positive class, and vice versa
Survival analysis methods can be applied to a variety of domains:
Survival data is characterized by censoring, i.e. when an individual is lost to followup while still “alive”, and is assumed to be independent of “failure” (violations of this assumption requires different methods)
This plot looks at how long we followed each individual and whether they failed (red) or were censored (white) the last time we saw them
library(survival)
pbc$status2 = factor(ifelse(pbc$status==2, 'Dead', 'Alive'))
pbc <- pbc %>%
mutate(id = factor(id)) %>%
mutate(id = fct_reorder(id, time))
ggplot(pbc, aes(x = time,y=id))+
geom_point( aes(color=factor(status2)), size=2)+
geom_segment(aes(x=time, xend=0, y=id, yend=id, color=factor(status2)))+
scale_color_manual(values = c('white','red'))+
labs(x = 'Time (days)', title = 'PBC study', color='Status')+
theme(axis.text.y=element_blank(),
axis.title.y=element_blank())
The Kaplan-Meier estimator is used to estimate the survival function \(\Pr(T > t)\), where \(T\) is the time to failure.
The visual representation of this function is usually called the Kaplan-Meier curve, and it shows what the probability of failure is at a certain point of time.
From the plot, at 1000 days, about 75% of patients are expected to survive.
Note
Add KM curves from ggsurvfit, and Python examples from https://erdogant.github.io/kaplanmeier/pages/html/index.html and
