心理統計Rワークショップ
実験調査法・統計解析法・データ可視化法・Tidyverseを学べる
杉本海里 / Kairi Sugimoto
公開: 2024/7/19、最終更新: 2025-10-05
1 はじめに
1.1 この教材について
1.1.1 作成者
心理学分野の研究・教育に携わっている杉本海里（https://kairi-pf.com）が作成した教材です。所属ラボメンバーの利用が主な目的ですが、一般の方も利用できるように公開しました。
本教材で使用しているデータはすべて架空のもので、実際の事象を反映したデータではありません。本教材はR Markdownで作成しています。
本資料に関してお気づきの点やご感想などがありましたら、下記フォームからご連絡をお願いいたします。
1.1.2 教材の特徴
- 論文を読み書きできるようになることが目標です。
- 分かりやすい事例を通して、研究デザイン設計と「データの前処理→作図→仮説検定→結果の報告」の一連の解析フローを、Tidyverseパッケージの主要関数とともに学べます。
- アメリカ統計協会による「統計的有意性とP値に関するASA声明」での指摘を組み込んでおり、p値によるデータ解釈の注意点や、伝統的な検定手法の限界を理解できるように配慮しています。
- そのほか、下記の重要トピックに触れています。
- 効果量と信頼区間の報告
- p-hackingと事前登録
- サンプルサイズ設計
- 検定の多重性問題
- ノンパラメトリック検定の限界
- 線形混合モデル・一般化線形混合モデル
1.1.3 想定利用者
仮説検定の概念に触れたことがあり、プログラミング言語における変数・関数の概念を理解している大学学部生であれば理解できます。
1.1.4 本教材の動作環境
R.version.string## [1] "R version 4.5.1 (2025-06-13)"1.2 学べること
この教材は、計3コマ（100分 × 3コマ = 5時間）で学習することを想定しています。事前準備と任意学習を含めると、7時間程度です。上から順番に学ぶように設計していますので、事例1を飛ばして事例2を学ぶようなことはできません。
| 章 | 所要時間 | 学べる検定 | 学べるグラフ |
|---|---|---|---|
| 1〜3章（事前準備） | 60分 | - | - |
| 4章（事例1：無相関検定） | 100分 | Pearsonの相関係数を用いた無相関検定 Spearmanの相関係数を用いた無相関検定 |
ヒストグラム カーネル密度推定図 散布図 |
| 5章（事例2：t検定） | 100分 | 対応のあるt検定 対応のないWelchのt検定 Brunner-Munzel検定 （サンプルサイズ設計） |
箱ひげ図 ドットプロット 棒グラフ |
| 6章（事例3：分散分析） | 100分 | 参加者内1要因分散分析 混合2要因分散分析 参加者間2要因分散分析 |
折れ線グラフ ファセット機能 |
| 7章（任意学習） | 60分 | \(\chi\)二乗検定（1x2分割表、2x2分割表） 単回帰分析 重回帰分析 線形混合モデル（LMM） |
- |
記載されているRのコードを自身のPC上のRソフトウェアにコピー&ペーストして、コードが動作することを確かめながら理解を深めていく学習スタイルです。
Rによる統計解析は、実際に自分で考えて自分の手で解析をしてみないと全く身につかないため、事例ごとに練習問題を用意しています。
1.3 なぜR言語なのか
心理統計の現場では主に下記のような言語・統計ソフトウェアが使用される。心理学部で頻繁に使われるのはRとSPSSの2種類である。
- R：心理学分野のみならず、幅広い学問分野で広くよく使われる統計ソフトウェア。
- SPSS：大規模調査データを扱うのが得意。有料。
- Python：もとからPythonに慣れている人は使うことがある。Rは統計解析に特化しているが、Pythonは汎用性が高い。
- Matlab：生理データの分析・可視化が得意。有料。
- Stata：医療系などでよく使われる。有料。
- HAD：Excelだけで分析ができるのが特徴。
- JASP：GUIなので簡単に解析ができる。ベイズ統計に対応している点がポイント。
- Jamovi：JASPに似たソフトウェア。
- EZR：医療統計分野でよく使われる。
とくにRは、拡張性が極めて高いためどんな解析でも実現でき、作図機能が強力であり、さらに無料で使えるオープンソースであることから、国内外を問わず非常に有力な統計ソフトウェアである。
近年は、JASPなどのGUI方式の無料統計ソフトウェアも普及してきたが、再現性のある解析処理のために、Rを共通言語とするのが望ましい。
実際、統計学分野の学術論文では、提唱した解析法をR言語で実装したコードを付録として公開することが一般的である。また、統計理論の書籍でも、Rによる処理を行いながら説明するものが多く、R言語を習得していないとこれらの文献を通読できない。
2 事前準備
2.1 解析環境の整備
はじめに、以下の手順に従って、Rによる開発環境を用意すること。
PC上のわかりやすい場所に、本教材の処理を行うための作業用フォルダを1つ新規作成する。フォルダ名は何でもよいが、たとえば「r_practice」などと、半角英数字かつスペースを含まないように命名する。
- Windowsユーザーのうち、OneDriveオンラインストレージにPC上のファイルを同期させている場合は、OneDriveの同期が行われていないディレクトリにフォルダを作成すること。
「R」と「R Studio」を「https://posit.co/download/rstudio-desktop/」からインストールする。「R」のインストールを完了させてから、「R Studio」をインストールすること。
- 【Windowsユーザーが「R」をインストールするとき】Rのダウンロードサイトに行き、画面最上部のDownload and Install Rの中の「Download R for Windows」→「base」をクリックして、一番上に表示される最新版のインストーラー（exeファイル）をダウンロードし、ダウンロードしたインストーラーを開いてインストールを完了させる。インストーラーで色々とインストール設定を選べるが、基本的にすべて初期設定のまま進めばよい。
- 【Macユーザーが「R」をインストールするとき】Rのダウンロードサイトに行き、画面最上部のDownload and Install Rの中の「Download R for Mac」をクリックする。Latest releaseという欄でインストーラー（pkgファイル）をダウンロードできるが、「For Apple silicon Macs」と「For older Intel Macs」の2つの種類があるため、自身のMacのOSにマッチした方をダウンロードする。ダウンロードしたpkgファイルを開いてインストールを完了させる。インストーラーで色々とインストール設定を選べるが、基本的にすべて初期設定のまま進めばよい。
Macユーザーのみ、追加で「XQuartz」を「https://www.xquartz.org」からインストールする。
「R Studio」を開く。R Studioは、Rソフトウェアを便利に使うための統合開発環境である。すべての操作はR Studio上で行うため、Rソフトウェアは開くことはない。
R Studioの画面左上の「File」→「New File」→「R Script」を選択して、新しいRスクリプトファイルを作成する。スクリプトファイルは画面左上の領域に表示される。スクリプトファイルを作成したら、Ctrl+S（MacはCommand+S）で、先ほど作成した作業用フォルダにスクリプトファイルを保存する。ファイル名は何でもよいが、半角英数字で「script.R」などと命名するとよい。
この時点で、R Studioの画面は大きく4分割されているはずである。主に使用するのは左側2つである。それそれの画面は以下の役割を持っている。
- 左上（ソース画面）：スクリプトファイルを編集するところ。コードを実行する時は、実行したいコードを選択した上で、Ctrl+Enter（MacはCommand+Enter）をする。コードを実行すると、実行したコードと出力結果が左下のコンソール画面に表示される。
- 左下（コンソール画面）：コードを実行したり出力結果が表示されるところ。
- 右上（変数画面）：作成した変数の中身を確認できる。
- 右下（プロット画面）：作成したグラフを確認できる。
練習として、試しになにかコードを打ってみる。左上画面のスクリプトファイルに半角で「1+2」と入力し、打ち込んだ行をマウスドラッグで選択して、Ctrl+Enter（MacはCommand+Enter）を押す。そうすると、左下のコンソール画面に「[1] 3」と表示され、1+2が正しく計算されたことがわかる。
画面中央上端の「Session」→「Set Working Directory」→「Choose Directory」を選び、先ほど作成した作業用フォルダを選択する（Macの場合は、画面上端にマウスカーソルを移動すると表示される）。ここで選択したフォルダは作業ディレクトリとなり、この作業ディレクトリの階層でR Studioが動作するようになる。ディレクトリの概念がわからない場合は、下記説明を参考にすること。
- ディレクトリとは、コンピュータ内においてファイルが置かれている階層のことであり、フォルダとほぼ同義である。
- Rでデータを読み込む時は、作業ディレクトリに読み込むデータを置いておかなければならない。
- R上で図などのファイルを出力すると、作業ディレクトリにそのファイルが保存される。
作業ディレクトリの設定を行うと、画面左下のコンソール画面に「setwd(“/~~~/~~~/~~~…”)」というコードが表示されるので、このコードをスクリプトファイルの一番上に貼り付ける。本来であればR Studioを開くたびに作業ディレクトリはリセットされてしまうが、スクリプトファイルに「setwd(“/~~~/~~~/~~~…”)」と書いておけば、毎回このコードを実行するだけで自動的に作業ディレクトリが設定されるようになり便利である。
2.2 解析データの入手
下記リンク先より、6種類のデータ（data1.xlsx / data1p.xlsx / data2.xlsx / data2p.xlsx / data3.xlsx / data3p.xlsx）をダウンロードし、作業ディレクトリに格納すること。
2.3 パッケージの読み込み
Rに最初から用意されている関数だけでは、限られた統計解析しかできない。さまざまな関数を使用するために、外部パッケージを最初に読み込む必要がある。
まずはパッケージをコンピュータにインストールし、そのパッケージを読み込むという２段階の手順をとる。パッケージのインストールは初回のみ行えばよいが、パッケージの読み込みはR Studioを起動するたびに最初に行わなければならない。
まずはinstall.packages("パッケージ名")でパッケージのインストールが必要となる。以下に、本教材で使用するすべてのパッケージをまとめてインストールできるコードを示すので、R
Studioの左下のコンソール画面にコピー&ペーストしてエンターキーで実行すること。なお、スクリプトファイルではなくコンソール画面に直接入力するのは、初回しか実行しないためである。コンソール画面が落ち着き、プロンプト（>マーク）が表示されたらインストール完了である。もし「〜〜〜(yes/no)?」などとコンソール画面に表示されたら、yesと入力してエンターキーを押すとインストールが進む。なお、はじめてパッケージをインストールするときは、インストール元のサーバーの指定として、Primary
CRAN
repositoryを選ぶ画面が出てくることがある。その際は、Japanと書かれているサーバーを選ぶこと。
install.packages(c("tidyverse", "openxlsx", "plotrix", "TOSTER", "effsize", "pwr", "see", "ggridges"))例外として、分散分析用パッケージである「anovakun_489.txt」は、install.packages関数ではインストールできないため、開発者である井関先生のWEBサイト【井関-ANOVA君】より手動でダウンロードして、作業ディレクトリに保存すること。
すべてのパッケージのインストールが終わったら、パッケージの読み込みを行う。下記に使用するすべてのパッケージを読み込むコードを記載している。この読み込みコードは、R Studioを起動するたびに実行しなければならないため、コンソール画面ではなく、スクリプトファイルに記入して保存しておき、毎回実行するようにする。
library(tidyverse) #dplyr, tidyr, ggplot2を含む総合パッケージ
library(openxlsx) #Excel読み込み（read.xlsx関数）
library(plotrix) #標準誤差計算（std.error関数）
library(TOSTER) #Brunner-Munzel検定（brunner_munzel関数）
library(effectsize) #効果量算出（cohens_d関数）
library(pwr) #サンプルサイズ設計（pwr関数）
library(see) #Okabe-Itoカラーパレット（scale_colour_okabeito関数）
library(ggridges) #カーネル密度推定（geom_density_line関数）
source("anovakun_489.txt") #分散分析用パッケージ（anovakun関数）この時点で、スクリプトファイルは下図のようになっているはずである。1〜9行目のコードは、R studioを開くたびに最初に実行する必要がある。これ以降は、見やすさのために10行目は空けて、11行目からコードを打ち込んでいくこととなる。
また、作業ディレクトリには以下の8個のファイルが格納されているはずである。
3 基礎知識の確認
3.1 重要語句の意味
本教材を読み進めるために最低限理解が必要となる基本語句である。
| 語句 | 意味 |
|---|---|
| 母集団 | 推測の対象とする集団全体のこと。多くの研究は、理想的には全人類を母集団とおく。 |
| 標本（サンプル） | 母集団から抽出したデータのこと。たとえば、20人分の心理測定データを集めたら、標本サイズ（サンプルサイズ）は20となる。 |
| 記述統計 | 標本の特徴を記述する統計手法。たとえば、平均値・中央値などが該当する。記述統計では、標本の傾向を捉えることはできても、母集団における事象の推測は難しい。 |
| 推測統計 | 標本に基づいて母集団における事象を推測する統計手法。仮説検定は、推測統計の最も主要な手法である。 |
| 正規性 | 検定の仮定の1つ。母集団の確率分布が正規分布していること。なお、正規分布は左右対称の山なりの形状である。 |
| 等分散性 | 検定の仮定の1つ。比較する母集団の分散（ばらつき）が等しいこと。 |
| 独立性 | 検定の仮定の1つ。標本がそれぞれ母集団から無作為（ランダム）に抽出されたもので、互いに独立していること（≒サンプリングバイアスがないこと）。 |
| 要因と水準 | たとえば、「つぶあん派」と「こしあん派」の違いを調べるとき、要因は「あんこの状態」となり、水準は「つぶあん」と「こしあん」の2水準となる。 |
| 参加者間要因 | 水準間で参加者（標本）が異なる要因。たとえば、A組とB組のテストの得点を比較するとき、「クラス」が参加者間要因となる。 |
| 参加者内要因 | 水準間で参加者（標本）が同一の要因。たとえば、A組における5月と6月のテストの得点を比較するとき、「テストの実施時期」が参加者内要因となる。 |
3.2 仮説検定の考え方
推測統計では、我々が実際に取得した標本データに基づいて、母集団の事象を推測することとなる。その際に、母集団に関する仮説を調べるためのプロセスを検定と呼ぶ。たとえば以下のような問いを検証する際に使われる。
- 変数Aと変数Bに有意な相関関係があるかどうか。
- 変数Aと変数Bに有意な差があるかどうか。
検定では、はじめに帰無仮説と対立仮説の2つの仮説を立てる。原則的に、帰無仮説は有意性がないという仮説となり、対立仮説は有意性があるという仮説となる。たとえば、問いごとの帰無仮説と対立仮説は下記のようになる。
- 変数Aと変数Bに有意な相関関係があるかどうか。
- 帰無仮説：変数Aと変数Bの母集団における相関係数は0である（相関関係はない）。
- 対立仮説：変数Aと変数Bの母集団における相関係数は0ではない（相関関係がある）。
- 変数Aと変数Bに有意な差があるかどうか。
- 帰無仮説：変数Aと変数Bの母集団における平均値は同じである（差はない）。
- 対立仮説：変数Aと変数Bの母集団における平均値は同じではない（差がある）。
仮説を立てたら、帰無仮説が棄却されるかどうかを調べ、帰無仮説が棄却された場合に対立仮説を採択する。このような背理法的な仕組みを用いるのは、有意性を立証するときに、有意性があることを直接証明するよりも、有意性がないことを否定するほうが統計的に容易だからである。
注意点として、帰無仮説が棄却されれば対立仮説は採択されるが、帰無仮説が採択されることはない（ないことを証明することは困難であるため）。したがって、有意な相関（あるいは差）が見られなかったときに、相関（あるいは差）がないと主張することは誤りである。
- 帰無仮説が棄却された場合 → 対立仮説が支持される。
- 帰無仮説が棄却されなかった場合 → 帰無仮説と対立仮説の判定は保留となる（帰無仮説が支持されたり対立仮説が棄却されることはない）。
帰無仮説が正しいのに誤って棄却してしまうことをアルファエラー（第一種の過誤）と呼ぶ。当然アルファエラー率は小さいほうがよいが、小さくしようとすると、対立仮説が正しいのに帰無仮説を棄却できなくなるベータエラー（第二種の過誤）の危険性が高まる。検定では、これら2つのエラー率のバランスと、研究間での一貫性が重要である。そこで、アルファエラーを容認するパーセンテージ（有意水準α）を5%と設定するという科学界のコンセンサスが存在する。簡単に言えば、帰無仮説のもとで20回中1回しか起きないような珍しい現象が起きたら、それは帰無仮説では説明できずに対立仮説で説明すべき特異的現象と判断するということである。
帰無仮説の棄却を判断するためには、伝統的にp値と呼ばれる指標が用いられる。p値は下記のように説明される。（一部参考：アメリカ統計協会、2017、統計的有意性とP値に関するASA声明、佐藤訳）
- p値はデータと特定の統計モデルが矛盾する程度をしめす指標のひとつである。p値が大きいほど矛盾しないことを意味する。
- p値が0.05を下回る場合に帰無仮説を棄却するというコンセンサスがある。（ただし絶対的基準ではないため盲信してはいけない）
- サンプルサイズが大きいほどp値は小さくなる。サンプルサイズの過剰・不足を避けるために、事前のサンプルサイズ設計（最適なサンプルサイズを決める手続き）が推奨される。
- p値のみに基づいて効果の有無や強さなどを結論づけてはならない。p値による二分法的判断は論証材料の1つに過ぎず、研究デザインの妥当性や効果の大きさを示す効果量なども含めて総合的に判断する必要がある。
- よくある間違いとして、p値は調べている仮説が正しい確率や、データが偶然のみでえられた確率を測るものではない。
つまり、特定の統計モデルを仮定した母集団の確率分布に対して、標本データがどのくらい矛盾しているかをp値として算出し、しきい値として設定した有意水準5%よりもp値が小さい（＝矛盾度が大きい）場合に、帰無仮説を棄却するという仕組みである。
検定過程では、p値のほかにも以下のような数値が算出される。論文では、これらの情報を漏れなく報告することが、再現性のある実証科学のために重要となる。
- p値を算出するための「検定統計量」。検定統計量にはt値、F値、カイ二乗値などがあり、使用する検定手法によって異なる。
- サンプルサイズに基づく「自由度（df: degree of freedom）」。
- 対立仮説で仮定した効果の大きさを示す「効果量」。
- 母平均が含まれると推測される区間を示す「95%信頼区間（CI: confidence interval）」。なお「95%」とは、「信頼区間に95%の確率で母平均が含まれる」という意味ではなく、「標本抽出と信頼区間の計算を何度も繰り返したときに95%の頻度で信頼区間が母平均を含む」という意味である。
重要なこととして、母集団に対して仮定した統計モデルが誤っていると、検定は意味をなさない。統計モデルとは一般的に、確率分布の形状、正規性、等分散性、独立性などの仮定を指す。検定を行う前に、仮定した統計モデルが適切かどうかを十分に疑う必要がある。
なお、p値が0.05を下回るかどうかによって結果が明瞭に分かれる二分法的手法への批判が近年高まっている。それに対処するための方法として、効果の大きさと精度を同時に示すことができる信頼区間を明示することが求められている（Amrhein, Greenland, & McShane, 2019）。さらに、対立仮説の条件付き確率（≒対立仮説が正しい確率）を直接算出することによって事象を説明するベイズ統計が、検定とならぶ有力な手法として普及しつつある。すべてにおいてベイズが優位であるわけではなく、どちらの手法も長所短所があるため、目的に応じて仮説検定とベイズ統計を使い分けたり、あるいは併用することが重要である。
- 筆者は統計学の専門家はではないため、解説は必要最低限に留めた。数学的な議論は下記文献を参照したい。
- 久保川（2023）「データ解析のための数理統計入門」
- 南風原（2002）「心理統計学の基礎」
- 南風原（2014）「続・心理統計学の基礎」
- 奥村（2016）「Rで楽しむ統計」
3.2.1 理解度確認クイズ
とある掃除機メーカーの旧製品と新製品の満足度評価に有意な差があるかどうかを調べたい。20人にアンケートをとった。以下の記述のうち、下線部が正しいものをすべてチェックしなさい。
正解を表示する
正解は【1】【3】
3.3 読み込むデータのフォーマット
Rに読み込ませるデータをExcelなどで作成する場合が多いだろう。その時、正しいフォーマットにしていなければ、Rで読み込めなかったり、その後のデータ操作でエラーが出ることがある。以下のフォーマットルールに従って、実験データを整理したい。
- 一番上の行を列ラベル名とし、2行目から参加者のデータを記録する。
- 基本的には、1行あたり1人分のデータを記録する。したがって、要因構造は以下のようになる。
- 参加者間要因：縦方向に分割する。
- 参加者内要因：横方向に分割する。
- すべてのデータは以下のように命名する。
- 半角小文字の英数字にする。
- 文字型データの頭文字は数字にしない。
- スペースを含めない。単語を区切りたい場合は半角アンダーバーを使う。
- 要因内で、各水準をアルファベット順に並べる。
3.4 Tidyverseについて
Rでデータサイエンスを行う際には、Tidyverse（たいでぃばーす）パッケージをいかに活用できるかが重要となる。
Tidyverseパッケージとは、データの読み込みから視覚化までの一連のデータ処理をスムーズに行うためのR言語のパッケージ群である。Tidyverseが内包するパッケージのうち、とくに本資料において重要となるのは以下の4つのパッケージである。
- ggplot2: グラフを作図するためのパッケージ
- dplyr: データを操作、抽出、加工したりするためのパッケージ
- tidyr: 「tidy-data」を作成したりするためのパッケージ
- magrittr: パイプ演算子「%>%」の機能を提供するパッケージ
3.5 エラーを恐れない
コードを書いていると、正常に動作することよりもエラーが出ることのほうが圧倒的に多い。エラーが出ても、すぐに諦める必要はない。
エラーメッセージ（英文）を読めばエラーの原因が分かるし、自分ではわからなくてもGoogle検索にかければ世界中のだれかが解決策を掲載してくれている。さらに、ChatGPTなどのAIサービスに自身が書いたコードとそれによって生じたエラーメッセージを貼り付ければ、ほとんどの場合で的確な解決策を提示してくれる。
4 事例1｜無相関検定
4.1 データ読み込み
事例1では、気温とアイスクリーム販売個数に関するデータを扱う。
read.xlsx関数で、作業ディレクトリに保存されているExcelデータの「data1.xlsx」を、「df1」という名前のデータフレームに格納する。なお、データフレームとは、行列データを格納できる変数のことである。
# 事例1 ####
df1 = read.xlsx("data1.xlsx")変数名をそのまま実行すると、その変数の中身が表示される。
df1## temperature number_of_ice_creams_sold
## 1 22 23
## 2 26 30
## 3 32 40
## 4 28 39
## 5 35 50
## 6 27 28
## 7 30 20
## 8 38 50
## 9 37 34
## 10 13 241列目に気温、2列目にアイスクリーム販売個数が格納されている、10行2列のデータフレームであることがわかる。
なお、データフレーム名のあとにドルマーク（$）を挟んで列名を指定すると、その列だけをベクトルとして抽出できる。たとえば、データフレームdf1から気温データを抽出したい場合は、下記のように書く。
df1$temperature## [1] 22 26 32 28 35 27 30 38 37 13ベクトルを引数として関数による計算ができる。例えば、平均値を算出する際は、mean関数を用いる。下記コードは、気温データの平均値を算出している。
mean(df1$temperature)## [1] 28.8平均値だけでなく最小値・最大値・中央値なども知りたければ、summary関数で一括取得できる。
summary(df1$temperature)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 13.00 26.25 29.00 28.80 34.25 38.004.2 グラフ作成
4.2.1 ヒストグラム
仮説検定を行う際は、想定する統計モデルを決定するために、データの分布構造を確認することが重要である（より厳密には、データを取得する前に統計モデルを決定すべきであり、データ取得後にデータを見て統計モデルを決める行為は統計不正とみなされる恐れがある）。とくに、これから行う仮説検定では、使用条件として正規性（データが正規分布していること）を仮定しているものが多いため、正規性を満たすかどうかをチェックしたい（より厳密には、各変数の正規性ではなく残差の正規性が求められているのだが、ここでは説明を省略する）。
データの分布構造を知りたい場合、まずはヒストグラムを確認したい。下図は、気温データの分布を示したものである。
このグラフを自分の手で作ってみる。Rでのグラフの作図には、ggplot2（じーじーぷろっとつー）と呼ばれる非常に強力なパッケージを使用する。ggplot2では、プラスマークを用いてコードをつなげることで、次々と要素や設定を足していくことができる。
まず、ggplot関数内で、使用するデータフレーム（df1）を指定し、aes関数内でx軸に配置するデータの列ラベル名を指定する。今回は、x軸が気温、y軸が指定した範囲に該当する気温のデータ数となる。y軸のデータ数はソフトウェア側が自動的に計算してくれるため、今回はx軸の気温データのみを指定すればよい。
ggplot(df1, aes(x=temperature))
ggplot関数は、使用するデータを指定して、x軸になにを配置するかを決める関数なので、グラフを描くためのキャンバスを用意しただけである。ここから、具体的にどのようなグラフを描くかを指定しなければならない。今回は、ヒストグラムを描画するgeom_histogram（じおむ_ひすとぐらむ）関数を、先ほどのコードにプラスマークで連結する。なお、プラスマークでコードを連結するときは、見やすさのためにその都度改行する。
ggplot(df1, aes(x=temperature))+
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
このままだと、気温1℃ごとにデータ数が計算されてしまっており、さらに色のせいで隣接する棒の区別がつきにくい。そこで、ヒストグラムの描き方を具体的に指定するために、geom_histogram関数の丸括弧内にオプションを追加する。
今回は、x軸上のどの範囲で集計するかをbreaksオプションで決める。さらに、fillオプションで棒の塗りつぶしを灰色にし、colourオプションで棒の枠線を黒色に指定している。ggplot2では、fillが塗りつぶしの色、colourが線の色を意味する。
ggplot(df1, aes(x=temperature))+
geom_histogram(breaks=c(18,20,22,24,26,28,30,32,34,36,38,40), fill="gray", colour="black")
なお、breaksオプションの中にc(~,~,~,...)というc関数（combineに由来）が使われているが、これはベクトルを作るための関数である。集計範囲を決めるbreaksオプションに、ベクトルを渡しているわけである。
breaksオプションの指定方法は、今回のような手動指定の他にも色々とやり方がある（今回は触れない）。また、x軸の目盛は別章で扱うscale_x_continuous関数で調整可能である。
今回はデータ数が少ないため正規性の有無は判断しづらいが、教材の便宜上、正規性が認められるものと仮定して次に進む。
4.2.2 カーネル密度推定図
近年は、ヒストグラムの代わりに、カーネル密度推定図を使用することが多くなった。ggplot2パッケージにはすでにカーネル密度推定を描くためのgeom_density関数が存在するが、x軸上にも線が描かれてしまうなどの問題があるため、本教材ではggridgesパッケージのgeom_density_line関数を用いる。
曲線の滑らかさは、adjustオプションで変更可能できる。1.0を基準として、1.0より大きくするとより滑らかになり、1.0より小さくするとより細部が強調される。ぜひ色々な数値を試してみてほしい。
ggplot(df1, aes(x=temperature))+
geom_density_line(fill="skyblue", colour="black", adjust=0.5)
カーネル密度推定図は、サンプルサイズが大きい場合に適した可視化法であり、今回のようにサンプルサイズが少ないときに使用すると、実態から乖離した図になる懸念があるため、注意が必要である。
4.2.3 散布図
事例1では、気温とアイスクリーム販売個数の相関関係を調べる。有意な相関があるかどうかを調べる前に、散布図を描画したい。
ここでは、x軸が気温、y軸がアイスクリーム販売個数となる散布図を作る。なお、散布図では、x軸が原因、y軸が結果となるように変数を配置するのが原則である。今回は「気温が高くなる→アイスクリーム販売個数が増加する」という関係性がイメージしやすいため、x軸に気温を設定している。
ggplot関数内で、使用するデータフレーム（df1）、x軸のデータ（気温）、y軸のデータ（アイスクリーム販売個数）をそれぞれ指定する。
その後、点を描画するgeom_point関数をプラスマークで連結する。
ggplot(df1, aes(x=temperature, y=number_of_ice_creams_sold))+
geom_point()
散布図に回帰線を追加したい場合はgeom_smooth関数を追加し、直線回帰を指定するためにオプションとしてmethod="lm"を追加する。なお、グレーの範囲は95%信頼区間を示す。
なお、グラフデータは変数に格納することができる。以下のコードでは、グラフデータを「fig1_1」という変数に格納している。あとから同じグラフを何度も使う場合は、このように変数に格納しておくと便利である。
fig1_1 = ggplot(df1, aes(x=temperature, y=number_of_ice_creams_sold))+
geom_point()+
geom_smooth(method="lm")
fig1_1## `geom_smooth()` using formula = 'y ~ x'
作図したグラフを画像として保存する場合はggsave関数を用いる。実行すると、直前に生成されたグラフが作業ディレクトリに画像として保存される。
下記設定は、上の散布図を「fig1.png」というファイル名で、解像度300dpi、縦横1200pxで保存するという意味である。コードを実行したあとに、実際に作業ディレクトリを開いて、画像が正しく保存されているかを確かめるとよい。
ggsave("fig1.png", dpi=300, width=1200, height=1200, units="px", device="png")## `geom_smooth()` using formula = 'y ~ x'4.3 無相関検定
4.3.1 知りたいこと
気温とアイスクリーム販売個数が相関するか（暑い日ほどアイスクリームがよく売れているか）を調べたい。
4.3.2 記述統計量
記述統計量として、平均値と不偏標準偏差を算出する。これらの記述統計量は論文中で報告する。平均値はmean関数、不偏標準偏差はsd関数で算出できる。
mean(df1$temperature)## [1] 28.8sd(df1$temperature)## [1] 7.524774mean(df1$number_of_ice_creams_sold)## [1] 33.8sd(df1$number_of_ice_creams_sold)## [1] 10.757944.3.3 検定の実施
| 項目 | 説明 |
|---|---|
| 無相関検定とは | 2変数間の母相関係数（母集団における相関係数）が0かどうかを検定する。 |
| 帰無仮説 | 気温とアイスクリーム販売個数の母相関係数は0である。 |
| 対立仮説 | 気温とアイスクリーム販売個数の母相関係数は0ではない。 |
| 使う関数 | cor.test(変数1, 変数2, method=“手法”) ※今回はPearsonの方法を用いる。 |
| 結果の解釈：p < .05 | 帰無仮説が棄却されて対立仮説が採択され、母相関係数が有意に0を上回る（あるいは下回る）と判断される。 |
| 結果の解釈：p ≧ .05 | 帰無仮説は棄却されず、母相関係数が0ではないとはいえない（0かどうかの判定を保留する）。 |
| 報告する内容 | 記述統計量、検定手法、効果の方向性（正の相関／負の相関）、相関係数r（cor）、検定統計量t、自由度df、p値、95%信頼区間Confidence Interval |
cor.test(df1$temperature, df1$number_of_ice_creams_sold, method="pearson")##
## Pearson's product-moment correlation
##
## data: df1$temperature and df1$number_of_ice_creams_sold
## t = 2.6648, df = 8, p-value = 0.02859
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.09874335 0.91870485
## sample estimates:
## cor
## 0.6857364検定結果は、（r = 相関係数,
t(自由度) = t値, p =
p値, 95% CI [CI下限値,
CI上限値]）のように書く。
4.3.4 結果の報告例
気温（M = 28.8, SD = 7.5）とアイスクリーム販売個数（M = 33.8, SD = 10.8）の関連を調べるために、Pearsonの相関係数を算出して無相関検定を行ったところ、有意な正の相関が見られ、気温が高い日ほどアイスクリーム販売個数が多いことが明らかとなった（r = .69, t(8) = 2.66, p = .029, 95% CI [0.10, 0.92]）（Figure 00）。
Figure 00
気温とアイスクリーム販売個数の関連
【APAマニュアルに基づく検定結果のフォーマットルール】
- 左端と右端の丸括弧は全角で、丸括弧の内側はすべて半角英数字である。
- p値は小数点以下3桁、それ以外は小数点以下2桁まで表記する（記述統計量は指標によって調整する）。
- p値を不等号で示すのは、.001を下回る場合のみである（p < .001）。
- M、SD、t、r、p、F、dなどはイタリック体にする。ηやεなどのギリシャ文字や、CIはイタリック体にしない。
- rやpなど、絶対値が必ず1以下となる数値においては、小数点の左の0を省略する。
- 詳細は「新美（2023）心理学論文における数値と統計の書き方」を参照。
4.4 押さえておきたいポイント
4.4.1 相関関係と因果関係
相関係数は、あくまでも2変数間の関連性（相関関係）を示すものであって、因果関係を説明するものではない点に注意が必要である。今回の結果の解釈は以下のようになる。
- 【正しい】気温が高い日ほど、アイスクリーム販売個数も多い。
- 【間違い】気温の高さが、アイスクリーム販売個数に影響を与える。
研究法には、大きく実験法と調査法がある。実験とは、独立変数を操作して、その操作によって従属変数が変わるかを調べる手法である。調査とは、研究者が変数を操作せずに、現状を記録して分析する手法である。調査は、介入操作を行わないため因果関係の特定が難しい。
相関分析は、主に調査によって得られた指標に対して行う。因果関係を特定するためには、操作を行う群（実験群）と操作を行わない群（統制群）の差の有無を調べる実験的手法が適している。
4.4.2 効果量r
r値は、相関係数であり効果量でもある。効果量の目安は下記の通りである（水本・竹内、2008、研究論文における効果量の報告のために）。ただし、これはあくまでも目安であって絶対的な指標ではないため、自身がデータの傾向を把握するためだけに利用を留め、客観的考察の根拠にはしないほうがよい。
- 0.50：効果量-大
- 0.30：効果量-中
- 0.10：効果量-小
p値が小さいからといって、効果が大きいとは限らないため、効果量はp値と分離して考えなければならない。p値はサンプルサイズに依存するが、効果量はサンプルサイズに依存しない。たとえp値が.000001だとしても、相関係数rが.10程度であれば、議論するのが難しいほど弱い相関となる。結果の解釈の際には、p値による二分法的判断のみを頼りにするのではなく、効果の大きさの指標である効果量も十分に考慮する必要がある。
4.4.3 散布図を必ず確認
散布図を見ると、相関分析の手法が適切かがわかったり、相関の適切な解釈が可能になるため、相関分析を行った際は必ず散布図を描画する必要がある。
たとえば、下図のように標本グループが明らかに分離しているケースでも有意な相関は見られてしまうが、この相関は無意味である。このケースの場合、有意な相関が見られるのは小学5年生と小学6年生に明らかな能力差があるためであって、数学試験の得点と体力測定の得点が本当に相関するかは不明である。つまり、学年という剰余変数が両変数に影響しており、交絡（適切に相関・因果関係を検証できていない状態のこと）を引き起こしている。
このような層別データを扱う場合は全体の相関を見てはならず、グループごとに相関分析を行ったり、偏相関分析、一般化線形混合モデル、マルチレベルモデル、などの別の分析手法を用いたりする必要がある。
4.4.4 検定の使用条件を満たすかを確認
Pearsonの相関係数は、2つの変数間に線形関係があり、残差の正規性と等分散性が仮定できる場合に使える手法である。
正規性を仮定できない場合は、Pearsonの方法は適切ではない可能性がある。サンプルサイズが十分に大きい場合には、中心極限定理に基づき、正規性が仮定されなくてもPearsonの方法を使うことができる。必要サンプルサイズは、実験調査デザインや研究目的ならびに想定する統計モデルによって異なるものの、多くの文献では目安として30以上とされている（絶対的基準はない）。しかし、サンプルサイズが小さい場合や、外れ値の影響が強く懸念されるような場合は、別の方法を考える必要がある。
Pearsonの方法以外には、順位化したデータに基づく順位相関係数を算出するSpearmanの方法を利用できる。オプションにmethod="spearman"を追加すると、Spearmanの順位相関係数に基づく無相関検定を実施できる。
cor.test(df1$temperature, df1$number_of_ice_creams_sold, method="spearman")## Warning in cor.test.default(df1$temperature, df1$number_of_ice_creams_sold, :
## Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: df1$temperature and df1$number_of_ice_creams_sold
## S = 49.65, p-value = 0.02447
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6990914なお、Pearsonの相関係数の検定のように母集団の分布の影響を受ける検定手法のことをパラメトリック検定と呼ぶ。一方で、Spearmanの順位相関係数の検定のように母集団の分布の影響を受けない（より正確には影響を受けにくい）検定手法はノンパラメトリック検定と呼ばれる。
注意点として、Spearmanの相関係数のように元のデータを順位データに変換する手法は、変換過程で情報が大きく損失してしまうデメリットがある。外れ値の影響が強く懸念される場合などを除いて、順位データに変換する検定手法は極力避けたい。
4.4.5 検定の多重性 Part 1
目の前のデータが、PearsonとSpearmanのどちらの方法が適切かを判断したいケースは多いだろう。そのような場合に、無相関検定を行う前に、正規性の検定を別途行うことを推奨する書籍・WEBサイトがある。たとえば、統計ソフトウェアSPSSの公式サポートにも、正規性の検定が有意でない場合はPeasonの方法、正規性の検定が有意の場合はSpearmanの方法を使うように指示がある。
しかし、メインの検定の前に正規性検定や等分散性検定を行い、その事前検定結果に基づいてメインの検定の手法を選択する方法には、主に下記2つの問題がある。
- 正規性検定では、正規性がない（対立仮説）ことはわかっても、正規性がある（帰無仮説）ことはわからない。正規性検定でp値が0.05以上となっても、正規性があるということは主張できないため、正規性検定は無意味である。等分散性の検定も同様である。
- 同一データに対して仮説検定を繰り返すと、事前に設定したアルファエラー率とベータエラー率を維持できなくなる。
これに対処する方法はいくつかあるが、重要な指針として、仮説検定の際には正規性や等分散性などの検定は行わずに、ヒストグラムなどでデータの分布構造を確認する程度に留めた方がよい。もっと言えば、データを取得してから分布構造をヒストグラムで確認するのではなく、データを取る前からどのような分布構造になるかを想定できるように研究デザインを組み立てることが求められる。
4.4.6 外れ値を含むデータの解析
外れ値を含むデータに対して検定を行う際には、下記のような方法がよく用いられる。
- 平均値から大きく逸脱した値を外れ値として除外する。
- よく使われる基準は「平均値±2SD」あるいは「平均値±3SD」。標準正規分布において、データの約95%が平均値±2SDの範囲内に収まり、データの約99.7%が平均値±3SDの範囲内に収まる。
- たとえば、平均値10で標準偏差3のデータ群において、「平均値±2SD」を基準に除外する場合、（10-3×2）〜（10+3×2）の範囲、つまり4〜16の範囲に入らないデータを外れ値として扱い除外する。
- 順位化する。
- 順位化をイメージしやすいように、図で順位化前と順位化後のデータを比較してみる。左が順位化前のデータ、右が順位化後のデータである。順位化は
rank関数で容易にできる。
- 順位化をイメージしやすいように、図で順位化前と順位化後のデータを比較してみる。左が順位化前のデータ、右が順位化後のデータである。順位化は
ただし、外れ値が単なるエラーではなく重要な情報を含んでいる可能性は、常に考えておかなければならない。また、外れ値の除外基準をころころと変えながら分析を繰り返してはならない。データを取得する前に、除外基準を明瞭に定めておくべきである。
4.4.7 理解度確認クイズ
年収と血圧の関連を調べるために、20名分のデータを用いてSpearmanの相関係数の無相関検定を実施したところ、r = .50、p = .03と出力された。有意水準は5%とする。以下の記述のうち、下線部が正しいものをすべてチェックしなさい。
正解を表示する
正解は【2】
4.5 練習問題1
注意課題の反応時間と自閉症スペクトラム指数（AQ-J-10質問紙に基づく）の関連を調べてみましょう。注意課題の反応時間が短いほど注意能力が高く、自閉症スペクトラム指数が高いほど自閉症傾向が高いと解釈されます。
※AQ-J-10に関する詳細：Kurita, H., Koyama, T., & Osada, H. (2005). Autism‐Spectrum Quotient–Japanese version and its short forms for screening normally intelligent persons with pervasive developmental disorders. Psychiatry and clinical neurosciences, 59(4), 490-496.
STEP1〜5の手順を参考に解析を行い、下記フォームからグラフと解析結果を解答してみてください。大学アカウントでログインしているとGoogleフォームにアクセスできませんので、自身のGoogleアカウントでログインしてください。
STEP1
まず、作業ディレクトリに保存されているExcelデータの「data1p.xlsx」を、「df1p」という名前のデータフレームに格納します。
1列目（attention_task_reaction_time）に注意課題の反応時間の平均値、2列目（aq_score）に自閉症スペクトラム指数の得点が格納されている、20行2列のデータです。
## attention_task_reaction_time aq_score
## 1 334.6619 6
## 2 317.3599 3
## 3 319.2476 10
## 4 365.7743 7
## 5 271.2179 4
## 6 342.0688 6
## 7 339.9791 9
## 8 306.7227 2
## 9 300.4734 6
## 10 379.7396 5
## 11 349.5109 5
## 12 342.1110 7
## 13 328.3402 4
## 14 331.6492 3
## 15 329.1320 7
## 16 334.4032 7
## 17 311.9025 2
## 18 338.2931 5
## 19 317.1956 4
## 20 344.3361 1STEP2
続いて、以下のグラフと同じになるように、ggplot関数でグラフを作図してみましょう。
なお、出力されたグラフの大きさ（＝画面の大きさ）によって軸目盛りが変わるかもしれませんが、それは問題ありません。
STEP3
グラフが作図できたらggsave関数で図を保存しましょう。
保存したいグラフを作図した直後に下記コードを実行すると、「fig1p.png」というファイル名で画像データが作業ディレクトリに保存されます。
ここで作成したグラフを、Googleフォームにアップロードしてください。
ggsave("fig1p.png", dpi=300, width=1200, height=1200, units="px", device="png")STEP4
記述統計量として、変数ごとに平均値と不偏標準偏差を算出しましょう。
平均値はmean関数、不偏標準偏差はsd関数で算出できます。
STEP5
Pearsonの相関係数を用いた無相関検定を行いましょう。
検定結果の報告文章を、Googleフォームに入力してください。
5 事例2｜t検定
5.1 データ読み込み
事例2では、漢字テストの得点に関するデータを扱う。学習前と学習後の2種類のデータを用意した。
作業ディレクトリに保存されているExcelデータの「data2.xlsx」を、「df2」という名前のデータフレームに格納する。
# 事例2 ####
df2 = read.xlsx("data2.xlsx")
df2## sex score_before score_after
## 1 f 3 5
## 2 f 6 9
## 3 f 5 9
## 4 f 6 8
## 5 f 5 10
## 6 m 5 7
## 7 m 4 6
## 8 m 5 8
## 9 m 4 7
## 10 m 3 6
## 11 m 4 6
## 12 other 3 81列目に性別、2列目に学習前の得点、3列目に学習後の得点が格納されている、12行3列のデータフレームであることがわかる。
count関数を使えば、各列の質的カテゴリーデータの内訳を集計してくれる。たとえば、このデータの性別ごとの人数を集計したければ、下記のコードを実行する。
count(df2, sex)## sex n
## 1 f 5
## 2 m 6
## 3 other 15.2 データの前処理
学習によって漢字テストの得点が伸びた度合いを計算したい。得点の伸びは、学習後得点（score_after）から学習前得点（score_before）を減算すれば算出できる。
データの加工には、Tidyverseパッケージに内包されているdplyrパッケージやtidyrパッケージを用いる。
下記のコードは、学習後得点から学習前得点を減算した値をscore_diffという列に算出し、その列をデータフレームdf2に追加した上で、既存のデータフレームdf2を新しいデータフレームdf2に置き換える処理をしている。
%>%はパイプ演算子と呼ばれ、演算子の前のデータを演算子の後に引き渡す役割を担う。%>%の直後は、見やすさのために改行することが多い。
mutate関数は、新しい列を作成する関数である。イコールの左に作成したい列の列名を入力し、イコールの右にその計算式を書く。
df2 = df2 %>%
mutate(score_diff = score_after - score_before)
df2## sex score_before score_after score_diff
## 1 f 3 5 2
## 2 f 6 9 3
## 3 f 5 9 4
## 4 f 6 8 2
## 5 f 5 10 5
## 6 m 5 7 2
## 7 m 4 6 2
## 8 m 5 8 3
## 9 m 4 7 3
## 10 m 3 6 3
## 11 m 4 6 2
## 12 other 3 8 55.3 グラフ作成
5.3.1 ロング型データとは
心理実験データは、1行あたり1名の参加者のデータをまとめ、参加者内要因の水準が横方向に広がるワイド型データであることが多い。
しかし、ggplot関数で作図する際は、1つの列が1つの変数となるようなロング型データ（別名：Tidyデータ／整然データ）でなければならない。ワイド型データはロング型データに変換する必要がある。
以下に、ワイド型データとロング型データの変換例を挙げる。どちらも同じデータだが、データフレームの構造が異なる。
◆ワイド型データ
| 参加者ID | 年齢 | 1回目 | 2回目 | 3回目 |
|---|---|---|---|---|
| 1 | 18 | ~ | ~ | ~ |
| 2 | 20 | ~ | ~ | ~ |
| 3 | 22 | ~ | ~ | ~ |
◆ロング型データに変換した場合
| 参加者ID | 年齢 | 測定タイミング | スコア |
|---|---|---|---|
| 1 | 18 | 1回目 | ~ |
| 1 | 18 | 2回目 | ~ |
| 1 | 18 | 3回目 | ~ |
| 2 | 20 | 1回目 | ~ |
| 2 | 20 | 2回目 | ~ |
| 2 | 20 | 3回目 | ~ |
| 3 | 22 | 1回目 | ~ |
| 3 | 22 | 2回目 | ~ |
| 3 | 22 | 3回目 | ~ |
5.3.2 箱ひげ図
まずは、データフレームdf2をワイド型データからロング型データに変換する。
ロング型データへの変換には、ロング形式にピボット（回転変換）するという意味で、pivot_longer関数を使用する。丸括弧内では、下記3種類の値を指定する。
cols：ロング型データに変換したい参加者内要因の各水準の列名をすべて指定する。names_to：変換後の参加者内要因の列名を指定する。value_to：変換後の値の列名を指定する。
df2_plot = df2 %>%
pivot_longer(
cols = c(score_before, score_after),
names_to = "test_condition",
values_to = "score")
df2_plot## # A tibble: 24 × 4
## sex score_diff test_condition score
## <chr> <dbl> <chr> <dbl>
## 1 f 2 score_before 3
## 2 f 2 score_after 5
## 3 f 3 score_before 6
## 4 f 3 score_after 9
## 5 f 4 score_before 5
## 6 f 4 score_after 9
## 7 f 2 score_before 6
## 8 f 2 score_after 8
## 9 f 5 score_before 5
## 10 f 5 score_after 10
## # ℹ 14 more rowsggplot関数内で、使用するデータフレーム（df2_plot）、x軸のデータ（学習前か学習後か）、y軸のデータ（得点）をそれぞれ指定する。
その後、箱ひげ図を描画するgeom_boxplot関数をプラスマークで連結する。
labs関数では、x軸とy軸の軸ラベル名を指定できる。
ggplot(df2_plot, aes(x=test_condition, y=score))+
geom_boxplot()+
labs(x="Test Condition", y="Test Score")
現状では、左に学習後、右に学習前のデータが配置されてしまっていて直感的に分かりにくい。ggplotにおける配置順は、ファクターレベル（処理の順番を決める内部変数）が事前に設定されていない場合、自動的にアルファベット順となってしまう。
beforeとafterの位置を逆転させるために、scale_x_discrete関数を用いる。この関数は名前の通り、「離散値（discrete）であるx軸」について設定する関数である。limitsオプションで並べる順番を決定し、labelsオプションで目盛りの項目名を変更できる。
今回は、左に学習前のデータが配置されるようにし、学習前をbefore、学習後をafterと命名する。
ggplot(df2_plot, aes(x=test_condition, y=score))+
geom_boxplot()+
labs(x="Test Condition", y="Test Score")+
scale_x_discrete(limits=c("score_before","score_after"), labels=c("before","after"))
箱ひげ図では、データの分布をより詳細に可視化するために、箱とひげに加えて、データ点をプロットすることが多い。そこで、geom_jitter関数でデータ点をプロットする。
geom_jitter関数では、データ点の描画位置をランダムにばらつかせてくれる。同じ値のデータがあると重なって見えなくなってしまうが、ばらつかせることですべての点を可視化できる。縦方向にばらつかせてしまうと、計測値からずれてしまうため、ばらつかせるのは横方向のみとする。widthオプションで横方向のばらつかせ具合を調整し、heightオプションは0と指定する。さらに、点の重なりがわかるように、点の透過度をalphaオプションで50%に指定する。
fig2_1 = ggplot(df2_plot, aes(x=test_condition, y=score))+
geom_boxplot()+
geom_jitter(width=0.2, height=0, alpha=0.5)+
labs(x="Test Condition", y="Test Score")+
scale_x_discrete(limits=c("score_before","score_after"), labels=c("before","after"))
fig2_1
5.3.3 ドットプロット
ドットプロットでは、平均値と不偏標準誤差（SE: Standard Error）を描画することとなる。不偏標準偏差ではないので注意すること。したがって、今回のデータを可視化するためには、下表のような2行3列の集計データフレームを作成する必要がある。
| test_condition | mean | se |
|---|---|---|
| score_after | 7.42 | 0.43 |
| score_before | 4.42 | 0.31 |
集計データフレームを作成するために、まずは先ほどと同様に、pivot_longer関数でscore_beforeとscore_afterのワイド型データをロング型データに変換する。
df2_plot = df2 %>%
pivot_longer(
cols = c(score_before, score_after),
names_to = "test_condition",
values_to = "score")
df2_plot## # A tibble: 24 × 4
## sex score_diff test_condition score
## <chr> <dbl> <chr> <dbl>
## 1 f 2 score_before 3
## 2 f 2 score_after 5
## 3 f 3 score_before 6
## 4 f 3 score_after 9
## 5 f 4 score_before 5
## 6 f 4 score_after 9
## 7 f 2 score_before 6
## 8 f 2 score_after 8
## 9 f 5 score_before 5
## 10 f 5 score_after 10
## # ℹ 14 more rows続いて、平均値と不偏標準誤差の集計をするために、group_by関数とsummarize関数を組み合わせる。まず、group_by関数でtest_conditionに基づいてグルーピングし、summarize関数でそれぞれのグループにおける平均値と不偏標準誤差を集計する。平均値はmean、不偏標準誤差はseという列名を指定する。
df2_plot = df2_plot %>%
group_by(test_condition) %>%
summarize(mean = mean(score),
se = std.error(score))
df2_plot## # A tibble: 2 × 3
## test_condition mean se
## <chr> <dbl> <dbl>
## 1 score_after 7.42 0.434
## 2 score_before 4.42 0.313ggplot関数内で、使用するデータフレーム（df2_plot）、x軸のデータ（学習前か学習後か）、y軸のデータ（得点）をそれぞれ指定する。
ドットプロットを描画する、geom_point関数を追加する。
さらに、不偏標準誤差をエラーバーとして表示するために、geom_errorbar関数を追加する。丸括弧内の中身は、aes関数ではエラーバーの縦の長さ、widthオプションではエラーバーの横幅を指定している。
ggplot(df2_plot, aes(x=test_condition, y=mean))+
geom_point()+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
labs(x="Test Condition", y="Test Score")+
scale_x_discrete(limits=c("score_before","score_after"), labels=c("before","after"))
仮にこれら2つの間に統計的な有意差が見られた場合を想定して、統計的有意性を示すアスタリスクをつけることを考えたい。
現状の図では、アスタリスクをつける領域がないことから、y軸の表示範囲を調整して、グラフ中の上側に余白を作りたい。
y軸の調整は、scale_y_continuous関数を用いる。この関数は名前の通り、「連続値（continuous）であるy軸」について設定する関数である。
limitsオプションでは、y軸を表示する範囲を設定する。下記コードでは、下限3.5から上限9.5の範囲だけy軸を表示するようにしている。
breaksオプションでは、目盛りの表示の仕方を設定する。下記コードでは4から9までの区間を1で区切るように目盛りをつけるようにしている。つまり、4、5、6、7、8、9の箇所に目盛りを描画する。
ggplot(df2_plot, aes(x=test_condition, y=mean))+
geom_point()+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
labs(x="Test Condition", y="Test Score")+
scale_x_discrete(limits=c("score_before","score_after"), labels=c("before","after"))+
scale_y_continuous(limits=c(3.5,9.5), breaks=seq(4,9,1))
このグラフに、統計的有意性を示すアスタリスクを描画する。任意の要素をグラフ中に描画するには、anotate関数を用いる。
アスタリスクの描画のためのコード（“text”）と、横線の描画のためのコード（“segment”）の2つに分かれている。xオプションはx軸のどこに描画するか、yオプションはy軸のどこに描画するかを示す。注意点として、線を描画する際は、x軸とy軸それぞれで始点と終点を指定しなければならない。labelオプションには描画したいテキスト、sizeオプションにはテキストの大きさを入力する。
fig2_2 = ggplot(df2_plot, aes(x=test_condition, y=mean))+
geom_point()+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
labs(x="Test Condition", y="Test Score")+
scale_x_discrete(limits=c("score_before","score_after"), labels=c("before","after"))+
scale_y_continuous(limits=c(3.5,9.5), breaks=seq(4,9,1))+
annotate("text", x=1.5, y=9, label="***", size=5)+
annotate("segment", x=1, xend=2, y=9, yend=9)
fig2_2
5.3.4 棒グラフ
別のグラフとして、漢字テストの得点の伸びの程度を男女別に示した棒グラフも作図したい。その場合は、下表のような2行3列の集計データフレームを作成する必要がある。
| sex | mean | se |
|---|---|---|
| f | 3.2 | 0.58 |
| m | 2.5 | 0.22 |
先ほどと同様の手順で集計データを作成する。なお、今回はscore_beforeとscore_afterのデータは使用せず、score_diffのデータを使用するため、ワイド型データからロング型データに変換する必要はない。
今回は、出生時に割り当てられた性別をx軸に配置するため、自身の性別を「other」と回答した参加者のデータを除外したい。条件に当てはまる行を抽出する場合は、filter関数を用いる。下記コードは、sex列の値がmあるいはfの行のみを抽出するという意味である。なお、「|」はプログラミングにおいて「または」を意味する。
先ほどと同様に、group_by関数で性別に基づいてグルーピングし、summarize関数で性別ごとの平均値と不偏標準誤差を集計する。
df2_plot = df2 %>%
filter(sex=="m" | sex=="f") %>%
group_by(sex) %>%
summarize(mean = mean(score_diff),
se = std.error(score_diff))
df2_plot## # A tibble: 2 × 3
## sex mean se
## <chr> <dbl> <dbl>
## 1 f 3.2 0.583
## 2 m 2.5 0.224ggplot関数内で、使用するデータフレーム（df2_plot）、x軸のデータ（性別）、y軸のデータ（得点）をそれぞれ指定する。
棒グラフはgeom_col関数で描画できる。先ほどと同様に、geom_errorbar関数も追加する。
男女の順番についてはジェンダーバイアスが表れるが、データ数が多い方を先とすれば悩まないだろう。
fig2_3 = ggplot(df2_plot, aes(x=sex, y=mean))+
geom_col()+
geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0.1)+
labs(x="Sex", y="Test Score Difference")+
scale_x_discrete(limits=c("m","f"), labels=c("male","female"))
fig2_3
5.3.5 グラフの使い分け
これまでに作成した箱ひげ図、ドットプロット、棒グラフは、いずれもx軸に離散値、y軸に連続値を設定する似たようなデータ可視化法である。棒グラフに慣れ親しんでいる人が多いと考えられるが、箱ひげ図やドットプロットのほうが適しているケースも多く、データ構造に合った可視化法を採用する必要がある。それぞれの特徴は以下の通りである。
- 箱ひげ図：データの分布構造を分かりやすく示すことができる。近年のデータサイエンスでは、記述統計量である平均値と標準誤差を示すドットプロットや棒グラフよりも、箱ひげ図のようなデータの分布構造が見える可視化法のほうが重視されている。
- ドットプロット：棒グラフと同じ情報を示しているが、間隔尺度（0を絶対的基準としない）を示すときによく使われる。たとえば、心理測定尺度の値などはドットプロットが適する場合が多い。
- 棒グラフ：ドットプロットと同じ情報を示しているが、比例尺度（0を絶対的基準とする）を示すときによく使われる。たとえば、課題の正答数などは棒グラフが適する場合が多い。
5.4 2つの値の比較手法の選び方
データAとデータBの差が有意かどうかを調べるときに、比較手法を選択しなければならない。目的はシンプルだが手法は多岐にわたる。最適な手法の選択は研究者の中でも意見が分かれることがあり、統計専門書でも記載内容にばらつきがある（統計学の専門家による啓蒙の効果がなかなか表れていない）。
2つの値を比較する場面は主に2種類である。それぞれの場面における手法を確認する。
- 場面例1：「A組の5月と6月のテストの得点に差はあるか？」
- 1つの標本から2つのデータ群を抽出しており、このような反復測定データに対しては「対応のあるt検定」を用いる。
- ほとんどの場合で、差データの正規性と分布対称性を仮定する「対応のあるスチューデントのt検定」を用いて平均値の比較を行う。正規性を仮定しているためパラメトリック検定である。本やネットで対応のあるt検定と見かけたら、とくに明記されていなくてもスチューデントの方法を使った検定である。
- 正規性が仮定できない場合は、一般的にノンパラメトリック検定であるウィルコクソンの符号順位検定を用いて中央値を比較する、と説明する文献が多い。しかし、ウィルコクソンの符号順位検定も一切の仮定を持たないわけではなく分布対称性を仮定している。現実世界において、分布対称性を完璧に仮定できる状況は極めて少ないため、ウィルコクソンの符号順位検定は使いにくい。さらに、分布対称性が仮定されたとしても不適切な結果を出力することがあるため、使用を避けるべきである。
- そもそも、中心極限定理に基づき、サンプルサイズが十分に大きい場合には、正規性が仮定されなくてもパラメトリック検定を使うことができる場合がある。必要サンプルサイズは、実験調査デザインや研究目的ならびに想定する統計モデルによって異なるものの、多くの文献では目安として30以上とされている。しかし、サンプルサイズが小さい場合は、別の方法を考える必要がでてくる。
- 近年では、t検定が使えない場合に、対応のあるBrunner-Munzel（ブルンナー=ムンツェル）検定が有用であるとする指摘がある。対応のあるBrunner-Munzel検定は、正規性と分布対称性を仮定しない、確率的等質性に関する検定である。サンプルサイズが小さいときでも、対応のある並べ替えBrunner-Munzel検定が使える。ただし、t検定と根本的に原理が異なり、t検定とは異なる指標の比較を行っているため、研究目的に基づいて最適な手法を選択すべきである。
- 場面例2：「A組とB組のテストの得点に差はあるか？」
- 2つの標本からそれぞれ1つのデータ群を抽出しており、このような群間データに対しては「対応のないt検定（2標本のt検定）」を用いる。各群のデータが対応しておらず独立しているから、このような名前となっている。
- 残差の正規性と等分散性が仮定できる場合は、「対応のないスチューデントのt検定」を用いて平均値の比較を行う。多くの書籍やネットには、事前に等分散性の検定を行い、等分散性が確認できたらスチューデントのt検定、等分散性が確認できなかったら代わりに「welchのt検定」を用いると記載がある。
- しかし、等分散性を完璧に仮定できる状況は極めて少なく、welchのt検定は等分散でも非等分散でも使える万能な手法であることから、事前の等分散性検定による検定の多重性問題を避けるためにも、常にwelchのt検定を用いることが強く推奨される。なお、Rの「t.test」関数は初期状態でwelchのt検定を実施するようになっている。つまり、スチューデントのt検定を使うことはない。
- ただし、welchのt検定は正規性を仮定するパラメトリック検定である。サンプルサイズが十分に大きいときは、正規性が仮定できなくても中心極限定理に基づき、使うことができる場合がある。必要サンプルサイズは、実験調査デザインや研究目的ならびに想定する統計モデルによって異なるものの、多くの文献では目安として30以上とされている。しかし、サンプルサイズが小さい場合は、別の方法を考える必要がでてくる。
- そのような場合は、一般的にノンパラメトリック検定であるウィルコクソンの順位和検定（マンホイットニーのU検定と同義）を用いて順位平均を比較する、と説明する文献が多い。しかし、ウィルコクソンの順位和検定も一切の仮定を持たないわけではなく等分散性を仮定しており使いにくい。さらに、等分散性が仮定されたとしても不適切な結果を出力することがあるため、使用を避けるべきである。
- 近年では、ウィルコクソンの順位和検定の代わりに、Brunner-Munzel（ブルンナー=ムンツェル）検定が有用であるとする指摘がある。Brunner-Munzel検定は、正規性と等分散性を仮定しない、確率的等質性に関する検定である。サンプルサイズが小さいときでも、並べ替えBrunner-Munzel検定が使用できる。ただし、t検定と根本的に原理が異なり、t検定とは異なる指標の比較を行っているため、研究目的に基づいて最適な手法を選択すべきである。
- まとめ
- 対応のある場合：「対応のある（スチューデントの）t検定」を基本とする。差データの正規性や分布対称性が仮定されず、サンプルサイズが中心極限定理を適用できないほど小さく、さらに研究目的に合致すれば、「対応のあるBrunner-Munzel検定」あるいは「対応のある並べ替えBrunner-Munzel検定」を使う。ウィルコクソンの符号順位検定は基本的に使わない。
- 対応のない場合：「welchのt検定」を基本とする。残差の正規性や等分散性が仮定されず、サンプルサイズが中心極限定理を適用できないほど小さく、さらに研究目的に合致すれば、「Brunner-Munzel検定」あるいは「並べ替えBrunner-Munzel検定」を使う。スチューデントのt検定やウィルコクソンの順位和検定（マンホイットニーのU検定）は基本的に使わない。
- 参考文献
- 井口「Wilcoxon符号順位検定の対称性条件と変数変換」
- 井口「t検定の正規性も残差を調べる：検定の多重性問題」
- 清水「分析のいろんな仮定と，それに対する頑健さ・対処法」
- 黒木「Twitter」午後11:21 · 2023年11月5日
- 奥村「Brunner-Munzel検定」
- hoxo_m「マイナーだけど最強の統計的検定 Brunner-Munzel 検定」
- 名取「マン・ホイットニーのU検定と不等分散時における代表値の検定法」
- ねこすたっと「対応のない2群の代表値を比較するときに用いるべき検定方法」
- Psychologists Should Use Brunner-Munzel’s Instead of Mann-Whitney’s U Test as the Default Nonparametric Procedure
5.5 対応のあるt検定
5.5.1 知りたいこと
学習前後で得点が変わったかを調べたい。
5.5.2 記述統計量
記述統計量として、平均値と不偏標準偏差を算出する。これらの記述統計量は論文中で報告する。
#学習前の得点の平均値
mean(df2$score_before)## [1] 4.416667#学習前の得点の不偏標準偏差
sd(df2$score_before)## [1] 1.083625#学習後の得点の平均値
mean(df2$score_after)## [1] 7.416667#学習後の得点の不偏標準偏差
sd(df2$score_after)## [1] 1.5050425.5.3 検定の実施
| 項目 | 説明 |
|---|---|
| 対応のあるt検定とは | 1つの標本による2つの母平均（母集団における平均値）が異なるかどうかを検定する。参加者内要因の反復測定データに対して使用する手法。 |
| 帰無仮説 | 学習前と学習後の漢字テストの得点の母平均に差はない。 |
| 対立仮説 | 学習前と学習後の漢字テストの得点の母平均に差がある。 |
| 使う関数1（t検定） | t.test(変数1, 変数2, paired=TRUE) ※左に大きい変数、右に小さい変数を配置する。 |
| 使う関数2（効果量） | cohens_d(変数1, 変数2, paired=TRUE) ※effectsizeパッケージに含まれる関数。左に大きい変数、右に小さい変数を配置する。 |
| 結果の解釈：p < .05 | 帰無仮説が棄却されて対立仮説が採択され、2つの母平均の差が有意に0を上回る（あるいは下回る）と判断される。 |
| 結果の解釈：p ≧ .05 | 帰無仮説は棄却されず、2つの母平均の差が0ではないとはいえない（0かどうかの判定を保留する）。 |
| 報告する内容 | 記述統計量、検定手法、効果の方向性、検定統計量t、自由度df、p値、効果量d、95%信頼区間Confidence Interval |
t.test(df2$score_after, df2$score_before, paired=TRUE)##
## Paired t-test
##
## data: df2$score_after and df2$score_before
## t = 9.2118, df = 11, p-value = 1.669e-06
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 2.283206 3.716794
## sample estimates:
## mean difference
## 3cohens_d(df2$score_after, df2$score_before, paired=TRUE) #効果量dを算出する関数## For paired samples, 'repeated_measures_d()' provides more options.## Cohen's d | 95% CI
## ------------------------
## 2.66 | [1.41, 3.88]5.5.4 結果の報告例
学習前の漢字テストの得点（M = 4.42, SD = 1.08）と学習後の漢字テストの得点（M = 7.42, SD = 1.51）に対して対応のあるt検定を実施したところ、学習後の得点のほうが有意に高いことが明らかとなった（t(11) = 9.21, p < .001, d = 2.66, 95% CI [2.28, 3.72]）（Figure 00）。
Figure 00
学習前後の漢字テストの得点
注）エラーバーは標準誤差を示す。***p < .001。
【APAマニュアルに基づく検定結果のフォーマットルール】
- 左端と右端の丸括弧は全角で、丸括弧の内側はすべて半角英数字である。
- p値は小数点以下3桁、それ以外は小数点以下2桁まで表記する（記述統計量は指標によって調整する）。
- p値を不等号で示すのは、.001を下回る場合のみである（p < .001）。
- M、SD、t、r、p、F、dなどはイタリック体にする。ηやεなどのギリシャ文字や、CIはイタリック体にしない。
- rやpなど、絶対値が必ず1以下となる数値においては、小数点の左の0を省略する。
- 詳細は「新美（2023）心理学論文における数値と統計の書き方」を参照。
5.6 対応のないt検定
5.6.1 知りたいこと
漢字テストの得点の学習前後の変化量に性差があるかを調べたい。
5.6.2 記述統計量
記述統計量の算出のために、男性のデータと女性のデータをそれぞれ抽出しなければならない。以下のコードは、男性のscore_diffのデータをdf2_m、女性のscore_diffのデータをdf2_fに格納している。
filter関数でsex列の性別を指定したあとに、pull関数でscore_diff列のデータのみを抽出している。pull関数は、選択した列のデータをベクトルとして抽出する関数である。
df2_m = df2 %>% filter(sex=="m") %>% pull(score_diff)
df2_f = df2 %>% filter(sex=="f") %>% pull(score_diff)記述統計量として、平均値と不偏標準偏差を算出する。これらの記述統計量は論文中で報告する。
#男性の平均値
mean(df2_m)## [1] 2.5#男性の不偏標準偏差
sd(df2_m)## [1] 0.5477226#女性の平均値
mean(df2_f)## [1] 3.2#女性の不偏標準偏差
sd(df2_f)## [1] 1.303845.6.3 検定の実施
| 項目 | 説明 |
|---|---|
| 対応のないt検定とは | 2つの標本におけるそれぞれの母平均（母集団における平均値）が異なるかどうかを検定する。参加者間要因のデータに対して使用する手法。 |
| 帰無仮説 | 学習前後の漢字テストの得点の変化量について、男性と女性の母平均に差はない。 |
| 対立仮説 | 学習前後の漢字テストの得点の変化量について、男性と女性の母平均に差がある。 |
| 使う関数1（t検定） | t.test(変数1, 変数2) ※welchの方法によるt検定が標準で行われる。左に大きい変数、右に小さい変数を配置する。 |
| 使う関数2（効果量） | cohens_d(変数1, 変数2) ※effectsizeパッケージに含まれる関数。左に大きい変数、右に小さい変数を配置する。 |
| 結果の解釈：p < .05 | 帰無仮説が棄却されて対立仮説が採択され、2つの母平均の差が有意に0を上回る（あるいは下回る）と判断される。 |
| 結果の解釈：p ≧ .05 | 帰無仮説は棄却されず、2つの母平均の差が0ではないとはいえない（0かどうかの判定を保留する）。 |
| 報告する内容 | 記述統計量、検定手法、効果の方向性、検定統計量t、自由度df、p値、効果量d、95%信頼区間Confidence Interval |
t.test(df2_f, df2_m)##
## Welch Two Sample t-test
##
## data: df2_f and df2_m
## t = 1.1209, df = 5.1735, p-value = 0.3116
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.889272 2.289272
## sample estimates:
## mean of x mean of y
## 3.2 2.5cohens_d(df2_f, df2_m) #効果量dの算出## Warning: 'y' is numeric but has only 2 unique values.
## If this is a grouping variable, convert it to a factor.## Cohen's d | 95% CI
## -------------------------
## 0.73 | [-0.52, 1.94]
##
## - Estimated using pooled SD.5.6.4 結果の報告例
漢字テストの得点の学習前後の変化量に性差があるかを調べるために、男性（M = 2.50, SD = 0.55）と女性（M = 3.20, SD = 1.30）のデータに対して対応のないwelchのt検定を実施したところ、有意な性差はみられなかった（t(5.17) = 1.12, p = .312, d = 0.73, 95% CI [-0.89, 2.29]）（Figure 00）。
Figure 00
性別ごとの漢字テストの得点の変化量
注）エラーバーは標準誤差を示す。
【APAマニュアルに基づく検定結果のフォーマットルール】
- 左端と右端の丸括弧は全角で、丸括弧の内側はすべて半角英数字である。
- p値は小数点以下3桁、それ以外は小数点以下2桁まで表記する（記述統計量は指標によって調整する）。
- p値を不等号で示すのは、.001を下回る場合のみである（p < .001）。
- M、SD、t、r、p、F、dなどはイタリック体にする。ηやεなどのギリシャ文字や、CIはイタリック体にしない。
- rやpなど、絶対値が必ず1以下となる数値においては、小数点の左の0を省略する。
- 詳細は「新美（2023）心理学論文における数値と統計の書き方」を参照。
5.7 押さえておきたいポイント
5.7.1 効果量d
t検定の効果量d（Cohen’s d）の大きさの目安は下記の通りである（水本・竹内、2008、研究論文における効果量の報告のために）。ただし、これはあくまでも目安であって絶対的な指標ではないため、自身がデータの傾向を把握するためだけに利用を留め、客観的考察の根拠にはしないほうがよい。
- 0.80：効果量-大
- 0.50：効果量-中
- 0.20：効果量-小
p値が小さいからといって、効果が大きいとは限らないため、効果量はp値と分離して考えなければならない。p値はサンプルサイズに依存するが、効果量はサンプルサイズに依存しない。結果の解釈の際には、p値による二分法的判断のみを頼りにするのではなく、効果の大きさの指標である効果量の値も十分に考慮する必要がある。
5.7.2 p-hackingを避ける
p-hackingとは、有意な結果を出すための、意識的あるいは無意識的な操作のことで、研究不正の一種である。たとえば、以下のような行為が挙げられる。
- n増し：50人のデータをとってみたところ、有意な傾向が見られなかったから100人に増やす。
- サンプルサイズを増やすとp値は小さくなっていく。これを避けるためには、データを取得し始める前に、サンプルサイズ設計によって必要サンプルサイズを決めておく必要がある。
- 数撃ちゃ当たる戦法：20種類の検定ができるデータを取っておけば、そのすべてが本当は有意ではないのに、平均1回は有意と判断する。
- 検定回数が多ければ多いほど、アルファエラー（間違って有意と判断してしまうこと）の可能性が高まる。実験調査デザインはできる限りシンプルにすべきである。事前に明瞭に仮説を立てておけば、必然的にシンプルになるはずである。
- データをこねくり回す：試しに検定にかけたところ、有意な傾向が見られなかったから、外れ値を除外して再度検定にかける。
- データ整理の方法は、検定を行う前に決めておかなければならない。データを調整して検定にかけることを繰り返していれば、どこかのタイミングで偶然p値が0.05を下回ってしまう。
- いろんな検定手法を試す：手法Aで検定にかけたところ、有意な傾向が見られなかったから、手法Bでも検定を行ってみた。
- 検定手法も、実験調査デザインの段階で事前に決めて置かなければならない。
これらを避けるために、事前登録（プレレジ：プレレジストレーション）が普及しつつある。事前登録とは、データを集める前に、第三者機関に下記のような情報を登録しておくシステムである。
- 背景・目的・仮説
- 研究デザイン・必要サンプルサイズ・手続き
- データの除外基準・欠測値の処理
- 解析手法
事前登録はしないとしても、p値の恣意的な操作が生じないように、実験調査デザインは事前に十分に練り上げる必要がある。とくに、解析手法を事前に決めずにデータを取る行為は絶対に避けたい。
この問題が生じる背景にはp値の過信があり、p値による二分法的思考に囚われすぎないようにすれば、もっと自由な研究活動ができると考えられる。たとえば、数撃ちゃ当たる戦法は、有意性を見出すという意味では問題があるが、仮説探索という観点では重要な試みであり、科学の発展に寄与する。また、データをいろいろと調整してみたり、検定手法を工夫することも、最適な統計モデルを探すために必要となることがある。
より詳しくは「池田・平石（2016）心理学における再現可能性危機：問題の構造と解決策」や「平石・中村（2022）心理学における再現性危機の10年」を参照したい。
5.7.3 サンプルサイズ設計
サンプルサイズの過剰・不足を回避し、適切にp値による判断を行うために、データを取得しはじめる前にサンプルサイズ設計によって最適なサンプルサイズを事前に決定しておくことが強く推奨される。サンプルサイズ設計の実施内容は、必ず論文中で報告する。
サンプルサイズ設計の理論、ならびにサンプルサイズ設計用フリーソフトウェア「G*power」の使い方は、「水本・竹内（2011）効果量と検定力分析入門: 統計的検定を正しく使うために」が詳しい。
G*powerではなく、Rでもpwrパッケージを用いることでサンプルサイズ設計ができる。推定効果量、有意水準（sig.level）、検出力（power
= 1 -
ベータエラー率）を指定すれば、必要サンプルサイズが導出できる。有意水準は5%、検出力は0.8で基本的に固定なので、効果量を適切な値に設定すればよい。なお「検定力
=
0.8」とは、実際に有意差があるときには80％の確率でそれを検出できることを意味する（水本,
2009）。
以下に、各検定におけるサンプルサイズ設計のコードを示す。
Pearsonの相関係数を用いた無相関検定を両側検定で行う場合のサンプルサイズ設計は、pwr.r.test関数を用いる。推定効果量rを中程度の0.3とすると、必要サンプルサイズnが85人であることがわかる（小数点以下切り上げ）。
pwr.r.test(r=0.3, sig.level=0.05, power=0.8)##
## approximate correlation power calculation (arctangh transformation)
##
## n = 84.07364
## r = 0.3
## sig.level = 0.05
## power = 0.8
## alternative = two.sided対応のあるt検定を両側検定で行う場合のサンプルサイズ設計は、pwr.t.test関数を用いる。対応がある場合、typeオプションにpairedを指定する。推定効果量dを中程度の0.5とすると、必要サンプルサイズnが34人であることがわかる（小数点以下切り上げ）。
pwr.t.test(d=0.5, sig.level=0.05, power=0.8, type="paired")##
## Paired t test power calculation
##
## n = 33.36713
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number of *pairs*対応のないt検定を両側検定で行う場合のサンプルサイズ設計は、pwr.t.test関数を用いる。対応がない場合、typeオプションにtwo.sampleを指定する。推定効果量dを中程度の0.5とすると、必要サンプルサイズnが128人（各群64人ずつ）であることがわかる（小数点以下切り上げ）。表示される人数は各群の人数なので、2倍する必要がある。
pwr.t.test(d=0.5, sig.level=0.05, power=0.8, type="two.sample")##
## Two-sample t test power calculation
##
## n = 63.76561
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* groupなお、効果量をどのように設定するかは、研究内容によって変わる。先行研究の効果量を使うのが最も望ましいが、新規性の高い研究で先行研究が存在しないような場合は、中程度の効果量を使う事が多い。
5.7.4 検定の多重性 Part 2
無相関検定と同様に、t検定の前に正規性の検定や等分散性の検定を行う必要はない。これらの検定を行ったところで、正規性があることや等分散であることは立証できないし、検定の多重性問題が生じてアルファエラー率とベータエラー率を維持できなくなる。
さらに、今回は学習のために、1つのデータセットに対して2回の検定（対応のあるt検定と対応のないt検定）を実施した。このように、同一データに対して複数回の検定を行う場合でも、その回数だけ抽選を行っていることと同義であり、検定の多重性問題が生じるため、本来であれば補正が必要となる。
今回のデータ構造の場合、本来であれば、学習前後の比較と性別の比較を別々に行うのではなく、後述する分散分析で検討すべきである。そうすれば、t検定で調べることのできることを超える事象を、1回の検定で調べることができる。
5.7.5 両側検定と片側検定
AとBに差がある、という対立仮説を設定する場合は両側検定を用いるが、AよりもBのほうが大きい（あるいはBよりもAのほうが大きい）という対立仮説を設定する場合は片側検定を用いる。
しかし、実際の状況において、片側検定を用いることはほとんどなく、よほどの例外的な仮説を扱わない限りは両側検定を使うことが推奨される。
たとえば、新薬開発の現場における、新薬と偽薬の効用を比較する実験を想定する。研究者の期待は、新薬の効用が偽薬の効用を上回ることなので、片側検定を使いたくなるかもしれない。ただ、有害事象などの影響によって新薬の効用が偽薬の効用を下回る可能性もあるため、片側検定ではなく両側検定を使うのが基本となるだろう。
なお、片側検定を用いた場合は、有意水準を5%から2.5%に引き下げることが推奨されており、片側検定にすれば有意になりやすいというわけでもない（奥村「検定と区間推定」）。
5.7.6 Brunner-Munzel検定
Brunner-Munzel検定は、TOSTERパッケージに含まれるbrunner_munzel関数で実施できる。
対応のあるBrunner-Munzel検定は下記のコードで解析できる。
brunner_munzel(df2$score_after, df2$score_before, paired=TRUE)## Sample size in at least one group is small. Permutation test (perm = TRUE) is highly recommended.##
## exact paired Brunner-Munzel test
##
## data: df2$score_after and df2$score_before
## t = 15.478, df = 11, p-value = 8.189e-09
## alternative hypothesis: true relative effect is not equal to 0.5
## 95 percent confidence interval:
## 0.8872001 1.0155777
## sample estimates:
## p(X>Y) + .5*P(X=Y)
## 0.9513889サンプルサイズが10あるいは15未満の場合は、対応のある並べ替えBrunner-Munzel検定を行う。オプションにperm=TRUEをつければよい。
brunner_munzel(df2$score_after, df2$score_before, paired=TRUE, perm=TRUE)## NOTE: Confidence intervals derived from permutation tests may differ from conclusions drawn from p-values. When discrepancies occur, consider additional diagnostics or alternative inference methods.##
## Paired Brunner-Munzel permutation test
##
## data: df2$score_after and df2$score_before
## t-observed = 15.478, df = 11, p-value < 2.2e-16
## alternative hypothesis: true relative effect is not equal to 0.5
## 95 percent confidence interval:
## 0.8844945 1.0182832
## sample estimates:
## p(X>Y) + .5*P(X=Y)
## 0.9513889対応のないBrunner-Munzel検定は下記のコードで解析できる。
brunner_munzel(df2_f, df2_m)## Sample size in at least one group is small. Permutation test (perm = TRUE) is highly recommended.##
## two-sample Brunner-Munzel test
##
## data: df2_f and df2_m
## t = 0.82261, df = 5.2473, p-value = 0.4465
## alternative hypothesis: true relative effect is not equal to 0.5
## 95 percent confidence interval:
## 0.187833 1.000000
## sample estimates:
## p(X>Y) + .5*P(X=Y)
## 0.65サンプルサイズが10あるいは15未満の場合は、対応のない並べ替えBrunner-Munzel検定を行う。オプションにperm=TRUEをつければよい。
brunner_munzel(df2_f, df2_m, perm=TRUE)## NOTE: Confidence intervals derived from permutation tests may differ from conclusions drawn from p-values. When discrepancies occur, consider additional diagnostics or alternative inference methods.##
## two-sample Brunner-Munzel permutation test
##
## data: df2_f and df2_m
## t-observed = 0.82261, df = 5.2473, p-value = 0.4008
## alternative hypothesis: true relative effect is not equal to 0.5
## 95 percent confidence interval:
## 0.08905257 1.00000000
## sample estimates:
## p(X>Y) + .5*P(X=Y)
## 0.65Brunner-Munzel検定の効果量は、出力結果の最下部のestimatesの値である。この効果量の大きさの目安は下記の通りである。（井口、Brunner-Munzel 検定の効果量と Cohen の d）
- 0.71：効果量-大
- 0.64：効果量-中
- 0.56：効果量-小
理解度確認クイズ
以下の記述のうち、下線部が正しいものをすべてチェックしなさい。
正解を表示する
正解は【3】
5.8 練習問題2
通常の机を使う群（統制条件）とスタンディングデスクを使う群（実験条件）で、計算課題の成績に差があるかを調べてみましょう。計算課題の成績は高いほどよく、課題に集中できていると解釈されます。
STEP1〜5のヒントを参考に解析を行い、下記フォームからグラフと解析結果を解答してみてください。大学アカウントでログインしているとGoogleフォームにアクセスできませんので、自身のGoogleアカウントでログインしてください。
STEP1
まず、作業ディレクトリに保存されているExcelデータの「data2p.xlsx」を、「df2p」という名前のデータフレームに格納します。
1列目（desk_condition）に使用した机の種類（normal / standing）、2列目（calculation_task_score）に計算課題の得点が格納されている、20行2列のデータです。
## desk_condition calculation_task_score
## 1 normal 80
## 2 normal 74
## 3 normal 81
## 4 normal 70
## 5 normal 73
## 6 normal 73
## 7 normal 60
## 8 normal 80
## 9 normal 70
## 10 normal 80
## 11 standing 70
## 12 standing 70
## 13 standing 77
## 14 standing 56
## 15 standing 58
## 16 standing 69
## 17 standing 65
## 18 standing 78
## 19 standing 66
## 20 standing 61STEP2
続いて、以下のグラフと同じになるように、ggplot関数でグラフを作図してみましょう。
なお、出力されたグラフの大きさ（＝画面の大きさ）によって軸目盛りが変わるかもしれませんが、それは問題ありません。
STEP3
グラフが作図できたらggsave関数で図を保存しましょう。
保存したいグラフを作図した直後に下記コードを実行すると、「fig2p.png」というファイル名で画像データが作業ディレクトリに保存されます。
ここで作成したグラフを、Googleフォームにアップロードしてください。
ggsave("fig2p.png", dpi=300, width=1200, height=1200, units="px", device="png")STEP4
記述統計量として、変数ごとの平均値と不偏標準偏差を算出しましょう。
平均値はmean関数、不偏標準偏差はsd関数で算出できます。
STEP5
対応のないwelchのt検定を行いましょう。
このとき、左側の引数で平均値が大きい変数を指定します。
検定結果の報告文章を、Googleフォームに入力してください。
6 事例3｜分散分析・多重比較
6.1 データ読み込み
事例3では、3つのセッションがある記憶課題の得点に関するデータを扱う。
作業ディレクトリに保存されているExcelデータの「data3.xlsx」を、「df3」という名前のデータフレームに格納する。
# 事例3 ####
df3 = read.xlsx("data3.xlsx")
df3## age_group session1 session2 session3
## 1 older 30 34 35
## 2 older 27 31 30
## 3 older 33 30 34
## 4 older 34 35 32
## 5 older 23 25 23
## 6 older 26 24 28
## 7 young 30 52 54
## 8 young 32 40 65
## 9 young 25 45 50
## 10 young 28 33 55
## 11 young 31 55 53
## 12 young 35 30 401列目に年齢層、2列目にセッション1の記憶課題得点、3列目にセッション2の記憶課題得点、4列目にセッション3の記憶課題得点が格納されている、12行4列のデータフレームであることがわかる。
6.2 グラフによる可視化
6.2.1 折れ線グラフ
分散分析の結果を示す際は、棒グラフあるいは折れ線グラフがよく使用される。今回は折れ線グラフを作成する。
折れ線グラフの土台となるのはドットプロットであるため、まずはドットプロットを描画する。ドットプロットには、平均値と標準誤差の集計値が必要となる。
ドットプロットを描画するために、pivot_longer関数（ワイド型データをロング型データに変換）、group_by関数（グルーピング）、summarize関数（集計）によって、下表のような集計データフレームを作成する。
| session | mean | se |
|---|---|---|
| session1 | 29.5 | 1.08 |
| session2 | 36.2 | 2.87 |
| session3 | 41.6 | 3.84 |
df3_plot = df3 %>%
pivot_longer(
cols = c(session1, session2, session3),
names_to = "session",
values_to = "score"
) %>%
group_by(session) %>%
summarize(mean = mean(score),
se = std.error(score))
df3_plot## # A tibble: 3 × 3
## session mean se
## <chr> <dbl> <dbl>
## 1 session1 29.5 1.08
## 2 session2 36.2 2.87
## 3 session3 41.6 3.84ggplot関数内で、使用するデータフレーム（df3_plot）、x軸（セッション）、y軸（記憶課題得点）をそれぞれ指定する。
ドットプロットを描くためにgeom_point関数を追加し、さらに不偏標準誤差をエラーバーとして表示するために、geom_errorbar関数を追加する。
ggplot(df3_plot, aes(x=session, y=mean))+
geom_point()+
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), width=0.1)+
labs(x="Session", y="Score of memory task")
折れ線グラフにするために、geom_line関数を追加して点と点を線で結ぶ。
また、ggplot関数内のaes関数内に、group=1オプションを追加する。これは3つのセッションを1つのグループとして認識するという意味である。このオプションを追加しないと線を引くことができない。
fig3_1 = ggplot(df3_plot, aes(x=session, y=mean, group=1))+
geom_point()+
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), width=0.1)+
geom_line()+
labs(x="Session", y="Score of memory task")
fig3_1
さらに、年齢層ごとに折れ線グラフを分けて描画してみたい。そのためには、下表のような集計データフレームを作成する必要がある。
| age_group | session | mean | se |
|---|---|---|---|
| older | session1 | 28.8 | 1.74 |
| older | session2 | 29.8 | 1.85 |
| older | session3 | 30.3 | 1.80 |
| young | session1 | 30.2 | 1.40 |
| young | session2 | 42.5 | 4.10 |
| young | session3 | 52.8 | 3.30 |
group_by関数内で、セッションと年齢層を同時にグルーピング指定することで実現できる。
df3_plot = df3 %>%
pivot_longer(
cols = c(session1, session2, session3),
names_to = "session",
values_to = "score"
) %>%
group_by(age_group, session) %>%
summarize(mean = mean(score),
se = std.error(score))## `summarise()` has grouped output by 'age_group'. You can override using the
## `.groups` argument.df3_plot## # A tibble: 6 × 4
## # Groups: age_group [2]
## age_group session mean se
## <chr> <chr> <dbl> <dbl>
## 1 older session1 28.8 1.74
## 2 older session2 29.8 1.85
## 3 older session3 30.3 1.80
## 4 young session1 30.2 1.40
## 5 young session2 42.5 4.10
## 6 young session3 52.8 3.30年齢層ごとにグルーピングして、それぞれのグループで折れ線グラフを描くため、aes関数内でgroup=age_groupと指定する。このようにすると、olderグループとyoungグループの2つのグループが認識され、それぞれのグループで線が引かれるようになる。
ggplot(df3_plot, aes(x=session, y=mean, group=age_group))+
geom_point()+
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), width=0.1)+
geom_line()+
labs(x="Session", y="Score of memory task")
このままだと、どちらの線が若年層でどちらの線が高齢層なのか区別がつかない。区別をつけるためには、たとえば線の色（colour）、ドットの形状（shape）、線の種類（linetype）をグループごとに変える必要がある。
今回は線の色のみを変えてみる。一番上の行のaes関数内にcolour=age_groupオプションを追加する。年齢層グループに応じて線の色を変えなさいという指示である。
ggplot(df3_plot, aes(x=session, y=mean, colour=age_group, group=age_group))+
geom_point()+
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), width=0.1)+
geom_line()+
labs(x="Session", y="Score of memory task")
ggplot2の標準色はカラーユニバーサル（色覚異常者でも明瞭に区別ができること）に対応していない。おすすめとして、カラーユニバーサルの専門家が制作した「岡部・伊藤カラーパレット」を使うとよい。
塗りつぶし色はscale_fill_okabeito関数、線の色はscale_colour_okabeito関数を追加することで、岡部・伊藤カラーパレットに準じた色に変更できる。これらの関数はseeパッケージに含まれている。
今回は線しか使っていないため、scale_colour_okabeito関数のみを追加する。
※稀にエラーとなってしまう人がいるようである。原因不明のため、エラーとなってしまう場合は本コードはスキップして進めてほしい。
fig3_2 = ggplot(df3_plot, aes(x=session, y=mean, colour=age_group, group=age_group))+
geom_point()+
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), width=0.1)+
geom_line()+
labs(x="Session", y="Score of memory task")+
scale_colour_okabeito()
fig3_2
6.2.2 ファセット
特定の変数に基づいて、図を分割表示する機能をファセットと呼ぶ。
ファセット機能は、facet_grid関数で実現できる。丸括弧内は、波線の左が縦方向への分割、波線の右が横方向への分割を意味し、ドットは変数がないことを意味する。
今回は、横方向に年齢層で分割するため、波線の左にドット、波線の右にage_groupを指定する。
ggplot(df3_plot, aes(x=session, y=mean, group=1))+
geom_point()+
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), width=0.1)+
geom_line()+
labs(x="Session", y="Score of memory task")+
facet_grid(. ~ age_group)
6.3 参加者内1要因分散分析
6.3.1 知りたいこと
セッションによって記憶課題得点に違いがあるかを調べたい（分散分析）。また、違いがあるとしたら、どのセッションとどのセッションの間に差があるかを調べたい（多重比較検定）。
6.3.2 検定の流れ
分散分析は、「3つ以上の水準をもつ要因」や「2つ以上の要因の組み合わせ」を扱うときに利用する検定である。今回のデータは1要因3水準なので、1要因分散分析を実施する。なお、参加者間要因の場合は「参加者間1要因分散分析」、参加者内要因の場合は「参加者内1要因分散分析」などと呼ぶ。
一般的に分散分析というと、分散分析に加えて、多重比較検定や単純主効果検定などを内包する段階的な検定プロセスを指すケースが多い。今回の場合は、以下の考え方で分析を進めていく。この分析手順は、分散分析のために使用するANOVA君パッケージの処理に基づいている。
- 今回分析するデータにおける要因と水準を確認する。要因が参加者内要因か参加者間要因かによって、分析過程が変わる。
- 要因：セッション（参加者内要因）
- 水準：session1、session2、session3（3つの水準）
- Step1：分散分析
- まず、分散分析によって、要因（セッション）が単独でデータに影響を与えているかどうかを検定にかける。この影響のことを主効果と呼ぶ。今回の場合は、セッションという要因が記憶課題得点に与える影響のことを、セッションの主効果と表現する。分散分析で検定にかけるのは主効果があるかどうかであって、水準間の差はこの時点ではわからない。
- 分散分析は、t検定などと同様に、正規性と等分散性を仮定している。参加者内要因では、球面性という水準間の差の等分散性に関する仮定も存在する（詳細は「ねこすたっと、分散分析（ANOVA）：球面性仮定」や「ANOVA君/球面性検定の出力」を参照）。
- 球面性仮定を確かめるために、ANOVA君では分散分析の前にMendozaの球面性検定が行われる。ただし、現実的に球面性仮定が満たされる状況は少ないため、最初から球面性仮定が満たされないと仮定して、常にε（球面性が成立する度合いの指標）による自由度補正を行う。球面性検定の省略は、検定の多重性問題への対処にもなる。
- 自由度補正法としては、「Greenhouse-Geisserのε」や「Huynh-Feldtのε」が有名だが、どちらもεの大きさに基づく使用条件があるため、筆者はより汎用的に使える手法として、「Chi-Muller（チ・ミュラー）法のε」を用いた補正を用いる。ただし、Chi-Muller法のεが1を超えた場合は自由度調整が行われず、εが下限値（LB: lower.bound）を下回った場合は下限値による自由度調整が行われるため、常にChi-Muller法のεが使われるわけではない点に注意したい。
- 有意な主効果が見られなかった場合、要因の主効果の有無は不明であり、今回は判断を保留することとなる。
- Step2：多重比較
- 3つ以上の水準をもつ要因において有意な主効果が見られた場合、どの水準間に差があるかを調べるために多重比較検定を行う。多重比較検定は、t検定を水準ペアごとに繰り返し、さらに検定の多重性問題を回避するために補正を行う検定である。今回の場合は、「session1とsession2の比較」「session1とsession3の比較」「session2とsession3の比較」と、計3種類のペアが存在しているので、t検定を3回行うこととなる。
- 多重比較の補正法はさまざまで、ANOVA君では「Holm法」や「Shaffer法」などが選択できる。今回は、仮定が少ないため汎用性が高く、さらに検出力が過度に下がりすぎない「Holm法」を用いて補正を行う。多重比較法の詳細は「ANOVA君/多重比較の方法」を参照。
- 注意点として、分散分析と多重比較検定は全く別の検定であるため、これらが常にセットであると考える必要はない。多くの文献では、「分散分析で主効果が有意だったら多重比較検定をする」と記述があるが、多重比較が当初からの研究目的であるならば、分散分析を行わずに多重比較検定をしてもよい。分散分析→多重比較という検定の流れだと、検定を2回繰り返すこととなり、検定全体でのアルファエラー率とベータエラー率を維持できない。（参考：井口「分散分析の下位に多重検定を置くな」）
6.3.3 多重比較の補正とは？
検定では、実際には一切差がなかったとしても、およそ20回に1回はアルファエラーが生じてしまう。検定を複数回行えば、その1/20の抽選を何度も行うことを意味し、検定全体としてアルファエラー率が著しく上昇してしまう。そこで、有意水準やp値を補正することによって、事前に設定したアルファエラー許容率5%を維持しようと考える。
様々な補正法が存在するが、ここでは分析者自身に厳しい制約を課す最も保守的な補正法である「Bonferroni法」の補正について簡単に説明する。なぜ厳しい制約を課すのかというと、保守的な補正法のもとで有意性がみられれば、だれも結果に文句は言わないだろうと考えられるからである。
たとえば以下のように3つの水準があるデータのうち、どの群間に差があるかを調べる場合、3C2 = 計3回のt検定を繰り返すこととなる。この場合、有意水準5%を1/3（つまり1.66…%）にしたり、p値を3倍に補正したりすることによって、3回のt検定全体の有意水準を5%に維持しようとする。したがって、3水準の多重比較検定で有意な差を認めるためには、補正前のp値が0.0166..を下回らなければならない。
なお、本教材でこれから使用するHolm法は、Bonferroni法よりも制約が緩く有意性がみられやすい。Bonferroni法が厳しすぎてベータエラー率が過度に高まってしまうことを避けるために、Bonferroni法を改良して作成された手法である。
6.3.4 ANOVA君の解説
Rでの分散分析は、井関が開発したANOVA君パッケージに含まれているanovakun関数が大変便利である。たとえばanovakun関数では、分散分析を下記のようなコードで実施する。
anovakun(データフレーム名, "AsB", 要因Aの名前=c("水準名1","水準名2","水準名3"), 要因Bの名前=c("水準名1","水準名2"), holm=T, cm=T, geta=T)括弧内のオプションはそれぞれ下記のような意味である。holm法、Chi-Mullerのε、一般化イータ二乗の利用は筆者の方針なので、適宜変更すること。
| オプション | 意味 |
|---|---|
| データフレーム名 | データフレームの名前を記入する |
| “AsB” | 要因の個数と種別を示している。小文字sの左には参加者間要因の数、小文字sの右には参加者内要因の数だけ、大文字アルファベットを左から順番に割り当てる。 |
| 要因名=c(“水準名1”, …) | 要因と水準の自分にとって分かりやすい名前をつける。ここでの名前が出力結果に表示される。この要因＆水準の名前セットは、要因の数だけ記入する。 |
| holm=T | 多重比較検定の補正法を指定するオプション。holm法による多重比較補正を行う。 |
| cm=T | 分散分析の球面性仮定不成立に対処するための、参加者内要因の主効果に対する自由度調整法を指定するオプション。Chi-Mullerのεによる自由度調整を行う（1を超えた場合は1に修正され、下限値を下回った場合は下限値に修正される）。 |
| geta=T | 分散分析の効果量に関するオプション。一般化イータ二乗の値を出力する。 |
要因を指定する”AsB”の部分は、たとえばケースごとに以下のようになる。
- 参加者間要因が1つの場合：“As”
- 参加者内要因が3つの場合：“sABC”
- 参加者間要因が2つと参加者内要因が1つの場合：“ABsC”
6.3.5 検定の実施
まずは、anovakun関数に読み込ませるデータフレームを作成する。読み込ませるデータは、分析に使うデータのみを残し、ほかの不要なデータを排除しておかなければならない。また、参加者内要因の水準は横方向、参加者間要因の水準は縦方向に並べる必要がある（つまりワイド型データ）。
df_anova = df3 %>% select(session1:session3)
df_anova## session1 session2 session3
## 1 30 34 35
## 2 27 31 30
## 3 33 30 34
## 4 34 35 32
## 5 23 25 23
## 6 26 24 28
## 7 30 52 54
## 8 32 40 65
## 9 25 45 50
## 10 28 33 55
## 11 31 55 53
## 12 35 30 40今回は参加者内1要因なので、“sA”と入力する（sの右側に1つアルファベットを割り当て）。また、要因の名前（Session）と各水準の名前（session1~3）を命名する。holm、cm、getaオプションは、要因水準計画によらず固定でよい。
anovakun(df_anova, "sA", Session=c("session1","session2","session3"), holm=T, cm=T, geta=T)##
## [ sA-Type Design ]
##
## This output was generated by anovakun 4.8.9 under R version 4.5.1.
## It was executed on Sun Oct 5 04:18:44 2025.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## ----------------------------------
## Session n Mean S.D.
## ----------------------------------
## session1 12 29.5000 3.7538
## session2 12 36.1667 9.9529
## session3 12 41.5833 13.2902
## ----------------------------------
##
##
## << SPHERICITY INDICES >>
##
## == Mendoza's Multisample Sphericity Test and Epsilons ==
##
## -------------------------------------------------------------------------
## Effect Lambda approx.Chi df p LB GG HF CM
## -------------------------------------------------------------------------
## Session 0.3356 1.9850 2 0.3706 ns 0.5000 0.8474 0.9854 0.8963
## -------------------------------------------------------------------------
## LB = lower.bound, GG = Greenhouse-Geisser
## HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller
##
##
## << ANOVA TABLE >>
##
## == Adjusted by Chi-Muller's Epsilon ==
##
## ----------------------------------------------------------------------
## Source SS df MS F-ratio p-value G.eta^2
## ----------------------------------------------------------------------
## s 1905.4167 11 173.2197
## ----------------------------------------------------------------------
## Session 879.1667 1.79 490.4476 7.5426 0.0046 ** 0.2162
## s x Session 1282.1667 19.72 65.0239
## ----------------------------------------------------------------------
## Total 4066.7500 35 116.1929
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## << POST ANALYSES >>
##
## < MULTIPLE COMPARISON for "Session" >
##
## == Holm's Sequentially Rejective Bonferroni Procedure ==
## == The factor < Session > is analysed as dependent means. ==
## == Alpha level is 0.05. ==
##
## ----------------------------------
## Session n Mean S.D.
## ----------------------------------
## session1 12 29.5000 3.7538
## session2 12 36.1667 9.9529
## session3 12 41.5833 13.2902
## ----------------------------------
##
## -----------------------------------------------------------------------------------
## Pair Diff t-value df p adj.p
## -----------------------------------------------------------------------------------
## session1-session3 -12.0833 3.2551 11 0.0077 0.0230 session1 < session3 *
## session1-session2 -6.6667 2.3193 11 0.0406 0.0813 session1 = session2
## session2-session3 -5.4167 2.0331 11 0.0669 0.0813 session2 = session3
## -----------------------------------------------------------------------------------
##
##
## output is over --------------------///6.3.6 検定結果の解釈
出力結果を上から順に解説する。
- DESCRIPTIVE STATISTICS（記述統計量）
- 記述統計量として、各水準の平均値と不偏標準偏差が算出されている。
- 各水準の平均値と不偏標準偏差は、論文中で報告する。
- SPHERICITY INDICES（球面性検定）
- 参加者内要因を含む場合のみ、Mendoza法による球面性の検定が実施される。
- p値が0.05を下回った場合は、球面性が仮定されないと判断する。一般論として、球面性が仮定されない場合に、LB（ε下限値）、GG（Greenhouse-Geisser法のε）、HF（Huynh-Feldt-Lecoutre法のε）、CM（Chi-Muller法のε）のいずれかのεの値を用いて、分散分析の自由度を調整する。
- しかし、ほとんどの状況で球面性が仮定されないと見なし、球面性検定の結果によらず、本教材では常にεによる自由度調整を行う。球面性検定を用いないのは、検定の多重性を回避するためでもある。
- 今回は、CM（Chi-Mullet）のεが0.90であり、1.00を超えておらず、ε下限値（LB）の0.50を上回っているため、Chi-Mullet法のεを用いて自由度調整している。
- Chi-Mullet法のε、あるいはε下限値による自由度調整を行った場合は、その旨とεの値を論文で報告する。
- ANOVA TABLE（分散分析）
- 各行は下記の意味である。中央横の破線の上側が参加者内要因、下側が参加者間要因を示す。
- s：参加者内での変動（つまり今回は3つのセッション間の個人内変動）。subjects（参加者）の略語。
- Session：自身が設定したセッション要因の主効果。
- s × Session：誤差項と呼ばれ、s要因とSession要因では説明できない残差を示している。
- Total：全体の変動。
- 各列は下記の意味である。
- 平方和SS：データの変動量。
- 自由度df
- 平均平方MS：変動の平均量。SSをdfで割った値。
- 検定統計量F：その行のMSを誤差項のMSで割った値。
- p値
- 効果量G.eta^2：一般化イータ二乗の値。
- 要因の行で、p値が0.05を下回っていれば、その要因における「すべての水準の母平均は等しい」という帰無仮説が棄却され、「すべての水準の母平均は等しいとはいえない」という対立仮説が採択される。
- 今回は、Session要因のp値が0.05を下回っているため、Session要因の記憶課題得点への主効果が確認された。
- 論文では、要因における、F値、自由度df、誤差項自由度df、p値、一般化イータ二乗を報告する。さらに、Chi-Muller法による自由度調整を行った旨と、Chi-Muller法のεを報告する。
- 各行は下記の意味である。中央横の破線の上側が参加者内要因、下側が参加者間要因を示す。
- POST ANALYSIS（下位検定）
- 下位検定とは、分散分析で要因の主効果や交互作用（後述）が有意となった場合に行われる検定のことである。ここからはすべて下位検定結果となる。
- MULTIPLE COMPARISON for “Session”（Session要因における多重比較検定）
Holm=Tをオプションに指定したため、冒頭にHolm法を用いたと書かれている。- まずは各水準における記述統計量として、平均値と不偏標準偏差が書かれている。この内容は、1要因分散分析の場合は、分析冒頭の「DESCRIPTIVE STATISTICS」の欄と同じ値である。
- その下には、ペアごとのt検定の結果が書かれている。各行には、Diff（水準間の平均値の差）、t値、自由度df、p値、調整済みp値（adj.p：adjusted p-valueの略）が記載されている。
- 検定結果の解釈には、単なるp値ではなく、調整済みp値を用いる。調整済みp値はHolm法によって調整された値であるため、検定の多重性問題に対応している。
- 今回は、session1とsession3の水準間でのみ、有意な差が見られた。
- 論文では、t値、自由度df、調整済みp値を報告する。
6.3.7 結果の報告例
まず、各セッションにおける記憶課題得点の平均値を計算したところ、セッション1は29.5点（SD = 3.8）、セッション2は36.2点（SD = 10.0）、セッション3は41.6点（SD = 13.3）であった（Figure 00）。
セッション間に差があるかを調べるために、セッションを参加者内要因とする1要因分散分析を実施したところ、セッションの有意な主効果が見られた（F(1.79, 19.72) = 7.54, p = .005, ηG2 = 0.22, εCM= 0.90, Chi-Muller法のεによる自由度補正あり）。そこで、セッションの要因においてHolm法を用いた多重比較検定を行ったところ、セッション1よりもセッション3のほうが有意に大きいことが明らかとなった（t(11) = 3.26, pholm = .023）。なお、セッション1とセッション2の間（t(11) = 2.32, pholm = .081）と、セッション2とセッション3の間（t(11) = 2.03, pholm = .081）には有意な差は見られなかった。
Figure 00
セッションごとの記憶課題得点
注）エラーバーは標準誤差を示す。
【APAマニュアルに基づく検定結果のフォーマットルール】
- 左端と右端の丸括弧は全角で、丸括弧の内側はすべて半角英数字である。
- p値は小数点以下3桁、それ以外は小数点以下2桁まで表記する（記述統計量は指標によって調整する）。
- p値を不等号で示すのは、.001を下回る場合のみである（p < .001）。
- M、SD、t、r、p、F、dなどはイタリック体にする。ηやεなどのギリシャ文字や、CIはイタリック体にしない。
- rやpなど、絶対値が必ず1以下となる数値においては、小数点の左の0を省略する。
- 詳細は「新美（2023）心理学論文における数値と統計の書き方」を参照。
6.4 混合2要因分散分析
6.4.1 知りたいこと
年齢層とセッションによって記憶課題得点に違いがあるかを調べたい（分散分析）。とくに、年齢層によって記憶課題得点の変化に違いがあるかを調べたい（交互作用）。
6.4.2 検定の流れ
つづいて、1つの参加者間要因（2水準）と1つの参加者内要因（3水準）による混合2要因分散分析について考えたい。混合とは、参加者内要因と参加者間要因が混ざっていることを意味する。
2要因以上の分散分析の場合、主効果に加えて、「交互作用」と「単純主効果」という概念が登場する。それぞれの意味を、事例に照らし合わせて確認する。
- 主効果（要因の影響）
- 年齢の主効果：年齢層によって記憶課題得点が異なる。
- セッションの主効果：セッションによって記憶課題得点が異なる。
- 交互作用（複数の要因の組み合わせ作用）
- 年齢とセッションの交互作用：年齢層とセッションの組み合わせ効果がある。より分かりやすく言うと、グラフ中の2本の折れ線が平行ではなく傾きが異なっている。たとえば、若年層はセッションを重ねるごとに記憶課題得点が高くなるが高齢層は変わらない、など。
- セッション要因の単純主効果（年齢層要因の1つの水準におけるセッションの主効果）
- 若年層におけるセッションの単純主効果：若年層において、セッションによって記憶課題得点が異なるかどうか。
- 高齢層におけるセッションの単純主効果：高齢層において、セッションによって記憶課題得点が異なるかどうか。
- 年齢層要因の単純主効果（セッション要因の1つの水準における年齢層の主効果）
- セッション1における年齢層の単純主効果：セッション1において、年齢層によって記憶課題得点が異なるかどうか。
- セッション2における年齢層の単純主効果：セッション2において、年齢層によって記憶課題得点が異なるかどうか。
- セッション3における年齢層の単純主効果：セッション2において、年齢層によって記憶課題得点が異なるかどうか。
主効果・交互作用・単純主効果は、実際に問題を解いてみないと理解が極めて困難である。下記リンク先にこれらのワードを理解するためのワークシートを用意しているので、ぜひ取り組んでみてほしい。
今回の場合は、以下の考え方で分散分析を進めていく。この分析手順は、分散分析のために使用するANOVA君パッケージの処理に基づいている。
- 今回分析するデータにおける要因と水準を確認する。
- 要因1：年齢層（参加者間要因）
- 要因1の水準：若年層、高齢層（2水準）
- 要因2：セッション（参加者内要因）
- 要因2の水準：session1、session2、session3（3水準）
- Step1：分散分析
- まず、分散分析によって、要因1（年齢層）と要因2（セッション）がそれぞれデータに影響を与えているかどうかを検定にかける。この影響のことを主効果と呼ぶ。分散分析で検定にかけるのは主効果があるかどうかであって、水準間の差はわからない。
- さらに交互作用も調べる。交互作用が見られた場合は、主効果は解釈せずに、単純主効果検定の結果を重視するのが基本となる（単純主効果のほうが主効果よりも高い解像度でデータを説明するため）。
- 球面性仮定が満たされないことに対する自由度調整法として、Chi-Mullerのεを用いた補正を用いる。ただし、Chi-Muller法のεが1を超えた場合は自由度調整が行われず、εが下限値（LB: lower.bound）を下回った場合は下限値による自由度調整が行われるため、常にChi-Muller法のεが使われるわけではない点に注意したい。
- Step2：多重比較
- 3つ以上の水準をもつ要因の有意な主効果が見られた場合、どの水準間に差があるかを調べるために、多重比較検定を行う。多重比較検定は、t検定を水準ペアごとに繰り返し、さらに検定の多重性問題を回避するために補正を行う検定である。今回の場合は、「session1とsession2の比較」「session1とsession3の比較」「session2とsession3の比較」と、計3種類のペアが存在している。
- 有意な交互作用が見られた場合は、主効果に対する多重比較は、主効果と同様に解釈しない方針となる。
- Step3：単純主効果検定
- 有意な交互作用が見られた場合に、単純主効果検定が行われる。単純主効果とは、片方の要因の1つの水準に対する、もう片方の要因の主効果のことである。したがって、この検定は水準の個数分だけ行われることになる。
- Step4：単純主効果の多重比較
- 3つ以上の水準をもつ要因の単純主効果が有意となった場合、多重比較検定が行われる。
6.4.3 検定の実施
まずは、anovakun関数に読み込ませるデータフレームを作成する。読み込ませるデータは、分析に使うデータのみを残し、他の不要なデータを排除しておかなければならない。また、参加者内要因の水準は横方向、参加者間要因の水準は縦方向に並べる必要がある（つまりワイド型データ）。
df_anova2 = df3
df_anova2## age_group session1 session2 session3
## 1 older 30 34 35
## 2 older 27 31 30
## 3 older 33 30 34
## 4 older 34 35 32
## 5 older 23 25 23
## 6 older 26 24 28
## 7 young 30 52 54
## 8 young 32 40 65
## 9 young 25 45 50
## 10 young 28 33 55
## 11 young 31 55 53
## 12 young 35 30 40今回は参加者間要因が1つと参加者内要因が1つなので、AsBと入力する。また、要因の名前（AgeGroup、Session）と各水準の名前（young-older、session1-3）を、要因ごとに命名する。holm、cm、getaオプションは、要因水準計画によらず固定でよい。
anovakun(df_anova2, "AsB", AgeGroup=c("older", "young"), Session=c("session1","session2","session3"), holm=T, cm=T, geta=T)##
## [ AsB-Type Design ]
##
## This output was generated by anovakun 4.8.9 under R version 4.5.1.
## It was executed on Sun Oct 5 04:18:45 2025.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## --------------------------------------------
## AgeGroup Session n Mean S.D.
## --------------------------------------------
## older session1 6 28.8333 4.2622
## older session2 6 29.8333 4.5350
## older session3 6 30.3333 4.4121
## young session1 6 30.1667 3.4303
## young session2 6 42.5000 10.0548
## young session3 6 52.8333 8.0850
## --------------------------------------------
##
##
## << SPHERICITY INDICES >>
##
## == Mendoza's Multisample Sphericity Test and Epsilons ==
##
## -------------------------------------------------------------------------
## Effect Lambda approx.Chi df p LB GG HF CM
## -------------------------------------------------------------------------
## Session 0.0006 12.2246 5 0.0318 * 0.5000 0.9509 1.1682 1.0416
## -------------------------------------------------------------------------
## LB = lower.bound, GG = Greenhouse-Geisser
## HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller
##
##
## << ANOVA TABLE >>
##
## == Adjusted by Chi-Muller's Epsilon ==
##
## --------------------------------------------------------------------------------
## Source SS df MS F-ratio p-value G.eta^2
## --------------------------------------------------------------------------------
## AgeGroup 1332.2500 1 1332.2500 23.2437 0.0007 *** 0.5298
## s x AgeGroup 573.1667 10 57.3167
## --------------------------------------------------------------------------------
## Session 879.1667 2 439.5833 14.4362 0.0001 *** 0.4265
## AgeGroup x Session 673.1667 2 336.5833 11.0536 0.0006 *** 0.3628
## s x AgeGroup x Session 609.0000 20 30.4500
## --------------------------------------------------------------------------------
## Total 4066.7500 35 116.1929
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## << POST ANALYSES >>
##
## < MULTIPLE COMPARISON for "Session" >
##
## == Holm's Sequentially Rejective Bonferroni Procedure ==
## == The factor < Session > is analysed as dependent means. ==
## == Alpha level is 0.05. ==
##
## ----------------------------------
## Session n Mean S.D.
## ----------------------------------
## session1 12 29.5000 3.7538
## session2 12 36.1667 9.9529
## session3 12 41.5833 13.2902
## ----------------------------------
##
## -----------------------------------------------------------------------------------
## Pair Diff t-value df p adj.p
## -----------------------------------------------------------------------------------
## session1-session3 -12.0833 6.0744 10 0.0001 0.0004 session1 < session3 *
## session1-session2 -6.6667 2.7498 10 0.0205 0.0410 session1 < session2 *
## session2-session3 -5.4167 2.3331 10 0.0418 0.0418 session2 < session3 *
## -----------------------------------------------------------------------------------
##
##
## < SIMPLE EFFECTS for "AgeGroup x Session" INTERACTION >
##
## -----------------------------------------------------------------------------------
## Effect Lambda approx.Chi df p LB GG HF CM
## -----------------------------------------------------------------------------------
## Session at older 0.8427 0.2738 2 0.8720 ns 0.5000 0.9379 1.4813 0.9069
## Session at young 0.8770 0.2100 2 0.9003 ns 0.5000 0.9513 1.5201 0.9307
## -----------------------------------------------------------------------------------
## LB = lower.bound, GG = Greenhouse-Geisser
## HF = Huynh-Feldt-Lecoutre, CM = Chi-Muller
##
## -------------------------------------------------------------------------------
## Source SS df MS F-ratio p-value G.eta^2
## -------------------------------------------------------------------------------
## AgeGroup at session1 5.3333 1 5.3333 0.3563 0.5638 ns 0.0344
## Er at session1 149.6667 10 14.9667
## -------------------------------------------------------------------------------
## AgeGroup at session2 481.3333 1 481.3333 7.9123 0.0184 * 0.4417
## Er at session2 608.3333 10 60.8333
## -------------------------------------------------------------------------------
## AgeGroup at session3 1518.7500 1 1518.7500 35.8055 0.0001 *** 0.7817
## Er at session3 424.1667 10 42.4167
## -------------------------------------------------------------------------------
## Session at older 7.0000 1.81 3.8593 0.8824 0.4365 ns 0.0235
## s x Session at older 39.6667 9.07 4.3739
## -------------------------------------------------------------------------------
## Session at young 1545.3333 1.86 830.2425 13.5714 0.0019 ** 0.6342
## s x Session at young 569.3333 9.31 61.1758
## -------------------------------------------------------------------------------
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## < MULTIPLE COMPARISON for "Session at young" >
##
## == Holm's Sequentially Rejective Bonferroni Procedure ==
## == The factor < Session at young > is analysed as dependent means. ==
## == Alpha level is 0.05. ==
##
## -----------------------------------------------------------------------------------
## Pair Diff t-value df p adj.p
## -----------------------------------------------------------------------------------
## session1-session3 -22.6667 5.8831 5 0.0020 0.0060 session1 < session3 *
## session1-session2 -12.3333 2.6268 5 0.0467 0.0934 session1 = session2
## session2-session3 -10.3333 2.3080 5 0.0691 0.0934 session2 = session3
## -----------------------------------------------------------------------------------
##
## output is over --------------------///6.4.4 検定結果の解釈
出力結果を上から順に解説する。
- DESCRIPTIVE STATISTICS（記述統計量）
- 記述統計量として、各水準の平均値と不偏標準偏差が算出されている。
- 各水準の平均値と不偏標準偏差は、論文中で報告する。
- SPHERICITY INDICES（球面性検定）
- 参加者内要因を含む場合のみ、Mendoza法による球面性の検定が実施される。
- ほとんどの状況で球面性が仮定されないと見なし、球面性検定の結果によらず、本教材では常にεによる自由度調整を行う。球面性検定を用いないのは、検定の多重性を回避するためでもある。ただし今回は、CM（Chi-Mullet）のεが1.04であり、1.00を超えているため、εによる自由度調整を行わない。
- ANOVA TABLE（分散分析）
- 各行は下記の意味である。中央横の破線の上側が参加者内要因、下側が参加者間要因を示す。
- AgeGroup：自身が設定した年齢層要因の主効果。
- s x AgeGroup：参加者内要因の誤差項。
- Session：自身が設定したセッション要因の主効果。
- AgeGroup x Session：年齢層要因とセッション要因の交互作用。
- s x AgeGroup x Session：参加者間要因の誤差項。
- Total：全体の変動。
- 各列の意味は、前と同じである。
- 要因の行で、p値が0.05を下回っていれば、その要因における「すべての水準の母平均は等しい」という帰無仮説が棄却され、「すべての水準の母平均は等しいとはいえない」という対立仮説が採択される。
- 交互作用の行で、p値が0.05を下回っていれば、その要因における「交互作用はない」という帰無仮説が棄却され、「交互作用がある」という対立仮説が採択される。
- 今回は、検討したすべての効果（年齢層の主効果、セッションの主効果、年齢層とセッションの交互作用）が有意となった。交互作用が有意となった場合は、主効果ではなく単純主効果に重きを置いて解釈する。
- 論文では、要因と交互作用における、F値、自由度df、誤差項自由度df、p値、一般化イータ二乗を報告する。
- 各行は下記の意味である。中央横の破線の上側が参加者内要因、下側が参加者間要因を示す。
- POST ANALYSIS（下位検定）
- 下位検定とは、分散分析で要因の主効果や要因間の交互作用が有意となった場合に行われる検定のことである。ここからは、すべて下位検定結果となる。
- MULTIPLE COMPARISON for “Session”（Session要因における多重比較検定）
- この多重比較検定は、交互作用が有意となった場合は論文では報告しない。単純主効果をベースに考える必要があるからである。
- やっていることは、以前に実施した参加者内1要因分散分析と同じである（ただし検定結果は変わることがある）。
- SIMPLE EFFECTS for “AgeGroup x Session”
INTERACTION（単純主効果検定）
- 冒頭に、球面性検定の結果が書かれている。今回も、分散分析のときと同じように、Chi-Muller法のεによる自由度調整を行う。いずれもCM（Chi-Mullet）のεが下限値（LB）よりも大きく1.00よりも小さいため、Chi-Mullet法のεによる自由度調整を行う。
- 球面性検定の下に単純主効果検定の結果が表示されている。Source列の意味は、たとえばAgeGroup at session1と表記されていたら、session1におけるAgeGroupの単純主効果という意味である。
- 今回は、「AgeGroup at session2」「AgeGroup at session3」「Session at young」の3つの単純主効果が有意となった。それぞれの意味は、ぜひグラフで確認したい。
- AgeGroupの単純主効果は自由度が調整されていないが、Sessionの単純主効果は自由度が調整されている（整数ではないことから判断）。これは、Chi-Muller法による自由度調整は、球面性仮定が求められる参加者内要因にしか実施されないためである（参加者間要因しかなければ自由度調整は行われない）。
- 論文では、各単純主効果における、F値、自由度df、誤差項自由度df、p値、一般化イータ二乗を報告する。さらに、参加者内要因に対してChi-Muller法による自由度調整を行った旨と、Chi-Muller法のεを報告する。
- MULTIPLE COMPARISON for “Session at
young”（若年層におけるセッション要因の多重比較検定）
- 単純主効果検定にて、「Session at young」（若年層におけるセッション要因）で有意な単純主効果が見られたことから、若年層におけるセッション要因の3つの水準の多重比較が行われる。
- Holm=Tをオプションに指定したため、冒頭にHolm法を用いたと書かれている。
- ペアごとのt検定の結果が書かれている。各行には、Diff（水準間の平均値の差）、t値、自由度df、p値、調整済みp値（adj.p：adjusted p-valueの略）が記載されている。
- 検定結果の解釈には、単なるp値ではなく、調整済みp値を用いる。調整済みp値はHolm法によって調整された値であるため、検定の多重性問題に対応している。
- 今回は、session1とsession3の水準間でのみ、有意な差が見られた。
- 論文では、t値、自由度df、調整済みp値を報告する。
6.4.5 結果の報告例
まず、各年齢層のセッションごとの記憶課題得点の平均値を計算したところ、若年層のセッション1は30.2点（SD = 3.4）、若年層のセッション2は42.5点（SD = 10.1）、若年層のセッション3は52.8点（SD = 8.1）、高齢層のセッション1は28.8点（SD = 4.3）、高齢層のセッション2は29.8点（SD = 4.5）、高齢層のセッション3は30.3点（SD = 4.4）であった（Figure 00）。
年齢層を参加者間要因、セッションを参加者内要因とする混合2要因分散分析を実施したところ、年齢層の主効果（F(1, 10) = 23.24, p < .001, ηG2 = 0.53）、セッションの主効果（F(2, 20) = 14.44, p < .001, ηG2 = 0.43）、年齢層とセッションの交互作用（F(2, 20) = 11.05, p < .001, ηG2 = 0.36）がそれぞれ有意となった。
有意な交互作用が見られたため、単純主効果検定を行った。その結果、セッション2における年齢層の単純主効果（F(1, 10) = 7.91, p = .018, ηG2 = 0.44）、セッション3における年齢層の単純主効果（F(1, 10) = 35.81, p < .001, ηG2 = 0.78）、若年層におけるセッションの単純主効果（F(1.86, 9.31) = 13.57, p = .002, ηG2 = 0.63, εCM= 0.93）がそれぞれ有意となった。セッション1における年齢層の単純主効果（F(1, 10) = 0.36, p = .564, ηG2 = 0.03）と、高齢層におけるセッションの単純主効果（F(1.81, 9.07) = 0.88, p = .437, ηG2 = 0.02, εCM= 0.91）は有意性はみられなかった。なお、 セッション要因に対しては、Chi-Muller法のεによる自由度の補正が行われた。
若年層におけるセッションの有意な単純主効果がみられたため、若年層における各セッションに対してHolm法を用いた多重比較検定を行ったところ、セッション1よりもセッション3のほうが有意に大きいことが明らかとなった（t(5) = 5.88, pholm = .006）。なお、セッション1とセッション2の間（t(5) = 2.63, pholm = .093）と、セッション2とセッション3の間（t(5) = 2.31, pholm = .093）には有意な差は見られなかった。
Figure 00
各年齢層におけるセッションごとの記憶課題得点
注）エラーバーは標準誤差を示す。
【APAマニュアルに基づく検定結果のフォーマットルール】
- 左端と右端の丸括弧は全角で、丸括弧の内側はすべて半角英数字である。
- p値は小数点以下3桁、それ以外は小数点以下2桁まで表記する（記述統計量は指標によって調整する）。
- p値を不等号で示すのは、.001を下回る場合のみである（p < .001）。
- M、SD、t、r、p、F、dなどはイタリック体にする。ηやεなどのギリシャ文字や、CIはイタリック体にしない。
- rやpなど、絶対値が必ず1以下となる数値においては、小数点の左の0を省略する。
- 詳細は「新美（2023）心理学論文における数値と統計の書き方」を参照。
6.5 押さえておきたいポイント
6.5.1 研究デザインをシンプルにする大切さ
多重比較補正では、水準数が多ければ多いほど、強い補正がかかることとなる。たとえば、3水準の多重比較を行った場合には3回検定を繰り返すため、Bonferroni法ではp値を3倍する（あるいは有意水準5%を3で割る）ことで調整する。そうなると、p値が0.016…を下回らないと有意性を認められなくなる。
このように、たくさんのデータを取得すると、その分だけ個々の検定で有意性がみられにくくなってしまうため、研究デザインはシンプルにするのが重要である。データはたくさん取れば取るほど、その分だけ意義があると思ってしまいやすいが、検定戦略をシンプルにすることが特に実験計画法では求められる。
6.5.2 多重比較における補正は本当に正しい？
多重比較において、Holm法やBonferroni法による補正を行うのは、検定全体のアルファエラー率を維持するためである。このような補正を行うと、たしかに読者に納得してもらいやすいかもしれないが、個別の検定のベータエラー率がかなり上がってしまう。
一部の研究者は、多重比較時に盲信的に補正を行うのは避け、補正前のp値を報告し、多重検定の問題への対処は読者に委ねることを推奨している。
特に、仮説探索型の調査研究は、多重比較補正との相性が悪い。研究目的に応じた柔軟な対応が求められる。
6.5.3 分散分析のノンパラメトリック検定
分散分析には様々な仮定があり、それらの仮定を満たさない場合も多いだろう。
たとえば、正規性が仮定されない場合には、クラスカル・ウォリス検定（参加者間要因用）やフリードマン検定（参加者内要因用）などのノンパラメトリック検定を使うように一般的に説明されることが多い。しかし、これらの検定は、t検定におけるウィルコクソンの符号順位検定や符号和検定と同様に扱いが難しく、万能に使える手法ではない（使うべきではないかもしれない）。
また、等分散性が仮定されない場合の対処としてwelchの分散分析なども存在するが、あまり使われていないようである。
仮定を満たさない場合は（というより仮定を満たす場合であっても）、後述する一般化線形混合モデルを用いたほうがよいかもしれない。
6.5.4 一般化混合モデルの利用の推奨
実は、これまで見てきた無相関検定、t検定、分散分析は、いずれも等分散・正規性を仮定する線形モデルに基づく検定手法である。しかし、線形モデルは極めて限定的なモデルであり、柔軟性に欠けていて、実際の観測データに合致しないことも多い。代替案としてしばしば提案されるノンパラメトリック検定は、一切の仮定なしに使えると勘違いされやすいが、実際には等分布性などを仮定しており、結果が歪みやすく大変使いにくい手法である。
そこで、正規分布以外の確率分布を扱える一般化線形モデル（GLM）、ランダム効果を扱える一般化線形混合モデル（GLMM）、さらに自由度の高い階層ベイズモデルなどを積極的に活用していく必要がある（久保、2012、データ解析のための統計モデリング入門）。名前だけ聞くと難しそうだが、少なくとも一般化線形混合モデルまでは比較的簡単に実装できる。
詳細は下記書籍を参照したい。特に緑本と呼ばれる久保（2012）の書籍はたいへん分かりやすい。
- 久保（2012）、データ解析のための統計モデリング入門
- 馬場（2019）、RとStanではじめる ベイズ統計モデリングによるデータ分析入門
- 松浦（2016）、StanとRでベイズ統計モデリング
- 豊田編（2018）、たのしいベイズモデリング
6.6 練習問題3
授業の実施形態と授業へのモチベーションを要因として、期末テストの成績の違いを調べてみましょう。つまり、異なる実施形態で学習を進め、どのくらい学んだかを期末テストで測定した形となります。要因水準の詳細は以下の通りです。今回は、2つの要因があり、2つとも参加者間要因ですので、参加者間2要因分散分析を行うこととなります。
- 1つ目の参加者間要因：授業の実施形態（対面／オンライン）
- 2つ目の参加者間要因：授業へのモチベーション（高／低）
STEP1〜5のヒントを参考に解析を行い、下記フォームからグラフと解析結果を解答してみてください。大学アカウントでログインしているとGoogleフォームにアクセスできませんので、自身のGoogleアカウントでログインしてください。
なお、分散分析の結果を書くのは大変なので、練習問題3では結果からどういうことが言えそうかを簡潔に書くだけで大丈夫です（統計量の記載は不要）。
STEP1
まず、作業ディレクトリに保存されているExcelデータの「data3p.xlsx」を、「df3p」という名前のデータフレームに格納します。
1列目（class_style）に授業の実施形態（in_person / online）、2列目（motivation）に授業へのモチベーション（high / low）、3列目にテストの得点（0 ~ 100点）が格納されている、30行3列のデータです。
## class_style motivation test_score
## 1 in_person high 76
## 2 in_person high 70
## 3 in_person high 94
## 4 in_person high 81
## 5 in_person high 73
## 6 in_person high 86
## 7 in_person high 74
## 8 in_person high 75
## 9 in_person low 62
## 10 in_person low 64
## 11 in_person low 76
## 12 in_person low 61
## 13 in_person low 60
## 14 in_person low 65
## 15 in_person low 75
## 16 online high 73
## 17 online high 76
## 18 online high 80
## 19 online high 71
## 20 online high 80
## 21 online high 70
## 22 online high 75
## 23 online high 72
## 24 online high 88
## 25 online low 41
## 26 online low 44
## 27 online low 45
## 28 online low 59
## 29 online low 48
## 30 online low 43STEP2
続いて、以下のグラフと同じになるように、ggplot関数でグラフを作図してみましょう。
なお、出力されたグラフの大きさ（＝画面の大きさ）によって軸目盛りが変わるかもしれませんが、それは問題ありません。
ggplot関数に読み込ませるデータフレームを作成するコードは、自力で書くのが難しいので、公開しておきます。
df3p_plot = df3p %>%
group_by(class_style, motivation) %>%
summarize(mean = mean(test_score),
se = std.error(test_score))
df3p_plot## # A tibble: 4 × 4
## # Groups: class_style [2]
## class_style motivation mean se
## <chr> <chr> <dbl> <dbl>
## 1 in_person high 78.6 2.82
## 2 in_person low 66.1 2.50
## 3 online high 76.1 1.91
## 4 online low 46.7 2.64
STEP3
グラフが作図できたらggsave関数で図を保存しましょう。
保存したいグラフを作図した直後に下記コードを実行すると、「fig3p.png」というファイル名で画像データが作業ディレクトリに保存されます。
ここで作成したグラフを、Googleフォームにアップロードしてください。
ggsave("fig3p.png", dpi=300, width=1200, height=1200, units="px", device="png")STEP4
今回、mean関数とsd関数で記述統計量を求める必要はありません。anovakun関数の実行結果内に平均値と標準偏差が出力されるためです。
STEP5
anovakun関数に読み込ませるデータフレームは、df3pがそのまま使えます。データの前処理は不要です。以下のコードを実行すれば解析できます。
今回は、参加者間要因が2つなので、anovakun関数内の要因設定は「“ABs”」となります。
また、参加者内要因が存在しないため、cmオプションは不要です。さらに、どちらの要因も2水準なので多重比較が行われないため、holmオプションも不要です。
なお、今回の分析は参加者内要因を含まないため、誤差項の表記は「Error」となりますが、考え方はこれまでと一緒です。
検定結果の解釈を、Googleフォームに入力してください。統計量などの詳細な検定結果を入力する必要はなく、得られた検定結果からどのようなことが言えそうかを日常言語でわかりやすくまとめてみてください。
anovakun(df3p, "ABs", Style=c("in_person","online"), Motivation=c("high","low"), geta=T)7 時間が余った方へ
以下の分析手法を手軽に体験できる章である。なお、本章はあくまでも補足的内容であり、第三者の論文で登場する解析手法を理解できるようになることが目的なので、詳細な説明はしない。
- \(\chi\)二乗検定（1x2分割表、2x2分割表）
- 単回帰分析
- 重回帰分析
- 線形混合モデル（LMM）
7.1 \(\chi\)二乗検定
7.1.1 1×2分割表
まずは以下のコードを実行して、分析するデータフレーム「data_extra1」を作成する。
name = c("たこ", "いか")
count = c(60, 40)
data_extra1 = data.frame(name=name, count=count)
data_extra1## name count
## 1 たこ 60
## 2 いか 40このデータは、100名に「たこ」と「いか」のどちらがより好きかを聞いた結果を示している。
たこのほうが好まれているように見えるが、統計的に有意な差があるかどうかを調べたい。
度数データ（観測頻度のこと）の比較には、\(\chi\)二乗検定を使用する。\(\chi\)は「エックス」ではなく「カイ」と読む。\(\chi\)二乗検定では、このような度数分布表（つまり集計表）を用いて分析を行う。
\(\chi\)二乗検定はchisq.test関数を用いる。
chisq.test(data_extra1$count)##
## Chi-squared test for given probabilities
##
## data: data_extra1$count
## X-squared = 4, df = 1, p-value = 0.0455結果はシンプルで、これまで同様にp値により判断する。p値が0.05未満のときに、たこといかに有意な差があると判断される。
ちなみに、割合は6:4で同じまま、人数を変えて再度分析をしてみたい。100人ではなく10人分の集計結果にしてみる。
name = c("たこ", "いか")
count = c(6, 4)
data_extra2 = data.frame(name=name, count=count)
data_extra2## name count
## 1 たこ 6
## 2 いか 4chisq.test(data_extra2$count)##
## Chi-squared test for given probabilities
##
## data: data_extra2$count
## X-squared = 0.4, df = 1, p-value = 0.5271\(\chi^2\)値とp値が変わったことからも分かる通り、割合が同じでも人数によって結果が変わる。
7.1.2 2×2分割表
先ほどは、1行2列の度数分布表の分析だったが、次は2行2列の度数分布表の分析をしてみたい。
まずは以下のコードを実行して、分析するデータフレーム「data_extra3」を作成する。
smoking = c("喫煙あり", "喫煙なし")
cancer = c("肺がん", "健康")
count_matrix = matrix(c(5, 10, 5, 80), nrow=2)
data_extra3 = as.table(count_matrix)
dimnames(data_extra3) = list("行"=smoking, "列"=cancer)
data_extra3## 列
## 行 肺がん 健康
## 喫煙あり 5 5
## 喫煙なし 10 802✕2の度数分布として有名な、喫煙と肺がんの例である。100人を対象に、喫煙の有無と肺がん罹患の有無をまとめており、喫煙者は1/2、非喫煙者は1/9の割合で肺がんになっていると仮定したデータである。
喫煙者のほうが肺がんに罹患しやすいように見えるが、統計的に有意な差があるかどうかを調べたい。
chisq.test(data_extra3)## Warning in chisq.test(data_extra3): Chi-squared approximation may be incorrect##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: data_extra3
## X-squared = 7.8431, df = 1, p-value = 0.005101結果の解釈は先ほどと同様である。
なお、chisq.testで2✕2分割表の検定を行う場合は、オプションに「correct=FALSE」をつけずに実行すると、イェーツの連続性補正が標準で行われ、検定結果が過大に見積もられるのを防ぐ。
また、上記コードでエラーメッセージが出ているように、度数が少ない（5以下が目安）セルが存在する場合は、\(\chi\)二乗検定ではなくFisherの正確確率検定がより適している場合がある。
7.2 単回帰分析
相関分析は、2変数間の関連を調べるものだったが、相関分析に似た解析として回帰分析がある。
回帰分析では、変数を目的変数と説明変数の2つに割り当てて、説明変数によって目的変数を予測することを試みる。説明変数が1つのみの場合、単回帰分析を行う。回帰式は以下のようになる。
目的変数 = \(\beta_0\) + (\(\beta_1\) × 説明変数) + 誤差項
それぞれの意味は下記の通りである。
- \(\beta_0\)：切片（Intercept）。説明変数が0のときの目的変数の値。
- \(\beta_1\)：回帰係数。回帰分析では、この回帰係数が有意かどうかを調べることとなる。
- 誤差項：そのほかの未知の要因や偶然のばらつきなど、切片と回帰係数だけではモデルを説明できない部分。残差とも呼ぶ。
第3章の無相関検定の練習で扱った、気温とアイス販売個数のデータ（df1）を改めて用いて、回帰分析を行ってみる。念の為、データの中身を再度確認する。1列目に気温、2列目にアイス販売個数が格納されている。
df1## temperature number_of_ice_creams_sold
## 1 22 23
## 2 26 30
## 3 32 40
## 4 28 39
## 5 35 50
## 6 27 28
## 7 30 20
## 8 38 50
## 9 37 34
## 10 13 24今回は、目的変数がアイス販売個数、説明変数が気温となるため、回帰式は以下のようになる。単回帰分析を行うことによって、\(\beta_0\)と\(\beta_1\)の値が算出される。
アイス販売個数 = \(\beta_0\) + \(\beta_1\) 気温 + 誤差項
早速、単回帰分析をしてみる。lm関数を用い、波線の左に目的変数（アイス販売個数）、波線の右に説明変数（気温）を設置する。最後に使用するデータフレームを指定する。
回帰分析の結果は、summary関数によってわかりやすい形式にまとめられる。
model1 = lm(number_of_ice_creams_sold ~ temperature, data=df1)
summary(model1)##
## Call:
## lm(formula = number_of_ice_creams_sold ~ temperature, data = df1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.976 -4.109 1.004 5.911 10.122
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.5651 10.9162 0.510 0.6240
## temperature 0.9804 0.3679 2.665 0.0286 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.305 on 8 degrees of freedom
## Multiple R-squared: 0.4702, Adjusted R-squared: 0.404
## F-statistic: 7.101 on 1 and 8 DF, p-value: 0.02859一番上の「Residuals」が誤差項（残差）を意味する。「Intercept」の行のEstimateが切片、「temperature」の行のEstimateが回帰係数である。どちらの値も、t値とp値が算出されている。回帰係数のp値を見てみると、0.05を下回っていることから、有意な回帰係数が得られたことがわかる。
切片が「5.5651」なので、気温が0度のときにアイス販売個数が約5.6個であると予測されたことがわかる。また、回帰係数は「0.9804」であることから、気温が1度上昇するごとにアイス販売個数が約0.98個増えると予測されたことがわかる。
グラフで視覚的に確認してみよう。グラフ中のどの部分が、切片\(\beta_0\)・回帰係数\(\beta_1\)・誤差項にそれぞれ対応しているだろうか？
ggplot(df1, aes(x=temperature, y=number_of_ice_creams_sold))+
geom_point()+
geom_smooth(method="lm", se=FALSE, fullrange=TRUE)+
scale_x_continuous(limits=c(-10,50))+
scale_y_continuous(limits=c(-10,60))## `geom_smooth()` using formula = 'y ~ x'
Multiple R-squared（重決定係数R2）とAdjusted R-squared（自由度調整済み決定係数R2adj）は、いずれも説明変数が目的変数の変動をどの程度予測できるかを示す指標である。単回帰分析では前者のR2を使用すればよい。後述する重回帰分析では、R2adjを用いる。1に近いほど説明変数が目的変数の変動をよく予測できることを意味する。決定係数の値が低い場合は、説明変数として投入されていない未知の要因の影響（誤差項の部分）が大きいことを意味する。絶対的な基準はないものの、決定係数が0.5以上のモデルは、投入した説明変数によって目的変数をよく予測できていると考えられるかもしれない。
なお、F値とp値は、説明変数全体が目的変数に有意に影響しているかを検定した結果である。ただし、単回帰分析では説明変数を1つしか投入しないため、説明変数の回帰係数のt検定と結果は変わらないことから、t検定の結果を報告しておけば、F検定については言及不要と考えられる。説明変数を2つ以上投入する重回帰分析では、F検定の結果が重要な意味を持つ。
ちなみに、目的変数と説明変数を逆にしてみたらどうなるかも試してみよう。
model2 = lm(temperature ~ number_of_ice_creams_sold, data=df1)
summary(model2)##
## Call:
## lm(formula = temperature ~ number_of_ice_creams_sold, data = df1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.0995 -1.6074 -0.3756 1.3178 8.1041
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.5879 6.3551 1.981 0.0830 .
## number_of_ice_creams_sold 0.4796 0.1800 2.665 0.0286 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.809 on 8 degrees of freedom
## Multiple R-squared: 0.4702, Adjusted R-squared: 0.404
## F-statistic: 7.101 on 1 and 8 DF, p-value: 0.02859目的変数と説明変数を逆にしても、t値とp値は変わらない（もちろん切片や回帰係数などは変わる）。これは、2変数の相関係数に基づいてこれらの計算を行っているためである。さらにいうと、これらのt値とp値は、Pearsonの無相関検定の結果とも一致する。
無相関検定と単回帰分析の使い分けは、分析目的が「関連の解明（＝無相関検定）」と「モデル作成による予測（＝単回帰分析）」のどちらなのかによる。とくに、因果関係を想定しやすい場合は、単回帰分析のほうが理にかなっているかもしれない。
7.3 重回帰分析
単回帰分析は説明変数が1つのみの回帰モデルだったが、説明変数が2つ以上存在する場合は重回帰分析を行う。
まずは以下のコードを実行してデータフレーム「data_extra4」を作成する。
set.seed(0)
distance = runif(100, min=1, max=10) #駅からの所要時間（分）
size = rnorm(100, mean=50, sd=10) #店の広さ（平方メートル）
kids_menu = sample(c(0, 1), 100, replace=TRUE) #お子様ランチの有無
sales = (-5)*distance+15*kids_menu+rnorm(100, mean=200, sd=20) #売上（架空の計算式）
data_extra4 = data.frame(distance=distance, size=size, kids_menu=kids_menu, sales=sales)
head(data_extra4, 10)## distance size kids_menu sales
## 1 9.070275 52.66137 1 148.7343
## 2 3.389578 46.23297 0 165.1279
## 3 4.349115 74.41365 0 203.6422
## 4 6.155680 42.04661 0 181.0984
## 5 9.173870 49.45123 0 169.6433
## 6 2.815137 52.50141 0 217.0717
## 7 9.085507 56.18243 0 147.2644
## 8 9.502077 48.27376 1 183.8207
## 9 6.947180 27.76100 1 179.0514
## 10 6.662026 37.36386 1 171.6623このデータフレームには、とある地域に存在する100個のレストランのデータが記録されている（head関数によって冒頭10行分のデータのみが表示されている）。
1列目は駅からの徒歩での所要時間、2列目は店の広さ、3列目はお子様ランチメニューの有無（0＝なし／1＝あり）、4列目は1ヶ月の売上が格納されている。
ここでは、目的変数に売上、説明変数にそのほかの3つの変数を設定し、重回帰分析を行いたい。以下のような式となる。
売上 = \(\beta_0\) + \(\beta_1\) 駅からの所要時間 + \(\beta_2\) 店の広さ + \(\beta_3\) お子様ランチの有無 + 誤差項
回帰分析の注目ポイントは、距離や広さなどの量的変数のみならず、お子様ランチメニューの有無という質的変数についても説明変数に組み込める点である。このような2値の質的変数は、0と1という仮の数値（ダミー変数と呼ぶ）を割り当てる。
コードの基本構造は単回帰分析と一緒で、波線の右側に3つの変数をプラスマークで連結して並べるだけで分析できる。
model3 = lm(sales ~ distance + size + kids_menu, data=data_extra4)
summary(model3)##
## Call:
## lm(formula = sales ~ distance + size + kids_menu, data = data_extra4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -60.250 -15.157 -0.454 10.882 45.648
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 182.3267 13.3242 13.684 < 2e-16 ***
## distance -3.6708 0.8718 -4.211 5.74e-05 ***
## size 0.2561 0.2406 1.064 0.2898
## kids_menu 12.5492 4.2901 2.925 0.0043 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 20.96 on 96 degrees of freedom
## Multiple R-squared: 0.2027, Adjusted R-squared: 0.1778
## F-statistic: 8.138 on 3 and 96 DF, p-value: 6.966e-05この結果より、距離とお子様ランチの有無の回帰係数が有意であることがわかる。
また、重回帰分析では、交互作用項の説明変数への投入も可能である。たとえば、お子様ランチの有無が売上に影響するのが広い店に限定されるという仮説が立てられるのであれば、交互作用項を投入すべきかもしれない。交互作用項はコロンを用いて表現する。たとえば「店の広さとお子様ランチの交互作用」をモデルに組み込みたいなら、size:kids_menuという変数を追加すればよい。
model4 = lm(sales ~ distance + size + kids_menu + size:kids_menu, data=data_extra4)
summary(model4)##
## Call:
## lm(formula = sales ~ distance + size + kids_menu + size:kids_menu,
## data = data_extra4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -58.827 -15.121 -0.671 10.987 47.779
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 176.1370 17.2073 10.236 < 2e-16 ***
## distance -3.5981 0.8840 -4.070 9.72e-05 ***
## size 0.3713 0.3145 1.180 0.241
## kids_menu 26.3274 24.4902 1.075 0.285
## size:kids_menu -0.2840 0.4969 -0.571 0.569
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.03 on 95 degrees of freedom
## Multiple R-squared: 0.2055, Adjusted R-squared: 0.172
## F-statistic: 6.142 on 4 and 95 DF, p-value: 0.0001935交互作用項を入れたところ、先ほどまでは見られていたお子様ランチの有無の有意性がみられなくなった。このように、変数を追加したり減らしたりすると、その他の説明変数の回帰係数や検定結果が変わるため、変数選択は非常に重要である。
変数をあれこれ調整しながら検定を繰り返すと統計不正につながるため、データを分析するときは事前に十分に解析方略を練っておく必要がある。
7.4 線形混合モデル（LMM）
本項で紹介する線形混合モデルは、それ単体で長い記事が書けるほど多くの知識が必要となるが、本教材ではあくまでも補足的紹介に留めるため、詳細な説明は第三者の書籍やWEBサイトなどを参照されたい。本項をなぞっただけでは理解は難しいだろう。とくに「久保（2012）データ解析のための統計モデリング入門」は必読である。
ここまでの分析は、参加者ごとのばらつきや刺激ごとの違いなどを無視した線形モデルであった（たとえば回帰分析では誤差項として括ってしまっていた）。実際には、参加者・刺激・観測値ごとに、ばらつき度合い・切片・傾きが異なるはずである。これらをランダム変数と呼び、ランダム変数を無視して解析することは、アルファエラー・ベータエラーが崩れる要因となる。
そこで、単なる線形モデルから発展して、線形混合モデル（LMM：Liner Mixed Model）を使用してモデリングすることを考える。線形混合モデルでは、目的変数を変動させる要因として、固定効果とランダム効果の2種類を設定する。それぞれの意味は、以下の通りである（一部、久保（2012）より引用）。重回帰分析などの線形モデルでは、固定効果のみを説明変数に組み込んでいた。線形混合モデルは、固定効果とランダム効果が混合しているという意味で、混合という名前が使われている。
| 意味 | 例 | |
|---|---|---|
| 固定効果 | サンプル全体の特徴≒データのかなり広い範囲を説明する大域的な効果 | 実験条件 |
| ランダム効果 | 対象ごとの特徴≒データのごく一部だけを担当している局所的な効果 | 参加者の個人差／刺激ごとの特性／観測値ごとの特性 |
まずは、線形混合モデリングを行うために、「lme4」パッケージと「lmerTest」パッケージをインストールする。下記のコードをコンソール画面（左下の画面）で実行すること。lme4パッケージは、線形混合モデルやこの後登場する一般化線形混合モデルを実行するための関数を含んでいる。lmerTestパッケージは、lme4パッケージで出力した回帰係数のp値を表示するために使用される。
install.packages("lme4")
install.packages("lmerTest")パッケージをインストールできたら、library関数で2つのパッケージを読み込む。
library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## step続いて、今回分析するデータフレーム「data_extra5」を作成する。
participant_id = rep(1:20, each=6)
stimulus_id = rep(c("m1", "m2", "m3", "o1", "o2", "o3"), times=20)
university = rep(c("keio", "kyoto", "tokyo", "waseda"), each=6, length.out=length(participant_id))
condition = ifelse(stimulus_id %in% c("m1", "m2", "m3"), "m", "o")
set.seed(123)
reaction_time = rnorm(120, mean=500, sd=100)
data_extra5 = data.frame(
participant_id = participant_id,
stimulus_id = stimulus_id,
university = university,
condition = condition,
reaction_time = reaction_time
)
head(data_extra5, 20)## participant_id stimulus_id university condition reaction_time
## 1 1 m1 keio m 443.9524
## 2 1 m2 keio m 476.9823
## 3 1 m3 keio m 655.8708
## 4 1 o1 keio o 507.0508
## 5 1 o2 keio o 512.9288
## 6 1 o3 keio o 671.5065
## 7 2 m1 kyoto m 546.0916
## 8 2 m2 kyoto m 373.4939
## 9 2 m3 kyoto m 431.3147
## 10 2 o1 kyoto o 455.4338
## 11 2 o2 kyoto o 622.4082
## 12 2 o3 kyoto o 535.9814
## 13 3 m1 tokyo m 540.0771
## 14 3 m2 tokyo m 511.0683
## 15 3 m3 tokyo m 444.4159
## 16 3 o1 tokyo o 678.6913
## 17 3 o2 tokyo o 549.7850
## 18 3 o3 tokyo o 303.3383
## 19 4 m1 waseda m 570.1356
## 20 4 m2 waseda m 452.7209本データは、視覚探索課題の反応時間を記録した20名分のデータフレームである。ターゲット刺激がMの条件とターゲット刺激がOの条件の間で、反応時間に違いがあるかどうかを検証することを目的としている。
- 1列目：参加者ID
- 2列目：刺激番号（m1, m2, m3, o1, o2, o3）
- 3列目：データ収集をした大学（慶應大、京大、東大、早稲田大）
- 4列目：実験条件（m, o）
- 5列目：反応時間（単位=ms）
このようなデータを分析するとき、これまでの考え方であれば、次のようなデータ構造に集計をして、対応のあるt検定を行っていたはずである。このデータ構造では、1行あたり1人分の参加者のデータを含んでおり、m条件とo条件のそれぞれの条件における反応時間平均値を集計している。
data_extra6 = data_extra5 %>%
group_by(participant_id, condition) %>%
summarize(mean_reaction_time = mean(reaction_time), .groups = "drop") %>%
pivot_wider(names_from = condition, values_from = mean_reaction_time, names_prefix = "mean_")
data_extra6## # A tibble: 20 × 3
## participant_id mean_m mean_o
## <int> <dbl> <dbl>
## 1 1 526. 564.
## 2 2 450. 538.
## 3 3 499. 511.
## 4 4 472. 434.
## 5 5 451. 509.
## 6 6 534. 580.
## 7 7 506. 457.
## 8 8 570. 434.
## 9 9 532. 543.
## 10 10 491. 531.
## 11 11 485. 440.
## 12 12 547. 475.
## 13 13 487. 484.
## 14 14 502. 522.
## 15 15 540. 542.
## 16 16 559. 504.
## 17 17 616. 451.
## 18 18 448. 417.
## 19 19 499. 464.
## 20 20 531. 416.実際に対応のあるt検定を行うとすれば、下記のようになる。
t.test(data_extra6$mean_m, data_extra6$mean_o, paired=TRUE)##
## Paired t-test
##
## data: data_extra6$mean_m and data_extra6$mean_o
## t = 1.454, df = 19, p-value = 0.1623
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -9.449576 52.449792
## sample estimates:
## mean difference
## 21.50011しかし、このように集計をして分析をすることは、以下の要素を潰してしまうことを意味する。これらの要素を潰した集計データを用いて検定を行うことは、アルファエラーの増大に繋がることが知られている。
- 参加者ごとの、ベースライン（切片）と実験条件の影響の受けやすさ（傾き）
- 刺激ごとの、ベースライン（切片）と実験条件の影響の受けやすさ（傾き）
- データ収集地ごとの、ベースライン（切片）と実験条件の影響の受けやすさ（傾き）
たとえば参加者でいえば、AさんとBさんは2人とも平均反応時間は500 msだが、Aさんの分散は300 msで、Bさんの分散は10 msかもしれない。このような分散情報が、集計することによって消失してしまう。このように、参加者や刺激やデータ収集地にはそれぞれ特性があり、その特性を潰さずに、できる限りリッチな情報量を維持したまま分析したい。
線形混合モデリングは、固定効果である実験条件（M条件とO条件）に加えて、これらの要素をランダム効果として投入することで、より正確なモデルを作成しようとする試みである。
実際に、参加者の個人差、刺激ごとの特性、観測値（大学）ごとの特性について、ランダム効果としてモデルに組み込んでみたい。線形混合モデリングには、lmer関数を用いる。基本的な変数投入方法は重回帰分析と同じで、波線の左で目的変数、波線の右で説明変数を指定する。ランダム効果は、たとえば参加者の個人差をモデルに組み込む場合、以下のように記述する。よく使うのは上の2つである。
- ランダム切片効果とランダム傾き効果の両方：
(condition|participant_id) - ランダム切片効果のみ：
(1|participant_id) - ランダム傾き効果：
(0 + condition|participant_id)
今回は、参加者の個人差と刺激ごとの特性は切片と傾きの両方を投入し、観測値ごとの特性は切片のみ投入してみる。
model5 = lmer(reaction_time ~ condition + (condition|participant_id) + (condition|stimulus_id) + (1|university), data=data_extra5)## boundary (singular) fit: see help('isSingular')## Warning: Model failed to converge with 1 negative eigenvalue: -1.6e-03summary(model5)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## reaction_time ~ condition + (condition | participant_id) + (condition |
## stimulus_id) + (1 | university)
## Data: data_extra5
##
## REML criterion at convergence: 1402.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.48642 -0.67304 -0.04892 0.65399 2.39487
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## participant_id (Intercept) 0.000e+00 0.000000
## conditiono 9.869e-06 0.003142 NaN
## stimulus_id (Intercept) 0.000e+00 0.000000
## conditiono 1.434e+01 3.786216 NaN
## university (Intercept) 0.000e+00 0.000000
## Residual 7.945e+03 89.137136
## Number of obs: 120, groups: participant_id, 20; stimulus_id, 6; university, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 512.29 11.51 116.00 44.518 <2e-16 ***
## conditiono -21.50 16.42 7.60 -1.309 0.229
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## conditiono -0.701
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')結果の見方は重回帰分析と基本的には同じで、主要な相違点はランダム効果に関する欄が追加されていることくらいである。
なお、ここまでの分析はすべて正規性を仮定してきたが、正規性が仮定されないデータを分析する場合は、一般化線形混合モデル（GLMM：Generalized Linear Mixed Model）を使用する。GLMMについては、本教材の域を超えるため説明しない。
- 参考文献
- 久保（2012）データ解析のための統計モデリング入門
- 井関「心理学者は反応時間をどう分析するか」
- 清水「社会心理学とGlmm」
8 練習問題の解答
8.1 練習問題1の解答
自力で解いてから、本欄で答え合わせをしてみてください。
#データフレーム読み込み
df1p = read.xlsx("data1p.xlsx")
df1p## attention_task_reaction_time aq_score
## 1 334.6619 6
## 2 317.3599 3
## 3 319.2476 10
## 4 365.7743 7
## 5 271.2179 4
## 6 342.0688 6
## 7 339.9791 9
## 8 306.7227 2
## 9 300.4734 6
## 10 379.7396 5
## 11 349.5109 5
## 12 342.1110 7
## 13 328.3402 4
## 14 331.6492 3
## 15 329.1320 7
## 16 334.4032 7
## 17 311.9025 2
## 18 338.2931 5
## 19 317.1956 4
## 20 344.3361 1#グラフによる可視化
ggplot(df1p, aes(x=aq_score, y=attention_task_reaction_time))+
geom_point()+
geom_smooth(method="lm")## `geom_smooth()` using formula = 'y ~ x'
#グラフの保存
ggsave("fig1p.png", dpi=300, width=1200, height=1200, units="px", device="png")## `geom_smooth()` using formula = 'y ~ x'#記述統計量
mean(df1p$attention_task_reaction_time)## [1] 330.206sd(df1p$attention_task_reaction_time)## [1] 23.48916mean(df1p$aq_score)## [1] 5.15sd(df1p$aq_score)## [1] 2.345769#無相関検定
cor.test(df1p$attention_task_reaction_time, df1p$aq_score, method="pearson")##
## Pearson's product-moment correlation
##
## data: df1p$attention_task_reaction_time and df1p$aq_score
## t = 0.98864, df = 18, p-value = 0.3359
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2396423 0.6083688
## sample estimates:
## cor
## 0.2269453注意課題の反応時間（M = 330.2, SD = 23.5）と自閉症スペクトラム指数（M = 5.15, SD = 2.35）の関連を調べるために、Pearsonの相関係数を算出して無相関検定を行ったところ、有意な相関はみられなかった（r = .23, t(18) = 0.99, p = .336, 95% CI [-0.24, 0.61]）。
8.2 練習問題2の解答
自力で解いてから、本欄で答え合わせをしてみてください。
#データフレーム読み込み
df2p = read.xlsx("data2p.xlsx")
df2p## desk_condition calculation_task_score
## 1 normal 80
## 2 normal 74
## 3 normal 81
## 4 normal 70
## 5 normal 73
## 6 normal 73
## 7 normal 60
## 8 normal 80
## 9 normal 70
## 10 normal 80
## 11 standing 70
## 12 standing 70
## 13 standing 77
## 14 standing 56
## 15 standing 58
## 16 standing 69
## 17 standing 65
## 18 standing 78
## 19 standing 66
## 20 standing 61#グラフによる可視化
ggplot(df2p, aes(x=desk_condition, y=calculation_task_score))+
geom_boxplot()+
geom_jitter(width=0.2, height=0, alpha=0.5)+
labs(x="Desk Condition", y="Score of Calculation Task")
#グラフの保存
ggsave("fig2p.png", dpi=300, width=1200, height=1200, units="px", device="png")
#データ抽出
df2p_normal = df2p %>% filter(desk_condition=="normal") %>% pull(calculation_task_score)
df2p_standing = df2p %>% filter(desk_condition=="standing") %>% pull(calculation_task_score)
#記述統計量
mean(df2p_normal)## [1] 74.1sd(df2p_normal)## [1] 6.556591mean(df2p_standing)## [1] 67sd(df2p_standing)## [1] 7.348469#t検定
t.test(df2p_normal, df2p_standing)##
## Welch Two Sample t-test
##
## data: df2p_normal and df2p_standing
## t = 2.2798, df = 17.771, p-value = 0.0352
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.5510402 13.6489598
## sample estimates:
## mean of x mean of y
## 74.1 67.0#効果量d
cohens_d(df2p_normal, df2p_standing)## Cohen's d | 95% CI
## ------------------------
## 1.02 | [0.07, 1.94]
##
## - Estimated using pooled SD.机の種類によって計算課題の得点に違いがあるかを調べるために、普通の机（M = 74.1, SD = 6.6）とスタンディングデスク（M = 67.0, SD = 7.3）のデータに対して対応のないwelchのt検定を実施したところ、普通の机を使用したほうがスタンディングデスクを使用したときよりも計算課題の得点が有意に高いことが明らかとなった（t(17.77) = 2.28, p = .035, d = 1.02, 95% CI [0.55, 13.65]）。
8.3 練習問題3の解答
自力で解いてから、本欄で答え合わせをしてみてください。
#データフレーム読み込み
df3p = read.xlsx("data3p.xlsx")
df3p## class_style motivation test_score
## 1 in_person high 76
## 2 in_person high 70
## 3 in_person high 94
## 4 in_person high 81
## 5 in_person high 73
## 6 in_person high 86
## 7 in_person high 74
## 8 in_person high 75
## 9 in_person low 62
## 10 in_person low 64
## 11 in_person low 76
## 12 in_person low 61
## 13 in_person low 60
## 14 in_person low 65
## 15 in_person low 75
## 16 online high 73
## 17 online high 76
## 18 online high 80
## 19 online high 71
## 20 online high 80
## 21 online high 70
## 22 online high 75
## 23 online high 72
## 24 online high 88
## 25 online low 41
## 26 online low 44
## 27 online low 45
## 28 online low 59
## 29 online low 48
## 30 online low 43#グラフを作図するためのデータフレーム作成
df3p_plot = df3p %>%
group_by(class_style, motivation) %>%
summarize(mean = mean(test_score),
se = std.error(test_score))## `summarise()` has grouped output by 'class_style'. You can override using the
## `.groups` argument.df3p_plot## # A tibble: 4 × 4
## # Groups: class_style [2]
## class_style motivation mean se
## <chr> <chr> <dbl> <dbl>
## 1 in_person high 78.6 2.82
## 2 in_person low 66.1 2.50
## 3 online high 76.1 1.91
## 4 online low 46.7 2.64#グラフによる可視化
ggplot(df3p_plot, aes(x=class_style, y=mean, colour=motivation, group=motivation))+
geom_point()+
geom_errorbar(aes(ymin=mean-se,ymax=mean+se), width=0.1)+
geom_line()+
labs(x="Style", y="Test Score")
#グラフの保存
ggsave("fig3p.png", dpi=300, width=1200, height=1200, units="px", device="png")
#分散分析
anovakun(df3p, "ABs", Style=c("in_person","online"), Motivation=c("high","low"), geta=T)##
## [ ABs-Type Design ]
##
## This output was generated by anovakun 4.8.9 under R version 4.5.1.
## It was executed on Sun Oct 5 04:18:46 2025.
##
##
## << DESCRIPTIVE STATISTICS >>
##
## ----------------------------------------------
## Style Motivation n Mean S.D.
## ----------------------------------------------
## in_person high 8 78.6250 7.9631
## in_person low 7 66.1429 6.6189
## online high 9 76.1111 5.7325
## online low 6 46.6667 6.4704
## ----------------------------------------------
##
##
## << ANOVA TABLE >>
##
## == This data is UNBALANCED!! ==
## == Type III SS is applied. ==
##
## ----------------------------------------------------------------------------
## Source SS df MS F-ratio p-value G.eta^2
## ----------------------------------------------------------------------------
## Style 886.2402 1 886.2402 19.5446 0.0002 *** 0.4291
## Motivation 3221.6390 1 3221.6390 71.0482 0.0000 *** 0.7321
## Style x Motivation 527.3117 1 527.3117 11.6290 0.0021 ** 0.3090
## Error 1178.9544 26 45.3444
## ----------------------------------------------------------------------------
## Total 5419.3667 29 186.8747
## +p < .10, *p < .05, **p < .01, ***p < .001
##
##
## << POST ANALYSES >>
##
## < SIMPLE EFFECTS for "Style x Motivation" INTERACTION >
##
## ---------------------------------------------------------------------------------
## Source SS df MS F-ratio p-value G.eta^2
## ---------------------------------------------------------------------------------
## Style at high 26.7655 1 26.7655 0.5903 0.4492 ns 0.0222
## Style at low 1225.5018 1 1225.5018 27.0265 0.0000 *** 0.5097
## Motivation at in_person 581.6679 1 581.6679 12.8278 0.0014 ** 0.3304
## Motivation at online 3121.1111 1 3121.1111 68.8312 0.0000 *** 0.7258
## Error 1178.9544 26 45.3444
## ---------------------------------------------------------------------------------
## +p < .10, *p < .05, **p < .01, ***p < .001
##
## output is over --------------------///まず、各実施形態のモチベーションの高さごとのテストの得点の平均値を計算したところ、対面方式においてモチベーションが高い人は78.6点（SD = 8.0）、対面方式においてモチベーションが低い人は66.1点（SD = 6.6）、オンライン方式においてモチベーションが高い人は76.1点（SD = 5.7）、オンライン方式においてモチベーションが低い人は46.7点（SD = 6.5）であった。
実施形態とモチベーションの高さを要因とする参加者間2要因分散分析を実施したところ、実施形態の主効果（F(1, 26) = 19.54, p < .001, ηG2 = 0.43）、モチベーションの高さの主効果（F(1, 26) = 71.05, p < .001, ηG2 = 0.73）、実施形態とモチベーションの高さの交互作用（F(1, 26) = 11.63, p = .002, ηG2 = 0.31）がそれぞれ有意となった。
有意な交互作用が見られたため、単純主効果検定を行った。その結果、モチベーションが低い人における実施形態の単純主効果（F(1, 26) = 27.03, p < .001, ηG2 = 0.51）、対面方式におけるモチベーションの高さの単純主効果（F(1, 26) = 12.83, p = .001, ηG2 = 0.33）、オンライン方式におけるモチベーションの高さの単純主効果（F(1, 26) = 68.83, p < .001, ηG2 = 0.73）がそれぞれ有意となった。モチベーションが高い人における実施形態の有意な単純主効果はみられなかった（F(1, 26) = 0.59, p = .449, ηG2 = 0.02）。
【この結果から言えそうなことを日常言語で説明】モチベーションが高い人は実施形態の影響を受けにくく、対面方式でもオンライン方式でも十分に学べる。しかし、モチベーションが低い人は対面方式よりもオンライン方式のほうがテストの成績が悪く、オンライン方式だと十分に学べないことが示唆される。なお、実施形態によらず、モチベーションが高い人が低い人よりもテストの成績がよい。