|
Mahalanobis distance
The following code illustrates the calculation of Mahalanobis distances in a "climate space"
described by two climate variables.
# Mahalonobis distances
# July and annual temperature
X <- as.matrix(midwtf2[,4:5])
X.cov <- var(scale(X)) # standardize first
X.cov <- var(X)
X.mean <- apply(X,2,mean)
X.mah <- mahalanobis(X, X.mean, X.cov)
# plot the variables, conf. ellipses and M.D.
library(car)
plot(tmeanjul, tmeanyr)
ellipse(X.mean, X.cov, radius=1, add=TRUE)
ellipse(X.mean, X.cov, radius=2, add=TRUE)
points(tmeanjul, tmeanyr, cex=10*(X.mah-min(X.mah))/(max(X.mah)-
min(X.mah)))
The following
script finds the Mahalanobis distances for each pollen site, where distances are
now in "pollen space.".
# Mahalonobis distances for individual
spectra, Midwestern pollen data
X <- as.matrix(midwtf2[,7:21])
#X.cov <- var(scale(X)) # standardize first
X.cov <- var(X)
X.mean <- apply(X,2,mean)
X.mah <- mahalanobis(X, X.mean, X.cov, tol.inv=1e-16)
The Mahalanobis distance can be graphically displayed using a
2-D bubble plot.
# Map the study area
library(maptools)
library(RColorBrewer)
# read and plot shape file
midwotl.shp <- readShapeLines(file.choose(),
proj4string=CRS("+proj=longlat"))
plot(midwotl.shp)
points(longitud, latitude, cex=10*(X.mah-min(X.mah))/(max(X.mah)-
min(X.mah)))
|