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

1.

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

2.

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

3.

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

4.

3.1 欠測値の処理 1 [1] TRUE TRUE TRUE 2 [12] TRUE TRUE TRUE FALSE 3 [23] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE 4 [34] TRUE FALSE FALSE FALSE FALSE FALSE 5 [45] FALSE FALSE FALSE FALSE FALSE FALSE 6 TRUE TRUE FALSE TRUE FALSE FALSE TRUE FALSE TRUE TRUE TRUE TRUE FALSE FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE FALSE FALSE [ reached getOption("max.print") -- omitted 2750 entries ] is.na() を強制的に数値に変換してみると 1 as.numeric(is.na(dat[, "education"])) 1 [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 2 [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 3 [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 4 [97] 0 0 0 0 5 [ reached getOption("max.print") -- omitted 2700 entries ] そのため，is.na() の出力は何も考えずにそのまま mean() に突っ込むことで，TRUE の割合がカウントできます。同様に，sum() をとれば TRUE の人数（この場合欠測の人数）をカウントすることも可能です。１項目の欠測（無回答）率と人数を出す 1 mean(is.na(dat[, "education"])) 1 [1] 0.07964286 1 sum(is.na(dat[, "education"])) 1 [1] 223 ということで，これを全項目に対して適用していけば良いのですが，下のように全部書いていてはキリがありません。今回は ID を除いて 28 項目なのでまだなんとかなりますが，数百項目になったらどうしましょう。こんなふうに全部書くのかい？ 4 1 mean(is.na(dat[, "Q1_1"])) 2 mean(is.na(dat[, "Q1_2"])) 3 mean(is.na(dat[, "Q1_3"])) 4 mean(is.na(dat[, "Q1_4"])) 5 ############# (中略) ############## 6 mean(is.na(dat[, "Q1_25"])) 7 mean(is.na(dat[, "age"])) 8 mean(is.na(dat[, "gender"]))

5.

3.1 欠測値の処理 9 mean(is.na(dat[, "education"])) このような場合に備えて「一気に複数の列（または行）に同じ処理を施す方法」を紹介します。apply() 関数は，データフレーム・行列・多次元配列について，指定した次元（行または列）にそれぞれ同じ関数を適用するという関数です。…よくわからないと思うので実際にやってみましょう。まずは以下のコードと出力を見てみてください。 apply のキホン 1 apply(dat, 2, head) 1 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 2 [1,] 1 2 4 3 4 4 2 3 3 4 4 3 3 3 [2,] 2 2 4 5 2 5 5 4 4 3 4 1 1 4 [3,] 3 5 4 5 4 4 4 5 4 2 5 2 4 5 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 Q1_23 6 [1,] 3 4 4 3 4 2 2 3 3 6 3 7 [2,] 6 4 3 3 3 3 5 5 4 2 4 8 [3,] 4 4 5 4 5 4 2 3 4 2 5 9 Q1_24 Q1_25 gender education age 10 [1,] 4 3 1 NA 16 11 [2,] 3 3 2 NA 18 12 [3,] 5 2 2 NA 17 13 [ reached getOption("max.print") -- omitted 3 rows ] head(dat) のときと同じような結果になりました。この時，裏では図3.1 のような処理が行われています。つまり head(dat[,1])，head(dat[,2])，…をという感じで列ごとに同じ関数を適用し，その結果をまとめているのです。 apply() は基本的に 3 つの引数からなります。apply(X, MARGIN, FUN) 関数を適用する元データです。 X MARGIN どの次元に関して，それぞれ関数を適用するかを指します。数字は [ , ] のときの順番と同じです。したがって MARGIN = 1 の場合は行ごと，MARGIN = 2 の場合は列ごとに関数を適用します。 FUN どの関数を適用するかです。もしその関数が別の引数を取る場合は，この後に引数を続けて書いていきます。では，当初の目的である「dat の列ごとに NA の割合を出したい」場合にはどうしたら良いでしょうか。引数を一つずつあてはめていくと，X には is.na(dat) を，MARGIN には 2(列ごと) を，FUN には mean が入れば良さそうです。データフレームに is.na() 1 5 is.na(dat)

6.

3.1 欠測値の処理 ID Q1̲1 Q1̲2 … gender education age 1 2 4 2 3 4 2 5 3 4 4 4 2799 2800 2 5 2 3 … 1 NA 16 2 2 2 NA NA NA 18 17 17 1 2 4 4 31 50 ⁝ … apply(dat,2,head) head(ID) ID Q1̲1 Q1̲2 gender education age 1 2 3 4 ⁝ 5 5 2 3 ⁝ 4 4 4 4 ⁝ 1 2 2 2 ⁝ NA NA NA NA ⁝ 16 18 17 17 ⁝ 2799 2800 2 5 2 3 1 2 4 4 31 50 head(Q1̲1) … head(Q2̲2) ID Q1̲1 Q1̲2 1 2 5 5 4 4 3 2 4 ... head(gender) … head(education) head(age) gender education age 1 2 NA NA 16 18 2 NA 17 ID Q1̲1 Q1̲2 … gender education age 1 2 3 2 2 5 4 4 4 … 1 2 2 NA NA NA 16 18 17 図3.1 apply(dat, 2, head) 関数の挙動 1 ID Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 2 [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 3 [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 4 [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 5 Q1_10 Q1_11 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 6 [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 7 [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 8 [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 9 Q1_20 Q1_21 Q1_22 Q1_23 Q1_24 Q1_25 gender education age 10 [1,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE 11 [2,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE [3,] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE 12 13 [ reached getOption("max.print") -- omitted 2797 rows ] 各項目（列）の無回答（NA）率を確認 1 1 2 3 6 apply(is.na(dat), 2, mean) ID Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 0.000000000 0.005714286 0.009642857 0.009285714 0.006785714 0.005714286 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11

7.

3.1 欠測値の処理 4 5 6 7 8 9 10 0.007500000 0.008571429 0.007142857 0.009285714 0.005714286 0.008214286 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 0.005714286 0.008928571 0.003214286 0.007500000 0.007857143 0.007500000 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 Q1_23 0.003928571 0.012857143 0.010357143 0.007857143 0.000000000 0.010000000 Q1_24 Q1_25 gender education age 0.005000000 0.007142857 0.000000000 0.079642857 0.000000000 education の無回答率が 7.96% とやや高いですが，25 項目中には無回答率の高い項目はなかったので，予定通り進めていきます。欠測値があるレコードを取り除く関数としては na.omit() というものがあります。ですがこの関数は一つでも欠測値がある行は全て除外するため，na.omit(dat) としてしまうと education が NA の人もごっそりと抜け落ちてしまいます。今は 25 項目の中に欠測がある人だけ取り除きたいので，ここでは subset() に適切な条件式を与えることで目的を達成していきましょう。材料は以下の関数たち＋ is.na() ＋ apply() です。 all() 関数 1 # 与えたものが全てTRUEであればTRUE，一つでもFALSEがあればFALSEを返します。 2 all(c(TRUE, TRUE, TRUE)) 1 [1] TRUE 1 all(c(TRUE, FALSE, TRUE)) 1 [1] FALSE paste0() 関数 1 # 与えた文字列を結合した文字列を返します。 2 paste0("This ", "is ", "a pen.") 1 [1] "This is a pen." 1 # ベクトルを与えた場合，ベクトルの各要素に対して結合を行います。 2 # 大量の列名の指定に使えます。 3 paste0("No.", 1:10) 1 [1] "No.1" "No.2" 2 [9] "No.9" "No.10" "No.3" "No.4" "No.5" "No.6" "No.7" "No.8" これらを組み合わせると，以下のようにして目的を達成できます。最後の結果を dat に代入するのを忘れないように気をつけてください。 7

8.

[beta]

3.1 欠測値の処理
部分的に欠測が無いデータを残す
1

# まず"Q1_1", "Q1_2", ...という文字列のベクトルを作る

2

cols <- paste0("Q1_", 1:25)

3

cols

1

[1] "Q1_1"

"Q1_2"

2

[9] "Q1_9"

"Q1_10" "Q1_11" "Q1_12" "Q1_13" "Q1_14" "Q1_15" "Q1_16"

"Q1_4"

"Q1_5"

"Q1_6"

"Q1_7"

"Q1_8"

3

[17] "Q1_17" "Q1_18" "Q1_19" "Q1_20" "Q1_21" "Q1_22" "Q1_23" "Q1_24"

4

[25] "Q1_25"

1

# 25列について，各要素がNAかどうかを判定した（is.na）データフレームを作る

2

# !を使っているので，NAならばFALSE，値が入っていればTRUEになっている点に注意

3

dat_not_na <- !is.na(dat[, cols])

4

dat_not_na

1

Q1_1

Q1_2

Q1_3

Q1_4

Q1_5

Q1_6

Q1_7

Q1_8

Q1_9 Q1_10

2

[1,]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

3

[2,]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

4

[3,]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

5

[4,]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

6

Q1_11 Q1_12 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20

7

[1,]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

8

[2,]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

9

[3,]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

10

[4,]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

11

Q1_21 Q1_22 Q1_23 Q1_24 Q1_25

12

[1,]

TRUE

TRUE

TRUE

TRUE

TRUE

13

[2,]

TRUE

TRUE

TRUE

TRUE

TRUE

14

[3,]

TRUE

TRUE

TRUE

TRUE

TRUE

15

[4,]

TRUE

TRUE

TRUE

TRUE

TRUE

16

[ reached getOption("max.print") -- omitted 2796 rows ]

1

# 行方向にapplyを適用，一つでもNAがある場合はFALSEを返している

2

no_missing <- apply(dat_not_na, 1, all)

3

no_missing

1

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE FALSE

TRUE

TRUE

2

[12] FALSE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

3

[23]

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

4

[34]

TRUE FALSE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE FALSE

TRUE

TRUE

5

[45]

TRUE

TRUE

TRUE

TRUE

TRUE

6

8

"Q1_3"

[1]

TRUE

[ reached getOption("max.print") -- omitted 2750 entries ]

9.

[beta]

3.2 適当な回答を除外する

1

# no_missingがTRUEの人だけをsubsetしている

2

dat <- subset(dat, no_missing)

もちろん，
1 行にまとめて dat <- subset(dat,apply(!is.na(dat[,paste0("Q1_",1:25)]),1,all))
と書いても OK です。でも読みにくいですねェ*4 。
一応，欠測値を除外した後のデータフレームのサイズを確認しておきましょう。上で紹介し
た str() でも行数は出してくれるのですが，ここではまた別の方法で確認してみます。
データフレームのサイズを確認
1

dim(dat)

1

[1] 2436

29

行数だけ
1

# 列数を見たい場合はncol()

2

nrow(dat)

1

[1] 2436

ちなみに ncol() は，一番右の列を取り出したい時なんかにも使えます。というように，

nrow() および ncol() はところどころ使うことがあると思うので覚えておきましょう。
3.2 適当な回答を除外する
調査においては，全ての回答者が必ず誠実に回答してくれるとは限りません。特にウェブ調
査では，驚くほど不真面目な回答が見られるものです。人間は基本的にラクをしたい生き物
なので，一部の人は調査に回答する際に「なるべく認知的なリソースを使わずに回答してし
まおう」と考えることがあり，これを努力の最小限化 (satisficing; Krosnick, 1991) と呼び
ます。
調査を行う人達は，長い間この「努力の最小限化」に対抗するために様々な方法を考えてき
ました (e.g., Curran, 2016)。ここでは，明らかな適当回答を除外するためのいくつかの方法
を紹介し，最後に実際に除外を行ってみます。

3.2.1 指示項目
質問紙を作る段階の話ですが，
「この項目は一番右を選んでください」など，尺度の内容とは
無関係に回答を指示する項目や，「私はここ 1 ヶ月寝ていない」など回答が絶対に一つに決
まるはずの項目を入れておく，というのが最も簡便に利用可能な検出方法だと思います。た
*4

9

カッコがたくさんあると読みにくくなってしまいます。だからといって，途中のオブジェクトをいちいち変
数に代入する場合，変数の名前を（ムダに）考える必要が出てきてしまいます。こういった問題をクリアする
のがパイプです。R ネイティブでも最近ようやく実装されました。例えば head(dat) は dat |> head() と
いう書き方が出来ます。パイプを使えば，処理の順番がわかりやすくなるかもしれません。

10.

3.2 適当な回答を除外するだ，1 問だけの場合，たまたま適当に指示項目を通過してしまう可能性があるため，鬱陶しくならない程度に複数個入れておくと良いかもしれません。また，Web 調査で項目をランダマイズする場合には，指示項目が一番上に来ないように，指示項目だけ表示順を固定するなど気をつけると良さそうです。指示項目を使用した場合には，その項目について subset() を使えば良いので除外も簡単ですね。除外する際には，論文に記載するために何％もしくは何人がそこで脱落したかを数えておきましょう。R でカウントする際には，以下のような方法が考えられます。特定の条件を満たす人の数を数える 1 # 仮にQ1_1が指示項目で，右から2番め(5)を選ぶように指示していたとしたら 2 # subsetによる方法 3 nrow(subset(dat, Q1_1 != 5)) 1 [1] 2244 1 # もっとシンプルな方法 2 sum(dat[, "Q1_1"] != 5) 1 [1] 2244 また，IMC (Instructional manipulation check, Oppenheimer et al., 2009) と呼ばれる方法も同じような効果を持ちます。この方法では，項目ではなく教示文に特別な指示を入れます。図 3.2 では，項目自体は「よくやる運動」を尋ねているのですが，IMC では「タイトル（Sport Participation）をクリックしてね」という指示を与えています。その結果，元論文では 46% もの人が引っかかってしまい，スポーツ名や下の”Continue” をクリックしてしまったそうです。IMC は心理尺度の中に混ぜると言うよりは，例えば教示として状況や特定の物事に関する説明文をおく必要がある場合に使えそうな方法と言えます。長い説明文の最後か途中に「この内容を読んで正しく理解できた方は，この後の『正しく理解できた』の項目で『いいえ』を選択してください」や「次の段落の最初の文字をマルで囲んでください」のような指示を入れておく，といったやり方が考えられます。ただ，IMC は場合によっては半数以上，多い時で８割以上も引っかかってしまうことがある (三浦・小林, 2016) ので，最終的なサンプルサイズが小さくなりすぎないように気をつける必要があるかもしれません。裏を返せば，それだけウェブ調査はまともに回答されていない，ということなのかもしれませんが…。 10

11.

3.2 適当な回答を除外する図3.2 IMC の例 3.2.2 回答時間 Web 調査の場合，回答時間のデータも取得可能です*5 。回答時間を用いて適当回答を除外する場合，基本的には閾値を設けて，早すぎる・遅すぎる回答者を除外するという方法がとられます (e.g., Huang et al., 2012)。一般的に回答時間を用いる場合，一項目ごとの時間を収集することが望ましいとされています。これは，例えば後半で疲れてしまいその後の回答が適当になってしまった場合を検出できるなど，回答行動内での動きをきめ細かくチェックできるようになるためです。実際にどの程度の閾値を用いるかはケースバイケースです。単純に設問文が長い項目だったり，内容的に複雑な思考を必要とするような項目では所要時間は必然的に長くなるでしょう。一部の研究では，内容的に絶対にまともに回答することができないであろう時間（3 秒など）を閾値とする考え方 (Kong et al., 2007) も提案されていますが，多くの先行研究ではデータから閾値を動的に決めています。Höhne and Schlosser (2018) では，先行研究で用いられた閾値のいくつかを表3.1にまとめてくれています。ということで，回答時間をもとに適当回答を検出する場合には，分位点を計算する必要があります。R ではお好きな分位点を以下のように求めることが出来ます。 *5 11 最も手軽に Web 調査フォームを作成する場合は Google Form が選択肢に上がるかと思いますが，どうやら Google Form には回答時間を記録する方法はなさそうです…同じ Web 調査作成プラットフォームの qualtrics であれば，回答時間の記録も可能です。あるいは javascript が使えるならば jsPsych といったフレームワークもあります。こちらは心理学系の実験が色々出来ますが，その一環として普通のアンケートも実装可能です。

12.

3.2 適当な回答を除外する表3.1 回答時間の閾値 (Höhne and Schlosser, 2018, Table 1) 下側の閾値上側の閾値平均値 −2 × SD 平均値 +2 × SD 中央値 −1.5× 四分位範囲中央値 +1.5× 四分位範囲中央値 −1.5×（中央値 − 第 1 四分位点）中央値 +1.5×（第 3 四分位点 − 中央値）中央値 −3×（中央値 − 第 1 四分位点）中央値 +3×（第 3 四分位点 − 中央値）下側 1% 上側 1% 分位点の求め方 1 # データに回答時間が無いので，今回はageの分位点を求めます 2 # 引数probsに設定した値の分位点が求められる 3 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 歳未満の大卒」は日本のルール的には不可能，ということがわかります。このようなケースでは，少なくとも年齢または学歴の項目について適当に回答または記入ミスがあるわけですが，一方で心理尺度の部分は正しく回答している可能性も十分にあります。そのため，該当者の全ての回答を除外するかは要検討です。例えば，「学歴によって因子構造が変わるか」を検証したい場合であれば，学歴が正しく記入されていない可能性のあるデータは使用できないでしょう。というように，目的に応じて使うかどうかを決めるのが良さそうです。他にも，成人の ADHD を診断するための尺度である CAARS (Conners’ Adult ADHD Rating Scales; Conners et al., 1999) には矛盾指標というものが用意されています。これは，例えば「1 つの場所に長時間いることが難しい」と「じっと座っているときでも，内心は落ち着かない」のように，普通に考えたらだいたい同じ回答が得られそうな項目ペアについ 12

13.

3.2 適当な回答を除外するて差得点を出すことで，誠実に回答しているかを判断する方法です。特定の構成概念を測定するための一般的な心理尺度では，CAARS の矛盾指標のように確立された項目対が用意されているわけではないですが，通常全ての項目が同じような内容を尋ねているわけですから，同じような考え方で，あまりにも回答が不安定な人は怪しい，と判断できるかもしれません。ということでここからは，不適切回答を検出するためのいくつかの手法を提供している careless パッケージを用いたやり方を紹介していきます。 1 # install.packages("careless") 2 library(careless) 心理尺度の場合，内容をきちんと読んで回答しているならば，基本的には同じ選択肢を選ぶはずです。言い換えると，そのような場合には回答の個人内標準偏差が小さくなっていると考えられます。…という考え方に基づいているのが Intra-individual Response Variability (IRV; Dunn et al., 2018) です。大層な名前がついていますが，やっていることは単純な標準偏差です。careless パッケージには irv() という関数が用意されています。 irv(): 尺度内標準偏差を求める 1 # 最初の5項目（Q1_1-Q1_5）は同じ特性を測定しているので，ばらつきは小さいはず 2 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 項目のように複数の異なる特性に関する項目が並んでいる場合には役に立つかもしれません。 irv() で分割 13 1 # split = TRUEにしてnum.splitで分割数を決める 2 # ただし等分割しかしてくれない 3 irv(dat[, paste0("Q1_", 1:25)], split = TRUE, num.split = 5)

14.

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

15.

3.2 適当な回答を除外する long string index のヒストグラムを見る 1 lsi <- longstring(dat[, paste0("Q1_", 1:25)]) 2 # hist()関数でヒストグラムを作れる 3 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 の値だけで機械的に除外するのではなく，可能な限り回答の内容も踏まえて除外するかを決めるようにしましょう。ちなみに，今回のデータでは項目はランダムに出題されていたと思われます。このような場合，ストレートラインによる除外は意味を成さないので気をつけてください*7 。また，ストレートラインと似た適当回答としては，複数個でパターンを繰り返すもの (e.g., 123451234512345…) や，階段状に回答するもの (e.g., 2345432345432…) などもあります。本当はこれらも全部検出してあげるのがベストなのですが，残念ながら longstring() 関数では対応できないので，自分で関数を作成する必要があります。これはめちゃめちゃ面倒なので今回は省略します*8 。【おまけ】以下では，longstring() 関数を使わずに最長ストレートラインを計算するコードの実装例を紹介します。一応各ステップでの状態も載せておくので，各行で何が行われているか確かめながら読み進めてください。また，やり方は他にも色々あるので，もしよければ自分で考えて実装してみてください。 *7 ランダムに出題されていた場合には，左から実際の出題順に回答を並べたデータを用意できれば同様にストレートラインの計算が可能です。 *8 気になる人は個人的に教えます。ただ相当大変なので覚悟してください。 15

16.

3.2 適当な回答を除外するストレートラインの検出 1 # ただし一人分だけ 2 vec <- dat[3, 2:26] # 面倒なので列番号で指定してしまいます 3 vec 1 2 Q1_1 Q1_2 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 Q1_9 Q1_10 Q1_11 Q1_12 Q1_13 3 3 4 4 5 4 4 4 5 4 2 5 2 4 5 4 5 4 2 3 4 2 3 1 is_same <- vec[, 2:25] == vec[, 1:24] # 前の項目と同じならTRUE 2 is_same # 一旦表示してみる 1 Q1_3 Q1_4 Q1_5 Q1_6 Q1_7 Q1_8 5 Q1_9 Q1_10 Q1_11 Q1_12 Q1_13 3 FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 Q1_23 Q1_24 3 TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 5 Q1_25 6 3 FALSE 1 is_same <- as.numeric(is_same) # 数字に変換 2 is_same # 一旦表示してみる 1 5 2 Q1_2 3 4 4 Q1_25 6 2 4 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 Q1_23 Q1_24 3 5 16 5 TRUE [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 is_same <- paste(is_same, collapse = "") # 全部くっつける 2 is_same # 一旦表示してみる 1 [1] "000110000001100000000010" 1 # 連続する1の最長を数えて1を足せば良いので 2 straight_line <- strsplit(is_same, split = "0") # 文字列を0の位置で切断 3 straight_line # 一旦表示してみる 1 [[1]] 2 [1] "" "" "" "11" "" "" 3 [14] "" "" "" "" "1" "" "" "" "" "11" "" "" ""

17.

3.2 適当な回答を除外する 1 # 文字数に変換 2 # strsplit()はリストを返すため，リストの1番目を選択 3 # リストについてはそのうち説明します 4 straight_line <- nchar(straight_line[[1]]) 5 straight_line # 一旦表示してみる 1 [1] 0 0 0 2 0 0 0 0 0 2 0 0 0 0 0 0 0 0 1 1 # 最大値に1を足せば最長が分かる 2 max(straight_line) + 1 1 [1] 3 本当はこれを全員に対して適用した上で，一定数より大きい人を除外するか検討する，といった手続きを取る必要があるわけです。 3.2.5 人と違う回答心理尺度の場合，同じ構成概念を尋ねている各項目はそれなりに相関しているはずです。もちろん項目によってそもそも当てはまりにくい・当てはまりやすい項目はあるのですが，総合すると項目への回答には「平均的なパターン」と言えるものがありそうです。「人と違う」を検出する方法は，異常検知の知見が使えたりします。したがって機械学習を利用して検出する方法なども研究されていたりしますが，ここでは伝統的なマハラノビス距離を用いた検出方法を紹介＆実践してみたいと思います。図3.4は，２項目の連続的な項目（６段階評価ではなく，スライダーみたいなもので回答したと想像してください）への回答をプロットした仮想的なものです。赤い点がそれぞれの軸の平均値を表しています。青い点 (1) と (2) は，それぞれ赤い点からの直線距離はほぼ同じなのですが，2 変数の相関を考えると (2) の方はさほどおかしな値ではなさそうです。正の相関があるため，平均値から離れていたとしても 2 項目とも一貫して高い or 低い値ならばありえる，ということですね。 1 Warning: The `size` argument of `element_line()` is deprecated as of 2 3.4.0. 3 i Please use the `linewidth` argument instead. 4 This warning is displayed once every 8 hours. 5 Call `lifecycle::last_lifecycle_warnings()` to see where this warning 6 was generated. ggplot2 このような時に使えるのがマハラノビス距離です。変数の相関関係を考慮した上でベクトル間の距離を出してくれます。回答者 𝑖 の回答ベクトルを 𝐱𝑖 ，全回答者の平均ベクトルを 𝐱̄ ， 17

18.

3.2 適当な回答を除外する 4 (1) 3 (2) 2 X2 1 0 -1 -2 -3 -4 -4 -3 -2 -1 0 X1 1 2 3 4 5 図3.4 仮想的な 2 項目変数の分散共分散行列を 𝚺𝐱 とすると，マハラノビス距離は ̄ −1 ̄ 𝑇 𝐷 = √(𝐱𝑖 − 𝐱)𝚺 𝐱 (𝐱𝑖 − 𝐱) (3.1) で計算することが出来ます。マハラノビス距離は，多変量正規分布と密接な関係を持っています。というのも，𝑘 個の変数に関する多変量正規分布の確率密度関数は 𝑓(𝐱; 𝛍, 𝚺) = 1 √(2𝜋)𝑘 |𝚺| 1 𝑇 exp (− (𝐱𝑖 − 𝛍)𝚺−1 𝐱 (𝐱𝑖 − 𝛍) ) 2 (3.2) と表されますが，よく見るとこの式の中には (3.1) 式のルートの中身によく似たものが入っています。実際に，データから多変量正規分布の平均ベクトル 𝛍 と分散共分散行列 𝚺 を最尤推定する場合，最尤推定量はそれぞれ標本平均ベクトル 𝐱̄ と標本共分散行列 𝚺𝐱 となることが知られています。したがって，マハラノビス距離はデータから推定された多変量正規分布の確率密度に反比例するということがわかります。言ってしまえば，マハラノビス距離は個人レベルでみたときの尤度みたいなもの，というわけですね。図3.4で言えば，塗りつぶされている楕円が，データから推定されたニ変量正規分布における大まかな確率密度を表していますが，マハラノビス距離はこの確率密度に応じて計算されているわけです。そう考えると，確かに (2) の点のほうが (1) の点よりも中心からの距離は短い，と判断できそうです。 R にはマハラノビス距離（の二乗）を計算するための関数 mahalanobis() が用意されているため，これを使って図3.4の 2 つの青い点のマハラノビス距離を計算してみます。 mahalanobis() では，基本的に 3 つの引数をとります。第一引数 x には x𝑖 を，第二引数 center*9 には x̄ を，第三引数 cov には 𝚺𝐱 に相当するものを入れてあげましょう。 *9 18 ふつうマハラノビス距離は「2 つのベクトルの距離」なので 2 つの引数は x と y のほうがしっくり来るのですが，この関数の第二引数が center なのは，この関数が「データフレームの各行について，平均ベクトルからのマハラノビス距離を計算する」という，まさにいま行っている分析の目的にピッタリの関数だからです。

19.

3.2 適当な回答を除外するマハラノビス距離の計算 1 # ここのコードは見ているだけでOKです。 2 mean_vec <- c(0, 0) # 平均ベクトル 3 cov_matrix <- matrix(c(1, 0.9, 0.9, 1), nrow = 2) # 分散共分散行列 4 # 各変数の分散が１，相関係数が0.9という想定です。 5 cov_matrix # 一応表示してみる 1 [,1] [,2] 2 [1,] 1.0 0.9 3 [2,] 0.9 1.0 1 vec1 <- c(-3, 3) # 青い点(1)の座標 2 # 青い点(1)のマハラノビス距離 3 mahalanobis(vec1, center = mean_vec, cov = cov_matrix) 1 [1] 180 1 vec2 <- c(3, 3) # 青い点(2)の座標 2 # 青い点(2)のマハラノビス距離 3 mahalanobis(vec2, center = mean_vec, cov = cov_matrix) 1 [1] 9.473684 ということで，青い点 (1)，(2) と平均ベクトルとのマハラノビス距離（の二乗）はそれぞれ 180，9.47 となり，たしかに (1) のほうが大きな外れ値であるということがわかりました。ありがたいことに careless() パッケージには，データフレームからマハラノビス距離を計算するための関数も用意されています。先程使用した mahalanobis() 関数では，自分でデータから平均ベクトルと分散共分散行列を計算する必要があったのですが，careless パッケージにある mahad() 関数ではデータフレームを一つ与えるだけで各行のマハラノビス距離をまとめて計算してくれます*10 。ありがたい。 mahad(): マハラノビス距離を計算する 1 1 2 19 mahad(dat[, paste0("Q1_", 1:25)]) [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 3 [13] 34.052697 38.710674 18.843160 18.056494 26.822273 38.451815 4 [19] 47.752717 30.685610 24.097729 27.202401 10.453943 30.186948 5 [25] 20.393634 12.667231 45.024926 23.116322 26.807564 12.173940 *10 内部の計算は psych::outlier() 関数を使っています。なのでこちらを使用しても OK です。

20.

[beta]

3.2 適当な回答を除外する

6

[31] 56.951093 20.122730 26.353218 30.358047 17.898820 28.093556

7

[37] 14.465344 10.483335 11.724812 27.950268

8

[43] 23.731495 20.956441 42.950337 14.636132 32.755405 52.321123

9

7.693067 57.603758

[ reached getOption("max.print") -- omitted 2388 entries ]

これで無事，全員分のマハラノビス距離を一気に計算できました。あとはこれを利用して怪
しい人を検出するだけです。この場合でも閾値の問題が登場します。これまで通り，ヒスト
グラムを描いてみたり，quantile() を利用してみても良いのですが，mahad() 関数にはも
う少しテクニカルな方法が実装されています。

𝑝 次元多変量正規分布に従うデータから計算されるマハラノビス距離の二乗は，理論的には
自由度 𝑝 の 𝜒2 分布に従うことが知られています。これを利用すると，「自由度 𝑝 の 𝜒2 分布
における上位 𝑘 パーセンタイル点より大きい人を検出しよう」という考えに至ります。
マハラノビス距離に基づく検出
1

# flag = TRUEにすると，フラグを立ててくれる

2

# confidenceによって，閾値のパーセンタイル点を決める

3

m_dist <- mahad(dat[, paste0("Q1_", 1:25)], flag = TRUE, confidence =

4

m_dist

0.99)

1

d_sq flagged

2

1

13.468425

FALSE

3

2

25.677291

FALSE

4

3

13.490783

FALSE

5

4

28.823203

FALSE

6

5

18.946470

FALSE

7

6

25.901522

FALSE

8

7

9.260945

FALSE

9

8

48.474541

TRUE

10

9

12.794634

FALSE

11

10 19.207672

FALSE

12

11 19.386422

FALSE

13

12 14.485101

FALSE

14

13 34.052697

FALSE

15

14 38.710674

FALSE

16

15 18.843160

FALSE

17

16 18.056494

FALSE

18

17 26.822273

FALSE

19

18 38.451815

FALSE

20

19 47.752717

TRUE

21

20 30.685610

FALSE

22

[ reached 'max' / getOption("max.print") -- omitted 2416 rows ]

ここで flagged 列が TRUE になっている人が，マハラノビス距離が閾値より大きい，すなわ

20

21.

3.2 適当な回答を除外するち異質な回答パターンを示した人だということです。ちなみに，マハラノビス距離が最大の人の回答ベクトルを見てみると，たしかにこの人は６か１かばかりで回答しており，どうやら適当に回答しているような気がします。ただ実際には日頃から極端な選択肢を選ぶ傾向がある人かもしれないので，これだけで機械的に除外するのではなく，回答ベクトルを見てストレートラインの傾向が無いかなども確認して除外の判断をしたほうが良いでしょう。 1 # which.max()はベクトルの中で最大値の「位置」を返す (arg max 的な) 2 # max()は最大値そのものを返す 3 # which.min()もあるよ 4 dat[which.max(m_dist$d_sq), ] 1 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 2 867 867 3 6 6 6 6 6 6 5 1 1 1 6 6 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 Q1_23 4 867 5 1 6 6 1 6 1 6 1 6 6 1 Q1_24 Q1_25 gender education age 6 867 6 1 1 1 31 マハラノビス距離による検出法の弱点マハラノビス距離は平均ベクトルからの距離を求めています。そのため，同じストレートラインでも 1 と 6 ばかりなど極端な選択肢を選んでいる人は引っかかりやすい一方で，3 と 4（どちらでもない周辺）ばかりを選んでいる人はどうしてもマハラノビス距離は小さくなってしまいます。同様に，どちらでもない周辺で適当な回答をしている人は検出しにくいという問題点があります。これを克服するためには，マハラノビス距離に加えて IRV やストレートラインなどの指標，あるいは可能であれば回答データ以外（回答時間など）の複数の指標を組み合わせて検出してあげるのが良さそうです。ここまで適当回答を検出するためのいくつかの方法を紹介してきました。実際のデータにおいて検出＆除外を行うときにはいろいろと考えるべきことはあるのですが，とりあえず今回は機械的に「long string index (lsi) が 8 以上」かつ「マハラノビス距離が上位 5%*11 」の人を除外することにします。ちなみにマハラノビス距離の上位 5% は，図3.5の分布でいえば赤い線より右側の人であり，極端に大きい人たちだろうということが言えそうです。では，実際に 2 つの基準をともに満たす人を除外してみましょう。ここでは，複数の条件式を結合するための論理演算子を使用します。論理演算とは，複数の TRUE/FALSE をもとに 1 つの TRUE/FALSE を返す演算のことです。何はともあれ，以下の例を見てください。 *11 21 勘違いしないでいただきたいのですが，「マハラノビス距離の上位 5% を除外するのが定番」ということではありません。あくまでも授業の演習としてとりあえず除外を体験してみよう，というくらいの意味合いです。

22.

3.2 適当な回答を除外する 600 400 0 200 Frequency 800 Histogram of m_dist$d_sq 0 20 40 60 80 m_dist$d_sq 図3.5 マハラノビス距離の二乗の分布 &: 2 つとも TRUE のときだけ TRUE を返す (AND) 1 TRUE & TRUE 1 [1] TRUE 1 TRUE & FALSE 1 [1] FALSE 1 FALSE & FALSE 1 [1] FALSE |: いずれかが TRUE であれば TRUE を返す (OR) 22 1 TRUE | TRUE 1 [1] TRUE 1 TRUE | FALSE 1 [1] TRUE 1 FALSE | FALSE 1 [1] FALSE 100

23.

3.2 適当な回答を除外する今回は，AND 演算子 (&) を用いて「long string index (lsi) が 8 以上」かつ「マハラノビス距離が上位 5%」の人を探し出したいと思います。複数条件を満たす人を見つける 1 # マハラノビス距離の閾値 2 thres <- quantile(m_dist$d_sq, prob = 0.95) 3 thres 1 95% 2 48.87949 1 # 除外対象を決める 2 is_remove <- (m_dist$d_sq >= thres) & (lsi >= 8) 3 is_remove 1 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 2 [12] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 3 [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 4 [34] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE 5 [45] FALSE FALSE FALSE FALSE FALSE FALSE 6 [ reached getOption("max.print") -- omitted 2386 entries ] 1 # 何行目の人を除外するか 2 which(is_remove) 1 [1] 1243 1359 1775 2018 1 # 一応除外対象の人の回答をチェック 2 subset(dat, is_remove) 1 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 2 1430 1430 1 1 1 1 1 1 1 1 1 1 1 1 3 1560 1560 1 6 6 1 1 5 6 6 4 6 1 6 4 2043 2043 1 1 1 1 1 1 1 1 1 1 1 1 5 2313 2313 6 1 1 6 1 6 5 6 1 1 6 6 6 Q1_13 Q1_14 Q1_15 Q1_16 Q1_17 Q1_18 Q1_19 Q1_20 Q1_21 Q1_22 Q1_23 7 1430 1 1 1 1 1 1 1 1 1 1 1 8 1560 6 1 6 6 6 6 6 6 6 6 4 9 2043 1 1 1 1 1 1 1 1 1 1 1 10 2313 1 1 1 1 1 1 1 1 6 1 1 11 23 Q1_24 Q1_25 gender education age 12 1430 1 1 1 3 23 13 1560 6 1 2 5 31

24.

[beta]

3.3 逆転項目の処理

14

2043

1

1

1

3

19

15

2313

6

1

1

3

23

1

# is_removeがFALSEの人だけ残せば良いので

2

dat <- subset(dat, !is_remove)

ということで，残ったデータは
1

nrow(dat)

1

[1] 2432

となりました。以後の分析はこの 2432 人分のデータで進めていきます。

3.3 逆転項目の処理
逆転項目がある場合には，得点をひっくり返す必要があります。psych::bfi では*12 1, 9,

10, 11, 12, 22, 25 番目の項目が逆転項目となっています。リッカート尺度の場合，逆転項目
に対する処理は「
（最大値 +1）から回答の値を引く」だけで OK です。ということで，一つ
の列について逆転項目の処理を行うならば以下のコードで可能となります。
逆転項目の処理：1 項目編
1

dat[, "Q1_1"] # 処理前

1

[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

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

3

[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

4

[97] 1 6 2 3

5

[ reached getOption("max.print") -- omitted 2332 entries ]

1

dat[, "Q1_1"] <- 7 - dat[, "Q1_1"]

2

dat[, "Q1_1"] # 処理後

1

[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

2

[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

3

[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

4

[97] 6 1 5 4

5

[ reached getOption("max.print") -- omitted 2332 entries ]

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

apply() の力を借りましょう。
*12

24

psych::bfi.keys で確認可能です。変数名の頭にマイナスがついているものが逆転項目です。

25.

[beta]

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

# function()関数は「関数を定義する」関数

2

# ()の中には新しく作る関数の引数を

3

# {}の中には関数の動作を記入する

4

tashizan <- function(x, y) {

5

x + y

6

}

7

tashizan(3, 2)

1

[1] 5

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

subtract_7 <- function(x) {

2

7 - x

3

}

4

# 挙動の確認

5

head(dat$Q1_9) # 処理前

1

[1] 4 3 2 5 3 1

1

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

1

[1] 3 4 5 2 4 6

これを apply() 関数の引数 FUN に指定してあげれば，すべての列に同じ処理を適用するこ
とができるわけです。
逆転項目の処理：複数項目編
1

# 処理する列名を指定

2

# 上のコードでQ1_1だけ逆転処理済みなのでここでは除外

3

cols <- paste0("Q1_", c(9, 10, 11, 12, 22, 25))

4

# cols <- c("Q1_9", "Q1_10", "Q1_11","Q1_12","Q1_22","Q1_25")

5

# ↑これでもよいが，列数が増えた時には大変

6
7

25

head(dat[, cols]) # 処理前

26.

3.3 逆転項目の処理 1 Q1_9 Q1_10 Q1_11 Q1_12 Q1_22 Q1_25 2 1 4 4 3 3 6 3 3 2 3 4 1 1 2 3 4 3 2 5 2 4 2 2 5 4 5 5 5 3 3 5 6 5 3 2 2 2 3 3 7 6 1 3 2 1 3 1 1 dat[, cols] <- apply(dat[, cols], 2, subtract_7) 2 head(dat[, cols]) # 処理後 1 Q1_9 Q1_10 Q1_11 Q1_12 Q1_22 Q1_25 2 1 3 3 4 4 1 4 3 2 4 3 6 6 5 4 4 3 5 2 5 3 5 5 5 4 2 2 2 4 4 2 6 5 4 5 5 5 4 4 7 6 6 4 5 6 4 6 これで，ようやく全ての逆転項目の処理が完了しました*13 。最後に保存 26 1 saveRDS(dat, "chapter03.rds") *13 実を言うと，apply() なんか使わなくても 7 - dat[,cols] で複数行まとめて処理は可能です。今回は自作関数の導入を兼ねてちょっと回りくどいやり方にしました。

27.

参考文献参考文献 Conners, C. K., Erhardt, D., and 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.jesp.2015.07.006 Dunn, A. M., Heggestad, E. D., Shanock, L. R., and Theilgard, N. (2018). Intraindividual response variability as an indicator of insuﬀicient 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 Huang, J. L., Curran, P. G., Keeney, J., Poposki, E. M., and DeShon, R. P. (2012). Detecting and deterring insuﬀicient effort responding to surveys. Journal of Business and Psychology, 27 (1), 99–114. https://doi.org/10.1007/s10869-011-9231-8 Höhne, J. K. and 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 Kong, X. J., Wise, S. L., and 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/00131644 06294779 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. URL： https://repository.tku.ac.jp/dspace/bitstream/11150/10813/1/JMIC01-03.pdf Oppenheimer, D. M., Meyvis, T., and Davidenko, N. (2009). Instructional manipulation 27

28.

参考文献 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. and 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 28

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

Kyosuke Bunji

関連スライド

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

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

ベイズ統計_03_尤度

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

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

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

各ページのテキスト