--- title: "Correlation and scagnostic guided parallel coordinate plots" author: "Catherine B. Hurley and R.W. Oldford" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Correlation and scagnostic guided parallel coordinate plots} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: chunk_output_type: console --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) options(digits=3) ``` PairViz is an R package which provides orderings of objects for visualisation purposes. This vignette demonstrates the use of PairViz constructing parallel coordinate plots, which place emphasis on different data patterns. ## Correlation guided plots of mtcars using Hamiltonians First, we select six of the variables in `mtcars`, and make a parallel coordinate plot. We will use the variant on a PCP provided by PairViz, which is the function `pcp`. ```{r, fig.align="center",fig.width=4, fig.height=2.5} suppressPackageStartupMessages(library(PairViz)) data <- mtcars[,c(1,3:7)] pcp(data,horizontal=TRUE,lwd=2, col="grey50", main = "Standard PCP") ``` The standard PCP will plot the variables in order of appearance in the data frame. It is obvious from this plot for instance that `mpg` and `disp` are negatively correlated. By following the line segments you might be able to see that `mpg` and `hp` are also negatively correlated, but but associations with other variables are difficult to ascertain. We could use one of the eulerian functions of PairViz to produce orderings of variables where all pairs of variables are adjacent. In the next display, we use `hpaths`, which gives three hamiltonians where each pair of variables is adjacent at least once. ```{r , fig.align="center",fig.width=7, fig.height=2.5} o <- hpaths(1:ncol(data),matrix=FALSE) par(cex.axis=.7) pcp(data,order=o,horizontal=TRUE,lwd=2, col="grey50", main = "Hamiltonian decomposition") ``` From this we see that `mpg` is also strongly negatively correlated with `wt`, `hp` and `drat`, and that there is not a strong association between `mpg` and `drat`. To make the patterns in the PCP a bit clearer, we will construct an eulerian where high correlation variables appear early on in the sequence. The code below makes a `dist` with the correlation between variables. ```{r} corw <- as.dist(cor(data)) ``` As the function `eulerian` constructs a path visiting lower weight edges first (and we want to visit high correlation variables first), we form the path as ```{r} o <- eulerian(-corw) o ``` Forming the PCP based on this order we get ```{r , fig.align="center",fig.width=7, fig.height=2.5} par(cex.axis=.7) pcp(data,order=o,horizontal=TRUE,lwd=2, col="grey50", main = "Weighted eulerian") ``` The first three panels on the left hand side have many parallel line segments indicative of positively correlated variables. On the right hand side of the PCP, the panels have many high-low line segments which means variables are negatively correlated. PairViz has an augmented version of a PCP which shows summary measures for each pair of variables, to assist in interpretation. The code below constructs the summary measures, one column for positive correlations, the second for negative correlations. ```{r , fig.align="center",fig.width=7, fig.height=3} corw <- dist2edge(corw) edgew <- cbind(corw*(corw>0), corw*(corw<0)) par(cex.axis=.7) guided_pcp(data,edgew, path=o,pcp.col="grey50" ,lwd=2, main="Weighted eulerian with correlation guide", bar.col = c("lightpink1","lightcyan1"), bar.ylim=c(-1,1),bar.axes=TRUE) ``` The above plot shows clearly that correlation generally decreases from left to right. If you prefer, omit the guides, but instead colour the panels according to the sign of the correlation. ```{r , fig.align="center",fig.width=7, fig.height=2.5} pathw <- path_weights(corw,o) corcols <- ifelse(pathw>0, "lightpink1", "lightcyan") par(cex.axis=.7) pcp(data,order=o,col="grey50" ,lwd=2, main="Weighted eulerian with correlation guide", panel.colors = corcols) ``` The guided_pcp function also has a panel.colors argument, so it is possible to keep the guides and also colour the panels. ## Scagnostic guided plots of sleep data using Hamiltonians We access the data, remove NAs, transform two highly skewed variables, give variables shortnames, and set up a colour vector for plots. ```{r} if (!requireNamespace("alr4", quietly = TRUE)){ install.packages("alr4") } suppressPackageStartupMessages(library(alr4)) data(sleep1) data <- na.omit(sleep1) # these vars changed to factors in alr4, change from alr3 data$D <- as.numeric(data$D) data$P <- as.numeric(data$P) data$SE <- as.numeric(data$SE) # logging the brain and body weights data[,4:5] <- log(data[,4:5]) # short variable names colnames(data) <- c("SW","PS" ,"TS" ,"Bd", "Br","L","GP","P" ,"SE" , "D" ) # colours for cases, split Life values into 3 equal sized groups cols1 <- scales::alpha(c("red","navy","lightblue3" ),.6) cols <- cols1[cut(rank(data$L),3,labels=FALSE)] ``` Calculate scagnostics for the data. `sc` is a matrix whose rows are values for 9 scagnostics, ```{r} library(scagnostics) library(RColorBrewer) sc <- scagnostics(data) scags <- rownames(sc) scags ``` As we will make plots involving different scagnostics, we will assign colours to scagnostics, for consistency across plots. ```{r} scag_cols <- rev(brewer.pal(9, "Pastel1")) names(scag_cols) <- scags ``` Define a utility function which returns a subset of scagnostics, retaining the class attribute. ```{r} select_scagnostics <- function(sc,names){ sc1 <- sc[names,] class(sc1) <- class(sc) return(sc1) } ``` Consider the outlying scagnostic. Suppose we want to construct a PCP, each variable appearing once, where pairs of variables with a high outlier score appear adjacently in the sequence. ```{r} scOut <- select_scagnostics(sc,"Outlying") dOut <- edge2dist(scOut) # dOut is a dist dOut <- as.matrix(dOut) rownames(dOut) <- colnames(dOut)<- names(data) ``` `dOut` is a symmetric matrix where each entry gives the outlying score for the scatterplot labelled by the row and column names. Note that the function `edge2dist` relies on the fact that objects of class scagnostics are in order (1,2), (1,3), (2,3),(1,4) etc To find the variable ordering with the highest (or nearly highest) total outlier score use one of ```{r eval=F} o <- order_best(-dOut, maxexact=10) o <-order_best(-dOut) o <-order_tsp(-dOut) ``` The first call to `order_best` above finds the best ordering by a brute force evaluation of all 10! factorial permutations, and so this is a bit slow. Without the `maxexact=10` input, order_best evaluates only a sample of permulations for sequences of length above 9. order_tsp uses a TSP solver from package TSP. We will use the result of the first call to order_best, which gives ```{r} o <-c( 2 , 4 , 1, 5 , 6, 7 , 3 , 8, 9, 10) ``` The guided PCP based on this ordering is ```{r , fig.align="center",fig.width=7, fig.height=3} par(tcl = -.2, cex.axis=.8,mgp=c(3,.3,0)) guided_pcp(data,scOut, path=o,pcp.col=cols,lwd=1.4, main= "Best Hamiltonian for outliers",bar.col = scag_cols["Outlying"],legend=FALSE,bar.axes=TRUE,bar.ylim=c(0,max(scOut))) ``` Notice that panels involving discrete-valued variables P,SE and D score zero on the outlying index. The L-GP pair of variables has the highest outlier score and two outliers are evident. The outlier cases are the two species (Human and Asian Elephant) with the highest life expectancy (L). Asian Elephant also has the highest value of gestation time (GP) ```{r} outliers <- order(data$L, decreasing=T)[1:2] rownames(data)[outliers] ``` For future use we will construct a colour vector marking these two outliers. ```{r} colOut <- rep("grey50", nrow(data)) colOut[outliers[1]] <- "red" # Human colOut[outliers[2]] <- "blue" ``` Suppose next we want to pick the ordering of PCP axes where high scores on the two measures Striated and Sparse are obtained. ```{r} scSS <- t(select_scagnostics(sc,c("Striated", "Sparse"))) ``` As we did with the calculation for Outliers, we can turn `scSS` into a distance matrix and then use one of `order_best` or `order_tsp` to produces the best ordering. ```{r eval=F} dSS <- edge2dist(scSS[,1]) + edge2dist(scSS[,2]) # You might think edge2dist(scSS[,1]+ scSS[,2]) would work, but as scSS[,1]+ scSS[,2] is # not of class scagnostics, edge2dist will not fill the dist in the correct order dSS <- as.matrix(dSS) rownames(dSS) <- colnames(dSS)<- names(data) order_best(-dSS,maxexact=10) ``` A shortcut calculation is ```{r eval=F} find_path(-scSS, order_best,maxexact=10) # for the best path # or find_path(-scSS, order_best) # for a nearly "best" path ``` The best path gives ```{r} o <- c(4, 10 , 2 , 9 , 1, 7 , 8, 6 , 5 ,3) ``` The guided PCP based on this ordering is ```{r , fig.align="center",fig.width=7, fig.height=3} par(tcl = -.2, cex.axis=.8,mgp=c(3,.3,0)) guided_pcp(data,scSS, path=o,pcp.col=cols,lwd=1.4, main= "Best Hamiltonian for Striated + Sparse", bar.col = scag_cols[c("Striated", "Sparse")], legend=FALSE,bar.axes=TRUE,bar.ylim=c(0,.6)) ``` The combined Striated and Sparse index rewards panels with discrete variables, as all panels with an index level above 0.4 involve a discrete variable. ## Scagnostic guided plots of sleep data using Eulerians For the sleep data there are 10 variables so an Eulerian path (which places all pairs of variables adjacently) will be rather long. Here is the Eulerian driven by the Outlyingness measure. ```{r} o <- find_path(-scOut, eulerian) o ``` The Outlying scores for the pairs of variables specified by this path are ```{r} pathw <- path_weights(scOut,o) head(pathw) ``` ```{r , fig.align="center",fig.width=7, fig.height=3,fig.show='hold'} par(tcl = -.2, cex.axis=.6,mgp=c(3,.3,0)) guided_pcp(data,scOut, path=o[1:25],pcp.col=cols,lwd=1.4, main= "First 25: Eulerian for Outlying", bar.col = scag_cols["Outlying"], legend=FALSE,bar.axes=TRUE,bar.ylim=c(0,max(scOut))) guided_pcp(data,scOut, path=o[25:50],pcp.col=cols,lwd=1.4, main= "Last 26: Eulerian for Outlying", bar.col = scag_cols["Outlying"], legend=FALSE,bar.axes=TRUE,bar.ylim=c(0,max(scOut))) ``` Many of the Outlying scores are zero, in fact the last 25 scores are zero. We see that GP-L has the highest overall Outlying score, and the higher Outlying scores involve the variables L,GP, Bd, SW,Br only. ## Focus on interesting plots The Eulerian path visiting all pairs of variables for the sleep data is of length 10. One way to focus on "interesting"" pairs (measuring interesting by outlyingness here) is the weighted Eulerian, where interesting pairs are visited early in the sequence, as above. Another approach is to drop all pairs of variables which are below some interestingness threshold, and construct a plot showing the rest. Here we will construct a PCP showing pairs of variables whose outlier score is .2 or higher. We will do this by constructing a graph containing variable pairs whose outlier score is .2 or higher, and then forming an Eulerian of this graph. One way to form the graph is to build the complete graph based on `dOut`, the matrix of outlier scores, and then remove low weight edges. ```{r} g <- mk_complete_graph(dOut) ``` Then form a graph from `g` keeping edges whose weight is .2 or above. The function `dn_graph` also deletes nodes with no edges. ```{r} g1 <- dn_graph(g,.2 , test=`>=`) ``` A utility function for plotting a graph: ```{r } requireNamespace("igraph") igplot <- function(g,weights=FALSE,layout=igraph::layout_in_circle, vertex.size=60, vertex.color="lightblue",...){ g <- igraph::graph_from_graphnel(as(g, "graphNEL")) op <- par(mar=c(1,1,1,1)) if (weights){ ew <- round(igraph::get.edge.attribute(g,"weight"),3) igraph::plot.igraph(g,layout=layout,edge.label=ew,vertex.size=vertex.size,vertex.color=vertex.color,...) } else igraph::plot.igraph(g,layout=layout,vertex.size=vertex.size,vertex.color=vertex.color,...) par(op) } ``` and then the plot ```{r fig.align="center", fig.width=3.5, fig.height=3.5} igplot(g1, weights=TRUE, layout=igraph::layout_as_tree) ``` For the weight-ordered Eulerian which follows high outlying edges first, we need to construct the graph whose edges are the negative of those in `g1`. It is not so straigtforward to do this using the graph package, ```{r eval=FALSE} e <- edgeMatrix(g1,duplicates=FALSE) ew <- eWV(g1,e) e <- matrix(nodes(g1)[e],ncol=2,byrow=TRUE) g2 <- ftM2graphNEL(e,-ew,edgemode="undirected") ``` and perhaps the following is easier: ```{r} g2 <- dn_graph(mk_complete_graph(-dOut),-.2 , test=`<=`) ``` Either way, the Eulerian and corresponding plot is ```{r fig.align="center",fig.width=5, fig.height=3} o <- eulerian(g2) o o <- match(o, names(data)) par(tcl = -.2, cex.axis=.86,mgp=c(3,.3,0)) guided_pcp(data,scOut, path=o,pcp.col=colOut,lwd=1.4, main= "Eulerian for high Outlying graph", bar.col = scag_cols["Outlying"], legend=FALSE,bar.axes=TRUE,bar.ylim=c(0,max(scOut))) ``` There is one panel in the above plot which has an outlier score below 0.2. This is SW-L. Looking back at the result of `eulerian(g2)`, the edge SW-L does not belong to `g2`, it was added simply to make the graph Eulerian. The case colours used in the PCP marked the outliers in the GP-L plot found previously, Humans are in red and the Asian Elephant in blue. It looks like these two cases are responsible for the high outlier scores in all panels except Bd-SW.