多変量解析_03_データの前処理

1.

Chapter 3 データの前処理 1 Chapter 3 データの前処理データを分析できる状態にするための前処理（欠測値の削除・適当な回答の検出・逆転項目の処理）の考え方および R での実行方法を紹介しています。本資料は，神戸大学経営学研究科で 2022 年度より担当している「統計的方法論特殊研究（多変量解析）」の講義資料です。CC BY-NC 4.0 ライセンスの下に提供されています。作成者連絡先神戸大学大学院経営学研究科分寺杏介（ぶんじ・きょうすけ） mail: [email protected] website: https://www2.kobe-u.ac.jp/~bunji/ Chapter 2 では，データを読み込み，元データのレベルでのチェックを行いました。ただ，まだデータを確認しただけで，このまま分析しても思ったような結果が得られない可能性があります。というのもデータには多かれ少なかれノイズ（ゴミ）が混ざっており，そのまま残しておくと今後の分析に良くない影響を与えかねません。ということで今回は主にデータの前処理を行っていきたいと思います。データの読み込みまずは前回保存したデータを読み込みましょう。 .rds ファイルを読み込む場合には， readRDS() 関数を使います。 .rds ファイルを読み込む 1 2 # 保存したときのオブジェクト名(dat)と同じである必要はありません dat <- readRDS("chapter02.rds") もし save() や save.image() 関数で保存した .rdata ファイルを読み込みたい場合には，

2.

3.1 欠測値の処理 2 load() 関数を使います。まとめて環境を読み込む 1 load("chapter02_all.rdata") load() 関数では，代入の記号（ <-）を使用していません。これは，save() 関数で保存した .rdata ファイルにはオブジェクトが中身から名前まで丸ごと保存されており，load() 関数はそれを丸ごと読み込んでいるためです。一方 saveRDS() 関数で保存した場合にはオブジェクトの中身だけが保存されるため，保存前と別の名前をあてても問題ないわけです。そういった理由から，保存する際の拡張子を .rds と .rdata と使い分けています（と私は思っています）。 3.1 欠測値の処理はじめに無回答の処理の仕方について少し説明します。一般的に分析の際には，無回答は欠測値 (missing value) として扱われます。テクニカルな話をすると欠測値には以下の 3 つのタイプがあり，そのタイプに応じた対処が求められます（Schafer & Graham, 2002）。ここでは，「体重」を尋ねた質問への無回答を例に考えてみましょう。 MCAR Missing Completely At Random とは，欠測値が完全にランダムに発生している状態です。簡単にいえばたまたま体重を記入し忘れた，という場合や明らかに異常値であるためにデータが使えない場合には MCAR として扱うことが出来ます。 MAR Missing At Random とは，欠測値がその変数以外に関連して発生している状態です。体重で言えば，女性の方が体重を答えるのをためらう傾向が強いと思われるので，「女性の方が体重の欠測が多い」と考えられる場合は MAR として扱います。名前に”At Random” と入っているのが少し違和感かもしれませんが，MAR では一応欠測の発生は確率的に発生する（ただし発生確率が他の変数によって変わる）と考えています。 MNAR Missing Not At Random とは，欠測値がその変数自体に関連して発生している状態です。体重で言えば，肥満傾向の人のほうが体重を答えるのをためらう傾向が強いと思われるので，「肥満の人の方が体重の欠測が多い」と考えられる場合や，体重計のスペック的に「測定不能」となってしまった場合は MNAR として扱います。具体的な対処法まで話すと長くなってしまうので省略しますが*1 ，簡単に言うと MNAR でなければ統計的に欠測値を埋める方法があると考えて良さそうです*2 。したがって，特定の属性について欠測値が多くなっていてもバイアスを補正した推定が可能となるわけです。ということで，質問紙調査を行う段階で，MNAR が発生しそうなセンシティブな項目などには特段の注意が必要となります。 *1 欠測データの処理についての日本語の書籍としては，例えば高橋・渡辺（2017）には R のコードも掲載されているので，興味があれば見てみると良いかもしれません。 *2 MCAR の場合には，最悪欠測値がある人をリストワイズ削除しても問題ありません（サンプルサイズが小さくなるだけでバイアスにはならない）。MAR の場合には，多重代入法 (MI) や完全情報最尤推定法 (FIML) を用いると，うまいこと欠測値を補完することができるとされています。

#p31
#p2

3.

3.1 欠測値の処理実際にデータで見られた欠測がどのタイプなのかについては，ドメイン知識やデータの状況などによって総合的に判断する必要があります。先程説明したように，同じ「体重」という変数をとっても，その背後に異なるメカニズムを想像するだけでどのタイプにも該当しそうでした。というわけで，欠測値のタイプを判断するための統計的な指標は今のところほぼ存在しないと思います*3 。ちなみに今回のデータにおける欠測値を考えてみると， • 心理尺度の 25 項目に関する欠測値は MCAR と考えてもよさそう • education の欠測は MNAR のような気もするけど… という感じでしょうか。とりあえず心理尺度の 25 項目に関しては MCAR とみなして，以後の分析では「25 項目の中には欠測値がない人」だけを残して使用したいと思います。今回は全体で 2800 人いるので，多少減っても問題ないでしょう。ただし，一つでも欠測値がある人を除外する場合は，事前に項目ごとの無回答率を確認しておいたほうが良いです。もしかしたら特定の項目が回答しにくくて無回答を誘発していたり，後半の設問ほど無回答率が高くなっていたり…ということがあるかもしれません。こうした状態を無視して機械的に除外してしまうと，「特定の項目に回答できる人」や「飽きっぽくない・回答が早い人」が多めに残ってしまい，何らかの偏りが発生してしまう可能性もあります。 R では欠測は NA として扱われています。そのため，列ごとに NA となっているデータの数（割合）を計算することで無回答率が算出できそうです。ちょっと厄介なのは，NA は比較演算では特殊な挙動をする，という点です。「 NA と一致するか」はシンプルに考えると等号 ( ==) を用いて以下のように書けそうですが，実際にはそうはいきません。等号の演算子を使って NA かどうかを判定すると 1 c(1, NA, 3, NA, 5, 6) == NA 1 [1] NA NA NA NA NA NA このように，NA との比較は「判断不能」的な感じで結果も NA になってしまいます。代わりに R には，NA かどうかを判断する関数として，is.na() という関数が用意されているのでこれを使いましょう*4 。is.na() 関数は，引数に与えたものに関して，NA であれば TRUE，そうでなければ FALSE を返します。ベクトルや行列・データフレームを与えた場合はその要素ごとに TRUE/FALSE を出してくれます。 *3 一応 MCAR かどうかを判別するための統計的検定法として Little’s MCAR test（Little, 1988）というものがあるようです。R には BaylorEdPsych::LittleMCAR() という関数が作られているようです。使ったことはありません。 *4 R では，「あるはずなのに無い（欠測）」は NA (Not Available)，「数字ではないもの (e.g., 0 で割ったもの)」は NaN (Not A Number)，「もともと何もない」は NULL，無限大は Inf という特殊な非数値が用意されており，それぞれ is.***() 関数で判定可能です。 3

#p3
#p30

4.

