Rで交互作用項を含めた相関係数のカラーマップを描く

はじめに

このブログのアイコンにも使っている相関係数のカラーマップ（上図）をRを使って描く方法について、読者の方から質問を受けましたので、説明します。

表1のような実験計画を作成した際、交互作用項や2乗項まで含めて相関係数のカラーマップを作成します。

コード

まず、以下の二つのコードをファイルで保存してください。

①要因A,B,C,...からなる計画の交互作用項と2乗項を計算するファイル#add_model_terms.R

②渡されたdataframeから相関係数のカラーマップを出力するファイル#my_corr.R

①の交互作用項と2乗項を追加するファイルは以下の通りです。

#add_model_terms.Rという名前で保存

add_quadratic <- function(df, only_additional=FALSE, only_quad=FALSE){
  k <- ncol(df)
  #2-factor interaction
  cnt <- 0
  if (!only_quad) {
    for (i in 1:k) {
      if (i<k) {
        for(j in (i+1):k){
          a <- colnames(df)[i]
          b <- colnames(df)[j]
          cnt <- cnt + 1
          df <- data.frame(df, df[a]*df[b])
          colnames(df)[k+cnt] <- paste0(a,b)
        } 
      }
    }
  }
  #quadratic
  for (i in 1:k) {
    a <- colnames(df)[i]
    cnt <- cnt + 1
    df <- data.frame(df, df[a]*df[a])
    colnames(df)[k+cnt] <- paste0(a,a)
  }
  if (only_additional) {
    return(df[(k+1):ncol(df)])
  } else{
    return(df)
  }
}


add_interaction <- function(df, only_additional=FALSE){
  k <- ncol(df)
  #2-factor interaction
  cnt <- 0
  for (i in 1:k) {
    if (i<k) {
      for(j in (i+1):k){
        a <- colnames(df)[i]
        b <- colnames(df)[j]
        cnt <- cnt + 1
        df <- data.frame(df, df[a]*df[b])
        colnames(df)[k+cnt] <- paste0(a,b)
      } 
    }
  }
  if (only_additional) {
    return(df[(k+1):ncol(df)])
  } else{
    return(df)
  }
}

②の相関係数のカラーマップを描くファイルは以下の通りです。

#my_corr.Rと言う名前で保存

library(RColorBrewer)
library(corrplot)


mycorr <- function(df, coef=FALSE){
  df.cor = cor(df, method = c("pearson"))
  palette = colorRampPalette(c("blue", "white", "black")) (20)
  
  acor <-  abs(df.cor)
  if (coef) {
    corrplot(acor, method="color", col=palette,  
             type="full", order="original", 
             addCoef.col = "black", # Add coefficient of correlation
             tl.col = "black", tl.srt=90, #Text label color and rotation
             cl.lim = c(min(acor),max(acor)),
             diag=TRUE 
    )
  }else{
    corrplot(acor, method="color", col=palette,  
             type="full", order="original", 
             #addCoef.col = "black", # Add coefficient of correlation
             tl.col = "black", tl.srt=90, #Text label color and rotation
             cl.lim = c(min(acor),max(acor)),
             diag=TRUE 
    )
  }
  return(acor)
}

①と②のファイルを保存したら、ある計画を作成した際は次のようなコードを実行します。ちなみに、別ファイルにしなくてもいいですが、私は毎回コピペするのではなく、別ファイルの関数として呼び出してます。

#計画作成（例は決定的スクリーニング計画）
library(daewr)
X <- DefScreen(m=8, c=0)

#交互作用項と2乗項の追加
source("../add_model_terms.R")
X <- add_quadratic(X)

#相関のヒートマップの作成
source("../mycorr.R")
corr <- mycorr(X, coef = FALSE)#相関係数を表示する場合はcoef=TRUE