########################################################################## ###################### # # R example code to calculate summary statistics for multivariate data: # ###################### # ########################################################################## # load 'HSAUR3' package: library(HSAUR3) # load USairpollution data set: data(USairpollution) # See first 5 data rows: head(USairpollution) # The observations are 41 cities (data were collected in 1981). # These are the 7 variables: # SO2: SO2 (sulphur dioxide) content of air in micrograms per cubic metre; # temp: average annual temperature in degrees Fahrenheit; # manu: number of manufacturing enterprises employing 20 or more workers; # popul: population size (1970 census) in thousands; # wind: average annual wind speed in miles per hour; # precip: average annual precipitation in inches; # predays: average number of days with precipitation per year. #attach(USairpollution) # attaching the data frame ############################################ ############################################ # Getting the sample mean vector: colMeans(USairpollution) # Getting the sample covariance matrix: var(USairpollution) # Rounding off: round(var(USairpollution),digits=2) # Note the diagonal elements of this are the sample variances for the 7 variables. ############################################ ############################################ # Getting the sample correlation matrix: cor(USairpollution) # Rounding off: round(cor(USairpollution),digits=2) ############################################ ############################################ # Calculations to show the relationship between the covariance and correlation matrices: # This is the long way... my.S <- var(USairpollution) D.minus.12 <- diag( 1/sqrt(diag(my.S) ) ) my.R <- D.minus.12 %*% my.S %*% D.minus.12 my.R # Getting the correlation matrix (if you are given the covariance matrix) the quick way: cov2cor(my.S) ## Be careful about the difference between the 'cor' and 'cov2cor' functions! # NOTE: If you have the DATA matrix and you want to get the sample correlation matrix, then use # the 'cor' function, e.g.: # cor(data_matrix) # If you have the COVARIANCE matrix and you want to calculate the correlation matrix from this, then use # the 'cov2cor' function, e.g.: # cov2cor(covar_matrix) ############################################ ############################################ # Getting distance matrix (Euclidean distances) between raw observations dis <- dist(USairpollution) dist2full<-function(ds){ n<-attr(ds,"Size") full<-matrix(0,n,n) full[lower.tri(full)]<-ds dist.mat<-full+t(full) rownames(dist.mat) <- colnames(dist.mat) <- attributes(ds)\$Labels return(dist.mat) } dis.matrix<-dist2full(dis) round(dis.matrix,digits=2) # Getting distance matrix (Euclidean distances) between SCALED observations std <- sapply(USairpollution, sd) # finding standard deviations of variables USairpollution.std <- sweep(USairpollution,2,std,FUN="/") # dividing each variable by its standard deviation dis <- dist(USairpollution.std) dis.matrix<-dist2full(dis) round(dis.matrix,digits=2) # Finding the smallest and largest distances within a distance matrix or distance object, # and also printing which pairs of objects these correspond to: # (Make sure dist2full function above has been copied into workspace first) # Copy and paste the following pick.smallest.largest function into R: ############ ## # pick.smallest.largest <- function(ds, howmany=1){ if (is.matrix(ds)==TRUE) { if (is.vector(rownames(ds))==TRUE) { ds.Labels <- rownames(ds) } else { ds.Labels <- 1:nrow(ds) } #print(ds.Labels) disLT <- ds[lower.tri(ds)] ds2 <- ds } else { if (is.vector(attributes(ds)\$Labels)==TRUE) { ds.Labels <- attributes(ds)\$Labels } else { ds.Labels <- 1:attributes(ds)\$Size } ds2 <- dist2full(ds) disLT <- ds2[lower.tri(ds2)] } smallest.few <- sort(disLT)[1:howmany] largest.few <- sort(disLT,decreasing=T)[1:howmany] print(paste("smallest", howmany, "distances are: ")); print(smallest.few) object.pairs <- object.pairs2 <- matrix(0,nrow=howmany,ncol=2) for (i in 1:howmany){ object.pairs[i,] <- which(ds2==smallest.few[i],arr.ind=T)[2,] } print("these smallest distances correspond to the pairs of objects given in the rows below:"); print(object.pairs) object.pairs.labels <- matrix(ds.Labels[object.pairs],ncol=2) print("the labels for these pairs of objects given in the rows below:"); print(object.pairs.labels) print(paste("largest", howmany, "distances are: ")); print(largest.few) for (i in 1:howmany){ object.pairs2[i,] <- which(ds2==largest.few[i],arr.ind=T)[2,] } print("these largest distances correspond to the pairs of objects given in the rows below:"); print(object.pairs2) object.pairs2.labels <- matrix(ds.Labels[object.pairs2],ncol=2) print("the labels for these pairs of objects are given in the rows below:"); print(object.pairs2.labels) } # ## ############ # Using the pick.smallest.largest function: pick.smallest.largest(dis,3) pick.smallest.largest(dis.matrix,3) # A measure of distance between the variables: 1 - abs(cor(USairpollution)) ############################################ ############################################ # Normal Q-Q plots for the first four variables in the data set, considered separately: # This checks normality of the marginal distributions par(mfrow=c(2,2)) # Setting up for a 2 by 2 array of plots qqnorm(USairpollution[,1], ylab="Ordered Observations", main = "Normal Q-Q Plot, Sulphur Dioxide") qqline(USairpollution[,1]) qqnorm(USairpollution[,2], ylab="Ordered Observations", main = "Normal Q-Q Plot, Temperature") qqline(USairpollution[,2]) qqnorm(USairpollution[,3], ylab="Ordered Observations", main = "Normal Q-Q Plot, Manufacturing") qqline(USairpollution[,3]) qqnorm(USairpollution[,4], ylab="Ordered Observations", main = "Normal Q-Q Plot, Population") qqline(USairpollution[,4]) par(mfrow=c(1,1)) # Using Everitt's chisplot function to check multivariate normality of entire data set: # Copy the chisplot function into R: ############### ### chisplot <- function(x) { if (!is.matrix(x)) stop("x is not a matrix") ### determine dimensions n <- nrow(x) p <- ncol(x) # xbar <- apply(x, 2, mean) S <- var(x) S <- solve(S,tol=1e-40) index <- (1:n)/(n+1) # xcent <- t(t(x) - xbar) di <- apply(xcent, 1, function(x,S) x %*% S %*% x,S) # quant <- qchisq(index,p) plot(quant, sort(di), ylab = "Ordered distances", xlab = "Chi-square quantile", lwd=2,pch=1) } ### ############### # Use the chisplot function: chisplot(as.matrix(USairpollution)) # Try it on a bigger data set (Fisher's iris data set, which is a built-in data set in R): fisher.iris <- iris[,1:4] chisplot(as.matrix(fisher.iris)) # For the three types of irises, separately: fisher.iris.setosa <- iris[1:50,1:4] chisplot(as.matrix(fisher.iris.setosa)) fisher.iris.versicolor <- iris[51:100,1:4] chisplot(as.matrix(fisher.iris.versicolor)) fisher.iris.virginica <- iris[101:150,1:4] chisplot(as.matrix(fisher.iris.virginica))