3.1 欠測値の処理 4 is.na(): NA であるかを判定する 1 is.na(c(1, NA, 3, NA, 5, 6)) 1 [1] FALSE TRUE FALSE TRUE FALSE FALSE これを使えば，それぞれの値が「 NA であるか」を判定できるようになります。あとは NA の数を数えるだけです。 is.na() 関数の出力は logical 型（ TRUE/FALSE しかとらない）のベクトルなのですが，実は内部的には TRUE は 1，FALSE は 0 としても扱われます。 1 列まるごと is.na() 1 is.na(dat[, "education"]) 1 2 3 4 5 6 [1] TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE [12] TRUE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE [23] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE [34] TRUE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE [45] FALSE FALSE FALSE FALSE FALSE FALSE [ reached getOption("max.print") -- omitted 2750 entries ] is.na() を強制的に数値に変換してみると 1 as.numeric(is.na(dat[, "education"])) 1 2 3 4 5 [1] 1 1 1 1 1 0 1 0 0 1 0 1 1 1 0 1 1 1 1 1 1 1 0 0 0 0 0 1 0 0 1 0 [33] 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [65] 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 [97] 0 0 0 0 [ reached getOption("max.print") -- omitted 2700 entries ] そのため，is.na() の出力は何も考えずにそのまま mean() に突っ込むことで，TRUE の割合（無回答率）がカウントできます。同様に，sum() をとれば TRUE の人数（この場合欠測の人数）をカウントすることも可能です。１項目の欠測（無回答）率を出す 1 mean(is.na(dat[, "education"])) 1 [1] 0.07964286

5.

3.1 欠測値の処理１項目の欠測人数を出す 1 sum(is.na(dat[, "education"])) 1 [1] 223 ということで，あとはこれを全項目に対して適用していけば良いのですが，下のように全部書いていてはキリがありません。今回は ID を除いて 28 項目なのでまだなんとかなりますが，数百項目になったらどうしましょう。変数ごとに愚直に書くなら 1 2 3 4 5 6 7 8 9 mean(is.na(dat[, "Q1_1"])) mean(is.na(dat[, "Q1_2"])) mean(is.na(dat[, "Q1_3"])) mean(is.na(dat[, "Q1_4"])) ############# (中略) ############## mean(is.na(dat[, "Q1_25"])) mean(is.na(dat[, "age"])) mean(is.na(dat[, "gender"])) mean(is.na(dat[, "education"])) このような場合に備えて「一気に複数の列（または行）に同じ処理を施す方法」を紹介します。 apply() 関数は，データフレーム・行列・多次元配列について，指定した次元（行または列）にそれぞれ同じ関数を適用するという関数です。…よくわからないと思うので実際にやってみましょう。まずは以下のコードと出力を見てみてください。 apply() のキホン 1 apply(dat, 2, head) 1 2 3 4 5 6 7 8 9 10 11 12 13 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 2 4 3 4 4 2 3 3 4 4 3 3 [2,] 2 2 4 5 2 5 5 4 4 3 4 1 1 [3,] 3 5 4 5 4 4 4 5 4 2 5 2 4 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 6 3 [2,] 6 4 3 3 3 3 5 5 4 2 4 [3,] 4 4 5 4 5 4 2 3 4 2 5 Q1_24 Q1_25 gender education age [1,] 4 3 1 NA 16 [2,] 3 3 2 NA 18 [3,] 5 2 2 NA 17 [ reached getOption("max.print") -- omitted 3 rows ] 5

6.

3.1 欠測値の処理 6 head(dat) のときと同じような結果になりました。この時，裏では図 3.1 のような処理が行われています。つまり head(dat[,1])，head(dat[,2])，…をという感じで列ごとに同じ関数を適用し，その結果をまとめているのです。図 3.1: apply(dat, 2, head) の挙動 apply() は基本的に 3 つの引数からなります。apply(X, MARGIN, FUN) 関数を適用する元データです。 X MARGIN どの次元に関して，それぞれ関数を適用するかを指します。数字はデータフレームにおける要素の指定 [ , ] のときの順番と同じです。したがって MARGIN = 1 の場合は行ごと，MARGIN = 2 の場合は列ごとに関数を適用します。 FUN どの関数を適用するかです。もしその関数が別の引数を取る場合は，この後に引数を続けて書いていきます。では，当初の目的である「 dat の列ごとに NA の割合を出したい」場合にはどうしたら良いでしょうか。引数を一つずつあてはめていくと，X には is.na(dat) を，MARGIN には 2(列ごと) を，FUN には mean が入れば良さそうです。

#p6

7.

