TypstとRでAPAスタイルの相関表を作成する

TypstとRの連携がどれぐらいできるかいろいろ試してるところです。

Rの分析結果から相関の表を作ろうと思ったところ、意外と苦戦したので、備忘録がてらその方法を記録しておきます。

Rでの準備

データの準備

まず、分析データはこちらの表をもとに架空のものを作成しました。最初は相関係数ができるだけ似たような数値になるようにと思ったんですが、わざと欠損値を入れたりしているせいもあって、必ずしもこの通りの値にはなっていません。ただ、値そのものはどうでもいいので、ここではこれで良いこととします。

また、この通りのデータである必要は全然ないので、手持ちのデータがあればそれを使うのがいいと思います。

R

library(MASS)     # 擬似データ生成用

set.seed(123)  # 再現性のため

# 変数名と相関行列の定義
# 相関行列は"https://apastyle.apa.org/style-grammar-guidelines/tables-figures/sample-tables"から
var_names <- c("IE", "MP", "Salary", "Promotion", "Tenure", "USP", "UFP")

# 相関行列の作成
cor_matrix <- matrix(c(
  1.00, -0.08,  0.45,  0.08, -0.29, -0.25,  0.00,
  -0.08,  1.00, -0.01,  0.07,  0.09, -0.39, -0.03,
  0.45, -0.01,  1.00,  0.04,  0.01, -0.24,  0.12,
  0.08,  0.07,  0.04,  1.00,  0.09,  0.08, -0.07,
  -0.29,  0.09,  0.01,  0.09,  1.00,  0.01, -0.02,
  -0.25, -0.39, -0.24,  0.08,  0.01,  1.00,  0.16,
  0.00, -0.03,  0.12, -0.07, -0.02,  0.16,  1.00
), nrow=7, byrow=TRUE)

colnames(cor_matrix) <- rownames(cor_matrix) <- var_names

# 各変数の平均と標準偏差
means <- c(0.43, 3.14, 1.01, 0.33, 6.45, 85.00, 42.61)
sds   <- c(0.49, 0.62, 0.27, 0.47, 6.62, 6.98, 5.86)

# 共分散行列を作成
cov_matrix <- diag(sds) %*% cor_matrix %*% diag(sds)

# データ数（nは最大の3697に合わせる）
n <- 3697

# 多変量正規分布に基づくデータ生成
data <- as.data.frame(mvrnorm(n = n, mu = means, Sigma = cov_matrix))
colnames(data) <- var_names

# 欠損を追加
data$MP[2135:n] <- NA
data$UFP[695:n] <- NA

# 変数の調整
# IE: 0 or 1 の二値変数に変換
data$IE <- ifelse(data$IE > 0.5, 1, 0)

# Promotion: 0 or 1 に変換（昇進の有無）
data$Promotion <- ifelse(data$Promotion > 0.5, 1, 0)

# 値の切り捨てや補正
data$Salary <- round(data$Salary, 2)
# パフォーマンススコア 1〜5に制限
data$MP <- pmin(pmax(round(data$MP, 2), 1), 5)

統計量の算出

まずは基本統計量です。

Rスクリプト

# 各変数の有効サンプルサイズ（NA除外数）
valid_n <- colSums(!is.na(data))

# 各変数の平均値（NA除外）
var_mean <- apply(data, 2, mean, na.rm = TRUE)

# 各変数の標準偏差（NA除外）
var_sd <- apply(data, 2, sd, na.rm = TRUE)

相関係数

次に、相関係数です。相関係数の算出はbaseモジュールでもできますが、ここではpsychパッケージのcor.test()を用いることにします。また、有意性のマークを処理するためにdplyrも使用します。

Rスクリプト

library(psych)    # 相関分析用
library(dplyr)    # データ整形用

# 相関分析
result <- corr.test(data)
cor_mat <- round(result$r, 3) # 小数第3位まで

表示用データの作成

分析そのものはこれで終わりですが、この結果をTypstで利用できるようにするため、結果の整形をここで行っておきます。Typstの方でも可能かもしれませんが、これらの処理はRでやっておく方がいいでしょう。

有意性マーカー

