4  予備的な分析と信頼性・妥当性

著者
所属

神戸大学大学院経営学研究科

初版公開日

June 7, 2022

最終更新日

December 12, 2024

abstract
収集したデータを問題なく分析にかけられるかを評価するための方法として,項目分析(I-T相関・トレースライン)を紹介しています。また,使用した尺度(項目)に問題がなかったかを判断するために,信頼性や妥当性の考え方について説明しています。
Keywords

R, 多変量解析, 項目分析, 信頼性, 妥当性

本資料は,神戸大学経営学研究科で2022年度より担当している「統計的方法論特殊研究(多変量解析)」の講義資料です。CC BY-NC 4.0ライセンスの下に提供されています。


前回までで,データを読み込み,転記ミスなどの異常が無いかを確認しました。 また,適当回答や欠測値の処理を行い,ようやく分析ができる「きれいな」データを用意することが出来ました。 今回は本格的な分析に入る前の前段階として,項目や尺度の性質・性能を簡単に確認していきます。 この段階で,明らかに良くない項目は事前に削除しておくことによって,その先の分析でより安定した結果が得られるようになるはずです。

コード 4.1: .rdsファイルを読み込む
dat <- readRDS("chapter03.rds")

4.1 分布の可視化

データ内のゴミをある程度取り除けたところで,各変数の分布を見てみましょう。 量的変数の場合,要約統計量だけではなく分布を可視化するのも大事です。 心理尺度の場合,合計点だけでなく各項目のレベルで,できるだけ回答が正規分布に近い形でばらついている方が都合が良い事が多いです。 例えば「毎日睡眠をとっている」や「ここ1ヶ月人と会話していない」のように,ほぼすべての人に当てはまりすぎる・全く当てはまらないような項目は,ほぼすべての人が同じ回答になるため,その項目に対する回答が個人の特性値に対して何の情報も与えない可能性があります。 分析の前に,このような天井効果床効果と呼ばれる状態が発生していないかを確認したり,分布のピークが一つになっているかを確認しておきましょう。 また余裕があれば,2変数の散布図を書いたりして,項目間の相関関係について想像したりするのも良いでしょう。

なお,ここから先では「とりあえず図を出す」レベルの話だけします。論文に載せるようなきれいなグラフを描きたい場合には色々引数をいじりまくる必要がありますが,それだけで1コマくらいは潰れると思うので…1

4.1.1 1変数の棒グラフ

リッカート尺度の項目反応など,取る値のパターンがせいぜい数個(離散変数や名義変数)の場合は,棒グラフが良いでしょう。棒グラフはbarplot()table()を組み合わせることで描けます。

コード 4.2: 離散変数の棒グラフ
barplot(table(dat$Q1_1))
図 4.1: 1変数の棒グラフ

barplot()は,与えられたベクトルの一つ一つの値を棒の長さで表現する関数です。そのため,barplot(dat$Q1_1)というように直接データを与えてしまうと, 図 4.2 のように左から長さが2, 2, 5, 4, 2, …という感じで2800本の棒グラフが誕生してしまいます。

コード 4.3: table()を噛ませないとダメ
barplot(dat$Q1_1)
図 4.2: table()を噛ませないと大変なことに

4.1.2 天井効果と床効果

天井効果(床効果)は,回答が一番上(下)の選択肢に偏っている傾向のことです。図 4.1 くらいの多少の効果は気にする必要がないのですが,極端な天井・床効果が生じている項目は除外したほうが良いとされます。 極端な天井・床効果が生じているということは,その項目ではほとんどの人が同じ選択肢を選んでいる,ということです。 通常,心理尺度の項目に対しては 図 4.3 の横軸のように,背後に連続量としての潜在特性を仮定し,各人の回答からその潜在特性の強さを推定します。 つまり各回答者(A,B,Cさん)は自分の潜在特性の強さがどの程度かに応じて,適切な選択肢(1から4)を選択しているはず,と考えます。 これに対して,例えば床効果が生じている状態は 図 4.4 のような状態です。 図 4.3 と比較して,各回答者の位置は変わっていないですが全員が1と回答しています。 実際に私達が観測できるのは回答結果としての「1」のみなので,この項目に「1」と回答したことはその回答者の潜在特性を推定する上でなんの情報も持っていないということになってしまいます。 もちろんごく少数の「極端な人」を識別し「どの程度極端か」を評価する上では有用な項目なので,項目数の制約がなければ入れておいても良いのですが,項目数をなるべく少なくしたいなら,こういった項目は削除したほうが良さそう,と判断できます。 Section 4.3 では,天井効果・床効果が生じている可能性のある項目に「あたり」をつける簡便な方法を紹介します。

図 4.3: ふつうの項目の背後に
図 4.4: 床効果が生じている項目の背後に

また,項目反応理論などの方法では,こういった極端な項目が含まれていると他の項目の特性の推定にもあまり良くない影響を及ぼしたり,推定自体が不安定になってしまうケースがあります。そういった観点からも,極端な天井効果や床効果が見られる項目は除外したほうが良いのです。

4.1.3 1変数のヒストグラム

年齢など,連続変数の場合は棒グラフよりもヒストグラムのほうが良いでしょう。ヒストグラムはhist()関数で描くことが出来ます。

コード 4.4: ヒストグラムを描く
hist(dat$age)
図 4.5: 1変数のヒストグラム

ヒストグラムでは,区切り幅が重要な役割を果たします。 デフォルトでは,スタージェスの公式 Sturges, 1926 と呼ばれるものに基づいて区切り幅が決定されています。 これは,データの総数をnnと置いた場合に,データの値が取る全範囲を1+log2n1+\log_2 n個に等分割する方法です。 したがって,2432人分のデータがあるdatでは1+log224321+\log_2 2432からおよそ12個に等分割された結果が 図 4.5 に表れています。 基本的にはデータの数にあわせていい感じに調整してくれるので良いのですが,場合によっては「年齢を10刻みで区切る」など,区切り幅をきっちり指定したいこともあるでしょう。 これを実現するためには,hist()の引数breaksを指定してあげます。

コード 4.5: 区切り幅を変える
hist(dat$age, breaks = 20) # 左上
hist(dat$age, breaks = c(0, 10, 20, 40, 50, 100)) # 左下
hist(dat$age, breaks = seq(0, 100, by = 10)) # 右上
hist(dat$age, breaks = seq(0, 100, length = 10)) # 右下
図 4.6: 区切り幅を変えたヒストグラム

hist()の引数breaksの指定方法は,大きく分けると単一の整数を与える方法と,区切り位置のベクトルを与える方法の2種類があります。 前者の場合は,データの値が取る全範囲をその個数で等分割してくれます。 後者の場合は,等分割以外の区切り方も指定可能( 図 4.6 の左下)ですが,上限と下限がデータの値が取る全範囲を含んでいないとエラーが返ってくるのでご注意ください。

ちなみにbreaksの指定で使用したseq()関数は,等差数列(ベクトル)を作成してくれる関数です。

コード 4.6: seq(): 等差数列を作る
# 最初の2つの引数が「どこからどこまで」
# 次に引数byを与えた場合,「いくつ刻み」
seq(0, 100, by = 10)
# 次に引数lengthを与えた場合,「長さいくつ」
seq(0, 100, length = 10)
 [1]   0  10  20  30  40  50  60  70  80  90 100
 [1]   0.00000  11.11111  22.22222  33.33333  44.44444  55.55556  66.66667
 [8]  77.77778  88.88889 100.00000