3.1 欠測値の処理まずはデータフレームに is.na() 1 is.na(dat) 1 2 3 4 5 6 7 8 9 10 11 12 13 ID Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE Q1_10 Q1_11 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE Q1_20 Q1_21 Q1_22 Q1_23 Q1_24 Q1_25 gender education age [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE [ reached getOption("max.print") -- omitted 2797 rows ] 各項目（列）の無回答（NA）率を確認 1 apply(is.na(dat), 2, mean) 1 2 3 4 5 6 7 8 9 10 11 12 ID Q1_1 Q1_2 Q1_3 Q1_4 0.000000000 0.005714286 0.009642857 0.009285714 0.006785714 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 0.005714286 0.007500000 0.008571429 0.007142857 0.009285714 Q1_10 Q1_11 Q1_12 Q1_13 Q1_14 0.005714286 0.008214286 0.005714286 0.008928571 0.003214286 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 0.007500000 0.007857143 0.007500000 0.003928571 0.012857143 Q1_20 Q1_21 Q1_22 Q1_23 Q1_24 0.010357143 0.007857143 0.000000000 0.010000000 0.005000000 Q1_25 gender education age 0.007142857 0.000000000 0.079642857 0.000000000 education の無回答率が 0.079642857，つまりおよそ 7.96% とやや高いですが，25 項目中には無回答率の高い項目はなかったので，予定通り進めていきます。欠測値があるレコードを取り除く関数としては na.omit() というものがあります。ですがこの関数は一つでも欠測値がある行は全て除外するため，na.omit(dat) としてしまうと education が NA の人（つまり最終学歴を答えることに抵抗がある層）もごっそりと抜け落ちてしまいます。今は 25 項目の中に欠測がある人だけ取り除きたいので，ここでは subset() に適切な条件式を与えることで目的を達成していきましょう。材料は以下の関数たち＋ is.na() ＋ apply() です。 7

8.

3.1 欠測値の処理 8 all(): 全て TRUE であるかどうかを判断する 1 2 # 与えたものが全てTRUEであればTRUEを返します。 all(c(TRUE, TRUE, TRUE)) 1 [1] TRUE 1 2 # 一つでもFALSEがあればFALSEを返します。 all(c(TRUE, FALSE, TRUE)) 1 [1] FALSE paste0(): 与えられた character 型を結合する 1 paste0("This ", "is ", "a pen.") 1 [1] "This is a pen." 1 2 3 # ベクトルを与えた場合，ベクトルの各要素に対して結合を行います。 # 大量の列名の指定に使えます。 paste0("No.", 1:10) 1 2 [1] "No.1" [9] "No.9" "No.2" "No.3" "No.10" "No.4" "No.5" "No.6" "No.7" "No.8" これらを組み合わせると，以下のようにして目的を達成できます。最後の結果を dat に代入するのを忘れないように気をつけてください。部分的に欠測が無いデータを残す手順 1 1 2 3 # まず"Q1_1", "Q1_2", ...という文字列のベクトルを作る cols <- paste0("Q1_", 1:25) cols # 結果の表示 1 2 3 4 [1] "Q1_1" "Q1_2" "Q1_3" "Q1_4" "Q1_5" "Q1_6" "Q1_7" "Q1_8" [9] "Q1_9" "Q1_10" "Q1_11" "Q1_12" "Q1_13" "Q1_14" "Q1_15" "Q1_16" [17] "Q1_17" "Q1_18" "Q1_19" "Q1_20" "Q1_21" "Q1_22" "Q1_23" "Q1_24" [25] "Q1_25"

9.

[beta]

3.1 欠測値の処理

9

部分的に欠測が無いデータを残す手順 2
1
2
3
4

# 25列について，各要素がNAかどうかを判定した（is.na）データフレームを作る
# !を使っているので，NAならばFALSE，値が入っていればTRUEになっている点に注意
dat_not_na <- !is.na(dat[, cols])
dat_not_na # 結果の表示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
Q1_11 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20
[1,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[2,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[3,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[4,] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
Q1_21 Q1_22 Q1_23 Q1_24 Q1_25
[1,] TRUE TRUE TRUE TRUE TRUE
[2,] TRUE TRUE TRUE TRUE TRUE
[3,] TRUE TRUE TRUE TRUE TRUE
[4,] TRUE TRUE TRUE TRUE TRUE
[ reached getOption("max.print") -- omitted 2796 rows ]
[1,]
[2,]
[3,]
[4,]

部分的に欠測が無いデータを残す手順 3
1
2
3
4
5

# 25列について，各要素がNAかどうかを判定した（is.na）データフレームを作る
# !を使っているので，NAならばFALSE，値が入っていればTRUEになっている点に注意
# 行方向にapplyを適用，一つでもNAがある場合はFALSEを返している
no_missing <- apply(dat_not_na, 1, all)
no_missing # 結果の表示

1
2
3
4
5
6

[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE
[12] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[23] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[34] TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE
[45] TRUE TRUE TRUE TRUE TRUE TRUE
[ reached getOption("max.print") -- omitted 2750 entries ]

TRUE
TRUE
TRUE
TRUE

部分的に欠測が無いデータを残す手順 4
1
2

# no_missingがTRUEの人だけをsubsetしている
dat <- subset(dat, no_missing)

もちろん，
上記の処理を全て 1 行にまとめて dat <- subset(dat,apply(!is.na(dat[,paste0("Q1_",1:25)]),1,all))

10.

3.1 欠測値の処理なんて書いても OK です。でも読みにくいですねェ*5 。きちんと NA がなくなっていることを確認するために，再度各列の NA の割合（または個数）をカウントしてみましょう。Q1_1 から Q1_25 が全て 0 になっていれば OK です。改めて無回答（NA）率を確認 1 apply(is.na(dat), 2, mean) 1 2 3 4 5 6 7 8 9 10 ID Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 Q1_23 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 Q1_24 Q1_25 gender education age 0.00000000 0.00000000 0.00000000 0.08210181 0.00000000 最後に，欠測値を除外した後のデータフレームのサイズを確認しておきましょう。上で紹介した str() でも行数は出してくれるのですが，ここではまた別の方法で確認してみます。 dim(): データフレームのサイズを確認 1 dim(dat) 1 [1] 2436 29 nrow(): データフレームの行数を確認 1 2 # 列数を見たい場合はncol() nrow(dat) 1 [1] 2436 nrow() 関数は，例えば「データフレームの各行に順番に何か処理をしたい」という場合に， 1:nrow(dat) という形で「行数だけの長さの整数ベクトル」を作ったりする場合に有用です。この方法ならば，データのサイズが変わったとしてもコードを修正する必要が無い点がおすすめポイントです。また，ncol() は，一番右の列を取り出したい時なんかにも使えます。というよう *5 カッコがたくさんあると読みにくくなってしまいます。だからといって，途中のオブジェクトをいちいち変数に代入する場合，変数の名前を（ムダに）考える必要が出てきてしまいます。こういった問題をクリアするのがパイプです。R ネイティブでも最近ようやく実装されました。例えば head(dat) は dat |> head() という書き方が出来ます。パイプを使えば，処理の順番がわかりやすくなるかもしれません。 10

#p10

11.

3.2 適当な回答を除外する 11 に，nrow() および ncol() はところどころ使うことがあると思うので覚えておきましょう。 3.2 適当な回答を除外する調査においては，全ての回答者が必ず誠実に回答してくれるとは限りません。特にウェブ調査では，驚くほど不真面目な回答が見られるものです。人間は基本的にラクをしたい生き物なので，一部の人は調査に回答する際に「なるべく認知的なリソースを使わずに回答してしまおう」と考えることがあり，これを努力の最小限化（satisficing: Krosnick, 1991）と呼びます。調査を行う人達は，長い間この「努力の最小限化」に対抗するために様々な方法を考えてきました（e.g., Curran, 2016）。ここでは，明らかな適当回答を除外するためのいくつかの方法を紹介し，最後に実際に除外を行ってみます。 3.2.1 指示項目質問紙を作る段階の話ですが，「この項目は一番右を選んでください」など，尺度の内容とは無関係に回答を指示する項目や，「私はここ 1 ヶ月寝ていない」など回答が絶対に一つに決まるはずの項目を入れておく，というのが最も簡便に利用可能な検出方法だと思います。ただ，1 問だけの場合，たまたま適当に指示項目を通過してしまう可能性があるため，鬱陶しくならない程度に複数個入れておくと良いかもしれません。また，Web 調査で項目をランダマイズする場合には，指示項目が一番上に来ないように，指示項目だけ表示順を固定するなど気をつけると良さそうです。指示項目を使用した場合には，その項目について subset() を使えば良いので除外も簡単ですね。除外する際には，論文に記載するために何％もしくは何人がそこで脱落したかを数えておきましょう。R でカウントする際には，以下のような方法が考えられます。特定の条件を満たさない人の数を数える 1 2 3 4 5 6 # 仮にQ1_1が指示項目で # 「この項目では，右から2番め(5)を選んでください」 # と指示していたとしたら 1 [1] 2244 1 2 # もっとシンプルな方法 sum(dat[, "Q1_1"] != 5) 1 [1] 2244 # subsetによる方法 nrow(subset(dat, Q1_1 != 5))

#p30

12.

3.2 適当な回答を除外する 12 また，IMC（Instructional manipulation check, Oppenheimer et al., 2009）と呼ばれる方法も同じような効果を持ちます。この方法では，項目ではなく教示文に特別な指示を入れます。図 3.2 では，項目自体は「よくやる運動」を尋ねているのですが，IMC では「タイトル（Sport Participation）をクリックしてね」という指示を与えています。その結果，元論文では 46% もの人が引っかかってしまい，スポーツ名や下の”Continue” をクリックしてしまったそうです。 IMC は心理尺度の中に混ぜると言うよりは，例えば教示として状況や特定の物事に関する説明文をおく必要がある場合に使えそうな方法と言えます。長い説明文の最後か途中に「この内容を読んで正しく理解できた方は，この後の『正しく理解できた』の項目で『いいえ』を選択してください」や「次の段落の最初の文字をマルで囲んでください」のような指示を入れておく，といったやり方が考えられます。図 3.2: IMC の例ただ，IMC は場合によっては半数以上，多い時で 8 割以上も引っかかってしまうことがある（三浦・小林，2016）ので，最終的なサンプルサイズが小さくなりすぎないように気をつける必要があるかもしれません。裏を返せば，それだけウェブ調査はまともに回答されていない，ということなのかもしれませんが…。

#p30
#p12

13.

3.2 適当な回答を除外する 13 3.2.2 回答時間 Web 調査の場合，回答時間のデータも取得可能です*6 。回答時間を用いて適当回答を除外する場合，基本的には閾値を設けて，早すぎる・遅すぎる回答者を除外するという方法がとられます（e.g., Huang et al., 2012）。一般的に回答時間を用いる場合，一項目ごとの時間を収集することが望ましい (OIOS: One-item-one-screen) とされています。これは，例えば後半で疲れてしまいその後の回答が適当になってしまった場合を検出できるなど，回答行動内での動きをきめ細かくチェックできるようになるためです。実際にどの程度の閾値を用いるかはケースバイケースです。単純に設問文が長い項目だったり，内容的に複雑な思考を必要とするような項目では所要時間は必然的に長くなるでしょう。一部の研究では，内容的に絶対にまともに回答することができないであろう時間（3 秒など）を閾値とする考え方（Kong et al., 2007）なども提案されていますが，多くの先行研究ではデータから閾値を動的に決めています。Höhne & Schlosser（2018）では，先行研究で用いられた閾値のいくつかを表 3.1 のようにまとめてくれています。表 3.1: 回答時間の閾値（Höhne & Schlosser, 2018, Table 1）下側の閾値上限の閾値平均値 −2 × SD 平均値 +2 × SD 中央値 −1.5× 四分位範囲中央値 +1.5× 四分位範囲中央値 −1.5×（中央値 − 第 1 四分位点）中央値 +1.5×（第 3 四分位点 − 中央値中央値 −3×（中央値 − 第 1 四分位点）中央値 +3×（第 3 四分位点 − 中央値）下側 1% 上側 1% ということで，回答時間をもとに適当回答を検出するためには，分位点を計算できるようになっておくと良さそうです。R ではお好きな分位点を以下のように求めることが出来ます。 quantile(): 分位点を求める 1 2 3 1 2 # データに回答時間が無いので，今回はageの分位点を求めます # 普通に実行すると，最大値・最小値と四分位数が出てくる quantile(dat[, "age"]) 0% 3 25% 20 50% 25 75% 100% 35 86 *6 最も手軽に Web 調査フォームを作成する場合は Google Form が選択肢に上がるかと思いますが，どうやら Google Form には回答時間を記録する方法はなさそうです…同じ Web 調査作成プラットフォームの中では， Questant や qualtrics では，項目ごとではなく回答全体での所要時間の記録が可能なようです（項目単位で記録できるかは不明）。あるいは javascript が使えるならば jsPsych といったフレームワークもあります。こちらは心理学系の実験が色々出来ますが，その一環として普通のアンケートも実装可能であり，javascript によって記録できる行動は（画面遷移やマウスカーソルの座標なども）理論上全て取得可能です。

14.

3.2 適当な回答を除外する 1 2 # 引数probsに0から1の間の値を指定すると，対応する分位点が求められる quantile(dat[, "age"], probs = c(0.01, 0.99)) 1 2 1% 99% 12.00 59.65 実際にデータを除外する場合には，1 項目でも基準に引っかかったらアウト，とはしないことが多いです。複数の項目の回答時間（引っかかった項目数）に加えて，この後紹介するような方法によって実際の回答内容を見ながら総合的に判断していきます。そして閾値の決め方に正解はありません。「疑わしきは罰せず」の姿勢でいった場合には適当回答の影響で回答の分布が歪むかもしれない一方で，「疑わしきは罰する」の姿勢の場合には一部の真面目な回答も除外してしまうために回答の分布が歪む可能性があります。余裕があれば，いろいろな閾値のもとで同じ分析を行って結果が変わらないことを確認（感度分析）すると良いでしょう。 3.2.3 回答内容の不一致項目の組み合わせや回答方法次第では，回答内容の不一致によって疑わしい人を検出できるかもしれません。例えば今回使用しているデータには age と education という項目があります。海外ならば飛び級制度があるので一概にルールを決めるのは難しいですが，日本であれば 1 年だけ「飛び入学」または「早期卒業」が可能なようです。ということは，「17 歳未満の高卒以上」や「20 歳未満の大卒」は日本のルール的には不可能なため，そのような回答をしている人は怪しいかもしれません。ただし，これは少なくとも年齢または学歴の項目について適当に回答または記入ミスがあるというだけであり，心理尺度の部分は正しく回答している可能性も十分にあります。そのため，上記のような「1 つのミス」によって該当者の全ての回答を除外するかは要検討です。例えば，「学歴によって因子構造が変わるか」を検証したい場合であれば，学歴が正しく記入されていない可能性のあるデータは使用できないでしょう。というように，目的に応じて使うかどうかを決めるのが良さそうです。他にも，成人の ADHD を診断するための尺度である CAARS（Conners’ Adult ADHD Rating Scales: Conners et al., 1999）には矛盾指標というものが用意されています。これは，例えば「1 つの場所に長時間いることが難しい」と「じっと座っているときでも，内心は落ち着かない」のように，普通に考えたらだいたい同じ回答が得られそうな項目ペアについて差得点を出すことで，誠実に回答しているかを判断する方法です。特定の構成概念を測定するための一般的な心理尺度では，CAARS の矛盾指標のように確立された項目対が用意されているわけではないですが，通常全ての項目が同じような内容を尋ねているわけですから，同じような考え方で，あまりにも回答が不安定な人は怪しい，と判断できるかもしれません。ということでここからは，不適切回答を検出するためのいくつかの手法を提供している careless パッケージを用いたやり方を紹介していきます。 14

#p30

15.

3.2 適当な回答を除外する 1 2 # install.packages("careless") library(careless) 同じ構成概念を測定しようとしている心理尺度の場合，内容をきちんと読んで回答しているならば，基本的には各項目に対して同じような選択肢（例えば全て「あてはまる」寄り）を選ぶはずです。言い換えると，そのような場合には回答の個人内標準偏差が小さくなっていると考えられます。…という考え方に基づいているのが Intra-individual Response Variability（IRV: Dunn et al., 2018）です。大層な名前がついていますが，やっていることは単純な標準偏差です。 careless パッケージには irv() という関数が用意されています。 irv(): 尺度内標準偏差を求める 1 2 # 最初の5項目（Q1_1-Q1_5）は同じ特性を測定しているので，ばらつきは小さいはず irv(dat[, paste0("Q1_", 1:5)]) 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 0.8944272 1.5165751 0.5477226 0.8366600 1.1401754 0.5477226 1.4142136 8 10 11 13 14 15 16 1.7888544 1.6431677 0.8366600 0.7071068 0.5477226 1.6431677 1.5165751 17 18 19 20 21 22 23 1.6733201 0.4472136 0.7071068 0.8366600 1.6431677 2.5884358 2.0736441 24 25 26 27 28 29 30 1.6431677 0.7071068 2.7386128 0.8944272 1.7320508 1.7888544 0.7071068 31 32 2.1679483 1.7320508 [ reached getOption("max.print") -- omitted 2406 entries ] irv() の基本的な挙動は，与えられたデータに対して行ごとに標準偏差を取っているだけなので，先程のコードを今までに登場した関数で表せば apply(dat[,paste0("Q1_",1:5)], 1, sd) と同じことです。一応 irv() 関数は他にも，引数次第で「前半の IRV」「後半の IRV」というようにデータを等分割してそれぞれのパートでの IRV を出してくれる機能もあります。調査で言えば，途中で疲れて後半が適当になっている人を検出したり，dat の 25 項目のように複数の異なる特性に関する項目が並んでいる場合には役に立つかもしれません。分割してから尺度内標準偏差を求める 1 2 3 # split = TRUEにしてnum.splitで分割数を決める # ただし等分割しかしてくれない irv(dat[, paste0("Q1_", 1:25)], split = TRUE, num.split = 5) 1 2 3 irvTotal irv1 irv2 irv3 irv4 irv5 1 0.900000 0.8944272 0.8366600 0.5477226 0.8366600 1.3038405 2 1.294862 1.5165751 0.7071068 2.1213203 1.0954451 0.8366600 15

#p30

16.

3.2 適当な回答を除外する 4 5 6 7 8 9 10 3 1.092398 0.5477226 1.2247449 1.0954451 1.1401754 1.5165751 4 1.166190 0.8366600 0.8366600 0.7071068 1.6431677 0.8944272 5 1.000000 1.1401754 1.1401754 1.5165751 0.8366600 0.4472136 6 1.863688 0.5477226 2.3021729 2.3452079 1.2247449 1.9235384 7 1.607275 1.4142136 1.1401754 0.8366600 0.5477226 2.1679483 8 1.519868 1.7888544 1.0000000 1.9235384 1.7888544 1.1401754 [ reached 'max' / getOption("max.print") -- omitted 2428 rows ] 3.2.4 ストレートライン上で説明したように，心理尺度では，内容をきちんと読んで回答しているならば，基本的には同じ選択肢を選ぶはずです。とはいえ，完全に同じ選択肢を選び続けると言うよりは，細かい内容の違い等によって，少しずつ異なる回答をすることが期待されています*7 。リッカート尺度は表形式で尋ねることが多く，そのようなフォーマットに対する努力の最小限化は，縦一列に同じ選択肢を選び続けるという形で表れると考えられます。特に逆転項目（「～ではない。」のような項目）があるにも関わらずそれを無視して同じ選択肢を選び続けている場合は要注意です。このように「ずっと同じ選択肢を選び続けている」状態をストレートラインや long string などと呼びます。一般的には，ずっと同じ選択肢を選び続けた最大の長さ (long string index) を計算して，これが閾値より大きい人は要注意，という感じで利用します。long string index を R で計算するのは結構骨が折れる作業なのですが，careless パッケージには longstring() という関数が用意されています。ありがたい。 longstring(): 最長の同一回答数を計算する 1 2 3 # ストレートラインは内容の違いに関わらず続くはず # ということで，25項目全部あたえてやる longstring(dat[, paste0("Q1_", 1:25)]) 1 2 3 4 5 [1] 3 4 3 3 3 3 2 1 5 2 4 3 3 3 5 4 3 3 2 2 6 3 3 4 3 3 2 4 3 3 3 3 [33] 3 5 2 7 2 3 4 3 4 4 2 2 4 4 3 4 3 3 5 4 4 4 7 3 2 2 2 4 3 3 3 3 [65] 3 3 3 3 2 4 3 3 4 3 2 2 3 3 3 4 3 3 5 3 3 3 2 3 4 3 3 3 2 3 3 2 [97] 2 3 3 3 [ reached getOption("max.print") -- omitted 2336 entries ] 実際に IRV や long string index を除外基準として利用する場合の閾値は，やはりケースバイケースです。一つの方法としては，ヒストグラムを作成して視覚的に決めるのが良いでしょう。 *7 本当に全ての項目に同じ回答を期待しているのであれば，複数項目用意する必要が無いですからね。複数項目用意しているのは，１項目では構成概念を十分に測定しきれないためです。 16

#p16

17.

3.2 適当な回答を除外する 17 long string index のヒストグラムを見る 1 2 3 lsi <- longstring(dat[, paste0("Q1_", 1:25)]) # hist()関数でヒストグラムを作れる hist(lsi) 600 400 0 200 Frequency 800 1000 Histogram of lsi 5 10 15 20 25 lsi 図 3.3: long string index のヒストグラム図 3.3 を見ると，大体 7 か 8 くらいより大きい人は怪しそうな気がします。ヒストグラムの他にも，quantile() 関数を使って上位何 % の人が怪しい，といった検出方法も考えられますが，いずれにしても long string index（や irv）の値だけで機械的に除外するのではなく，可能な限り回答の内容も踏まえて除外するかを決めるようにしましょう。ちなみに，今回のデータでは項目の出題順は個人ごとにランダムであったと思われます。このような場合，ストレートラインによる除外は意味を成さないので気をつけてください*8 。また，ストレートラインと似た適当回答としては，複数個でパターンを繰り返すもの (e.g., 123451234512345…) や，階段状に回答するもの (e.g., 2345432345432…) などもあります。本当はこれらも全部検出してあげるのがベストなのですが，残念ながら longstring() 関数では対応できないので，自分で関数を作成する必要があります。これはめちゃめちゃ面倒なので今回は省略します*9 。 *8 ランダムに出題されていた場合には，左から実際の出題順に回答を並べたデータを用意できれば同様にストレートラインの計算が可能です。 *9 気になる人は個人的に教えます。ただ相当大変なので覚悟してください。

#p17

18.

3.2 適当な回答を除外する LIGHTBULB longstring() 関数を使わずにストレートラインを計算する以下では，longstring() 関数を使わずに最長ストレートラインを計算するコードの実装例を紹介します。一応各ステップでの状態も載せておくので，各行で何が行われているか確かめながら読み進めてください。また，やり方は他にも色々あるので，もしよければ自分で考えて実装してみてください。 1 2 3 1 2 3 4 5 6 # 以下のデータ（3人目のQ1_1からQ1_25）について計算してみる vec <- dat[3, 2:26] # 面倒なので列番号で指定してしまいます vec Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 5 4 5 4 4 4 5 4 2 5 2 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 3 4 4 4 5 4 5 4 2 3 Q1_21 Q1_22 Q1_23 Q1_24 Q1_25 3 4 2 5 5 2 3 1 2 is_same <- vec[, 2:25] == vec[, 1:24] # 前の項目と同じならTRUE is_same # 一旦表示してみる 1 2 3 4 5 6 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 3 FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 3 FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE Q1_21 Q1_22 Q1_23 Q1_24 Q1_25 3 FALSE FALSE FALSE TRUE FALSE 1 2 is_same <- as.numeric(is_same) # 数字に変換 is_same # 一旦表示してみる 1 [1] 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 1 2 is_same <- paste(is_same, collapse = "") # 全部くっつける is_same # 一旦表示してみる 1 [1] "000110000001100000000010" 18

19.

3.2 適当な回答を除外する 19 1 2 3 # 連続する1の最長を数えて1を足せば良いので straight_line <- strsplit(is_same, split = "0") # 文字列を0の位置で切断 straight_line # 一旦表示してみる 1 2 3 [[1]] [1] "" [12] "" 1 2 3 4 5 # 文字数に変換 # strsplit()はリストを返すため，リストの1番目を選択 # リストについてはそのうち説明します straight_line <- nchar(straight_line[[1]]) straight_line # 一旦表示してみる 1 [1] 0 0 0 2 0 0 0 0 0 2 0 0 0 0 0 0 0 0 1 1 2 # 最大値に1を足せば最長が分かる max(straight_line) + 1 1 [1] 3 "" "" "" "" "11" "" "" "" "" "" "" "" "" "1" "" "11" "" これでようやく 1 人分の最長ストレートラインが計算できました。本当はこれを（ apply() などを使って）全員に対して計算した上で，一定数より大きい人を除外するか検討する，といった手続きを取る必要があるわけです。 3.2.5 人と違う回答心理尺度の場合，同じ構成概念を尋ねている各項目はそれなりに相関しているはずです。もちろん項目によってそもそも当てはまりにくい・当てはまりやすい項目はあるのですが，総合すると項目への回答には「平均的なパターン」と言えるものがありそうです。「人と違う」を検出する方法は，異常検知の知見が使えたりします。したがって機械学習を利用して検出する方法なども研究されていたりしますが，ここでは伝統的なマハラノビス距離を用いた検出方法を紹介＆実践してみたいと思います。図 3.4 は，２項目の連続的な項目（６段階評価ではなく，スライダーみたいなもので回答したと想像してください）への回答をプロットした仮想的なものです。赤い点がそれぞれの軸の平均値を表しています。青い点 (1) と (2) は，それぞれ赤い点からの直線距離はほぼ同じなのですが， 2 変数の相関を考えると (2) の方はさほどおかしな値ではなさそうです。正の相関があるため，平均値から離れていたとしても 2 項目とも一貫して高い or 低い値ならばありえる，ということ

#p20

20.

3.2 適当な回答を除外する 20 ですね。 4 (1) (2) 3 2 X2 1 0 -1 -2 -3 -4 -4 -3 -2 -1 0 1 2 3 4 5 X1 図 3.4: 仮想的な 2 項目への回答このような時に使えるのがマハラノビス距離です。変数の相関関係を考慮した上でベクトル間の距離を出してくれます。回答者 𝑖 の回答ベクトルを 𝐱𝑖 ，全回答者の平均ベクトルを 𝐱̄ ，変数の分散共分散行列を 𝚺𝐱 とすると，マハラノビス距離は 𝐷 = √(𝐱𝑖 − 𝐱)𝚺 ̄ −1 ̄ ⊤ 𝐱 (𝐱𝑖 − 𝐱) (3.1) で計算することが出来ます。マハラノビス距離は，多変量正規分布と密接な関係を持っています。というのも，𝑘 個の変数に関する多変量正規分布の確率密度関数は 𝑓(𝐱; 𝝁, 𝚺) = 1 √(2𝜋)𝑘 |𝚺| 1 ⊤ exp (− (𝐱𝑖 − 𝝁)𝚺−1 𝐱 (𝐱𝑖 − 𝝁) ) 2 (3.2) と表されますが，よく見るとこの式の中には (3.1) 式のルートの中身によく似たものが入っています。実際に，データから多変量正規分布の平均ベクトル 𝝁 と分散共分散行列 𝚺 を最尤推定する場合，最尤推定量はそれぞれ標本平均ベクトル 𝐱̄ と標本共分散行列 𝚺𝐱 となることが知られています。したがって，マハラノビス距離はデータから推定された多変量正規分布の確率密度に反比例するということがわかります。言ってしまえば，マハラノビス距離は個人レベルでみたときの尤度みたいなもの，というわけですね。図 3.4 で言えば，塗りつぶされている楕円が，データ

#p20

21.

3.2 適当な回答を除外するから推定されたニ変量正規分布における大まかな確率密度を表していますが，マハラノビス距離はこの確率密度に応じて計算されているわけです。そう考えると，確かに (2) の点のほうが (1) の点よりも中心からの距離は短い，と判断できそうです。 R にはマハラノビス距離（の二乗）を計算するための関数 mahalanobis() が用意されているため，これを使って図 3.4 の 2 つの青い点のマハラノビス距離を計算してみます。mahalanobis() では，基本的に 3 つの引数をとります。第一引数 x には x𝑖 を，第二引数 center*10 には x̄ を，第三引数 cov には 𝚺𝐱 に相当するものを入れてあげましょう。マハラノビス距離の計算 1 2 3 4 5 # ここのコードは見ているだけでOKです。 mean_vec <- c(0, 0) # 平均ベクトル cov_matrix <- matrix(c(1, 0.9, 0.9, 1), nrow = 2) # 分散共分散行列 # 各変数の分散が１，相関係数が0.9という想定です。 cov_matrix # 一応表示してみる 1 2 3 [,1] [,2] [1,] 1.0 0.9 [2,] 0.9 1.0 1 2 3 vec1 <- c(-3, 3) # 青い点(1)の座標 # 青い点(1)のマハラノビス距離 mahalanobis(vec1, center = mean_vec, cov = cov_matrix) 1 [1] 180 1 2 3 vec2 <- c(3, 3) # 青い点(2)の座標 # 青い点(2)のマハラノビス距離 mahalanobis(vec2, center = mean_vec, cov = cov_matrix) 1 [1] 9.473684 ということで，青い点 (1)，(2) と平均ベクトルとのマハラノビス距離（の二乗）はそれぞれ 180，9.47 となり，たしかに (1) のほうが大きな外れ値であるということがわかりました。ありがたいことに careless() パッケージには，データフレームからマハラノビス距離を計算するための関数も用意されています。先程使用した mahalanobis() 関数では，自分でデータから平均ベクトルと分散共分散行列を計算する必要があったのですが，careless パッケージにあ *10 ふつうマハラノビス距離は「2 つのベクトルの距離」なので 2 つの引数は x と y のほうがしっくり来るのですが，この関数の第二引数が center なのは，この関数が「データフレームの各行について，平均ベクトルからのマハラノビス距離を計算する」という，まさにいま行っている分析の目的にピッタリの関数だからです。 21

#p20
#p21

22.

3.2 適当な回答を除外するる mahad() 関数ではデータフレームを一つ与えるだけで各行のマハラノビス距離をまとめて計算してくれます*11 。ありがたい。 mahad(): データフレームの各行のマハラノビス距離を計算する 1 mahad(dat[, paste0("Q1_", 1:25)]) 1 2 3 4 5 6 7 8 9 [1] 13.468425 25.677291 13.490783 28.823203 18.946470 25.901522 [7] 9.260945 48.474541 12.794634 19.207672 19.386422 14.485101 [13] 34.052697 38.710674 18.843160 18.056494 26.822273 38.451815 [19] 47.752717 30.685610 24.097729 27.202401 10.453943 30.186948 [25] 20.393634 12.667231 45.024926 23.116322 26.807564 12.173940 [31] 56.951093 20.122730 26.353218 30.358047 17.898820 28.093556 [37] 14.465344 10.483335 11.724812 27.950268 7.693067 57.603758 [43] 23.731495 20.956441 42.950337 14.636132 32.755405 52.321123 [ reached getOption("max.print") -- omitted 2388 entries ] これで無事，全員分のマハラノビス距離を一気に計算できました。あとはこれを利用して怪しい人を検出するだけです。ですがこの場合でも閾値の問題が登場します。これまで通り，ヒストグラムを描いてみたり，quantile() を利用してみても良いのですが，mahad() 関数にはもう少しテクニカルな方法が実装されています。 𝑝 次元多変量正規分布に従うデータから計算されるマハラノビス距離の二乗は，理論的には自由度 𝑝 の 𝜒2 分布に従うことが知られています。これを利用すると，「自由度 𝑝 の 𝜒2 分布における上位 𝑘 パーセンタイル点より大きい人を検出しよう」という考えに至ります。マハラノビス距離に基づく検出 1 2 3 4 # flag = TRUEにすると，フラグを立ててくれる # confidenceによって，閾値のパーセンタイル点を決める m_dist <- mahad(dat[, paste0("Q1_", 1:25)], flag = TRUE, confidence = 0.99) m_dist 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 d_sq flagged 13.468425 FALSE 25.677291 FALSE 13.490783 FALSE 28.823203 FALSE 18.946470 FALSE 25.901522 FALSE 9.260945 FALSE 48.474541 TRUE 12.794634 FALSE *11 内部の計算は psych::outlier() 関数を使っています。なのでこちらを使用しても OK です。 22

#p22

23.

3.2 適当な回答を除外する 11 12 13 14 15 16 17 18 19 20 21 22 23 10 19.207672 FALSE 11 19.386422 FALSE 12 14.485101 FALSE 13 34.052697 FALSE 14 38.710674 FALSE 15 18.843160 FALSE 16 18.056494 FALSE 17 26.822273 FALSE 18 38.451815 FALSE 19 47.752717 TRUE 20 30.685610 FALSE [ reached 'max' / getOption("max.print") -- omitted 2416 rows ] ここで flagged 列が TRUE になっている人が，マハラノビス距離が閾値より大きい，すなわち異質な回答パターンを示した人だということです。ちなみに，マハラノビス距離が最大の人の回答ベクトルを見てみると，たしかにこの人は６か１かばかりで回答しており，どうやら適当に回答しているような気がします。ただ実際には日頃から極端な選択肢を選ぶ傾向がある人かもしれないので，これだけで機械的に除外するのではなく，回答ベクトルを見てストレートラインの傾向が無いかなども確認して除外の判断をしたほうが良いでしょう。マハラノビス距離が最大の人の回答を見てみる 1 2 3 4 # which.max()はベクトルの中で最大値の「位置」を返す (arg max 的な) # max()は最大値そのものを返す # which.min()もあるよ dat[which.max(m_dist$d_sq), ] 1 2 3 4 5 6 ID Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 867 867 6 6 6 6 6 6 5 1 1 1 6 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 867 6 1 6 6 1 6 1 6 1 6 6 Q1_23 Q1_24 Q1_25 gender education age 867 1 6 1 1 1 31 Exclamation-Triangle マハラノビス距離による検出法の弱点マハラノビス距離は平均ベクトルからの距離を求めています。そのため，同じストレートラインでも 1 と 6 ばかりなど極端な選択肢を選んでいる人は引っかかりやすい一方で，3 と 4（どちらでもない周辺）ばかりを選んでいる人はどうしてもマハラノビス距離は小さくなってしまいます。同様に，どちらでもない周辺で適当な回答をしている人は検出しにくいという問題点があります。これを克服するためには，マハラノビス距離に加えて IRV やストレートラインなどの指

24.

3.2 適当な回答を除外する 24 標，あるいは可能であれば回答データ以外（回答時間など）の複数の指標を組み合わせて検出してあげるのが良さそうです。また，そもそもマハラノビス距離は，多変量正規分布の確率密度に比例する値であるという点にも注意が必要です。天井・床効果程度なら多分問題ないのですが，例えば 6 件法で 1 と 6（極端な選択肢）に回答が集中するような場合には，回答の分布が二峰になっているので，そもそも正規分布を当てはめることが怪しくなってきます。そのような場合にはマハラノビス距離を使うこと自体もよく考えたほうが良さそうです。ここまで適当回答を検出するためのいくつかの方法を紹介してきました。実際のデータにおいて検出＆除外を行うときにはいろいろと考えるべきことはあるのですが，とりあえず今回は機械的に「long string index ( lsi) が 8 以上」かつ「マハラノビス距離が上位 5%*12 」の人を除外することにします。ちなみにマハラノビス距離の上位 5% は，図 3.5 の分布でいえば赤い線より右側の人であり，極端に大きい人たちだろうということが言えそうです。 400 200 0 Frequency Histogram of m_dist$d_sq 0 20 40 60 80 100 m_dist$d_sq 図 3.5: マハラノビス距離の二乗の分布では，実際に 2 つの基準をともに満たす人を除外してみましょう。ここでは，複数の条件式を結合するための論理演算子を使用します。論理演算とは，複数の TRUE/FALSE をもとに 1 つの TRUE/FALSE を返す演算のことです。何はともあれ，以下の例を見てください。 &演算子: 2 つとも TRUE のときだけ TRUE を返す (AND) 1 2 TRUE & TRUE TRUE & FALSE *12 勘違いしないでいただきたいのですが，「マハラノビス距離の上位 5% を除外するのが定番」ということではありません。例えば最初から「上位 5 ％の人にフラグを立てよう」と決め打ちにしていると，もし 100 人中 50 人が適当回答者だったとしても，そのうち 5 人しか検出できないことになります。なので今回の基準はあくまでも授業の演習としてとりあえず除外を体験してみよう，というくらいの意味合いです。

#p24

25.

[beta]

3.2 適当な回答を除外する

3

FALSE & FALSE

1
2
3

[1] TRUE
[1] FALSE
[1] FALSE

25

|演算子: いずれか一つでも TRUE であれば TRUE を返す (OR)
1
2
3

TRUE | TRUE
TRUE | FALSE
FALSE | FALSE

1
2
3

[1] TRUE
[1] TRUE
[1] FALSE

今回は，AND 演算子 ( &) を用いて「long string index ( lsi) が 8 以上」かつ「マハラノビス
距離が上位 5%」の人を探し出したいと思います。
複数条件を満たす人を見つける (1)
1
2
3

# まずはマハラノビス距離の閾値を求める
thres <- quantile(m_dist$d_sq, prob = 0.95)
thres

1
2

95%
48.87949























複数条件を満たす人を見つける (2)
1
2
3

# 除外対象となる人を判別する
is_remove <- (m_dist$d_sq >= thres) & (lsi >= 8)
is_remove

1
2
3
4
5
6

[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[45] FALSE FALSE FALSE FALSE FALSE FALSE
[ reached getOption("max.print") -- omitted 2386 entries ]

26.

3.2 適当な回答を除外する 26 1 2 # 何行目の人が除外されるかを一応チェック which(is_remove) 1 [1] 1243 1359 1775 2018 1 2 # 一応除外対象の人の回答をチェック subset(dat, is_remove) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 ID Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 1430 1430 1 1 1 1 1 1 1 1 1 1 1 1560 1560 1 6 6 1 1 5 6 6 4 6 1 2043 2043 1 1 1 1 1 1 1 1 1 1 1 2313 2313 6 1 1 6 1 6 5 6 1 1 6 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 Q1_21 1430 1 1 1 1 1 1 1 1 1 1 1560 6 6 1 6 6 6 6 6 6 6 2043 1 1 1 1 1 1 1 1 1 1 2313 6 1 1 1 1 1 1 1 1 6 Q1_22 Q1_23 Q1_24 Q1_25 gender education age 1430 1 1 1 1 1 3 23 1560 6 4 6 1 2 5 31 2043 1 1 1 1 1 3 19 2313 1 1 6 1 1 3 23 複数条件を満たす人を見つける (3) 1 2 # is_removeがFALSEの人だけ残せば良いので dat <- subset(dat, !is_remove) ということで，残ったデータは 1 nrow(dat) 1 [1] 2432 となりました。以後の分析はこの 2432 人分のデータで進めていきます。 LIGHTBULB AND/OR 演算子と NA 以前，論理演算子（ ==や >など）に NA を与えた場合には，「判断不可能」的な意味合いの NA が返ってきてしまう，という話がありました。では，AND/OR 演算子に NA が与えられた場合どうなるかというと，これはやや複雑な挙動となります。一言で言えば，NA に

27.

3.3 逆転項目の処理 27 TRUE と FALSE のどちらが入るのかによって結果が変わる場合には「判断不可能」の NA が返ってきます。あとは下の例を見て理解してください。 1 2 3 4 TRUE & NA FALSE & NA TRUE | NA FALSE | NA 1 2 3 4 [1] NA [1] FALSE [1] TRUE [1] NA 3.3 逆転項目の処理心理尺度では，内容を読まずに「慣れ」でストレートライン的回答を避けたい，といった目的から意図的に「＊＊ではない」のように否定語を使ったり，構成概念をきちんと反映した項目として，意味的に反対のことを尋ねる（例えば，社交性を測定する尺度において「休日は基本的に家から出ない」のような項目を入れておく）ことがあります。こうした項目は逆転項目と呼ばれ，分析の際には，多くのケースにおいて事前に得点をひっくり返す必要があります。psych::bfi では*13 1, 9, 10, 11, 12, 22, 25 番目の項目が逆転項目となっています。逆転項目の処理自体は何も難しいことはなく，リッカート尺度の場合は「（最大値 + 最小値）から回答の値を引く」だけで OK です。ということで，一つの列について逆転項目の処理を行うならば以下のコードで可能となります。逆転項目の処理（1 項目編） 1 dat[, "Q1_1"] # 処理前 1 2 3 4 5 [1] 2 2 5 4 2 6 2 4 2 4 5 5 4 4 4 5 4 4 5 1 1 2 4 1 2 2 2 4 1 2 1 2 [33] 5 1 1 1 1 1 1 5 2 1 5 2 1 1 1 1 3 4 1 1 1 3 1 3 2 2 4 2 1 4 2 1 [65] 4 4 1 3 2 2 5 1 2 1 1 1 1 1 4 1 1 1 1 1 2 1 2 2 3 1 4 2 3 2 1 1 [97] 1 6 2 3 [ reached getOption("max.print") -- omitted 2332 entries ] 1 2 dat[, "Q1_1"] <- 7 - dat[, "Q1_1"] dat[, "Q1_1"] # 処理後 *13 psych::bfi.keys で確認可能です。変数名の頭にマイナスがついているものが逆転項目です。

#p27

28.

[beta]

3.3 逆転項目の処理

1
2
3
4
5

28

[1] 5 5 2 3 5 1 5 3 5 3 2 2 3 3 3 2 3 3 2 6 6 5 3 6 5 5 5 3 6 5 6 5
[33] 2 6 6 6 6 6 6 2 5 6 2 5 6 6 6 6 4 3 6 6 6 4 6 4 5 5 3 5 6 3 5 6
[65] 3 3 6 4 5 5 2 6 5 6 6 6 6 6 3 6 6 6 6 6 5 6 5 5 4 6 3 5 4 5 6 6
[97] 6 1 5 4
[ reached getOption("max.print") -- omitted 2332 entries ]

今回は 7 項目だけなので一つ一つやっても良いのですが，もっと項目数が多いと大変なので

apply() の力を借りましょう。今回の場合，apply() 関数の引数 X には「 dat のうち逆転項目
を処理したい列だけ」，MARGIN は 2，FUN には「7 から引く」という関数を指定してあげたら良
いわけです。ですが困ったことに R には「7 から引く」なんていう specific な関数は存在してい
ません。こういう時には，関数を定義するという方法を使う必要があります。簡単な関数を定義
している例を挙げてみましょう。
関数を作ってみる
1
2
3
4
5
6
7

# function()関数は「関数を定義する」関数
# ()の中には新しく作る関数の引数を
# {}の中には関数の動作を記入する
tashizan <- function(x, y) {
x + y
}
tashizan(3, 2)

1

[1] 5

この関数 tashizan() はもちろん R にもともと存在しない，与えられた２つの引数 x と y の
和を返す関数です*14 。同じ要領で，「7 から与えられた値を引く」という関数を作りましょう。
「7 から引く」関数を作ってみる
1
2
3
4
5

subtract_7 <- function(x) {
7 - x
}
# 挙動の確認
head(dat$Q1_9) # 処理前

1

[1] 4 3 2 5 3 1

1

subtract_7(head(dat$Q1_9)) # 関数を適用









































*14 もちろん様々なパッケージにある関数も全て，どこかの誰かが必要にかられて function() によって定義して

くれたものなのです。

#p28

29.

3.3 逆転項目の処理 1 29 [1] 3 4 5 2 4 6 これを apply() 関数の引数 FUN に指定してあげれば，すべての列に同じ処理を適用することができるわけです。逆転項目の処理（複数項目編） 1 2 3 4 5 6 7 # 処理する列名を指定 # 上のコードでQ1_1だけ逆転処理済みなのでここでは除外 cols <- paste0("Q1_", c(9, 10, 11, 12, 22, 25)) # cols <- c("Q1_9", "Q1_10", "Q1_11","Q1_12","Q1_22","Q1_25") # ↑これでもよいが，列数が増えた時には大変 head(dat[, cols]) # 処理前 1 2 3 4 5 6 7 1 2 3 4 5 6 1 2 dat[, cols] <- apply(dat[, cols], 2, subtract_7) head(dat[, cols]) # 処理後 1 2 3 4 5 6 7 1 2 3 4 5 6 Q1_9 Q1_10 Q1_11 Q1_12 Q1_22 Q1_25 4 4 3 3 6 3 3 4 1 1 2 3 2 5 2 4 2 2 5 5 5 3 3 5 3 2 2 2 3 3 1 3 2 1 3 1 Q1_9 Q1_10 Q1_11 Q1_12 Q1_22 Q1_25 3 3 4 4 1 4 4 3 6 6 5 4 5 2 5 3 5 5 2 2 2 4 4 2 4 5 5 5 4 4 6 4 5 6 4 6 これで，ようやく全ての逆転項目の処理が完了しました*15 。次章では，データの性質を確認して，収集したデータが「本当に使えるものなのか」を見ていきます。最後に保存 1 saveRDS(dat, "chapter03.rds") *15 実を言うと，apply() なんか使わなくても 7 - dat[,cols] で複数行まとめて処理は可能です。今回は自作関数の導入を兼ねてちょっと回りくどいやり方にしました。

#p29

30.

30 参考文献参考文献 Conners, C. K., Erhardt, D., & Sparrow, E. P. (1999). Conners’ adult ADHD rating scales (CAARS): technical manual. Multi-Health Systems North Tonawanda, NY. Curran, P. G. (2016). Methods for the detection of carelessly invalid responses in survey data. Journal of Experimental Social Psychology, 66, 4–19. https://doi.org/10.1016/j.je sp.2015.07.006 Dunn, A. M., Heggestad, E. D., Shanock, L. R., & Theilgard, N. (2018). Intra-individual response variability as an indicator of insufficient effort responding: comparison to other indicators and relationships with individual differences. Journal of Business and Psychology, 33(1), 105–121. https://doi.org/10.1007/s10869-016-9479-0 Höhne, J. K., & Schlosser, S. (2018). Investigating the adequacy of response time outlier definitions in computer-based web surveys using paradata SurveyFocus. Social Science Computer Review, 36(3), 369–378. https://doi.org/10.1177/0894439317710450 Huang, J. L., Curran, P. G., Keeney, J., Poposki, E. M., & DeShon, R. P. (2012). Detecting and deterring insufficient effort responding to surveys. Journal of Business and Psychology, 27 (1), 99–114. https://doi.org/10.1007/s10869-011-9231-8 Kong, X. J., Wise, S. L., & Bhola, D. S. (2007). Setting the response time threshold parameter to differentiate solution behavior from rapid-guessing behavior. Educational and Psychological Measurement, 67 (4), 606–619. https://doi.org/10.1177/0013164406294779 Krosnick, J. A. (1991). Response strategies for coping with the cognitive demands of attitude measures in surveys. Applied Cognitive Psychology, 5(3), 213–236. https: //doi.org/10.1002/acp.2350050305 Little, R. J. A. (1988). A test of missing completely at random for multivariate data with missing values. Journal of the American Statistical Association, 83(404), 1198–1202. https://doi.org/10.1080/01621459.1988.10478722 三浦麻子・小林哲郎（2016）．オンライン調査における努力の最小限化（Satisfice）傾向の比較: IMC 違反率を指標としてメディア・情報・コミュニケーション研究，1, 27–42.

31.

31 参考文献 Oppenheimer, D. M., Meyvis, T., & Davidenko, N. (2009). Instructional manipulation checks: Detecting satisficing to increase statistical power. Journal of Experimental Social Psychology, 45(4), 867–872. https://doi.org/10.1016/j.jesp.2009.03.009 Schafer, J. L., & Graham, J. W. (2002). Missing data: Our view of the state of the art. Psychological Methods, 7 (2), 147–177. https://doi.org/10.1037/1082-989X.7.2.147 高橋将宜・渡辺美智子（2017）．欠測データ処理: R による単一代入法と多重代入法出版共立

多変量解析_03_データの前処理

Kyosuke Bunji

関連スライド

ベイズ統計_07_マルコフ連鎖モンテカルロ法(2)

ベイズ統計_02_確率の基本とベイズの定理

ベイズ統計_01_イントロダクション

ベイズ統計_03_尤度

ベイズ統計_05_基本的なベイズ推論(2)

ベイズ統計_04_事前分布・基本的なベイズ推論(1)

各ページのテキスト