File a10_dimension reduction By Daniel B. Carr Copyright 2005, 2006, 2007, 2008 class use permitted Sections 1. Discussion of variables and dimension reduction 2. Simplex cluster generation 3. Classical metric multidimensional scaling 4. Principal Components 5. Digression: Passing arguments to functions used in function arguments 6. Heptathlon data and principal components 6.1 Look at the data 6.2 Use identify()to name points in a plot 6.3 Prepare for principal components and easier visualization 6.4 Principal components 6.5 Assess the variation in the principal components 6.6 Relate the first principal component to official score 6.7 Revisions from removing an outlier Due 3. three plots 4. two plot 6. plots from 6.2, 6.3, 6.5 and 6.7 1. Discussion of variables and dimension reduction This disussion cover more than is addressed in class, but provide context. Dimension reduction serves at least two purposes: 1) Developing computationally feasible and useful models 2) Developing cognitively accessible graphics, data descriptions and models. Contexts for dimension reduction mentioned here 1) Data density characterizatio 2) Supervisd classification or regression models In supervised clasifiction the response or dependent variable is categorical In regression the response or dependent variable is continuous Explanatory or predictor variiables These are the targets of dimension reduction These are often call independent variables and sometimes controlled in experimental design setting. Here the label independent does not mean linear independent. Many regression models treat explantory variables as measured without error. In numerous circumstance this is not a good approximaation to reality. The names of the problem varies accross discipline. Literature is availalble under the labels errors in variables and total least squares. Singular value decomposition is helpful for some very special cases. In general good solution are lacking and modelers procede by making assumption. If the model is sensitive to predictor variable measurement the errors do not have to be very large to make a big difference. Some statistical models assume that some preditor variables are random variables. These are called typically mixed models. Labels suggest alternative roles for predictor variables can have various labels Causal variables: Establishing causes can be hard but this is important to our thinking Surrogates variables Used in place of the desired variables Happenstance variables Variables possibly related to the response but not causally Confounded variables I think of these causal variables being correlated with a happenstance variable. Correlated variable Context should indicate whether this is to the response or other predictor variables. Irrelevant or noise variables: These complicate, contaminate or totally spoil modeling efforts A surrogate and confound variable example__________________ In modeling cancer incidence, mortality, or survivial curves data for known or conjectured risk factor is often not available. Sometimes we seek available data to serve as surrogates. For example I was look at chronic obstructive pulmonary disease, and smoking was a know risk factor. Of course there is no good geographically detailed data available on the loccal populatin distribution of doses in term of packs per day. As a surrogate I chose lung cancer mortality. I also conjected that particulates in addition to cigarette smoke and coals dust would also be a risk factor. The available pariculate data is not national in scope and had other limitation. In wishful thinking I thought that precipitation might periodically clear the air. I picked preciptation as a primitive inverse surrogate for particulates. It turned that these two surrogate variables lung cancer mortality and precipitation were highly correlated. I think that precipitation is confounded with cigarette smoking rates in the US. Both are high in the southern US were the cigarette industry has had a strong influence for decades. In the absence of good cigarette smoking data it is not posssible to assess if precipitation might be a related factor, for example by influence cigarette smoking rates. If data were able for other nations, it would be possible to look for a similar correlation. Ecological studies like that above often seek to see what variables are related. Such studies are totally inadequate for establishin cause and effort. Rather, on lucky occasions they suggest relations that may warrant further study. Variables may be also described in many other ways Composite variables are functions other variables A composite variable might be a linear combinations of transformed variables. Dimension reduction techniques include the following list Some are only application in the context of a predictor variables 1) Selecting known or presumed important variables 2) Choosing a representative variable for a set of correlated variables 3) Computing a representative variable from related variables For example averaging radiosity values from adjacent spectral bands to several bands to spectrally local summary vector. 4) Multidimensional scaling (MDS) There are metric, non-metric, and constrained variants of MDS 5) Linear projection into lower dimensional subspaces principal components projection pursuit 6) Nonlinear manifold finding 7) Clustering can be applied to variables as well as cases This may be used in conjunction with 2)to select a representative variable for each cluster of variables 8) In clustering contexts search algorithms may identify variable subspaces likely to have clusters and other variables are dropped. 9) In supervised classification or regression contexts there are several approaches to variable selection Forward, backward and stepwise regression Lasso regression: less greedy forward selection Semi-supervised principal components All-subsets regression Various criteria: Cp statistics Smarter search algorithms: leaps and bounds The list is not complete. The number of cases and variables can strongly influence the approach used. Comments Given a response variable for 100 cases and can generate few thousand random variables and develop a model based on noise that fits the respone variable well. Cross validation or validation on new data are often crucial more model believability. New data provide a better test. We can often borrow strenght from related variables to help us see throught or reduce noise. Variables are often correlated in space and time so local smoothing to produce new variables can be helpful. Methods using multidimensional scaling break down with too many variables. Computational error comes into play. It is not just a matter of as computingtime. Semi-supervised principal components was developed to address a small number of cases and thousands of variables. It combines more than one method from the list above. In the computational context there can sometimes be more than a million explanatory variables. A paper from Andrew Moore's web site addressed logistic regression for such situations. Since collecting data is relatively easy in many contexts, it is often easy to put together data sets that contain many variables. A significant proportion of the variables available my be noise with respect to the task at hand. In such cases classical dimension reduction methods such as principal components will contain linear combinations of extraneous variables. If an analyst is really unlucky the extraneous variables will strongly influence the direction of greatest variability. The Random Forests methodology developed by Leo Brieman and Adele Cutler offers some hope. The methodology samples both variables and cases to build thousands of models that vote. The hope is that a sufficient number of models are educated voters. Educated models are largely influenced by the structure in the data that is relevant to the task at hand. Study of the variables used and not used in the trees of the forest provides clues to the variables irrelevant to the task a hand. Rulefit converts random forest trees into indicator functions that can be used in a regression model along with variables used to build the trees. Lasso regression shring many of the coeffients to zero thus removing many predictor variables. 2. Multivariate Data Generation: The example generates 9 Clusters centered on simplex vertices in 8 dimensions. One cluster is centered on rep(0,8)=(0,0,0,0,0,0,0,0). The remaining 8 cluster centers are 8 sigma away in one of their coordinates The non-origin centers are 8*sqrt(2) sigma apart from each other 8*sqrt(2) = sqrt( sum ((c(8,0,0,0,0,0,0,0)-c(0,8,0,0,0,0,0,0))^2 )) ##Run library(MASS) n = 50 m0 = rep(0,8) samp = mvrnorm(n,mu=m0 ,Sigma=diag(8)) for ( i in 1:8){ m = m0 # All zeros m[i] = 8 # ith component 8 samp = rbind(samp,mvrnorm(n,mu=m,diag(8))) } ##End 3. Classical multidimensional scaling and distant neighbors in low dimensions Metric multidimensional scaling is often used transform a data set into a new data set with few variables. Often we pick two or three variables for plotting purpose. Given the number of variables, The distance matrix for cases using MDS variables is the best Frobenius norm approximation to the distance matrix for case using the original variables. The norm is the square of the sum of squares matrix elements ans is applied to the difference of the two matrices. The lower dimensional approximation may not be very good. Historically when I constructed distinct clusters in high dimensions and saw them overplotted in low dimensional MDS plots I though of of space folding over. However,the full set of MDS variables preserves the distance between points. Picking only two variables and plotting is simple projection points points by dropping coordiantes. Thus overplotting is simply due to orthogonal projection. The simplex based construction in 2) above is such that the plotting the first two MDS variables should reveal 2 distinct clusters. and one overplotted cluster composed of the other 7 distinct clusters. The metric multidimensional scaling (when computable) can make cases that are far apart in high dimensional space become near neighbors when forced onto the a 2D plane. ## Run Setup # # The script below used colors 4:12 # 4:9 from ColorBrewer # I use the first three colors # often so include them for other # examples # You might choose modify this # for example replacing a lighter # color with black colors = matrix(c( 0,0,0, # black 255,255,255, # white 230,230,230, # light gray 228,26,28, # red 1 55,126,184, # green- blue 2 77,175,74, # green 3 152,78,163, # purple 4 255,127,0, # orange 5 255,255,51, # yellow 6 166,86,40, # brown 7 247,129,191, # pink 8 153,153,153), # gray 9 ncol=3,byrow=T) colorCodes=rgb(colors[,1],colors[,2],colors[,3],max=255) lazyPlot = function(dat,subs,extraLabel){ windows() dat =dat[,1:2] # Control plot shape and axes to preserve 2D distances datR = range(dat) par(pty='s') plot(datR,datR,type='n',main='Multidimensional Scaling', xlab="X",ylab="Y") points(dat,pch=16,col=colorCodes[subs],cex=1.1) mtext(side=3,line=.3,extraLabel,cex=1.1) } ##End 3.1 First plot ##Run dis = dist(samp) # data clusters from 2. mds = cmdscale(dis,eig=T,k=4)$points # could put k=2 subs = rep(4:12,rep(50,9)) # for the cluster colors extraLabel="Colors Codes 9 Distinct Cluster in Higher Dimensions" lazyPlot(mds[,1:2],subs,extraLabel) ##End Note: here sigma = 1. The cluster coded red was sampled from a 8 dimensional standard normal distribution centered at the origin. The origin is 8*sigma from the other cluster centers. The other 8 clusters centers are 8*sqrt(2)*sigma apart so further from each other. Multidimensional scaling is doing its best to keep all 9 clusters a part but this is possible. The red cluster at the origin is closer to the others we could expect it to be heavily overplotted. In my first example yellow and red clusters were substantially overplotted. The orientation of the other clusters with respect to red is likely a fragile computation influenced by the sampling. Your plots may be different. ##Run lazyPlot(mds[,c(1,3)],subs,"MDS Variables 1 and 3") lazyPlot(mds[,c(2,3)],subs,"MDS VAriables 2 and 3") lazyPlot(mds[,c(1,4)],subs,"MDS Variables 1 and 4") ##End Other projections overplot red clusters. Cluster overplotting the red are not all the same. 3.2 Adding a point at the simplex centroid. sampCenter = rbind(samp,rep(1/8,8)) disCenter = dist(sampCenter) mdsCenter = cmdscale(disCenter,eig=T,k=2)$points subsCenter=c(subs,1) extraLabel="A Single Simplex Center Point (Black) Claims the Center" lazyPlot(mdsCenter,subsCenter,extraLabel) 3.3 What to extreme outlier do ## Run sampOutlier = rbind(samp, c(12,0,0,0,0,0,0,0), c( 0,0,0,12,0,0,0,0), c( 0,0,0,0,0,0,0,12)) disOutlier = dist(sampOutlier) mdsOutlier = cmdscale(disOutlier,eig=T,k=2)$points subsOutlier=c(subs,1,1,1) extraLabel="Outliers (Black) can make some clusters relatively closer" lazyPlot(mdsOutlier,subsOutlier,extraLabel) ## End Note The blue cluster center has a large 1st coordinate The gray cluster center had a large 8th coordiante The orange cluster center had a large 4th coordinate 3.4 Other multidimensinal scaling approaches There are non-metric multidimensional scaling algorithms. The Kruskal algorithm attempts to preserve the ranks of the dissimilarities in the dissimilarity matrix and minimizes a particular stress function. Sum (Dij -Pij)^2 /Pij^2 ij where i != j, Dij are the original distances between case and Pij are the new 2D coordinate distances between cases. This is available is as isoMDS() in the MASS libary. In recent times various investigator have chose different stress functions to minimize. Another scaling function for 2D layouts in the MASS library in Sammon's non-linear mapping sammon() 4. Principal components and strange neighbors Principal components are often used to reduce the dimension of a dataset. The principal components transforms data matrix as follows. 1) It centers the variables at the origin 2*) By default prcomp bases the rotation on the covariance matrix of the data. If the variables are in roughly the same units there is merit in using the covariance matrix. If the variables are in substantially different units the common choice is to scale the data and use the correlation matrix. Setting the scale option to TRUE divides the variables by their standard deviation hence based the rotation on the correlation matrix. prcomp(x, scale=T) 3) Rotates the data so a projection into 1st dimension has the largest variance. This the first principal component and is often shown on the x axis. Then the transformation picks the orthonormal linear combination relative to the proceding linear combinations that gives the next largest variance and so on. Principal components or "scores" are the rotated variables. The rotation or "loading" matrix contains the linear combinations that when applied to the centered (and often scaled) data yields the principal components. The rotation values are linear weights that are needed to develop an interpretation of the principal components in terms of the centered original variables. Using the first few principal components reduces the dimensions. The covariance matrix for normalized variables is the correlation matrix. Specifying this causes this normalization to occur. 4.1 Check on the computations ##Run prin = prcomp(samp) samp.center = scale(samp,T,scale=F) # center the data samp.rotate = samp.center%*% prin$rotation all.equal(prin$x,samp.rotate) # compare the stored principal components # in prin$x to the calculation # 4.2 Plot the nine groups based on the first two principal components ##Run windows() par(pty='s') dat = prin$x[,1:2] datR = range(dat) plot(datR,datR,type='n',main="First Two Principal Components") points(dat,pch=16,col=colorCodes[subs]) # End # 4.3 The fraction of variability represented. ##Run windows() names(prin$sdev) = paste('Comp',1:length(prin$sdev),sep='') plot(prin) compVar = prin$sdev^2 # Fraction of total variance 100*compVar/sum(compVar) # Cumulative fraction of total variance 100*cumsum(compVar)/sum(compVar) ##End plot() or screeplot() shows the variance associated with the principal components. A drop to a much smaller variance is sometimes used as an indicator of how many principal components to keep. In the above example the big drop is after 7 components. Dropping just the last component does not provide much in the way of variable reduction. Fairly often people look at just the first 2 components because they are easy to plot. People like look at different pairs rather than looking at 3D or 4D plots. Some disciplines such as meteorology make heavy use of principal components for dimension reduction. Dianne Cook at Iowa State has examples using grand tour and projection pursuit where the crucial information turns up in the component with the least variance. Principal components are often used as a hammer but sometimes it is the wrong tool for the problem. 5. Digression on passing arguments to functions___________________ for use as arguments in functions passed. Some functions have arguments that are functions. In this context call such functions argFunctions. The question is how to pass arguments to argFunctions. 5.1 A single argfunction and the ... notation. In defining a function, the ... notation is often used to indicate one or more arguments to be passed into a function. where they can be used as argument to an argfunction. Consider apply(mat,1,sum,na.rm=T) This will apply the sum() function to each row of the matrix mat. However, some values in a row may be missing. The sum() function has an argmuent na.rm to address the treatment of missing vlues. If na.rm=T, sum() will omit missing values. The argument list of apply is apply = function (X, MARGIN, FUN, ...) In the body of the apply script there is something like FUN(stuff,...) For apply(mat,1,sum,na.rm=T) this will become sum(stuff,na.rm=T) 5.2 What about functions multiple argFunctions. The situation gets sticky since ... arguments many not be appropriate argument to all of the argfunctions. A work around is to define a new argfunctions with some of their arguments set to the desire values. For example pairs has five argFunctions panel = points upper.panel=panel lower.panel=panel diag.panel=NULL text.panel=textPanel The default functions absent the parentheses appear on the right. Rather than using the points() function we want to use panel.smoooth to plot points and smooth line. This has argument col.smooth that controls the color of the smooth line, but this cannot be passed via ... because it is invalid argment to some of the other argfunctions. If we want the line to be orange was can define a new function with the col.smooth argument already set to orange. ## Run myPanelSmooth = function(x,y,...)panel.smooth(x,y,col.smooth='orange',...) ## End The pairs function further below will myPanelSmooth() instead of panelSmooth() so line will orange. 6. First Principal component and official system scores______________ The example is adapted from "A Handbook of Statistical Analysis Using R" by Everitt and Hothorn and using in my previous class. Get the heptathlon data from the assigment 8 ## Run heptathlon = read.csv(file="heptathlon.csv",row.names=1) ##End 6.1 Look at the data__________________________________________________ ## Run #! uses myPanelSmooth defined in 5.) so lines will be orangea pairs(heptathlon,panel=myPanelSmooth,lwd=2,pch=21,cex=1.5, col="blue",bg='blue',gap=0, main="Original Data") ## End Note that some smooths fit the data pretty well. Several are are monotone increasing and several are monotone decrease Note that there is an outlier, for example in the lower left panel at the lower right side. 6.2 Use identify()to name points in a plot___________________________ Here is about mousing before making the plot. Left mouse click near a point for the name to appear. Clicking slightly on the left of the point will cause it to appear right adjusted on the left. Clicking slightly on the right of the point will cause it to appear left adjust on the right. Find the person's name for point at the bottom right of the next plot, and a few more names. Select stop from inside the stop menu at the top left of the window or Right mouse click to get a menu and select stop. ## Run plot(heptathlon[,1:2],pch=16,col='blue',cex=1.2) identify(heptathlon[,1],heptathlon[,2],lab=rownames(heptathlon)) # End 6.3 Prepare for principal components and easier visualization________ Remove the official score. We will compare the first pricipal against this score. Transform variables. Initially the values for hurdles, run200m and run800m are in seconds and the values are highjump, shot, long jump and javelin are meters. A good start to simplify appearanceis to make all the positive outcomes have larger values. Hence the following transformations modify timed outcomes so large values are positive. The outlier will mentioned about will be removed later to see the improved resolution of the values for the other athletes and to see if makes much difference in comparing the first principal component with the official scores. ## Run # Remove the official score hepDat = heptathlon[,-ncol(heptathlon)] # transform values for timed events hepDat$hurdles=max(hepDat$hurdles)-hepDat$hurdles hepDat$run200m=max(hepDat$run200m)-hepDat$run200m hepDat$run800m=max(hepDat$run800m)-hepDat$run800m pairs(hepDat,panel=myPanelSmooth,lwd=2,pch=21,cex=1.5, col="blue",bg='blue',gap=0, main="Transformed So Large Values Are Good") ## End 6.4 Principal components____________________________________________ Pricipal components centers the data by subtracting the average from each variable. When variable are not in the same units, the common approach is to base the rotation on the correlation matrix rather than the covariance matrix. This is equivalent to converting the variables to z scores and using the resulting covariance matrix. Dividing a centered variable's values by the it's standard deviation converts it to a score. ## Run # scale = TRUE converts centered columns to z scores heptPca = prcomp(hepDat,scale=TRUE) # centered=TRUE by default heptPca$center # averages that were subtracted heptPca$scale # standard deviation of the variables heptPca$x # principal components heptPca$sdev # the standand deviation of the principal components round(heptPca$rotation,2) ## End Look at the rotation columns. The javelin value in the second column, -.84, sticks out as having a large magnitude. Throwing the javelin has physical requirements less compatitible with the training for other events. 6.5 Assess the variation in the principal components______________ ## Run names(heptPca$sdev) = paste('Comp',1:length(heptPca$sdev),sep='') plot(heptPca) heptVar = heptPca$sdev**2 100*cumsum(heptVar)/sum(heptVar) ## End The first principal component accounts for about 64% of the variation. 6.6 Relate the first principal component to official score ## Run windows() plot(heptathlon$score,heptPca$x[,1], xlab="Official Score",ylab="First Principal Component", main="Comparing Official Scores and First Principal Component") correl = cor(heptPca$x[,1],heptathlon$score) dataLimits = par()$usr xloc = mean(dataLimits[1:2]) yloc = dataLimits[4]-.05*diff(dataLimits[3:4]) text(xloc,yloc,paste("Correlation =",round(correl,2)),adj=.5,cex=1.1) ## End 6.7 Revisions from removing an outlier # Add the names rownames(hepDat) = rownames(heptathlon) ## Run plot(hepDat[,1:2],pch=16,col='blue',cex=1.2, main="Transformed variables so large values are good") # Find the subcript for the outlier in the lower left # identify always returns subscripts # It it use subscript as laels if labels if they aren't given. outlierSub = identify(hepDat[,1],hepDat[,2]) outlierSub # should be 25 # End hepDatV2 = hepDat[-outlierSub,] hepV2Pca = prcomp(hepDatV2,scale=TRUE) # variables are in different units hepV2Pca$center hepV2Pca$scale # standard deviation of the variables hepV2Pca$x # principal components hepV2Pca$sdev # the standand deviation of the principal components round(hepV2Pca$rotation,2) names(hepV2Pca$sdev) = paste('Comp',1:length(hepV2Pca$sdev),sep='') plot(hepV2Pca) hepV2Var = heptPca$sdev**2 100*cumsum(hepV2Var)/sum(hepV2Var) windows() plot(heptathlon$score[-outlierSub],hepV2Pca$x[,1], xlab="Official Score",ylab="First Principal Component", main="Comparing Official Scores and First Principal Component Ver2") correl = cor(hepV2Pca$x[,1],heptathlon$score[-outlierSub]) dataLimits = par()$usr xloc = mean(dataLimits[1:2]) yloc = dataLimits[4]-.05*diff(dataLimits[3:4]) text(xloc,yloc,paste("Correlation =",round(correl,2)),adj=.5,cex=1.1) ## End In comparing principal components against the offical scores, removing the outliers change the ranges for the x and y about the same. A point for one person moves further from a straight line fit. There are some other slight changes. A better plot would be to remove the one case from the first principal component and look at the differnce directly.