## The partial regression-on-predecessors scatterplot matrix defined in section 4.3.4 "partial.r"<-function(x=matrix(rnorm(120),20,6),...){ # For testing purposes only (set to F if you don't need "2 vs. 1" labels): showlabels <- F # Make sure x is in matrix form x <- as.matrix(x) n <- ncol(x) # Make sure we clean up our mess (revert back to original "par" parameters) # when we're done oldpar <- par("pty", "oma", "mar", "cex", "tck", "mgp", "mex", "mfrow") on.exit({par(oldpar)}) # Set layout, size and margin parameters # make plots square par(pty="s") # determine new optimal plotting character size CEX <- par("cex") * max(7.7/(2*(n-1)+3), 0.6)*0.5 # specify grid of (n-1) by (n-1) plots par(mfcol = c(n-1,n-1)) # magically cause plots to pack together nicely par(oma = rep(3,4)) par(mar = rep(0.3,4)) dif <- diff(par("fin"))/2 if(dif > 0) par(omi = c(dif*(n-1), 0, dif*(n-1), 0) + par("omi")) else par(omi = c(0, (-dif)*(n-1), 0, (-dif)*(n-1)) + par("omi")) # re-specify plotting character size (since par(mfrow) screwed it up) par(cex = CEX) # specify the points type par(pch=16) dat <- x m <- n-1 mat.x <- NULL mat.y <- NULL for (k in (1:m)){ fres <- sapply(1:(m+1-k), function(i) if(k==1 && i==1) dat[,i] else residuals(lm(dat[,i] ~ dat[,c(1:(i+k-1))][,-c(i)]))) bres <- sapply(1:(m+1-k), function(i) if(k==1 && i==1) dat[,i+k] else residuals(lm(dat[,i+k] ~ dat[,c(1:(i+k-1))][,-c(i)]))) mat.x <- cbind(mat.x,fres) mat.y <- cbind(mat.y,bres) } locmat <- matrix(rep(0,m^2),nrow=m) count <- 0 for (i in (1:m)){ for (j in (1:(m+1-i))){ count <- count+1 locmat[i,j] <- count } } # Make plots (i from top to bottom, j from left to right) for (i in (2:n)){ for(j in (1:(n-1))){ if (j < i-1){ # prepare an empty plot with correct limits and no axes plot(range(x[!is.na(x[,j]),j]), range(x[!is.na(x[,i]),i]), type = "n", axes = F, ...) # create empty plot and put border around it box() # show points for j'th vs. i'th columns of x points(as.vector(x[,j]),as.vector(x[,i]), ...) # show which variables are being plotted if (showlabels) text(mean(range(x[!is.na(x[,j]),j])), mean(range(x[!is.na(x[,i]),i])), paste(j,"vs.",i),cex=1.5*CEX) } else{ # Get the information of which variables to be plotted against each other on this grid loc <- locmat[j-i+2,i-1] # prepare an empty plot with correct limits and no axes plot(range(mat.x[!is.na(mat.x[,loc]),loc]), range(mat.y[!is.na(mat.y[,loc]),loc]), type = "n", axes = F, ...) # create empty plot and put border around it box() # show points for j'th vs. i'th columns of x points(as.vector(mat.x[,loc]),as.vector(mat.y[,loc]), ...) } } } }