この例のようにサンプルサイズが大きい場合、相関係数の有意性検定に意味があるとは思えませんが、一般に相関係数が有意な場合には、数字の右にアスタリスクをつけて示すことが一般的なので、その処理をここで行います。Typstの方でエラーになるのを避けるために、有意でないところには半角空白を1つ入れるようにしてあります。

R

# 有意確率のマーカー
p_value <- as.numeric(result$r)
sigmark <- case_when(
  p_value < .001 ~ '#super("***")', # 表示を上付きにする
  p_value < .01  ~ '#super("**")',
  p_value < .05  ~ '#super("*")',
  TRUE     ~ " "
)
sig_table <- matrix(
  sigmark,
  nrow = nrow(result$r), ncol = ncol(result$r)
)

ABAの書式では、相関係数のように値が1.00を超えることがない数値については、整数位の0を表示しないのが一般的です。そのための処理をここで行います。また、相関係数が負の値の場合にマイナス記号が正しく表示されるようにもします。

R

# 相関係数の整数位の0を削除
cor_str <- sprintf("%.3f", result$r) # 小数第3位まで
cor_str <- sub("0.", ".", cor_str)   # 0を削除
cor_str <- sub("-", "#sym.minus;", cor_str)   # ハイフンをマイナス記号に置換
cor_table <- matrix(
  cor_str,
  nrow = nrow(result$r), ncol = ncol(result$r),
  dimnames = dimnames(result$r)
)

このままだと行列の上三角と下三角に同じ数値が表示されるので、行列の上三角を空欄にします。また対角成分は「—」に置き換えます。

R

cor_table[upper.tri(cor_table)] <- " " # 上三角を" "に
sig_table[upper.tri(sig_table)] <- " " # 上三角を" "に

diag(cor_table) <- rep("---", nrow(cor_table)) # 相関係数は対角成分を"---"に
diag(sig_table) <- rep(" ", nrow(sig_table)) # 有意確率は対角成分を" "に

変数名に簡略化した文字列を使っている場合は、本文中では正しい変数名で表示されるよう、変数名もデータ化しておきます。

R

# 変数名ラベル
var_lst <- list(
  "IE" = "社内外採用区分",
  "MP" = "管理職パフォーマンス", 
  "Salary" = "初任給", 
  "Promotion" ="昇進の有無", 
  "Tenure" = "在職年数", 
  "USP" = "部門サービス指標", 
  "UFP" = "部門財務指標"
)
# 相関行列の行名を変数名ラベルに変換
varnames <- unname(sapply(rownames(result$r), function(x) var_lst[[x]]))

これで分析に関する書類は終わりです。

JSONにエクスポート

分析結果をリストにまとめてJSON形式で書き出します。

R

# 必要なデータをJSONファイルに保存

library(jsonlite)

out_json <- list(
  names = varnames,      # 変数名のリスト
  n = valid_n,           # 各変数のサンプルサイズ
  mean = var_mean,       # 各変数の平均値
  sd = var_sd,           # 各変数の標準偏差
  cor_table = cor_table, # 表示用相関係数行列
  sig_table = sig_table, # 有意確率マーカー行列
  cor_mat = cor_mat      # 数値が必要になる場合に備えて念の為
)

write_json(out_json, paste0("rdata.json"),
           pretty = TRUE, auto_unbox = TRUE)

Typstの準備

データの整形

ここからはタイプストでの作業です。まず必要なパッケージを読み込み、文書の言語設定を行います。ここでは、数値表示形式の処理用にzeroパッケージ、表の見た目を整えるためにbooktabsパッケージを用いています。言語は日本語で、デフォルトのフォントはSource Han Serifにしてあります。

Typst

#import "@preview/zero:0.5.0": *
#import "@preview/booktabs:0.0.4": *
#show: booktabs-default-table-style

#let jfont = "Source Han Serif"
#set text(lang: "ja", font: jfont)

データの読み込みは次の1行だけです。

Typst

#let rdata = json("rdata.json")

次に、読み込んだデータを表で扱いやすいように整えます。

まず、相関行列の相関係数と有意性検定結果のマーカーをつなぎ合わせます。