Sturges (1926 の説明によれば,スタージェスの公式は「二項定理の係数の和」をもとに考えられた値のようです。ここでは,なぜ二項定理が出てくるのかを簡単に説明してみましょう。

まず,ヒストグラムを作りたい変数が正規分布に従っているとします。このとき,正規分布は二項分布で近似することが可能です(二項分布の正規近似の逆)。 とりあえず左右対称な二項分布を作るために,成功確率50%の二項分布Binom(ñ,12)Binom(\tilde{n}, \frac{1}{2})で近似してみましょう。 例えばñ=4\tilde{n}=4とした場合には,P(x=k|k=0,1,2,3,4)P(x=k|k=0,1,2,3,4)はそれぞれ,ñCk×(12)4=c(1,4,6,4,1)×(12)4{}_{\tilde{n}}\mathrm{C}_k\times(\frac{1}{2})^4=c(1,4,6,4,1)\times(\frac{1}{2})^4となり,左右対称な確率分布になります。 このñCk{}_{\tilde{n}}\mathrm{C}_kの総和が,「二項定理の係数の和」です。

上の(1,4,6,4,1)(1,4,6,4,1)は,もとのデータがきれいな正規分布に従うとき,それを5等分した各区切りの中に含まれるデータの割合が(1,4,6,4,1)(1,4,6,4,1)になることを意味しています。 つまり,もしデータが16個だった場合には,ヒストグラムは左から順に(1,4,6,4,1)(1,4,6,4,1)と並ぶだろう,ということです。 言い換えると,データが16個だった場合には,データを5等分してあげたらいい感じのヒストグラムになりそうだ,と言えるわけです。

ここで一つ重要なポイントがあります。それは,「二項定理の係数」は必ずñ+1\tilde{n}+1個になり,その総和が2ñ2^{\tilde{n}}になる,ということです。 したがって,データの数が2n2^n個のときには,ヒストグラムをn+1n+1個に区切ってあげたらよいだろう,という感じで提案されたのがスタージェスの公式です(たぶん)。

ちなみに,スタージェスの公式はヒストグラムの最適な区切り幅を計算するための公式です。実際のヒストグラムでは,区切り幅は意味のある(または解釈上都合の良い)値として,5や10などキリの良い値になることも多いです。そのため,スタージェスは「意味のある区切り幅のうち,この公式で算出された値を超えない最大の値」にすることを提案しています。

ただし心理尺度の得点などはそもそも値の単位に絶対的な意味がないため,こうした変数に対して用いるような場合にはスタージェスの公式にとりあえずそのまま乗っかってしまっても良い気がします。

4.1.4 2変数の散布図

2変数の散布図を描く場合には,plot()関数が使えます。なお,plot()関数は他にもデータの形式などに合わせて何かしらの図を出してくれたりするので,色んなところでお世話になると思います。

コード 4.7: 散布図を描く
# 1問目(Q1_1)と年齢(age)の散布図
# 2変数なので2列指定してあげる
plot(dat[, c("Q1_1", "age")])
図 4.7: 2変数の散布図

図 4.7 ではQ1_1の取りうる値が6通りしかないため散布図感が薄いですが,一つ一つのマルが各レコードを表しています。

また,2変数のtable()を用いるとクロス表が得られるのですが,これをplot()に与えると違った形の図が得られます。

コード 4.8: tableをplot
table(dat[, c("Q1_1", "Q1_2")])
    Q1_2
Q1_1   1   2   3   4   5   6
   1  10  13   5   7  15  22
   2   6  26  15  46  61  38
   3   6  21  33  97  86  49
   4   3  19  39 104 128  56
   5   7  25  39 152 333 163
   6   7   9   7  73 273 439
# そのままplot()に入れると
plot(table(dat[, c("Q1_1", "Q1_2")]))
図 4.8: いわゆるモザイク図

図 4.8 では,それぞれの四角形および行・列の幅が,各選択肢を選んだ人の割合をそれぞれ表しています。

変数の分布をプロットする方法は,他にも色々ありますが(例えば箱ひげ図やバイオリンプロットなど),必要に応じて適切な方法を選べるようになりましょう。 本講義では,あくまでも基本的なプロットの方法のみの紹介にとどめて次に進んでしまいます。

4.2 相関関係のチェック

Section 2.2 では1変数レベルでの記述統計量をチェックする関数としてsummary()などを紹介しました。 ここでは,同じようにして2変数の関係の記述統計量である相関係数を確認する方法を紹介しておきます。 相関係数はcor()に2つのベクトルを与えることによって計算可能です。これも欠測値がある場合NAを返してくるので,きちんと相関係数を出してほしい場合には適切な引数を与えます。ただし,mean()などの時に使用するna.rmではなくuseという引数を使用します。(説明は後ほど)

コード 4.9: 2変数の相関係数
# "pairwise"は"pairwise.complete.obs"の略
cor(dat[, "Q1_1"], dat[, "Q1_2"], use = "pairwise")
[1] 0.3527771

cor()には上のように2つのベクトルを与える他に,データフレームをまるごと与えることも可能です。データフレームを一つだけ与えた場合にはそのデータフレーム内の全ての列の相関行列が,データフレームを2つ与えた場合にはデータフレーム同士の全ての列の組み合わせの相関行列が得られます。

コード 4.10: データフレーム内のすべての相関係数
# 量が多いので省略しています
cor(dat, use = "pairwise")
              ID       Q1_1        Q1_2        Q1_3       Q1_4
ID    1.00000000 -0.0442558 -0.04853408 -0.01340771 0.01208236
Q1_1 -0.04425580  1.0000000  0.35277706  0.27444894 0.16169504
Q1_2 -0.04853408  0.3527771  1.00000000  0.49741109 0.35019071
Q1_3 -0.01340771  0.2744489  0.49741109  1.00000000 0.38480054
Q1_4  0.01208236  0.1616950  0.35019071  0.38480054 1.00000000
コード 4.11: 第1因子と第2因子の項目間相関行列
cor(dat[, 2:6], dat[, 7:11], use = "pairwise")
            Q1_6        Q1_7       Q1_8      Q1_9      Q1_10
Q1_1 -0.01132405 -0.01088135 0.02365191 0.1161066 0.03943949
Q1_2  0.09785479  0.12411971 0.18509733 0.1523665 0.12848202
Q1_3  0.10989033  0.14186345 0.12557243 0.1261299 0.16038337
Q1_4  0.08963104  0.22389063 0.12966842 0.1739516 0.25181181
Q1_5  0.13115144  0.11412746 0.12906468 0.1290221 0.17094545

引数useは,計算に使用するデータを決定します。デフォルトでは"everything",つまり全部使用します。"complete.obs"は,与えたデータフレーム内で欠測値が一つも無いレコードだけを使用して相関を計算します。"pairwise.complete.obs"では各相関を計算する際に,その2変数において欠測値が無いレコードを使用して相関を計算します。つまり相関行列の要素ごとに計算に使用するレコードの数・内容が異なる可能性がある,ということです。

4.2.1 (おまけ)相関関係も可視化してみる

心理尺度に対する多変量解析では,複数の質問項目が同じ構成概念を反映しているかが分析手法のキーとなります。 これを確認する方法の一つとして,後半では信頼性のお話が出てくるわけですが,ここでも少しだけ,相関関係を可視化することでその感覚を掴んでみましょう。

相関行列を可視化するためのパッケージとして,corrplotというパッケージが用意されています。

# install.packages("corrplot")
library(corrplot)

インストールできたら,早速プロットしてみましょう。

コード 4.12: 相関行列を可視化する
cols <- paste0("Q1_", 1:25) # 25項目だけ使う
corrplot(cor(dat[, cols]))
図 4.9: 相関行列の可視化

この図では,相関係数が高い変数の組のところほど濃い色(プラスなら青,マイナスなら赤)で表されています。 ざっと眺めるだけでも,いくつかのまとまりがあるのが見て取れる気がします(Q1_16からQ1_20など)。 質問項目間に正の相関が見られるということは,これらの項目の背後には,何か共通の要因があるのかもしれません。 後半で紹介する因子分析 Chapter 6 では,このような変数間の相関関係をもとに,その背後にあると考えられる要因に迫っていきます。

ちなみに,引数methodを指定することで,プロットの形を変えることも可能です。

コード 4.13: 相関行列を可視化する(2)
# 個人的にはこちらのほうがまとまりが見やすい気がします
corrplot(cor(dat[, cols]), method = "shade")
図 4.10: 相関行列の可視化(2)

4.3 項目分析

まずはじめに,項目単位での「変なやつ」がいないかをチェックしていきます。

4.3.1 天井効果・床効果について検討する

天井効果または床効果が発生しているかの判断について,一発で見抜ける方法というものはありません。 あくまでも一つの基準としてですが,よく言われているのは「平均値 ± 標準偏差が,取りうる値の範囲(今回の例では1から6)を超えている場合」というものらしいです(要出典)。特に変数の数が多くて全ての棒グラフを作って見るのも骨が折れる,という場合には,まずこの判別方法であたりを付けてみるのも良いかもしれません。

とりあえず上述の基準にのっとって天井効果または床効果を確認するためには,各変数(列)ごとに平均値や標準偏差を計算する必要があります。普通に考えると次のようなコードによって実現可能です。

コード 4.14: 1問目の平均値と標準偏差を足すと
# 欠測値があるため,na.rm=TRUEを引数に追加している
mean(dat[, "Q1_1"], na.rm = TRUE) + sd(dat[, "Q1_1"], na.rm = TRUE)
[1] 5.998918
コード 4.15: 1問目の平均値から標準偏差を引くと
# 欠測値があるため,na.rm=TRUEを引数に追加している
mean(dat[, "Q1_1"], na.rm = TRUE) - sd(dat[, "Q1_1"], na.rm = TRUE)
[1] 3.18776

ただ,今回は25項目もあります。上の計算をそれぞれの列にやっても良いのですが効率が悪いので,apply()を使って

コード 4.16: 列ごとの平均値
# 本当はQ1_1からQ1_25のところだけに適用したら良いが,
# 今回は全ての列がnumericかinteger型なので気にせず進めます。
apply(dat, 2, mean, na.rm = TRUE)
         ID        Q1_1        Q1_2        Q1_3        Q1_4 
1396.945312    4.593339    4.801398    4.602385    4.691612 
       Q1_5        Q1_6        Q1_7        Q1_8        Q1_9 
   4.548931    4.527138    4.374178    4.301398    4.449013 
      Q1_10       Q1_11       Q1_12       Q1_13       Q1_14 
   3.692434    4.020148    3.846217    3.987253    4.414474 
      Q1_15       Q1_16       Q1_17       Q1_18       Q1_19 
   4.394326    2.944901    3.519737    3.226151    3.203947 
      Q1_20       Q1_21       Q1_22       Q1_23       Q1_24 
   2.972451    4.814967    4.314556    4.454359    4.927632 
      Q1_25      gender   education         age 
   4.528783    1.670230    3.190860   28.582237 

とすればよいわけです2。したがって,各列の平均値±標準偏差をまとめて計算したい場合はそれぞれ以下のようになります。

コード 4.17: 列ごとの平均値+SD
apply(dat, 2, mean, na.rm = TRUE) + apply(dat, 2, sd, na.rm = TRUE)
         ID        Q1_1        Q1_2        Q1_3        Q1_4 
2207.424918    5.998918    5.974081    5.908401    6.172146 
       Q1_5        Q1_6        Q1_7        Q1_8        Q1_9 
   5.812634    5.758867    5.690390    5.589267    5.825444 
      Q1_10       Q1_11       Q1_12       Q1_13       Q1_14 
   5.323572    5.650285    5.458144    5.335440    5.876204 
      Q1_15       Q1_16       Q1_17       Q1_18       Q1_19 
   5.733057    4.519407    5.050856    4.819228    4.771938 
      Q1_20       Q1_21       Q1_22       Q1_23       Q1_24 
   4.594638    5.936672    5.866132    5.654409    6.116027 
      Q1_25      gender   education         age 
   5.852550    2.140456    4.302853   39.550527 
コード 4.18: 列ごとの平均値-SD
apply(dat, 2, mean, na.rm = TRUE) - apply(dat, 2, sd, na.rm = TRUE)
        ID       Q1_1       Q1_2       Q1_3       Q1_4 
586.465707   3.187760   3.628715   3.296368   3.211077 
      Q1_5       Q1_6       Q1_7       Q1_8       Q1_9 
  3.285227   3.295409   3.057966   3.013529   3.072583 
     Q1_10      Q1_11      Q1_12      Q1_13      Q1_14 
  2.061297   2.390011   2.234291   2.639066   2.952743 
     Q1_15      Q1_16      Q1_17      Q1_18      Q1_19 
  3.055594   1.370395   1.988618   1.633075   1.635957 
     Q1_20      Q1_21      Q1_22      Q1_23      Q1_24 
  1.350264   3.693263   2.762979   3.254308   3.739236 
     Q1_25     gender  education        age 
  3.205016   1.200004   2.078867  17.613946 

天井効果が出ている可能性がある(平均値+SDが6を超えている)変数はQ1_4Q1_24の2つです。また,床効果が出ている可能性のある項目はなさそうでした。

ということで,天井効果が出ている可能性のある変数についてプロットしたものが 図 4.11 です。

コード 4.19: 天井効果はあるのか?
barplot(table(dat[, "Q1_4"]))
barplot(table(dat[, "Q1_24"]))
(a) Q1_4
(b) Q1_24
図 4.11: 天井効果が疑われる項目の棒グラフ

実際に天井効果および床効果が生じていると思われる項目を分析から除外するかについては,定量的な基準(平均±標準偏差)だけで決めるべきではありません。実際には,この後紹介する妥当性などの観点から項目の内容を吟味したりして,構成概念の定義的に削除しても問題ないかを判断する必要があります。ちなみに今回の程度であれば,分布が多少偏っているとはいえ,除外するほどの大きな問題とはいえないでしょう。

私が 図 4.11 のようなレベルの天井効果について除外しなくても良いと考えている理由は,こうした項目は少なくとも特性値の低い人を区別する力は持っているためです。

心理尺度への回答は,その背後にある心理特性(社交性や誠実さなど)の強さを反映して決まっていると考えられます。 つまり天井効果(床効果)が起きているような項目では,「その心理特性が一定レベル以上(以下)の強さの人はみんな端の選択肢を選んでいる」状態なわけです。 その「一定レベル」の閾値を仮に30とすると,天井効果が起きている項目は,レベル80の人とレベル90の人はどちらも最大値の選択肢を選ぶために区別ができないわけですが,少なくともレベル10の人とレベル20の人は区別できるだろう,というわけです。 したがって,多少の天井/床効果が起きている程度で,その項目には全く意味がない,ということにはならないだろう,というのが私の考えです。

ただし,天井/床効果が起きていない項目と比べると,1項目が持つ「高低を区別する総合力」は低い可能性が高いので,例えば尺度を新規作成する際や,質問票にいれることのできる項目数に厳しい制限があるならば,どちらかといえば天井/床効果が起きていない項目を選んだほうが良いと思います。 一方で,上述の理由から,既存の尺度を使って,しかも既に取った後のデータの分析を行う段階ならば,そこでわざわざ多少の天井/床効果での項目削除はしなくても良いのではないか,と思っています。

また,除外するかを判断する場合に重要となるのは,天井/床効果が生じているかよりも,回答が一箇所に集中しすぎていないかだと思っています。もしも全ての人がちょうど中間の選択肢を選んでいた場合,天井/床効果はいずれも生じていなくても項目としては明らかにおかしなものになっているはずです。 そう考えると,単に標準偏差が異様に小さい項目が無いかを見ておけば良いような気がしています。

4.3.2 I-T相関

ここでは大前提として,心理尺度では複数の項目が同じ構成概念を測定しようとしているはずだというところから出発して考えます。 このとき,ある項目が,同じ尺度内の他の項目と確かに同じ構成概念を測定しているならば,他の項目の得点との相関は高いはずです。ということは当然,合計点との相関も高いだろうという予測が立ちます。 この考えに基づいて算出される指標が,各項目の値と,尺度の合計点の相関であるI-T相関 (Item-TestまたはItem-Total)です。 他の項目との一貫性という点では,後述する信頼性の分析の一部という見方もできます3

I-T相関を出すためには,まず合計点を出す必要があります。これはapply()を使えば簡単に出来ます。最も簡便な尺度得点としてそのまま他の変数との関係を見る分析にも使えたりするので,ここで各因子ごとに合計点を作成してdatに追加しておきましょう。 なおRではほとんどの場合,データフレームの存在しない列名を指定して代入すると,新しい列を追加してくれます4

コード 4.20: 列ごとの平均値+SD
# 列の指定の仕方はいろいろありました
# ↑復習を兼ねて,わざと異なる方法で列を指定しています
# Agreeableness
dat[, "Q1_A"] <- apply(dat[, paste0("Q1_", 1:5)], 1, sum)
# Conscientious
dat$Q1_C <- apply(dat[, paste0("Q1_", 6:10)], 1, sum)
# Extraversion
dat[, "Q1_E"] <- apply(dat[, paste0("Q1_", 11:15)], 1, sum)
# Neuroticism
dat["Q1_N"] <- apply(dat[, paste0("Q1_", 16:20)], 1, sum)
# Openness
dat$Q1_O <- apply(dat[, paste0("Q1_", 21:25)], 1, sum)

# 一旦確認
head(dat)
  ID Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 Q1_12
1  1    5    4    3    4    4    2    3    3    3     3     4     4
2  2    5    4    5    2    5    5    4    4    4     3     6     6
3  3    2    4    5    4    4    4    5    4    5     2     5     3
4  4    3    4    6    5    5    4    4    3    2     2     2     4
5  5    5    3    3    4    5    4    4    5    4     5     5     5
6  6    1    6    5    6    5    6    6    6    6     4     5     6
  Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 Q1_23
1     3     4     4     3     4     2     2     3     3     1     3
2     6     4     3     3     3     3     5     5     4     5     4
3     4     4     5     4     5     4     2     3     4     5     5
4     4     4     4     2     5     2     4     1     3     4     4
5     5     4     5     2     3     4     4     3     3     4     4
6     6     5     6     3     5     2     2     3     4     4     5
  Q1_24 Q1_25 gender education age Q1_A Q1_C Q1_E Q1_N Q1_O
1     4     4      1        NA  16   20   14   19   14   15
2     3     4      2        NA  18   21   20   25   19   20
3     5     5      2        NA  17   19   20   21   18   24
4     3     2      2        NA  17   23   15   18   14   16
5     3     4      1        NA  17   20   22   24   16   18
6     6     6      2         3  21   23   28   28   15   25

合計点が出来たら,後は相関係数を計算するだけです。ですが,ここで計算するべき相関係数は,普通の(ピアソンの積率)相関係数ではありません。

カテゴリカル変数の相関

各項目への回答は,本来連続量である潜在変数(因子)の強さによって規定されます。実際に観測できる項目反応(5件法なら1から5の数字)は,この潜在変数を適当なところで打ち切ってカテゴリ化された値とみなせるわけです。 図 4.12 はその関係を表した図です。潜在変数は通常,正規分布に従うと仮定されます。心理尺度への回答データを考える際には,正規分布の仮定のもとで,潜在変数の強さの程度に応じて,回答者は自分にあった強さの選択肢を選んでいる,というように考えます5

図 4.12: 5件法に対する潜在変数と回答の関係

ということは,合計点6との相関係数を計算する場合,本来は項目への回答そのものではなく,背後にある潜在変数の値で計算したほうが良いのです。むしろそうしないと,連続量をカテゴリに落とし込むことで失われる情報のせいで,本来の相関係数よりも低い値が算出されてしまいます。

今回はポリシリアル相関(polyserial correlation)というものを使います。これは,(多値)カテゴリカル変数の背後にあると想定される連続量と,連続変数との相関係数です。 他にも,同じような考え方の特殊な相関係数がいくつかあり,それらは 表 4.1 のようにまとめられます。

表 4.1: カテゴリカル変数との相関係数
連続 多値 二値
連続 ピアソンの積率相関 ポリシリアル相関 点双列相関
(バイシリアル相関)
多値 ポリシリアル相関 ポリコリック相関 ポリコリック相関
二値 点双列相関
(バイシリアル相関)
ポリコリック相関 テトラコリック相関

これらの特殊な相関係数を計算するためのパッケージとしては,polycorパッケージというものがあります。

# まだパッケージをインストールしていない場合は
# install.packages("polycor")
library(polycor)

今回はポリシリアル相関を計算するため,polyserial()を使用しましょう7。 少し注意が必要な点として,polyserial()関数では

  • 第1引数(x)に入れた変数は自動的に連続変数として扱われる
  • 第2引数(y)に入れた変数は自動的にカテゴリカル変数として扱われる

つまり引数の順番が変わると計算結果が変わってしまいます。 そして,このように引数の順番が重要な関数は,apply()関数と組み合わせて使用する場合にやや面倒なことが起こります。 図 4.13 に,簡単な概要をまとめたものを用意しました。 apply(X, MARGIN, FUN)では,MARGINの方向に分割されたXを,それぞれFUNで指定した関数の第一引数に与えるという挙動をします。 つまりFUNpolyserialを指定した場合,Xで与えたものは自動的に連続変数として扱われざるを得ないのです。

図 4.13: polyserial()関数をapply()関数の中で呼んだとき

こればっかりはapply()関数の仕様のため,たぶんどうしようもありません。 そこでやや面倒ですが,回避策としてpolyserial()関数の引数の順番を入れ替えたものを別の関数として先に定義しておくことにします。

コード 4.21: polyserial()関数の引数を入れ替えただけの関数を定義
# x: I-T相関を出したい項目
# total: 合計点
itcor <- function(x, total) {
  polycor::polyserial(total, x)
}
コード 4.22: I-T相関の計算(1)
# apply()の第一引数xは,FUNに与えた関数(itcor)の第一引数になる
# 第二引数以降は引数FUNより後ろに明示的に与える

# 因子ごとに,計算に使用する合計点(total)が変わるので,5項目ごとにそれぞれ書く
# Agreeableness
apply(dat[, paste0("Q1_", 1:5)], 2, itcor, dat[, "Q1_A"])
     Q1_1      Q1_2      Q1_3      Q1_4      Q1_5 
0.6343893 0.8006811 0.8320099 0.7379408 0.7392092 
# Conscientious
apply(dat[, paste0("Q1_", 6:10)], 2, itcor, dat[, "Q1_C"])
     Q1_6      Q1_7      Q1_8      Q1_9     Q1_10 
0.6980845 0.7432043 0.7076029 0.7993942 0.7615535 
# Extraversion
apply(dat[, paste0("Q1_", 11:15)], 2, itcor, dat[, "Q1_E"])
    Q1_11     Q1_12     Q1_13     Q1_14     Q1_15 
0.7650335 0.8265158 0.7100692 0.8071135 0.6899571 
# Neuroticism
apply(dat[, paste0("Q1_", 16:20)], 2, itcor, dat[, "Q1_N"])
    Q1_16     Q1_17     Q1_18     Q1_19     Q1_20 
0.8495472 0.8190185 0.8462378 0.7489533 0.7173007 
# Openness
apply(dat[, paste0("Q1_", 21:25)], 2, itcor, dat[, "Q1_O"])
    Q1_21     Q1_22     Q1_23     Q1_24     Q1_25 
0.6689789 0.7160443 0.7168503 0.5418518 0.7209368 

I-T相関という場合,実は2種類の相関係数が考えられます。一つはいま計算した「合計点との相関」です。しかし今回のように項目数が少ない場合,「合計点」はその項目自身の影響を結構強く受けていると考えられます。そこで,「合計点から自身の得点を引いたものとの相関」I-T相関として考えることもあります。例えばQ1_1のI-T相関であれば,Q1_2からQ1_5までの4項目の合計点との相関を見る,ということです。

「合計点から自身の得点を引いたものとの相関」を見る場合,先程よりはやや難しくなります。これもやり方はいくつかあると思いますが,ここでは関数を定義してやってみましょう。

コード 4.23: 関数の定義を少し変える
# 先程定義したitcor()関数をベースに
# 引数totalを「合計点から自身の得点を引いたもの」にしたら良いので
itcor2 <- function(x, sum) {
  polyserial(sum - x, x)
}
コード 4.24: I-T相関の計算(2)
# あとは先程と同じ

# 因子ごとに,計算に使用する合計点(total)が変わるので,5項目ごとにそれぞれ書く
# Agreeableness
apply(dat[, paste0("Q1_", 1:5)], 2, itcor2, dat[, "Q1_A"])
     Q1_1      Q1_2      Q1_3      Q1_4      Q1_5 
0.3497715 0.6254912 0.6506640 0.4597522 0.5309969 
# Conscientious
apply(dat[, paste0("Q1_", 6:10)], 2, itcor2, dat[, "Q1_C"])
     Q1_6      Q1_7      Q1_8      Q1_9     Q1_10 
0.4966869 0.5445442 0.5033232 0.6130345 0.5143994 
# Extraversion
apply(dat[, paste0("Q1_", 11:15)], 2, itcor2, dat[, "Q1_E"])
    Q1_11     Q1_12     Q1_13     Q1_14     Q1_15 
0.5473676 0.6487321 0.5260702 0.6280941 0.4922296 
# Neuroticism
apply(dat[, paste0("Q1_", 16:20)], 2, itcor2, dat[, "Q1_N"])
    Q1_16     Q1_17     Q1_18     Q1_19     Q1_20 
0.7127767 0.6793635 0.7088856 0.5709521 0.5129809 
# Openness
apply(dat[, paste0("Q1_", 21:25)], 2, itcor2, dat[, "Q1_O"])
    Q1_21     Q1_22     Q1_23     Q1_24     Q1_25 
0.4291648 0.3812691 0.4836023 0.2386248 0.4531874 

当然ですが,自身の得点を引いたものとのI-T相関のほうが低くなりました。ただ今回はいずれの項目も問題ないレベルのI-T相関を示しています。(当然と言えば当然ですが)サンプルデータなのでかなりきれいです8

ちなみにI-T相関が低いということは,それすなわち項目が悪い,というわけではありません。あくまでも他の項目との相関が低いだけです。ただ,I-T相関が低い項目が入っている場合,因子分析をしても思った通りに因子がまとまらなかったり,項目反応理論でもうまくパラメータ推定ができない可能性が高くなります。 そのため,内容次第ではありますがI-T相関が低い項目があればこの段階で分析から除外する,ということもあります。実際のところどの程度低かったら問題視されるかの明確な基準は多分無いのですが,経験的には自身の得点を引いたものとのI-T相関が0.2から0.1を下回ると,相当厳しいケースが多いです。

カテゴリカルな相関係数の計算方法の考え方

上述したようなカテゴリ変数に対する相関係数では, 図 4.12 に示したように,各項目への回答の背後に正規分布に従う連続量があると仮定します。 ここで,2つの二値項目(当てはまる/当てはまらない)への回答xi,xjx_i,x_jについて,(点双列)相関係数の出し方を説明します。 もちろん多値項目に対する相関係数も原理は同じです。

二値項目その背後にある連続量をそれぞれui,uju_i,u_jとしておきます。 カテゴリカルな相関係数では,連続量が閾値を超えていたら1, 超えなければ0と回答している,と考えるので,各項目についての閾値をそれぞれτi,τj\tau_i,\tau_jと置くと,項目への回答と連続量の関係性を以下のように表すことができます。 xi={0ui<τi1uiτixj={0uj<τj1ujτj \begin{aligned} x_i = \begin{cases} 0 & u_i < \tau_i \\ 1 & u_i \geq \tau_i \end{cases} \;\;\;\;\;\; x_j = \begin{cases} 0 & u_j < \tau_j \\ 1 & u_j \geq \tau_j \end{cases} \end{aligned}

uiu_iに関して,ここまでの関係を表したものが 図 4.14 です。 同様の関係がuju_jについても存在しているとします。

図 4.14: 2件法に対する潜在変数と回答と閾値の関係

いま知りたいのは,2つの回答そのものの相関係数rxi,xjr_{x_i,x_j}よりも,その背後にある連続量同士の相関係数rui,ujr_{u_i,u_j}です。 これを計算するためには,2つの項目についての 図 4.14 に2変量正規分布の図を組み合わせて,実際に得られる回答の割合を考えます。 図 4.15 の2つの図は,いずれも(τi,τj)=(0.2,0.8)(\tau_i,\tau_j)=(0.2,-0.8)として作図したものです。

(a) 相関係数0.3
(b) 相関係数0.8
図 4.15: 2つの潜在変数の相関関係と回答

周辺尤度(各図の上と右にある正規分布)にはいずれも標準正規分布を仮定して,それぞれτi=0.2\tau_i=0.2以上およびτj=0.8\tau_j=0.8のところが塗りつぶされています。 いま,各図の中の楕円は2変量正規分布の確率密度(左:相関弱め,右:相関強め)を表しています。 したがって楕円のうち塗りつぶされている箇所(右上)は,2変量正規分布においてτi=0.2\tau_i=0.2以上かつτj=0.8\tau_j=0.8となる確率=2項目の両方に「当てはまる」と回答する確率を表しています。 図からもわかるように,2変量正規分布の形は相関係数によって変化するため,2つの閾値(τi=0.2\tau_i=0.2τj=0.8\tau_j=0.8)で区切られた4つの区間の面積は相関係数によって変わります。

そして観測されたデータにも,2つの二値項目に対する回答の組み合わせが4パターン存在しうるので,それぞれの出現割合(0π10\leq \pi_{\cdot\cdot}\leq 1)を 表 4.2 のようにまとめておきます。

表 4.2: 2項目への回答のパターンとその割合
xj=0x_j=0 xj=1x_j=1
xi=0x_i=0 π00\pi_{00} π01\pi_{01}
xi=1x_i=1 π10\pi_{10} π11\pi_{11}

これで材料が出揃いました。 カテゴリ変数に対する相関係数では,図 4.15 の楕円の4つの領域の確率と 表 4.2 の4つのπ\piの割合が一致するように各項目の閾値τi,τj\tau_i,\tau_jおよび相関係数rui,ujr_{u_i,u_j}を求めます。 例えば「100人中25人が(1,1)(1,1)と回答した」という場合にはπ11=0.25\pi_{11}=0.25となるので, 図 4.15 の楕円のうち塗りつぶされている部分(右上)の確率も0.25になるように,(τi,τj,rui,uj)(\tau_i,\tau_j,r_{u_i,u_j})の3つを計算し,得られたrui,ujr_{u_i,u_j}の値が(点双列)相関係数となるわけです。

実際のアルゴリズムでは,これをもう少し簡便にした二段階推定法が用いられることが多いようです。 二段階推定法では,まず第一段階で閾値τi,τj\tau_i,\tau_j各回答の出現割合から決定してしまいます。 例えば「100人中30人がxi=1x_i=1と回答した」という場合には,標準正規分布のなかで上位30%の人がuiτiu_i\geq\tau_iであったと考えられます。 したがって,閾値τi\tau_iは,標準正規分布における上側30%の( 図 4.14 において点線より右側の面積が30%になる)点となり,これを計算するとおよそ0.524と求めることが出来ます。 そして第二段階では,第一段階で得られた閾値を固定して,相関係数rui,ujr_{u_i,u_j}のみを計算します。

4.3.3 G-P分析・トレースライン

Good-Poor分析とは,特に学力テストの文脈で用いられることのある分析方法です。テストの場合,合計点が高い人ほどGood,低い人ほどPoorとされるわけですが,仮にある項目がテスト全体と同じものを測定しているならば,合計点が高い人ほど正答率が高く,低い人ほど正答率が低い,ということが言えそうです。…この考え方,先程のI-T相関と非常に良く似ていますね。

G-P分析では,合計点によって全回答者を3群程度に分け,上位(Good)群と下位(Poor)群での正答率の差をチェックします。もしも正答率に差があれば,たしかにその項目はそのテスト全体が言うところの上位と下位を弁別できていそうだぞ,ということになるわけです。

トレースラインは,G-P分析のときと同じように,合計点によって全回答者をいくつかの群に分け,各群での各選択肢の選択率をプロットしていきます。一番左から順に,得点の低い群から並べていった時に,誤答選択肢の選択率は右肩下がりに,正答選択肢の選択率が右肩上がりになっていれば,たしかにその項目はそのテスト全体が言うところの上位と下位を弁別できていそうだぞ,ということになるわけです。

図 4.16 は,仮想的な5件法の項目に対するトレースラインの例です。 この例では,データを合計点で5等分しています。 したがって,各プロットの一番左は「合計点の下位20%の人たちの,この項目への回答の平均点」を表しているわけです。 では一つずつ図を見ていきましょう。

図 4.16: トレースラインの例

図 4.16 の左上の図は,合計点が高い人ほどこの項目にも高い点をつけているため,良い項目だと判断できます。 言い換えると,この項目に対して何と回答したかによって,その人の合計点(その背後にある潜在特性)の高低を正しく評価できている,ということです。 このような項目では,I-T相関も高く出ることでしょう。

図 4.16 の右上の図は,合計点の最上位層のみが高い点をつけていて,他の人達は軒並み低い点をつけている項目です。 このような項目は,床効果が出ている項目と見ることが出来ます。 一方で,この項目は「高得点か否か」を識別できる項目としてはよく機能しているため,この程度の床効果では削除する必要はない可能性が高いでしょう。 I-T相関はちょっと低くなるかもしれませんが,G-P分析を行えば問題ない結果が出るはずです。

図 4.16 の左下の図は,合計点の高低にかかわらず平均値が概ね同じ値になっています。 一般的には,このような項目が悪い項目と判断されます。 この項目に何と回答するかが,その人の合計点(その背後にある潜在特性)の高低と無関係になっているということは,この項目への回答が潜在特性の高低の推定に役立っていないということです。 このような項目ではI-T相関は低く,G-P分析もダメダメになると思われるので,事前に削除を検討しても良いと思います。

図 4.16 の右下の図は,これまでと違って右下りになっています。 ふつう心理尺度では,すべての項目は同じ特性を反映していると考えているため,単一項目への回答と合計点は正の相関になるはずなので,ちょっとおかしいです。 最もありえるのは「逆転項目の処理忘れ」でしょう。 一方で,もしも逆転項目のつもりではない項目でこのような関係が見られた場合は,ワーディングなどに問題を抱えている可能性が考えられるので,使用を慎重に検討したほうが良いかもしれません。

ここからは,実際に平均点のトレースラインを作成してみたいと思います。 トレースライン作成のために必要な工程は以下のとおりです。

  1. 合計点を用いてデータをいくつかのグループに分ける
  2. 各グループについて,特定の項目への回答の平均値を出す
  3. X軸を1. で求めたグループ(小さい順),Y軸を2. で求めた平均値としてプロットする

工程1. では新しい材料として,cut()関数を使用します。 この関数は,与えられたベクトルの中身を小さい順に並び替えてbreaks等分したときに,各レコードが下から何番目のグループに入るかを教えてくれる関数です。

コード 4.25: cut(): データを等分する
# まずはやってみる
# labelsにbreaksと同じ長さのベクトルを与えると表示がわかりやすい
cut(dat$Q1_A, breaks = 10, labels = 1:10)
  [1] 6  7  6  8  6  8  8  4  9  8  7  8  4  7  7  7  6  8  3  8  10 9 
 [23] 8  8  6  10 6  8  9  10 9  8  5  9  6  10 8  9  8  6  10 10 5  8 
 [45] 6  10 10 9  7  6  10 8  5  6  8  7  7  7  6  10 7  3  8  9  8  6 
 [67] 8  8  9  10 8  5  8  8  6  8  10 8  7  9  9  9  9  8  8  9  8  6 
 [89] 3  9  6  6  6  7  8  6  10 1  9  8 
 [ reached getOption("max.print") -- omitted 2332 entries ]
Levels: 1 2 3 4 5 6 7 8 9 10

試しに,89-100行目の人(上の出力で一番下の行に表示されている人たち)のdat$Q1_Aの値とcut()関数の返り値を見比べてみると,確かにdat$Q1_Aが大きいID=106番の人のcut10になっており,またその次の(dat$Q1_Aが小さい)108番の人のcutは確かに小さな値(1)になっていることがわかります( 表 4.3 )。

表 4.3: 89-100行目の人の合計値(Q1_A)とcut()の返り値
Error in dimnames(x) <- dn: length of 'dimnames' [2] not equal to array extent
Error in cat(tbl): argument 1 (type 'closure') cannot be handled by 'cat'
# これをデータにくっつけておく
dat$group_A <- cut(dat$Q1_A, breaks = 10, labels = 1:10)

工程2. において便利なのがaggregate()関数です。 この関数は,指定した変数の値ごとにグルーピングしたうえで,それぞれのグループに対して計算を行ってくれる関数です。 ちょっとわかりにくいと思うので,以下のコード例を見てみてください。

コード 4.26: aggregate(): 指定した変数の値でグルーピングしてから関数を適用する
aggregate(Q1_1 ~ education, dat, mean, na.rm = TRUE)
  education     Q1_1
1         1 4.419192
2         2 4.296000
3         3 4.581395
4         4 4.832370
5         5 4.953168

aggregate()関数には何通りかの記法があるのですが,ここではformulaを使った記法で紹介します9。 もともとformulaは回帰分析の式(y=b1x+b0y=b_1x+b_0)に合わせた形で,左辺に被説明変数,右辺に説明変数を書くというR特有の記法の一つです。 aggregate()関数では,この記法を利用して左辺に計算対象の変数名,右辺にグルーピング対象の変数名を書き,その後ろに「どのデータについて実行するか」「どの関数を実行するか」をそれぞれ書いていきます。 したがって先程のコードは,「datについて,educationの値ごとにグループを作り,Q1_1mean(na.rm=TRUE)を計算してくれ」という指示になっていたわけです。

ちなみにaggregate()では,複数の属性について同時にグルーピングを行い関数を適用することも出来ます。 この場合,回帰分析の説明変数を追加する要領で,右辺に変数名を足していくだけでOKです。 ここでは使いませんが,以下のように「gender×\timeseducationごとの平均点」といった計算も可能です。

コード 4.27: 複数属性でaggregate()
aggregate(Q1_1 ~ education + gender, dat, mean, na.rm = TRUE)
   education gender     Q1_1
1          1      1 4.146341
2          2      1 4.195652
3          3      1 4.184818
4          4      1 4.475000
5          5      1 4.577778
6          1      2 4.612069
7          2      2 4.354430
8          3      2 4.737047
9          4      2 5.022124
10         5      2 5.175439

話を戻すと,今やりたいことは「Q1_Aの値によって分割された各グループについて,特定の項目への回答の平均値を出す」でした。 aggregate()を使えば,以下のように実装することが出来ます。

コード 4.28: グループ平均値を出す
group_mean <- aggregate(Q1_1 ~ group_A, dat, mean, na.rm = TRUE)
group_mean
   group_A     Q1_1
1        1 1.250000
2        2 2.476190
3        3 2.937500
4        4 3.154545
5        5 3.412281
6        6 3.801858
7        7 4.271084
8        8 4.527216
9        9 5.262295
10      10 5.821596

あとはこれをプロットするだけです。

コード 4.29: 平均値のトレースラインの作成
plot(group_mean)
図 4.17: トレースラインに見える?

group_meanには2つの変数が入っているので,これで散布図が書けると思ったのですが,なんだか思ってたのとはちょっと違う図( 図 4.17 )が出てきました。 実はこれは,aggregate()関数によってグルーピングされた変数の列がnumeric型ではなくfactor型になってしまっているためです。

str(group_mean)
'data.frame':   10 obs. of  2 variables:
 $ group_A: Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10
 $ Q1_1   : num  1.25 2.48 2.94 3.15 3.41 ...

Rではこのように,関数の挙動によって想定していたものと異なる型の変数に勝手に置き換わってしまうことがちょくちょくあるので,ここは逃げずにきちんと変数の型を変更してから再度プロットしましょう10

コード 4.30: 型を修正してから再度プロット
group_mean$group_A <- as.numeric(group_mean$group_A)
plot(group_mean)
図 4.18: トレースラインっぽくなってきた

あとは 図 4.18 の点どうしを線で繋げば, 図 4.16 で見たようなトレースラインの完成です。 そのためには,以下のように引数typeを与えてあげると良いでしょう( 図 4.19 )。

コード 4.31: plot()の引数type
plot(group_mean, type = "b")
図 4.19: トレースラインの完成

引数typeとして良く用いられるのは,p(点をうつ),l(線で結ぶ),bplを両方),n(何も表示しない)あたりでしょうか。 他にもplot()には「線の太さ・色・形」「点の大きさ・色・形」「軸の名前」「図の上部に表示する名前」「目盛りの大きさ」を始めとして様々な設定が可能ですが,ここでは取り上げません11

完成したQ1_1についてのトレースライン( 図 4.19 )は, 図 4.16 でいえば左上のように,いい感じに右上がりの線になっています。 したがって,Q1_1の回答は,因子Aの他の項目Q1_2からQ1_5と同じ傾向にある,言い換えると,これらの項目と概ね同じような因子について尋ねている項目であると言えそうです。

…ということで余裕があれば,この調子で全ての項目についてトレースラインを描いてみてください。

dat <- dat[, colnames(dat) != "group_A"]

【おまけ1】こういうときこそ関数化

トレースラインを作成する方法はわかりましたが,実際にはこの作業を項目の数だけ行う必要があります。 一つ一つの工程はシンプルですが,いちいち複数の工程をこなす必要があり,結構面倒ですよね。 こういったときこそ,同じ作業を関数化することを考えてみましょう。 以下に,トレースライン作成のコードを関数化するための(私なりの)手順を紹介しておきます。 「関数化」というのは,基本的に同じ作業を何回も書くのが面倒だという面倒くさがりを救済するための技術なので,「われこそは面倒くさがりだ」という人は是非習得してください。

関数化する場合も,最初は普通にコードを作成します。 ということで,先程紹介したトレースライン作成のコードを以下にまとめて示しておきます。

コード 4.32: トレースラインの作り方まとめ
# 1. 合計点で10個にグルーピングしたものをdatに追加する
dat$group <- cut(dat$Q1_A, breaks = 10, labels = 1:10)
# 2. ある項目について,グループ平均値を出す
group_mean <- aggregate(Q1_1 ~ group, dat, mean, na.rm = TRUE)
# 3.1 factor型になっているところをnumericまたはinteger型に修正する
group_mean$group <- as.numeric(group_mean$group)
# 3.2 プロットする
plot(group_mean, type = "b")

きちんと目的を果たせるコードが書けたら,このコードのうち変わる部分はどこかを考えます12。 上記のコードでいえば,「どの項目についてトレースプロットを描くか」と「cut()でどの合計点(A,C,E,N,O)を使うか」の2点が変わるところです。 そこで,この「変わる部分」を一旦変数に代入します。

コード 4.33: 変数に置き換える
# データフレームの中の変数(列)名を指定したい場合はchr型が一般的
# トレースプロットを描きたい変数
x <- "Q1_1"
# 合計点
total <- "Q1_A"

そして,作成した変数x, totalを使って先程のコードを置き換えていきます。

コード 4.34: 変数x, totalに置き換えたバージョン
# 列をdat$の代わりにchr型で指定する
# 1. 合計点で10個にグルーピングしたものをdatに追加する
dat$group <- cut(dat[, total], breaks = 10, labels = 1:10)
# 2. ある項目について,グループ平均値を出す
group_mean <- aggregate(x ~ group, dat, mean, na.rm = TRUE)
# 3.1 factor型になっているところをnumericまたはinteger型に修正する
group_mean$group <- as.numeric(group_mean$group)
# 3.2 プロットする
plot(group_mean, type = "b")

ただ,実際にやってみるとわかりますが,実はこのままではうまくいきません。

Error in model.frame.default(formula = x ~ group, data = dat): variable lengths differ (found for 'group')

エラーメッセージを見てもわからないと思いますが,実はaggregate()関数の第一引数をformulaとして与える場合,「Q1_1 ~ group」のように" "でくくらない(character型ではない)必要があります13。 しかし,x <- "Q1_1"を用いてx ~ groupを書き換えると,「"Q1_1" ~ group」となってしまうために,正しくformulaとして解釈されないのです。

この場合の対処法は,formula()関数を用いることです。 例えば以下のように,第一引数にformula()関数を噛ませたcharacter型(文字列)として与えることで,文字列のまま正しく機能させることが出来ます。

コード 4.35: formula()関数を使う
# 同じ結果になる
aggregate(formula("Q1_1 ~ group"), dat, mean, na.rm = TRUE)
   group     Q1_1
1      1 1.250000
2      2 2.476190
3      3 2.937500
4      4 3.154545
5      5 3.412281
6      6 3.801858
7      7 4.271084
8      8 4.527216
9      9 5.262295
10    10 5.821596
# 元のコード
aggregate(Q1_1 ~ group, dat, mean, na.rm = TRUE)
   group     Q1_1
1      1 1.250000
2      2 2.476190
3      3 2.937500
4      4 3.154545
5      5 3.412281
6      6 3.801858
7      7 4.271084
8      8 4.527216
9      9 5.262295
10    10 5.821596

ということで,先程のコードを更に書き換えていきます。

コード 4.36: formula()関数を使って修正
# 列をdat$の代わりにchr型で指定する
# 1. 合計点で10個にグルーピングしたものをdatに追加する
dat$group <- cut(dat[, total], breaks = 10, labels = 1:10)
# 2.1 formulaの文字列を作る
fml <- paste0(x, " ~ group")
# 2.2 ある項目について,グループ平均値を出す
group_mean <- aggregate(formula(fml), dat, mean, na.rm = TRUE)
# 3.1 factor型になっているところをnumericまたはinteger型に修正する
group_mean$group <- as.numeric(group_mean$group)
# 3.2 プロットする
plot(group_mean, type = "b")

これでようやく完成です。 うまくできていれば, 図 4.19 と同じトレースラインができるはずです。

最後に,「変わる部分」を引数として,このコードをまるごと関数に入れてしまいましょう。 コメントは残しておくとあとあと見返したときにわかりやすいのでおすすめです。

コード 4.37: 全体を自作関数として定義
traceline <- function(x, total) {
  # 1. 合計点で10個にグルーピングしたものをdatに追加する
  dat$group <- cut(dat[, total], breaks = 10, labels = 1:10)
  # 2.1 formulaの文字列を作る
  fml <- paste0(x, " ~ group")
  # 2.2 ある項目について,グループ平均値を出す
  group_mean <- aggregate(formula(fml), dat, mean, na.rm = TRUE)
  # 3.1 factor型になっているところをnumericまたはinteger型に修正する
  group_mean$group <- as.numeric(group_mean$group)
  # 3.2 プロットする
  plot(group_mean, type = "b")
}
コード 4.38: 自作関数を使うぞ
traceline("Q1_16", "Q1_N")
図 4.20: traceline()関数の出力

関数化する必要があるのは,トレースラインの作成のように,同じような処理を繰り返す必要がある場合です。 同じ処理を何度も書いてしまうと,あとから修正を加えた際に一部だけ修正して,残りを修正し忘れる,といって面倒が発生してしまいます。 コードの保守や可読性の観点から,同じような処理はなるべく繰り返し記述しないようにすることが望ましいわけです。 このルールは,プログラミングの世界では,よくDRY(Don’t Repeat Yourself)原則と呼ばれたりしています。

ということで,関数化の手順をまとめると,以下のようになります。

関数化の手順(諸説あり)
  1. まずは愚直にコードを書く
  2. 繰り返しにおいて変わる場所を変数に置き換える
  3. うまく動くように必要に応じてコードを修正する
  4. function()関数を使って関数を作成する

特に3. の部分に関しては,毎回どのようなエラーが発生するかは分かりません。 したがって,Rの経験値が重要になってきます。 同じ結果をもたらす処理の書き方が複数通り存在するということを知っていることで,解決策が見いだせるかもしれないのです(traceline()の例の場合,データフレームの列名指定は$ではなくcharacter型を使った方法の可能であることや,formulaの指定時にcharacter型にformula()を噛ませる方法もあるいうこと)。

dat <- dat[, colnames(dat) != "group"]

【おまけ2】選択肢ごとのトレースライン

一応,以下に「選択肢ごとの選択率のトレースライン」のコード例を載せておきます14。ただこれは,多肢選択形式の(つまり誤答選択肢の間には順序関係が無い)問題について誤答選択肢を分析するような場合には有効ですが,心理尺度の項目分析の場合はここまでやらなくても良いかもしれません。

コード 4.39: 選択肢ごとのトレースラインの作成
# 今回は10群に分けてやってみる
group_A <- cut(dat[, "Q1_A"], breaks = 10)
# グループごとに合計点の平均値を出す(X軸用)
group_mean <- aggregate(Q1_A ~ group_A, dat, mean, na.rm = TRUE)
# 各選択肢のグループごとの選択率を出すために関数を作る
# ↓ベクトルの中で値がkの割合を計算する関数
prop_k <- function(vec, k) {
  mean(vec == k)
}
# 選択肢1のグループごとの選択率を出す
# 一時的な(無名)関数の定義を利用する
group_prop1_item <- aggregate(Q1_3 ~ group_A, dat, prop_k, 1)
# 他の選択肢についても
group_prop2_item <- aggregate(Q1_3 ~ group_A, dat, prop_k, 2)
group_prop3_item <- aggregate(Q1_3 ~ group_A, dat, prop_k, 3)
group_prop4_item <- aggregate(Q1_3 ~ group_A, dat, prop_k, 4)
group_prop5_item <- aggregate(Q1_3 ~ group_A, dat, prop_k, 5)
group_prop6_item <- aggregate(Q1_3 ~ group_A, dat, prop_k, 6)
# X軸はグループごとの「合計点」の平均値
# Y軸は各選択肢のグループごとの選択率
# type="n"とすると,何も表示しない=枠だけ作る
plot(x = c(min(group_mean$Q1_A), max(group_mean$Q1_A)), y = c(0, 1), type = "n", xlab = "グループ平均値", ylab = "選択率")
# 選択肢ごとに線と点を別々に描き足していく
# 選択肢1
lines(x = group_mean$Q1_A, y = group_prop1_item$Q1_3) # 線
points(x = group_mean$Q1_A, y = group_prop1_item$Q1_3, pch = "1") # 点
# 選択肢2
lines(x = group_mean$Q1_A, y = group_prop2_item$Q1_3) # 線
points(x = group_mean$Q1_A, y = group_prop2_item$Q1_3, pch = "2") # 点
# 選択肢3
lines(x = group_mean$Q1_A, y = group_prop3_item$Q1_3) # 線
points(x = group_mean$Q1_A, y = group_prop3_item$Q1_3, pch = "3") # 点
# 選択肢4
lines(x = group_mean$Q1_A, y = group_prop4_item$Q1_3) # 線
points(x = group_mean$Q1_A, y = group_prop4_item$Q1_3, pch = "4") # 点
# 選択肢5
lines(x = group_mean$Q1_A, y = group_prop5_item$Q1_3) # 線
points(x = group_mean$Q1_A, y = group_prop5_item$Q1_3, pch = "5") # 点
# 選択肢6
lines(x = group_mean$Q1_A, y = group_prop6_item$Q1_3) # 線
points(x = group_mean$Q1_A, y = group_prop6_item$Q1_3, pch = "6") # 点
図 4.21: 選択肢ごとのトレースラインの例(Q1_3

心理尺度項目では,各選択肢は順序変数になっています。ということは良い項目の場合, 図 4.21 のように,各選択肢のピークが順番に並ぶはずです。もしもX軸の値に関係なく「どちらでもない」の選択率が高い,といった事があれば,その項目では潜在特性以外の何かしらの理由で「どちらでもない」を選択している人が一定数いそうだ,ということを確認したりはできると思います。

学力テストの場合は,ふつう選択肢に順序関係はありません。したがって,「どの選択肢が誤答として良く働いているか」や「下位層が特に引っかかりやすい誤答選択肢がないか」など,選択肢レベルでの詳細な検討を行うことが出来るわけです。

4.4 信頼性と妥当性の考え方

デジタル体重計を開発する場合には,必ず出荷前に品質チェックを行うでしょう。その内容はきっと「10kgのおもりを載せた時にきちんと10kgと表示するか」や「様々な環境下で測定しても安定して同じ重さを表示できるか」といったものだと推測されます。

これと同じように心理尺度の開発においても品質チェックが重要になります。構成概念という目に見えないものを測定しようとする心理尺度の場合,そのチェックは難解を極めます。例えば「10kgのおもり」のように客観的にみて正しいと証明されている基準がないことや,通常複数の項目の集合で単一の構成概念を測定すること,測定の安定性を脅かす要因が多様であることなどがあるために,完全な形での品質チェックはほぼ不可能と言っても過言ではありません。

測定の質は,大きく分けると妥当性(validity)信頼性(reliability)という2つの概念のもとで議論されます。 簡単に言うと,妥当性は「測定値が,測定しようとしている構成概念を正しく表した値か」の程度です。例えば「デジタル体重計」という名前で身長を測定する機械があったとすると,これは明らかに妥当性がありません。さすがに体重計と身長計を間違えることは無いでしょうが,心理尺度の場合,なにせ構成概念が目に見えないので,「外向性」のつもりで作った尺度が実は「攻撃性」を測定していた,ということが起こる可能性はゼロとは言えないのです。

一方,信頼性は「測定値が,一貫して同じ値を表しているか」の程度です。例えば同じ人なのに乗るたびに表示される値が変わる体重計や,同じメーカーの同じ型番の体重計なのに物によって表示される値が変わる体重計があっても,その結果を信頼することは出来ません。心理尺度でも同じことが言えます。

妥当性と信頼性は完全に独立した概念ではありません。これらの関係性は,良く 図 4.22 のようなダーツのアナロジーによって説明されます。赤いマル一つ一つがダーツの刺さった場所を表していると考えてください。ダーツ本来の目的である「真ん中に当てる」ということが,心理尺度でいう「測定しようとしている構成概念を測定する」に相当します。真ん中に近いほど測定しようとしているものに近い,ということですね。

図 4.22: ダーツのアナロジー(信頼性と妥当性の関係)

図 4.22 の一番左は,信頼性も妥当性も低いパターンです。毎回あっちへ行ったりこっちへ行ったりで,測定が全く安定していません。毎回ぜんぜん違う結果が出てくるので,特定の回で見ればたまたま測定したい真の値に近い結果が出ることもあるかもしれませんが,結局いつ当たったかがわからない以上使い物にはならないでしょう。

図 4.22 の真ん中は,信頼性は高いが妥当性は低いパターンです。「デジタル体重計」という名前で身長を測定する機械があれば,これに該当するでしょう。毎回ほぼ同じ値を返してくるが,毎回「測定しようとしているもの」とは異なる値を返している,という状態です。この場合は,元々想定していた名前の尺度としては使えないかもしれませんが,もしこの尺度が本当は何を測定していたかがわかれば,(マト自体をうごかしてやることで)別の尺度としては使えるようになるかもしれません。

図 4.22 の一番右は,信頼性も妥当性も高いパターンです。この状態に持っていくことができれば,いつどこで測定しても安定して「測定しようとしている構成概念」を測定できるため,安心して心理尺度を利用することが出来ますね。

ということで,心理尺度を作成・利用する際には,信頼性も妥当性も高くないといけないということがわかりました。また,ダーツのアナロジーからは信頼性は妥当性の必要条件だということがわかります。つまり「妥当性が高いのに信頼性が低い」というケースはあり得ない,ということです。

また,妥当性・信頼性を考える際に注意しなければいけないのは,妥当性・信頼性は尺度自体が持つ性質ではなく,尺度得点に関する性質として評価されるという点です Streiner & Norman, 2008Appendix A では,心理尺度が翻訳や時代の変化によって変質してしまう危険性を指摘しました。 つまり,ある集団において「良い」と判断された尺度があったとしても,それが自分の集めたサンプルについても「良い」尺度であるかどうかはわからない,ということです。 したがって,先行研究ですでに妥当性や信頼性が検証済みの尺度を利用する場合でも,基本的には自分で集めたデータにおける妥当性・信頼性は改めて評価する必要があると言えるでしょう。 もちろん,先行研究で検証に使用されたサンプルの母集団と,自分が利用するデータの母集団が(測定したい内容に関して)概ね同じであることが示せるならばいちいち再評価する必要はありません。 まとめると,どこかから尺度を探して利用する場合には,妥当性・信頼性は改めて評価する必要があるかよく考える必要があるということです。

それでは,ここからは具体的に妥当性・信頼性を検証するための方法について見ていきましょう。

4.5 妥当性の検証

妥当性の考え方は,未だにちょくちょく変化しています。 その歴史的な変遷については例えば 村山 (2012Streiner & Norman (2008 などが詳しいですが,妥当性の考え方自体にトレンドのようなものがあり,「これぞ」という万人が納得できる答えはまだ見つかっていません。 そしてこのあたりの話には科学哲学も入り込んでくるかなり複雑な話になるため,ここでは深く取り上げないことにして,まずは現代の妥当性の主流となっていると思われる Messick (1995 の考え方を紹介します15

4.5.1 妥当性の考え方

もともと心理学界隈では,妥当性にはいくつかのタイプが存在すると考えられていました 村山,2012

1つ目は基準連関妥当性 (Criterion-related Validity)です。 これは,ある尺度によって測定された得点が,すでにその効果を認められている外的な基準とどの程度相関しているかによって評価されます。 何かを測定するための方法が確立しているが測定自体が大変(お金や時間がかかる,などの)ときに,それよりも簡便な方法で同じものが測定できたら嬉しい,という考えが背景にあります。例えば認知症の診断を行う場合には,MRIを撮って脳の萎縮度を見ればほぼ確実な診断ができますが,これには相当なお金がかかってしまいます。 もしもこれが,認知機能を測るいくつかの質問項目でほぼ同等の診断ができるとしたら,とても嬉しいわけです。 つまり,基準連関妥当性では,尺度を別の測定の代理変数 (proxy) として考えている節がありそうです。

2つ目は構成概念妥当性 (Construct Validity)です。 基準連関妥当性の考え方と比べると,構成概念妥当性では,ある尺度によって測定された得点を,その背後にある構成概念の顕在化した指標という感じで考えます。 そして,構成概念の世界において理論的に導出される仮説に合致する結果が得られた場合に,「確かにその尺度は想定していた構成概念を測定する尺度なんだ」と言えるわけです。 例えば「きょうだいの数」が多いほど大人になったときの「外向性」が高い,というナントカ理論が存在したとします。 このとき,「外向性」を測定するために作成した(がまだ本当に外向性を測定できているかわからない)尺度の得点と「きょうだいの数」に強い正の相関が見られたとしたら,この尺度はたしかに「外向性」を測定出来ている気がしてきます。 もちろん,このたった一つの結果だけで尺度の妥当性が示せるわけではないですが,このような要領で,構成概念の世界において理論的に導出される様々な仮説を一つ一つ検証していくことで,その尺度が「もともと測定したかったもの」を表している可能性が高められていくわけです。 俗に「収束的妥当性 (Convergent validity)」や「弁別的妥当性 (Discriminant validity)」と呼ばれるものは,「相関が高い(低い)構成概念の間では,それらを測定する尺度同士の得点も相関が高い(低い)はず」という構成概念妥当性の考え方にもとづいたものです。

3つ目は内容的妥当性 (Content Validity)です。 これは,構成概念妥当性と同様に,ある尺度は特定の構成概念を測定するためのものとして考えます。 ですが,内容的妥当性では「他の構成概念との関係」などではなく,その尺度自体の内容から妥当性を判断します。 具体的には尺度の項目セットが,構成概念が指し示すあらゆる現象の全体からまんべんなく・過不足なくサンプリングされたものであれば「内容的妥当性が高い」と判断することが出来ます。 例えば「高校数学の理解度」を評価するための(大学入試の)数学のテストでひたすら連立方程式問題ばかりが出題されたとしたら,これは「高校数学の理解度」と言うよりは「連立方程式の理解度」になってしまいます。 このような場合には,内容的妥当性が低いと判断されてしまうわけです。

以上のように,妥当性には基準連関妥当性,構成概念妥当性,内容的妥当性の3つのタイプがある,と長らく考えられてきました。 しかしこの考え方は一方で「この3つの妥当性をそれぞれクリアしたら良いんだな」というように,妥当性検証のプロセスをスタンプラリー化してしまった Landy, 1986といった批判を生み出してしまいます。

こうした批判を踏まえて Messick (1995 は,妥当性とは「構成概念妥当性」そのものであるという主張をしました。 つまり妥当性を「尺度得点の意味や解釈に関する様々な証拠を統合したもの」と定義したわけです。 ここで重要なのは,妥当性が尺度得点をつかって行う解釈や推論がどれだけ信頼できるものかを表すための指標として考えられている点です。 そしてこの定義に従うと,基準連関妥当性や内容的妥当性はいずれも,(構成概念)妥当性を示す傍証の一つとして考えることが出来ます。 その意味で,これらは妥当性の「基準連関的証拠」や「内容的証拠」と呼ばれることがあります。

Messickの考え方では,妥当性の証拠としてどのようなものが重要視されるかは,尺度の目的などによって変わります。 例えば「同じ構成概念を少ない項目数で測定したい」という短縮版尺度の目的からすると,元の尺度との相関によって示される「基準連関的証拠」が重要になるでしょう。 一方で,学力テストのように総合的な能力を評価したい尺度では,専門家などによって内容的側面の評価が重要になるはずです。 このように,本来妥当性の評価は「お決まりのパターン」があるわけではなく,目的と対象に応じて必要な証拠を収集していく作業です16

といった議論を踏まえて,ここからはアメリカ心理学会が発行している「テストスタンダード17American Educational Research Association et al., 2011に掲載されている妥当性の証拠のパターンについて紹介したいと思います。

内容について

内容に関する証拠(Evidence based on test content,内容的妥当性とも呼ばれる)は,データから確認するというよりも,尺度作成時の手続きに関する話です。構成概念および下位概念の定義を明確にした上で,その定義を満たすと考えられる内容に関する項目をきちんと用意する必要があり,その手続きに関して論文などの中での説明が求められます。いわゆる「複数の専門家によって議論した」であったり,項目案を先行研究から引用したり,という話は,この点をクリアするために必要な記述,というわけですね。

認知的なこと

回答過程に関する証拠(Evidence based on response processes)は,回答者が各項目に対して,作成者の想定していた通りに「測定しようとしている構成概念」のことを考えていたか,といった側面のことです。テストスタンダード内でもこの点に関する記述はかなり少なく,また実際の論文などでこの側面に関して妥当性検証を行っているものは見たことが無い(かなり少ない?)のですが,検証する方法としてCognitive interview (e.g., Beatty & Willis, 2007; Peterson et al., 2017といった方法があるようです。例えば「回答中に,(構成概念)についてどのようなことを考えていましたか?」と言った内容を事後的にインタビューする方法のようです。かなり大変そうですね…

内部構造

内部構造に関する証拠(Evidence based on internal structure,因子妥当性とも呼ばれる)は,いわゆる「因子構造」のことを指していると思われます。事前に因子の構造(下位概念)について仮定を置いて項目を作成した場合には,因子分析をした結果,項目が想定通りの分かれ方をしていれば良さそうな気がする,ということです。 この視点に関連する分析(因子分析)の結果を「妥当性の証拠」という位置づけで提示する研究は心理学ではあまり見られないですが,マーケティング系などでは割と用いられている気もします18。 「構成概念をこう定義して,下位概念としてこれとこれがあると考える」という仮定のもとで,項目が実際に下位概念と適合する形で分かれるようであれば,妥当性の傍証にはなりそうです19

また,もしも特定の属性の人では異なる挙動を示す項目(特異項目機能: DIF, Differential Item Functioning)や尺度があったとすると,これは妥当性を脅かす要因の一つになりえます。 DIFの対象となりうる属性の例としては,年齢や性別,人種や宗教など様々なものが考えられます。 DIFが発生しそうな項目の例も挙げておきましょう。例えば「他者に助けを求める傾向」という尺度の中に「苦しいときには神に祈る」という項目があったとします。項目作成者としては「他者に助けを求める傾向が強いほど,神にもすがるようになるだろう」という感覚で項目を作成したとしても,実際にはたぶん宗教観によってその回答は大きく変わるでしょう。 本来,尺度項目への回答は想定される構成概念(因子)の強さのみによって規定されるべきですが,DIFがある場合にはこれに加えて特定の属性が影響を与えていると考えられます。すなわち,その特定の属性が意図していなかった因子として因子構造に混ざってしまっている,という状態です。先程の例で言えば,その項目だけが「他者に助けを求める傾向」と「宗教観」という2つの因子の影響を受けてしまっている,ということになり,これは当初想定していた因子構造から逸脱してしまっている状態です。 こういったことがないかを確認するために,分析の際にはデータを属性ごとに分割してそれぞれで分析を行ったり,SEMの枠組みの中で多母集団同時分析を行ったりするわけですね。

他の変数との関係

他の変数との関係に関する証拠(Evidence based on relations to other variables,基準連関妥当性とも呼ばれる)は,たぶん妥当性検証の中では最もよく提示される証拠です。概念的に関連があると考えられる別の尺度や指標との相関が高いこと(収束的証拠)や,概念的に関連がないと考えられるものとの相関が低いこと(弁別的証拠)をもって,妥当性の証拠とします。

先程説明した通り,比較対象の変数は,心理尺度に限りません。例えば「外向性」の妥当性検証であれば,「SNSの友だちの数」や「直近1年で何回パーティーに行ったか」といった行動指標なども考えられるでしょう。 むしろ心理尺度どうしの相関では,両方の変数の誤差項には共通の「自己申告によるバイアス」が含まれることで,誤差項同士が相関してしまう可能性があります。その結果,本当の相関よりも過大に推定されてしまうかもしれません(common method bias, Podsakoff et al., 2003。そういった意味で,心理尺度以外の変数(特に行動指標)との関連を考えてみるのも良いでしょう。

また,「予測」という観点でもこの妥当性は重要な意味を持ちます(予測的妥当性)。例えば入社試験での性格検査は,「外向的な性格のひとだから営業に配属しよう」的な意思決定を行うために利用されているかもしれません。言い換えると,この性格検査には,最もパフォーマンスが良くなる配属先を予測する機能が求められていることになります。この場合,性格検査が「測定しようとしている構成概念を測定できている」ことよりも,「それによって将来のパフォーマンスが予測できている」ことのほうが重要になります。 この場合,やることとしては尺度得点と「1年後の営業成績」の相関を見るということで,やることは先程と同じですが妥当性の観点は異なっていたりします。

4.5.2 まとめ

妥当性を検証する最良の方法,というものは見つかっていません。というか多分一つに決めるのは不可能です。 したがって,上に紹介したような妥当性検証の方法についても,お決まりのルーティーン実行したら十分(スタンプラリー)というわけではなく,可能であれば複数の観点から検証した上で「これだけ証拠が出揃ってるんだからいいだろう」的なスタンスで行くことが望ましいですね。

4.6 信頼性の検証

信頼性は,妥当性よりは統計的に評価することが出来ます。そのために,古典的テスト理論 (Classical Test Theory [CTT]: Lord & Novick, 1968に基づく以下のモデルを考えます。

x=t+e(4.1) x = t + e \qquad(4.1) (4.1)式において,xxは実際に観測される得点を,eeは測定誤差を表します。そしてtt真値と呼ばれます。言い換えるとttこそが「測定しようとしている構成概念」を直接的に反映した値である,ということです。理想的には,真値ttのみによって回答値が決まれば良いのですが,それに対して何らかの測定誤差が混ざることで,実際に観測される得点が決まる,という状況を仮定しています。

このモデルにおける重要な(そして自然に受け入れられる)仮定として,真値と測定誤差は無相関ということが置かれます。したがって,観測値xxの分散は σx2=σt2+σe2(4.2) \sigma^2_x = \sigma^2_t + \sigma^2_e \qquad(4.2) というように,真値と測定誤差の分散の単純な和として表すことができます。このとき,信頼性係数ρ\rhoρ=σt2σx2=1σe2σx2(4.3) \rho=\frac{\sigma^2_t}{\sigma^2_x}=1-\frac{\sigma^2_e}{\sigma^2_x} \qquad(4.3) と表すことができます。つまり,観測変数の分散に占める真値の分散の割合ということです。もちろん実際に知ることができるのは観測値の分散σx2\sigma^2_xだけであり,真値の分散σt2\sigma^2_tを知ることはできません。

以下で紹介する信頼性係数の求め方は,このような仮定のもとで,信頼性係数ρ\rhoを直接求めようとしています。いくつかの方法があるのですが,測定誤差として何を取り上げているかが方法によって異なるので,その点には少し注意が必要です。なお,ここから先の説明については 岡田 (2015 を拝借しています。

4.2 式の要領で妥当性を定式化する方法も一応存在しています。 CTTでは,観測得点を「真値」と「誤差」に分解していましたが,このうち「誤差」をより詳細に見ていくと,「系統誤差 (systematic error)」と「偶然誤差 (random error)」に分けることが出来ます Maruyama & Ryan, 2014

系統誤差とは,繰り返し測定を行った場合に測定値を常に一定方向に動かすような誤差です。例えば体重計が壊れていて常に実際の体重よりも5%大きな値を表示してしまう場合を考えてみます。 この体重計では,同じ人が何回乗っても必ず実際の体重より少し大きな値が出てしまっているわけなので,これは系統誤差にあたります。

一方偶然誤差とは,繰り返し測定を行った場合に測定値が毎回ランダムに変動してしまうような誤差です。 体重計で言えば,「乗るたびに表示される値が変わる」程度を指しています。

ポイントは,真値と偶然誤差が無相関であるように,系統誤差と偶然誤差も無相関であるという点です。 これを踏まえると観測値の分散は,系統誤差による分散σSE2\sigma^2_{SE}と偶然誤差による分散σRE2\sigma^2_{RE}を用いて(4.2)式からさらに以下のように分解することが出来ます。 σx2=σt2+σe2=σt2+σSE2+σRE2 \begin{aligned} \sigma^2_x &= \sigma^2_t + \sigma^2_e \\ &=\sigma^2_t + \sigma^2_{SE} + \sigma^2_{RE} \end{aligned}

ここで信頼性の定義に立ち返ってみると,信頼性とは「測定値が,一貫して同じ値を表しているか」の程度でした。 この観点から言うと,常に同じだけ値をずらす系統誤差は信頼性に影響がないことがわかります。 したがって,信頼性は ρ=σt2+σSE2σx2=1σRE2σx2 \rho=\frac{\sigma^2_t + \sigma^2_{SE}}{\sigma^2_x}=1-\frac{\sigma^2_{RE}}{\sigma^2_x} と表すことが出来ます。 一方で,妥当性は「測定値が,測定しようとしている構成概念を正しく表した値か」の程度なので,系統誤差であってもそのズレは許容されません。 したがって,妥当性は ρ=σt2σx2=1σSE2+σRE2σx2 \rho=\frac{\sigma^2_t}{\sigma^2_x}=1-\frac{\sigma^2_{SE} + \sigma^2_{RE}}{\sigma^2_x} と表すことができるわけです。

4.3式では,信頼性を観測変数の分散に占める真値の分散の割合として表しました。 しかし実際には真値と系統誤差は区別できないことが多いのです。 例えば体重の真値がわからないときに,ある体重計が「乗るたびに値が変わる」ことはわかっても,「常に5%大きな値を示している」ことを証明するのはかなり難しいでしょう。 なのですが,信頼性の定義からすると,系統誤差の有無は実用上問題ではありません。 そのため,この後紹介する信頼性の評価方法においても,系統誤差は無視して扱っています。

4.6.1 (平行テスト法)

はじめに,真の信頼性係数ρ\rhoを求めるための「理想」のお話を少しだけしておきます。 真の信頼性係数ρ\rhoを求めるためには,真値と誤差分散が同じ2つのテスト(平行テストと呼ぶ)を独立に実施する必要があります。 平行テストがある場合,2つのテスト得点x1x_1およびx2x_2はそれぞれ x1=t+e1x2=t+e2(4.4) \begin{aligned} \begin{split} x_1 &= t + e_1 \\ x_2 &= t + e_2 \end{split} \end{aligned} \qquad(4.4) と置くことが出来ます。この2つのテスト得点の共分散σx1x2\sigma_{x_1x_2}は,真値と測定誤差は無相関であり,また独立に実施しているならば異なるテストの真値と測定誤差も無相関(σt,e2=σe1,t=σe1,e2=0\sigma_{t,e_2}=\sigma_{e_1,t}=\sigma_{e_1,e_2}=0)になるはずなので, σx1,x2=σ(t+e1),(t+e2)=σt,t+σt,e2+σe1,t+σe1,e2=σt2(4.5) \begin{aligned} \begin{split} \sigma_{x_1,x_2} &= \sigma_{(t+e_1),(t+e_2)} \\ &= \sigma_{t,t} + \sigma_{t,e_2} + \sigma_{e_1,t} + \sigma_{e_1,e_2} \\ &= \sigma_t^2 \end{split} \end{aligned} \qquad(4.5)

となり,結果的に真値の分散に一致します。また,誤差分散が同じであれば,2つのテストの得点の分散は一致する(σx12=σx22=σx2\sigma^2_{x_1}=\sigma^2_{x_2}=\sigma^2_x)ため,真の信頼性ρ\rhoρ=σt2σx2=σx1,x2σx1σx2=rx1,x2(4.6) \begin{aligned} \begin{split} \rho &= \frac{\sigma^2_t}{\sigma^2_x} \\ &= \frac{\sigma_{x_1,x_2}}{\sigma_{x_1}\sigma_{x_2}} \\ &= r_{x_1,x_2} \end{split} \end{aligned} \qquad(4.6) ということで,2つのテスト得点の相関係数に一致するため,理論上は計算可能となるわけです。

4.6.2 再テスト法

ただ,真値と誤差分散が同じ2つのテスト・心理尺度なんていうものは,言い換えると内容は異なるが平均点および分散が全く同じテスト・尺度ということであり,そう簡単には作れません。そこで,同じテスト・尺度を2回実施するという考えが浮かんできます。これならば,当然内容は全く同じなので平行テスト法の条件は満たされているといえます。

心理尺度の場合,この再テスト法でもそれなりに妥当に信頼性を測定できるような気がします。しいて言えば,2回間隔をあけて同じ人に実施することのコストがかかりますが,それさえクリアできれば…という感じです。

学力テストの場合,同じ問題のテストを2回実施するのはほぼ不可能です。出題がネタバレしているテストならば,事前準備が容易なため2回目の方がほぼ確実に得点が高くなってしまうでしょう。ということで再テスト法もできそうにない,ときました。さてどうしようか…

4.6.3 折半法

1回のテストで平行テストみたいなものを2つ用意するために,大昔には1つのテストを2つに分けるという考え方が用いられていました。1つの尺度ならば,うまいこと半分に分ければ平行テストに見えなくもない気がします。この折半法では,2つに分けたときのそれぞれの得点の相関rx1,x2r_{x_1,x_2}について,スピアマン・ブラウンの公式にもとづいて,信頼性を ρsh=2rx1,x21+rx1,x2(4.7) \rho_{sh} = \frac{2r_{x_1,x_2}}{1+r_{x_1,x_2}} \qquad(4.7) と計算します。

平行テスト法および再テスト法では,真値tt尺度全体の真値を,誤差分散eeには尺度の外側にある余計なものを想定していました。これに対して折半法での信頼性は,真値tt尺度の半分の真値(尺度の中で共通の成分)を,誤差分散eeには尺度の中の項目間(前半と後半)で異なる成分を表している,と見ることが出来ます20。その意味では真の信頼性ρ\rhoとはちょっと意味合いが異なるのですが,ひとまず折半法によって,1回のテストから信頼性係数が求められるようになった,ということでこの方法が一時期スタンダードだったわけです。 ちなみに,折半法的な意味での信頼性は内的一貫(整合)性などと呼ばれたりします。

ただ,当然ですが折半法では尺度の分け方によって信頼性係数の値が変わってしまいます。これはScienceとしてはあまりよろしくないですね。特に現代のコンピュータを使えば,全ての分け方を試した上で最も高い値を示す分け方で信頼性を主張してしまうことも可能でしょう。これは良くないということで,全ての分け方で算出されるρsh\rho_{sh}の平均値として,KR-20 Kuder & Richardson, 1937 と呼ばれる公式が提案されました。これならば分け方の恣意性は無視できるようになります。ただしKR-20は二値項目にしか適用できない公式だったため,このままではリッカート式の尺度などには使えません。

4.6.4 クロンバックのα\alpha係数

現在最も広く用いられているα\alpha係数 Cronbach, 1951 は,KR-20を多値項目でも使えるようにクロンバックが一般化したものです。nn項目の尺度に対するα\alpha係数の計算式は α=nn1(1i=1nσxi2σx2)(4.8) \alpha=\frac{n}{n-1}\left(1-\frac{\sum_{i=1}^{n}\sigma_{x_i}^2}{\sigma^2_{x}}\right) \qquad(4.8) となっており,これがKR-20と同様に全ての分け方で算出されるρsh\rho_{sh}の平均値になっています。

α\alpha係数は今でも様々な論文で利用されているものですが,様々な批判もあります。

  • ρ\rhoの下界の一つであるため,ほぼ確実にρ\rhoを過小推定する21
  • 項目数が多いほど大きな値になりやすい=項目数の影響を受ける
  • 一次元性の指標とはいえない側面がある

3点目については, 図 4.23 にあるように,実際には尺度が一次元ではなさそうな場合でも,α\alpha係数だけを見ていては判断できないケースがある,ということです。 現在でも多くの論文がα\alpha係数だけを見て「一次元性が確認された」などと供述していますが,実際にはα\alpha係数は一次元性の必要条件にすぎず,もうすこし詳細な確認をすべきだと言えます。

図 4.23: α\alpha係数と一次元性 岡田,2015, Table 1)

そんなα\alpha係数に変わって,現在ではω\omega係数などの信頼性係数や,SEMに基づく方法,一般化可能性理論と呼ばれる理論に基づく方法が提案されています。ただここで説明している余裕と事前知識が不足しているため,後ほど因子分析の回においてω\omega係数を紹介します。

とはいえα\alpha係数はまだ現役で使われている方法なので,ここでもα\alpha係数の使い方を紹介しておきます。Rではpsych::alpha()という関数があります。今回のように一因子ではなく複数の因子がある尺度の場合は,25項目全てではなく各因子ごとにalpha()関数を適用し,それぞれで出た値を報告するようにしましょう。

コード 4.40: α\alpha係数を計算する
# もしまだ読み込んでいなければ
# library(psych)

# 第一因子の5項目についてalpha()
# 残りの因子についても,コードは省略していますがチェックしておいてください。
alpha(dat[, paste0("Q1_", 1:5)])

Reliability analysis   
Call: alpha(x = dat[, paste0("Q1_", 1:5)])

  raw_alpha std.alpha G6(smc) average_r S/N    ase mean   sd median_r
      0.71      0.72    0.69      0.34 2.6 0.0092  4.6 0.91     0.35

    95% confidence boundaries 
         lower alpha upper
Feldt      0.7  0.71  0.73
Duhachek   0.7  0.71  0.73

 Reliability if an item is dropped:
     raw_alpha std.alpha G6(smc) average_r S/N alpha se  var.r med.r
Q1_1      0.73      0.74    0.68      0.41 2.8    0.009 0.0061  0.39
Q1_2      0.63      0.64    0.59      0.31 1.8    0.012 0.0167  0.30
Q1_3      0.61      0.63    0.57      0.30 1.7    0.013 0.0089  0.34
Q1_4      0.69      0.70    0.66      0.37 2.4    0.010 0.0155  0.37
Q1_5      0.66      0.67    0.62      0.34 2.0    0.011 0.0126  0.35

 Item statistics 
        n raw.r std.r r.cor r.drop mean  sd
Q1_1 2432  0.59  0.58  0.39   0.32  4.6 1.4
Q1_2 2432  0.73  0.75  0.67   0.57  4.8 1.2
Q1_3 2432  0.77  0.77  0.72   0.60  4.6 1.3
Q1_4 2432  0.67  0.64  0.49   0.41  4.7 1.5
Q1_5 2432  0.69  0.70  0.60   0.50  4.5 1.3

Non missing response frequency for each item
        1    2    3    4    5    6 miss
Q1_1 0.03 0.08 0.12 0.14 0.30 0.33    0
Q1_2 0.02 0.05 0.06 0.20 0.37 0.32    0
Q1_3 0.03 0.06 0.07 0.20 0.36 0.27    0
Q1_4 0.05 0.08 0.06 0.16 0.24 0.41    0
Q1_5 0.02 0.07 0.09 0.22 0.35 0.25    0

出力を上から見ていきましょう。

2-6行目 (Reliability analysis)

各種信頼性係数の指標です。基本的には最初の2つのalphaを見ておけばOKです。

raw_alpha
データの共分散行列から普通に計算したα\alpha係数です。本来的な定義( 4.8式)はこれなので,基本的にはこれを報告しておけばOKです。
std.alpha
共分散行列の代わりに相関行列を用いて計算したα\alpha係数です。例えば特に項目ごとに最大値・最小値が一定ではないような場合,取りうる値の幅が大きい項目のほうが分散も大きくなるため,α\alpha係数にかかる重みが極端に大きくなってしまいます。 このような場合には,項目の重みを揃えるため,相関行列を用いたこちらのα\alpha係数を報告するほうが良いかもしれません。 ただしリッカート尺度の場合はふつう最小値・最大値は同じなので,raw_alphaと極端に異なるケースは,天井・床効果が発生している場合を除けばほぼ無いと思います。

気になる人もいるかもしれないので解説しておきます。左から

G6(SMC)
俗にガットマンのλ6\lambda_6 Guttman, 1945 と呼ばれる指標です。各項目について「それ以外の項目で重回帰分析した際の重相関係数」をもとに計算する指標です。重相関係数が高いということは,その項目が「他の項目との相関が高い」ことを意味するため,これも内的一貫性の指標になっています。が,α\alpha係数を報告しておけば十分です。
average_r
項目間相関の平均値です。5項目ならば5C2=10{}_5C_2=10通りの相関があり,その平均値です。
s/n
Signal/Noise比です。上のaverage_rrrとすると,nr/(1r)nr/(1-r)で計算されています。
ase
α\alpha係数の漸近標準誤差(asymptotic standard error)だと思われます。
mean
項目レベルでの平均値です。
sd
項目レベルでの標準偏差です。
median_r
項目間相関の中央値です。

8-11行目 (95% confidence boundaries)

raw_alphaの95%信頼区間です。どちらの方法でも別にかまわないですが,余裕があればいずれかの信頼区間を報告するようにしましょう。

13-19行目Reliability if an item is dropped

1項目を除外した時に各種指標がどう変化するかを表しています。今回の場合,raw_alphaを見るとQ1_1だけは5行目に表示されているものよりも高い値を示しています。つまりこの場合,統計的にはQ1_1は除外したほうが内的一貫性は高まる,ということです。 ただ,実際にはこれを見て機械的に項目を除外するのではなく,内容面などを吟味した上で削除するかを決定しましょう。 図 4.24 のように,特定の項目が他の項目と相関が低いとしても,内容的にその項目が無いと構成概念の全要素をカバーしているとは言えなくなるような項目であれば,除外せず残したほうが妥当性は高い可能性があります22。 今回の場合は,raw_alphaも大した改善ではないので,このまま残して分析を進めることにします。

図 4.24: 項目と構成概念の関係

21行目以降Item statisticsおよびNon missing response frequency for each item

項目ごとの要約統計量および各カテゴリの選択率です。まあ,説明することも無いでしょう。ちなみに,23-27行目の一番左のraw.rはピアソンの積率相関係数として計算した場合のI-T相関を指しています。

4.6.5 信頼性係数の基準

信頼性係数は,その定義から0ρ10\leq\rho\leq1の値を取ります。直感的には,1に近いほどよい感じがしますね。実際に,平行テスト法や再テスト法的な意味でのρ\rhoは,1に近いほど良い指標です。一方で,内的一貫性としての信頼性係数は,1に近すぎるとかえって良くないとされています。というのも,「内的一貫性」という言葉から想像できるように,α\alpha係数が1になるのは全ての項目が全く同じものを測定しているような場合です。言い換えると,全ての項目間の相関係数が1の場合になります。 相関係数が1ということは,一方の値によってもう一方の値が完全に説明可能な状態にあるため,実質的に1項目の回答だけわかれば尺度内の残りの全項目の回答が分かってしまう状態です。つまり1項目だけあれば十分であり,尺度として何項目も用意した意味がなくなってしまうことになります。 本来,構成概念の多くは普通に1項目で聞いただけでは十分に測定できないために,その構成概念を表す様々な側面に関する質問を重ねることによって,その共通成分としてあぶり出そうとしています。 それなのに1項目分しか機能していないとしたら,それは尺度の失敗と言っても良いでしょう。

一方,最低でもどの程度の信頼性が欲しいかですが,これも一貫した答えはまだありません。α\alpha係数について言えば,項目数の影響を受けるため,一概にどれだけあれば良い,という話は出来ないはずです。世間では0.7だの0.8だの0.9だの言っている人がいるようですが,最終的には信頼性も連続的な指標のため,明確な閾値を考えるよりは,ありのままの値を報告するくらいの気持ちでいたほうが健全だと思います。…とはいいつつ「少なくとも0.7くらいは欲しいかなぁ」という個人的な気持ちもあったりしますが…

ということで,各項目の分析および尺度の信頼性(・妥当性)のチェックが無事完了したら,いよいよ分析の本編に進むことができます。 次のChapterでは,まずその先で紹介する因子分析・構造方程式モデリング・項目反応理論に通底する「回帰分析」の考え方についておさらいしながら,記法の導入などをしていきます。

コード 4.41: 最後に保存
saveRDS(dat, "chapter04.rds")

  1. 本気できれいな図を作りたい場合は,ggplot2というパッケージがおすすめです。ただ,この授業で紹介する方法とは設計の基本理念から異なっているので,紹介しているヒマがありません。論文用の図を作りたい時が来たら個別に質問してください。↩︎

  2. 実は,apply(dat, 2, mean)と同じ挙動をする,つまり列ごとの平均を出す関数としてcolMeans(dat)というものもあります。そしてrowMeans()colSums()という関数もあります。↩︎

  3. 一般的に信頼性(および妥当性)は尺度全体の性質(または尺度全体×\times母集団の組み合わせに対する性質)として考えられるものなので,項目ごとの性質の分析であるI-T相関の確認は信頼性とは別枠で紹介しています。↩︎

  4. ncol(dat)+1列目を指定して代入することも可能ですが,この場合列の名前を決められないのでやめたほうが良いでしょう。↩︎

  5. もちろん,回答者自身が自分の潜在特性の強さと項目の閾値の関係を考えて回答を選んでいるわけではありません。あくまでも「こう考えたらしっくり来るよね」というイメージ的なものです。↩︎

  6. それなりの項目数があれば,合計点は連続変数とみなしても良いのでそのまま使います。↩︎

  7. polyserial()関数は,「カテゴリカル変数と連続変数の相関を出す」関数です。したがって,ポリシリアル相関の特殊なケース(カテゴリカル変数が二値)である点双列相関もこの関数で対応可能です。 これとは別に,polychoric()関数というものも用意されています。こちらは「カテゴリカル変数同士の相関を出す」関数で,ポリコリック相関とテトラコリック相関係数を出すことができます。↩︎

  8. 現実のデータで,これだけきれいにI-T相関が高く出て(今後やりますが)きれいな因子構造が得られることはまずありません。この授業ではかなり理想に近いデータを扱っていたんだということは,みなさんもいずれ痛感するときが来るかもしれません。↩︎

  9. formula記法は回帰分析やSEMのところでガンガン使用するので,少しでも慣れておくと良いでしょう。↩︎

  10. そして関数の挙動が思った通りにならない場合は,変数の型の間違いを疑うのも一つの手だということです。↩︎

  11. 以前も紹介しましたが,発表や論文用に図を作る場合は,ggplot2というパッケージがおすすめです。plot()はあくまでもサッと図を作る用,と思っておくと良いと思います。とはいえもちろんplot()を駆使してきれいな図を作ることも可能です。お好きな方法を選んでください。↩︎

  12. もちろん,慣れてきたら「cut()でいくつのグループに分けるか」や「プロットの色・太さ」なども変えられるようにもできます。↩︎

  13. 一応,エラーメッセージにはvariable lengths differと表示されています。formula記法では,回帰分析のように左辺と右辺のデータの数は同じになっている必要があるのですが,"Q1_1" ~ groupとなった場合に,右辺はdat$groupつまりdatの行数の長さのデータがある一方で,左辺は(dat$Q1_1ではなく)"Q1_1"という長さ1の文字列になっているために,「変数の長さが違うよ」というエラーが出ているわけです。↩︎

  14. 【おまけ1】の考え方からすると,もっと関数化してシンプルに書くことも可能です。もし暇なら,どのくらいコード量で作成できるか挑戦してみてください。↩︎

  15. 妥当性に関する議論については他にも, Borsboom et al. (2004Borsboom (2005 仲嶺監訳 2022 なども参考になると思います。↩︎

  16. まあそんなものは理想論なのかもしれません。妥当性の定義についても未だに議論は尽きません。重要なのは,先行研究のやり方が常に正しいと思うな,もっと良い妥当性の証拠はあるかもしれないぞと常に考え続けることなのかもしれません…。↩︎

  17. 2024年現在,改訂作業が進んでいるとのことなので,数年後にはまた変わっているかもしれません。↩︎

  18. 因子分析やSEMの結果からCRやAVEと言われる指標を算出し,これが所定の値より高いことをもって「収束的妥当性が示された」という書き方をしているものがちょくちょく見られます。↩︎

  19. もちろん,下位概念の構造を含めてそっくりそのままコピーしたような別の構成概念が存在している可能性は否定できません。つまりCRやAVEだけで満足しているものは正直言うと「不十分」だと思っています。↩︎

  20. なおスピアマン・ブラウンの公式は,尺度の半分に対する信頼性rx1,x2r_{x_1,x_2}を補正することで,尺度全体での信頼性を求めようとしています。その点はちゃんと気をつけています。↩︎

  21. 一応,尺度全体において,項目レベルで(本質的)タウ等価(タウ,つまり真値が同じか,個人によらず一定のズレがある状態)であれば,α\alpha係数は真値ρ\rhoに一致するのですが,そうでない場合は,ρ\rhoを過小推定することが知られています。↩︎

  22. たぶんその場合は,下位概念を導入して複数因子の構造を考えるのが良いかもしれませんね。↩︎