散布図行列

順序カテゴリカルデータの視覚化方法について
スクリプトまで紹介しているすごい便利なホームページがあったので
使ってみた

Correlation scatter-plot matrix for ordered-categorical data | R-statistics blog
scratch-R: scatterplots

ver1:散布図と相関係数

#散布図行列の設定を記憶
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr)) 
  par(usr = c(0, 1, 0, 1)) 
  r <- round(cor(x, y),digit=2)
  txt <- format(c(r, 0.123456789), digits=digits)[1] 
  txt <- paste(prefix, txt, sep="") 
  if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
  test <- cor.test(x,y) 
  # borrowed from printCoefmat
  Signif <- symnum(test$p.value, corr = FALSE, na = FALSE, 
                   cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                   symbols = c("***", "**", "*", ".", " ")) 
  text(0.5, 0.5, txt, cex = cex*r ) 
  text(.8, .8, Signif, cex=cex*r , col=2)  
#  text(0.5, 0.5, txt, cex = 4 ) 
#  text(.8, .8, Signif, cex=3, col=2)
}

panel.hist <- function(x, ...) # ヒストグラムを描くための関数を定義
{ usr <- par("usr") # 現在のユーザー領域座標情報を得る
  on.exit(par(usr)) # 関数終了時に usr パラメータ復帰
  par(usr = c(usr[1:2], 0, 1.5) ) # ユーザー領域座標情報を変更
  h <- hist(x, breaks=30,plot = FALSE)
  breaks<-h$breaks
  nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], col = "#0050ff90", border = "#0000ff00", y, ...)
}

panel.lsfit <- function(x,y,...) {
  f <- lsfit(x,y)$coef
  xx<- c ( -99999 , 99999 )
  yy <- f["X"] * xx+ f["Intercept"]
  points(jitter(x),jitter(y), pch = ".", col = "#0000ff50",bg = "#0050ff50")
  lines(xx,yy,col="#00ff0090",lwd=2)
}

以下は使用例

library(psych)
data(bfi)
opar<-par()
par(family="serif")
par(bg="black",fg="white",col.axis="white",col.lab="white",col.main="white",col.sub="white",family="serif")
pairs(na.omit(bfi[1:5]),
      diag.panel=panel.hist,
      lower.panel=panel.cor, 
      upper.panel=panel.lsfit)
par(opar)

ver2:相関係数と偏相関係数

#散布図行列の設定を記憶
#散布図行列の設定を記憶
panel.cor <- function(x, y, digits=2, prefix="", cex.cor, ...)
{

  points(jitter(x),jitter(y), pch = ".", col = "#00ff0050",bg = "#00ff5050")
  
  r <- round(cor(x, y),digit=2)
  txt <- format(c(r, 0.123456789), digits=digits)[1] 
  txt <- paste(prefix, txt, sep="") 
  if(missing(cex.cor)) cex <- 0.8
  text( 3,3 , sprintf("cor=%s",txt), cex = cex*3 ) 
}

panel.hist <- function(x, ...) # ヒストグラムを描くための関数を定義
{ usr <- par("usr") # 現在のユーザー領域座標情報を得る
  on.exit(par(usr)) # 関数終了時に usr パラメータ復帰
  par(usr = c(usr[1:2], 0, 1.5) ) # ユーザー領域座標情報を変更
  h <- hist(x, breaks=30,plot = FALSE)
  breaks<-h$breaks
  nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], col = "#00ff5090", border = "#00ff0000", y, ...)
}


panel.pcor <- function(x, y, data=na.omit(bfi[1:5]), digits=2, prefix="", cex.cor, ...)
{
  X<-as.matrix(data)
  tmp<-X-x
  tmp<-tmp[,-which(colSums(tmp)==0)]
  tmp<-tmp+x
  tmp<-tmp-y
  tmp<-tmp[,-which(colSums(tmp)==0)]
  tmp<-tmp+y
  X<-tmp
  X<-cbind(1,X)
  fx<-(diag(nrow(data)) - (X %*% solve( t(X) %*% X ) %*% t(X)) ) %*% x
  fy<-(diag(nrow(data)) - (X %*% solve( t(X) %*% X ) %*% t(X)) ) %*% y
  
  points(fx+3,fy+3, pch = ".", col = "#00ff0050",bg = "#00ff5050")
  
  r <- round(cor(fx, fy),digit=2)
  txt <- format(c(r, 0.123456789), digits=digits)[1] 
  txt <- paste(prefix, txt, sep="") 
  if(missing(cex.cor)) cex <- 0.8
  text( 3,3 , sprintf("pcor=%s",txt), cex = cex*3 ) 
}

以下は使用例

library(psych)
data(bfi)
opar<-par()
par(family="serif")
par(bg="black",fg="white",col.axis="white",col.lab="white",col.main="white",col.sub="white",family="serif")
pairs(na.omit(bfi[1:5]),
      diag.panel=panel.hist,
      lower.panel=panel.cor, 
      upper.panel=panel.pcor)
par(opar)