Typst

#let soukan = rdata.cor_table.flatten().zip(rdata.sig_table.flatten()).chunks(rdata.n.len())

全部1行に書いているので、かなりごちゃごちゃして見えますが、ここでやっているのは次の処理です。

1. rdataのcor_table（相関係数）を取り出し、flatten()で1次元配列にする
2. それぞれの要素（相関係数）と、対応する有意性マーカー（これもflatten()しておく）をzip()でペアにする
3. (相関係数1, 有意性マーカー1), (相関係数2, 有意性マーカー2), ...という配列ができるので、chunk()で表の1行分ずつ折り返す

これで創刊の表の相関係数の部分のデータが整いました。このデータをサンプルサイズや平均値等の記述統計量のデータと連結します。ここではrange()で変数に番号をつけたうえで、zip()しています。なお、タイプストの配列は通常0から始まりますが、ここでは番号づけを1から始めるため、終了番号に+1しています。

Typst

#let table_data = range(1, rdata.names.len()+1).zip(rdata.names, rdata.n, rdata.mean, rdata.sd, soukan)

ここまでの処理で、このテーブルデータに含まれるデータは1行あたり次のような内容になります。

(変数番号, 変数名, サンプルサイズ, 平均値, 標準偏差, 変数1との相関係数, 有意性マーカー, 変数2との相関係数, 有意性マーカー,...)

表の整形

データが整ったので、今度は作成する表のフォーマットに関する設定を行います。

位置揃えの設定

表の中の要素がきれいに揃うよう、位置揃えの設定をしておきます。変数の個数が変わった場合でも対応できるようにしてみました。

Typst

#let table_align = (right, left, auto, auto, auto) // 変数番号、変数名、N、平均値、SD
#table_align.push( // 変数の数だけ（right, left）をつなげる
  rdata.names.map(x => (right, left))
)

ここでは相関係数を右寄せに、有意性マーカーを左寄せに設定してあります。

セル間のスペースの設定

そのままでは相関係数とマーカーの間がずいぶん離れてしまうので、セル間のスペースを調整します。ここでは、相関係数とマーカーの間を−9ptにして寄せました。

Typst

#let col_spaces = (auto, auto, auto, auto, auto) // 変数番号、変数名、N、平均値、SD
#col_spaces.push( // 変数の数だけ（-8pt, 0pt）をつなげる
  rdata.names.map(x => (-9pt, 0pt))
)

ヘッダ行の設定

相関係数の部分のヘッダに変数番号が表示されるようにします。ここではtable.cell(colspan:2, align: center)[#(x+1)]として、横2セル（相関係数と有意性マーカー）の部分を繋ぎ、中央寄せにしています。

Typst

#let header = range(0, rdata.names.len()).map(x => table.cell(colspan:2, align: center)[#(x+1)])

列数のカウント

Typstで表を作成する際に列数を指定する必要があるので、列数を数えておきます。

Typst

#let n_of_columns = table_data.at(0).flatten().len()

表の作成

これですべての準備が整ったので、後は表を作成します。3桁以上の数値で千の位に位取りのコンマを表示するために、zeroパッケージを用いています。

Typst

#set-group(separator: ",", threshold: 4) // 位取りコンマの設定（zeroパッケージ）
#block[
  #set text(7pt)
  #ztable(
    columns: n_of_columns,
    align: table_align.flatten(),
    column-gutter: col_spaces.flatten(),
    format:(none, none,
      (digits:0, math: false), 
      (digits:2, math: false), 
      (digits:2, math: false)),
      toprule(),
      table.header(
        table.cell(colspan:2, align: center)[変数],
        table.cell(colspan:1, align: center)[n],
        [平均値],[SD],
        ..header,
      ),
      midrule(),
      ..table_data.flatten().map(x => eval(str(x), mode: "markup")
    ),
    bottomrule()
  )
]

表に必要なヘッダやデータはすべて事前に整形してあるので、ここでは単に..でその内容を展開するだけで表が作成されます。万が一データに修正が入った場合でも、Rの分析スクリプトを実行し直してjsonを更新すれば、表の内容は自動で更新されます。変数の個数が変化してもきちんと対応できるはずです。

