Example R programs and commands 26. Multivariate analysis; linear discriminant analysis # All lines preceded by the "#" character are my comments. # All other left-justified lines are my input. # All other indented lines are the R program output. ################## Displaying multivariate data ########################## # Load a 4-variate data set on irises: data(iris) # ...now there is a data.frame named "iris" ### Basic data frame examination # See the names of the variables in "iris" names(iris) ## Output: ## [1] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width" "Species" # See what values the "Species" member takes levels(iris\$Species) ## Output: ## [1] "setosa" "versicolor" "virginica" # See how many data rows there are in "iris" nrow(iris) ## Output: ## [1] 150 ### Data visualizations ## Histograms # Plot the histogram of "Sepal.Length" for "virginica" hist(iris\$Sepal.Length[iris\$Species=="virginica"]) # Plot histograms of the 4 variables for "virginica" par(mfrow=c(2,2)); # arrange 4 plots in a 2x2 matrix iris1 <- iris[iris\$Species=="virginica",1:4]; # Use first 4 columns for( i in 1:4 ) hist( iris1[,i], xlab=NULL, main=names(iris1)[i] ); # ...emphasizing the variates in the plot titles # Plot histograms of the 4 variables for all 3 species par(mfrow=c(3,4)); # arrange 12 plots in 3 1x4 matrices for( j in levels(iris\$Species) ) { iris1 <- iris[iris\$Species==j, 1:4]; # First 4 columns of species "j" for( i in 1:4 ) { hist( iris1[,i], xlab=names(iris1)[i], main=j ); } } # ...emphasizing the species in the titles atop each row of 4 plots ## Box plots # Box plots comparing the 4 variables by species par(mfrow=c(2,2)) plot( Sepal.Length ~ Species, data=iris) plot( Sepal.Width ~ Species, data=iris) plot( Petal.Length ~ Species, data=iris) plot( Petal.Width ~ Species, data=iris) ## Scatter plots # Plot 2-D scatterplots of all (4 choose 2)=6 pairs of variables # Create point labels "s", "v", "c" for the 3 species point.labels <- rep("X", nrow(iris)); point.labels[iris\$Species=="setosa"] <-"s" point.labels[iris\$Species=="virginica"] <-"v" point.labels[iris\$Species=="versicolor"] <-"c" pairs( iris[,1:4], pch=point.labels); # Plot 3-D scatterplots of all (4 choose 3)=4 triples of variables library(scatterplot3d); # after install.packages("scatterplot3d") par(mfrow=c(2,2)); for(i in 1:4) scatterplot3d(iris[c(-i,-5)], pch=point.labels); # ...a loop which is a tricky way of executing the equivalent 4 commands: scatterplot3d(iris[c(1,2,3)], pch=point.labels); scatterplot3d(iris[c(1,2,4)], pch=point.labels); scatterplot3d(iris[c(1,3,4)], pch=point.labels); scatterplot3d(iris[c(2,3,4)], pch=point.labels); ################# Testing for outliers by Mahalanobis distance ############ # This is one way to detect misidentified samples. # Plot the Mahalanobis distance of each 4-variables sample dev.off() # resets par() and clears the previous plot. iris1 <- iris[iris\$Species=="setosa", 1:4]; # First 4 columns of "setosa" dist2 <- mahalanobis(iris1, center=colMeans(iris1), cov=var(iris1)); plot(sqrt(dist2), main="Setosa", ylab="Mahalanobis distance (std. devs)" ); # ...emphasizing the species in the titles atop each plot # Plot histograms of the 4 variables for all 3 species par(mfrow=c(1,3)); # arrange 3 species in adjacent plots for( j in levels(iris\$Species) ) { iris1 <- iris[iris\$Species==j, 1:4]; # First 4 columns of species "j" dist<-sqrt(mahalanobis(iris1, center=colMeans(iris1), cov=var(iris1))); plot(dist, main=j, ylab="Mahalanobis distance (standard deviations)" ); } # ...emphasizing the species in the titles atop each plot ####################### Mahalanobis classification ################### # Decide species from new measurements based on which species mean # is closest by Mahalanobis distance. # # First compute mean vectors and variance matrices for the 3 classes # # Species 1: setosa species1 <- iris[iris\$Species=="setosa",1:4]; # the 4 data for each setosa mean1 <- colMeans( species1 ); var1<-var(species1); # # Species 2: versicolor species2 <- iris[iris\$Species=="versicolor",1:4]; # 4 data for each versicolor mean2 <- colMeans( species2 ); var2<-var(species2); # # Species 3: virginica species3 <- iris[iris\$Species=="virginica",1:4]; # 4 data for each virginica mean3 <- colMeans( species3 ); var3<-var(species3); # # Predict the species from 4 measurements new<-data.frame(Sepal.Length=4.9,Sepal.Width=3.3,Petal.Length=1.2,Petal.Width=0.2) mahalanobis(new, mean1, var1) ## OUTPUT: 2.338388 mahalanobis(new, mean2, var2) ## OUTPUT: 106.9140 mahalanobis(new, mean3, var3) ## OUTPUT: 179.5317 # ==> probably species 1 (setosa) # # Second example: new<-data.frame(Sepal.Length=6.8,Sepal.Width=3.0,Petal.Length=6.1,Petal.Width=2.0) mahalanobis(new,mean1,var1) ## OUTPUT: 845.2354 mahalanobis(new,mean2,var2) ## OUTPUT: 22.36729 mahalanobis(new,mean3,var3) ## OUTPUT: 2.275968 # ==> probably species 3 (virginica) ##################### Principal Components Analysis ##################### # Identify how many variables are really needed to account for the variance. data(iris) iris1<-iris[iris\$Species=="setosa", 1:4] pca <- princomp(iris1); pca screeplot(pca) summary(pca) # ...which indicates that 3 components contain 97% of the variance ################### Linear Discriminant Analysis ################## # Find a linear function that predicts class membership library(MASS) # contains the lda() function lda(Species ~ ., data=iris) # determine "Species" from all other measurements ## Output: ## Call: ## lda(Species ~ ., data = iris) ## ## Prior probabilities of groups: ## setosa versicolor virginica ## 0.3333333 0.3333333 0.3333333 ## ## Group means: ## Sepal.Length Sepal.Width Petal.Length Petal.Width ## setosa 5.006 3.428 1.462 0.246 ## versicolor 5.936 2.770 4.260 1.326 ## virginica 6.588 2.974 5.552 2.026 ## ## Coefficients of linear discriminants: ## LD1 LD2 ## Sepal.Length 0.8293776 0.02410215 ## Sepal.Width 1.5344731 2.16452123 ## Petal.Length -2.2012117 -0.93192121 ## Petal.Width -2.8104603 2.83918785 ## ## Proportion of trace: ## LD1 LD2 ## 0.9912 0.0088 ## # Check how well the discriminator classifies the data that built it: disc<-lda(Species ~ ., data=iris) predict(disc) # ...which produces lots of output. # Check just the classifications: predict(disc)\$class # Find the indices of samples that were misclassified allindices <- 1:nrow(iris) misclassified <- allindices[ predict(disc)\$class != iris\$Species ] misclassified # Find how the species were misclassified right <- iris\$Species[misclassified]; wrong <- predict(disc)\$class[misclassified] # Find the misclassification table --- truth labels the rows freq <- table( iris\$Species, predict(disc)\$class ) freq # raw frequencies freq/rowSums(freq) # probabilities; Prob(class|truth) row-sums to 1 # Predict the species from 4 measurements # First example: new<-data.frame(Sepal.Length=4.9,Sepal.Width=3.3,Petal.Length=1.2,Petal.Width=0.2) predict(disc, new) ## OUTPUT: \$class: setosa # ==> setosa # # Second example: new<-data.frame(Sepal.Length=6.8,Sepal.Width=3.0,Petal.Length=6.1,Petal.Width=2.0) predict(disc, new) ## OUTPUT: \$class: virginica # ==> virginica ############### Cross-validation ####################### # Train a classifier on part of the data, use that to predict the rest train <- sample(1:150, 75) # Half the indices, no repeats table(iris\$Species[train]) ## Output [your answer may differ] ## setosa versicolor virginica ## 26 24 25 classifier <- lda(Species ~ ., data=iris, subset = train) class <- predict(classifier, iris[-train,])\$class truth <- iris\$Species[-train] table(truth,class) # Repeat 100 times to estimate the probability of correct classification pcc <- matrix(0,3,3); # There are 3 species, so 3 classes for( n in 1:100 ) { train <- sample(1:150, 75) # Train the classifier with half of the data... classifier <- lda(Species ~ ., data=iris, subset = train) class <- predict(classifier, iris[-train,])\$class # ...to predict the remainder truth <- iris\$Species[-train] # ...which will be compared with the true species pcc <- pcc + table(truth,class) } pcc pcc/rowSums(pcc) # Convert to probabilities