使用する表が1つしかない場合、ここまでやるのは手間がかかりすぎて正直ばかばかしく思えるかもしれませんが、表が複数ある場合には、こうしたスクリプトを1つ作成しておくと、使い回しができるので効率も上がるのではないかと思います。

なお、今のところzeroパッケージでは整数位の0を非表示にできないため、表示設定などでちょっとめんどくさいことになっていますが、今後0を非表示にできるオプションが追加された場合には、もっとシンプルに作成できるようになることでしょう。

サンプルスクリプト

最後に使用したRおよびTypstのファイルの内容を掲載しておきます。

R

library(MASS)     # 擬似データ生成用
library(jsonlite) # JSON用
library(psych)    # 相関分析用
library(dplyr)    # データ整形用

set.seed(123)  # 再現性のため

# 変数名と相関行列の定義
# 相関行列は"https://apastyle.apa.org/style-grammar-guidelines/tables-figures/sample-tables"から
var_names <- c("IE", "MP", "Salary", "Promotion", "Tenure", "USP", "UFP")

# 相関行列の作成
cor_matrix <- matrix(c(
  1.00, -0.08,  0.45,  0.08, -0.29, -0.25,  0.00,
  -0.08,  1.00, -0.01,  0.07,  0.09, -0.39, -0.03,
  0.45, -0.01,  1.00,  0.04,  0.01, -0.24,  0.12,
  0.08,  0.07,  0.04,  1.00,  0.09,  0.08, -0.07,
  -0.29,  0.09,  0.01,  0.09,  1.00,  0.01, -0.02,
  -0.25, -0.39, -0.24,  0.08,  0.01,  1.00,  0.16,
  0.00, -0.03,  0.12, -0.07, -0.02,  0.16,  1.00
), nrow=7, byrow=TRUE)

colnames(cor_matrix) <- rownames(cor_matrix) <- var_names

# 各変数の平均と標準偏差
means <- c(0.43, 3.14, 1.01, 0.33, 6.45, 85.00, 42.61)
sds   <- c(0.49, 0.62, 0.27, 0.47, 6.62, 6.98, 5.86)

# 共分散行列を作成
cov_matrix <- diag(sds) %*% cor_matrix %*% diag(sds)

# データ数（nは最大の3697に合わせる）
n <- 3697

# 多変量正規分布に基づくデータ生成
data <- as.data.frame(mvrnorm(n = n, mu = means, Sigma = cov_matrix))
colnames(data) <- var_names

# 欠損を追加
data$MP[2135:n] <- NA
data$UFP[695:n] <- NA

# 変数の調整
# IE: 0 or 1 の二値変数に変換
data$IE <- ifelse(data$IE > 0.5, 1, 0)

# Promotion: 0 or 1 に変換（昇進の有無）
data$Promotion <- ifelse(data$Promotion > 0.5, 1, 0)

# 値の切り捨てや補正
data$Salary <- round(data$Salary, 2)
# パフォーマンススコア 1〜5に制限
data$MP <- pmin(pmax(round(data$MP, 2), 1), 5)


# 分析はここから ====
# 各変数の有効サンプルサイズ（NA除外数）
valid_n <- colSums(!is.na(data))

# 各変数の平均値（NA除外）
var_mean <- apply(data, 2, mean, na.rm = TRUE)

# 各変数の標準偏差（NA除外）
var_sd <- apply(data, 2, sd, na.rm = TRUE)

# 相関分析
result <- corr.test(data)
cor_mat <- round(result$r, 3) # 小数第3位まで

# 表示用相関行列の作成

# 有意確率のマーカー
p_value <- as.numeric(result$r)
sigmark <- case_when(
  p_value < .001 ~ '#super("***")', # 表示を上付きにする
  p_value < .01  ~ '#super("**")',
  p_value < .05  ~ '#super("*")',
  TRUE     ~ " "
)
sig_table <- matrix(
  sigmark,
  nrow = nrow(result$r), ncol = ncol(result$r)
)

# 相関係数の整数位の0を削除
cor_str <- sprintf("%.3f", result$r) # 小数第3位まで
cor_str <- sub("0.", ".", cor_str)   # 0を削除
cor_str <- sub("-", "#sym.minus;", cor_str)   # マイナス記号の見た目を整える

cor_table <- matrix(
  cor_str,
  nrow = nrow(result$r), ncol = ncol(result$r),
  dimnames = dimnames(result$r)
)

cor_table[upper.tri(cor_table)] <- " " # 上三角を" "に
sig_table[upper.tri(sig_table)] <- " " # 上三角を" "に

diag(cor_table) <- rep("---", nrow(cor_table)) # 相関係数は対角成分を"---"に
diag(sig_table) <- rep(" ", nrow(sig_table)) # 有意確率は対角成分を" "に

# 変数名ラベル
var_lst <- list(
  "IE" = "社内外採用区分",
  "MP" = "管理職パフォーマンス", 
  "Salary" = "初任給", 
  "Promotion" ="昇進の有無", 
  "Tenure" = "在職年数", 
  "USP" = "部門サービス指標", 
  "UFP" = "部門財務指標"
)
# 相関行列の行名を変数名ラベルに変換
varnames <- unname(sapply(rownames(result$r), function(x) var_lst[[x]]))

# 必要なデータをJSONファイルに保存
out_json <- list(
  names = varnames,      # 変数名のリスト
  n = valid_n,           # 各変数のサンプルサイズ
  mean = var_mean,       # 各変数の平均値
  sd = var_sd,           # 各変数の標準偏差
  cor_table = cor_table, # 表示用相関係数行列
  sig_table = sig_table, # 有意確率マーカー行列
  cor_mat = cor_mat      # 数値が必要になる場合に備えて念の為
)

write_json(out_json, paste0("rdata.json"),
           pretty = TRUE, auto_unbox = TRUE)

Typst.typ

#import "@preview/zero:0.5.0": *
#import "@preview/booktabs:0.0.4": *
#show: booktabs-default-table-style

#let jfont = "Source Han Serif"
#set text(lang: "ja", font: jfont)


/////// データの処理
#let rdata = json("rdata.json")

// 相関係数行列を一度1次元配列にしてから、有意性マーカーとzipし、chunksで列を折り返す
#let soukan = rdata.cor_table.flatten().zip(rdata.sig_table.flatten()).chunks(rdata.n.len())

// 変数名、サンプルサイズ、平均値、標準偏差のデータと相関行列を連結
#let table_data = range(1, rdata.names.len()+1).zip(rdata.names, rdata.n, rdata.mean, rdata.sd, soukan)

/////// 表の準備
// 表の中の位置揃えの処理
#let table_align = (right, left, auto, auto, auto) // 変数番号、変数名、N、平均値、SD
#table_align.push( // 変数の数だけ（right, left）をつなげる
  rdata.names.map(x => (right, left))
)

// 列間のスペース調整
#let col_spaces = (auto, auto, auto, auto, auto) // 変数番号、変数名、N、平均値、SD
#col_spaces.push( // 変数の数だけ（-8pt, 0pt）をつなげる
  rdata.names.map(x => (-9pt, 0pt))
)

// 変数のヘッダ
#let header = range(0, rdata.names.len()).map(x => table.cell(colspan:2, align: center)[#(x+1)])

// 列数の計算
#let n_of_columns = table_data.at(0).flatten().len()



/////// 相関係数の表

#set-group(separator: ",", threshold: 4) // 位取りコンマの設定（zeroパッケージ）
#block[
  #set text(7pt)
  #ztable(
    columns: n_of_columns,
    align: table_align.flatten(),
    column-gutter: col_spaces.flatten(),
    format:(none, none,
      (digits:0, math: false), 
      (digits:2, math: false), 
      (digits:2, math: false)),
      toprule(),
      table.header(
        table.cell(colspan:2, align: center)[変数],
        table.cell(colspan:1, align: center)[n],
        [平均値],[SD],
        ..header,
      ),
      midrule(),
      ..table_data.flatten().map(x => eval(str(x), mode: "markup")
    ),
    bottomrule()
  )
]

作成された表