これから研究を始める高校生と指導教員のために 第2版
探究活動と課題研究の進め方･論文の書き方･口頭とポスター発表の仕方
(これ研)
副読文書
VIII-2. 一般化線形モデル入門
説明変数が範疇の場合の、各条件間での平均値の差の解析
確率分布が正規分布の場合
発射角度に依存した、ペットボトルロケットの平均飛行距離の差の解析
(｢これ研｣本文の第3部2.3節；p. 79)
｢これから研究を始める高校生と指導教員のために 第2版；探究活動と課題研究の進め方･論文の書き方･口頭とポスター発表の仕方｣(これ研)の第3部6.6節｢一般化線形モデルの奨め｣(p. 133)で紹介している一般化線形モデルを解説します。本章では、説明変数が範疇の場合の、各条件間での平均値の差の解析の仕方を説明します。VIII-3章での解析と異なり、確率分布が正規分布の場合です。発射角度に依存した、ペットボトルロケットの飛行距離の平均値の差の検定を例に用います(｢これ研｣の本文第3部2.3節；p. 79)。25°, 30°, 35°, 40°, 45°という発射角度間で飛行距離に差があるのかどうかを検定します。複数の飲料間でのカフェイン量の比較やメダカの体長･産卵数の地域間での比較や就寝前および起床後における単語記憶の比較；データ間に対応がない場合もまったく同じ方法で検定することできます。｢これ研｣の第3部第6章｢検定に挑戦しよう｣(p. 118)を読み、検定とは何かを理解しておいて下さい。
RStudioを起動して下さい。起動方法の詳しい説明は、RStusioの起動の仕方を参照して下さい。そして、作業ディレクトリの指定の説明に従って作業ディレクトリを指定します。
setwd("/Users/sakai/Documents/書籍等原稿/これ研2版/課題研究解析")
# パソコン内での作業ディレクトリの位置がわかっている場合はこの方法が便利です。作業ディレクトリの位置の知り方は作業ディレクトリの表示を参照して下さい。
# あなたの作業ディレクトリの位置を""で囲んで書きます。
# この命令文を実行しておきます。
解析と作図には、emmeansというものとtidyverseというものを用います。emmeans, tidyverseをインストールしてありますか？ まだならば、下記の説明に従ってインストールして下さい。
install.packages("emmeans")
install.packages("tidyverse")
# RStudioを起動しこの命令文を実行します。emmeans, tidyverseを""で囲みます。作図の準備も参照して下さい。
インストールしたら、emmeans, tidyverseをRに読み込みます。
library(emmeans)
library(tidyverse)
# この命令文を実行しておきます。Rstudio起動後、一度だけ実行すればよいです。Rstudioを終了して再び起動したときは再実行する必要があります。
解析では、データフレームの中の特定のデータ列を指定することを行います。指定の方法です。
データフレーム名$データ列名
# データフレーム名を書き、$を挟んで、指定したいデータ列名を書きます。
# 実行例
d$Angle.degree
d$Flying.distance.m
# データフレームdに入っているデータ列Angle.degree, flying.distance.mを指定します。
1. 解析に用いるcsvファイルのRへの読み込み
ペットボトルロケットの発射距離のデータがcsvファイル「Rocket.csv」(画面1)に入っています(架空データです)。発射角度（Angle.degree）25, 30, 35, 40, 45度におけるロケットの飛行距離(Flying.distance.m)のデータです。各角度で複数回の飛行実験を行っています。このファイルをダウンロードして下さい。
画面1
そして、ダウンロードしたファイルをあなたの作業ディレクトリに入れて下さい(｢作業ディレクトリの指定｣を参照)。
ご自身のデータを用いる場合は、そのデータが入ったcsvファイルを作業ディレクトリに入れて下さい。Excelで作ったファイルをcsvファイルに変換する方法は、Excelで作った解析用ファイルのcsv形式での保存を参照して下さい。
csvファイルをRに読み込み、データフレームに格納します。
d <- read.csv("Rocket.csv")
# csvファイル｢Rocket.csv｣を読み込んでデータフレームdに格納します。ファイル名を""で囲みます。ファイル名の拡張子｢.csv｣も忘れずに書きます。
# データフレームの名称(この例ではd)はお好みのものでよいです。
解析では、発射角度を範疇として扱います。｢25°｣｢30°｣｢35°｣ …… という発射角度を範疇とみなし、各発射角度の飛行距離を比較します。しかし、csvファイル｢Rocket.csv｣には、25, 30, 35, ……という数値で入っています。そのためRは、発射角度を数値として認識しています。以下の命令を実行し範疇として認識させます。
d$Angle.degree <- as.factor(d$Angle.degree)
# データフレームd中のデータ列Angle.degreeに発射角度が入っています。
# as.factorと命令して範疇と認識させ、d$Angle.degreeの中身を上書きするためにd$Angle.degreに再格納します。
扱うデータが、｢25°｣｢30°｣｢35°｣ (°が付いている)や｢Tohoku｣｢Kanto｣といった文字である場合はこの命令は不要です。何もしなくても範疇としてRは認識します。
d$Angle.degree <- factor(d$Angle.degree, levels=c("25", "30", "35", "40", "45"))
# 発射角度の並び順を指定します。データフレームd中のデータ列Angle.degreeに発射角度が入っています。並べたい順番に発射角度を書きます。数値なので、""で囲っても囲わなくてもよいです(文字の場合は必ず囲みます)。
# d$Angle.degreeの中身を上書きするためにd$Angle.degreeに再格納します。
# この命令文を実行しないと、数値は小さい順に、範疇(文字)はアルファベット順(数字の場合は小さい順)に並びます(数字を小さい順に並べているので、この命令文は実は不要です)。
2. 一般化線形モデルによる解析
一般化線形モデル実行の手順に従って一般化線形モデルを実行します。
2.1. 適用する確率分布の指定
適用する確率分布を指定します。確率分布とはと確率分布の説明を読んで下さい。ペットボトルロケットの飛行距離は連続する数値なので、正規分布かガンマ分布のどちらかです。下の説明に従い、各発射角度における飛行距離のヒストグラムを描いてみましょう。
画面2
左右対称の釣り鐘方に近いので、確率分布を正規分布にします。
# 命令文が続く場合には+で繋げ、最後の命令文の後には+を付けません。
# 必須命令文
ggplot(d, aes(x = Flying.distance.m)) +# データフレームdを指定します。
# x = Flying.distance.mとして、ヒストグラムを描く飛行距離のデータ列Flying.distance.mをx軸(横軸)に指定します。xは小文字です。
# 末尾の+を忘れないで下さい。
# ヒストグラムにおける、作図に用いるデータフレームとデータの指定
geom_histogram() +
# 描く図としてヒストグラムを指定します。
# 末尾の+を忘れないで下さい。
facet_wrap(~ Angle.degree, ncol = 1) +
# 複数の図を並べて描く場合はfacet_wrap()と命令します。
# データ列Angle.degreeに入っている発射角度ごとにヒストグラムを描きます。Angle.degreeの前の~を忘れずに書いて下さい。
# 図の並べ方をncol(またはnrow)で指定します。ncol=1にすると縦1列に並びます。
# 末尾の+を忘れないで下さい。ただし、以降の命令文を省略する場合は+を付けず、｢facet_wrap(~ Angle.degree, ncol = 1)｣とします。
# 省略してもよい命令文
# 以降の命令文を省略するとデフォルトで自動で描きます。
labs(x = "Flying distance (m)", y = "Frequency") +
# x軸(横軸)の名称をFlying distance (m)に、y軸(縦軸)の名称をFrequencyにします。名称を""で囲みます。x, yは小文字です。
scale_y_continuous(limits = c(0, 3)) +
# y軸(縦軸)を描く範囲を0から3に指定します。ヒストグラムの場合は、y軸の最小値を必ず0にします。
theme_classic() +
# 図の背景色を白にします。
theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 9),
axis.title.y = element_text(size = 18),
axis.text.y = element_text(size = 9)
)
# x軸(横軸)･y軸(縦軸)の名称(axis.title.x, axis.title.y)の文字の大きさを18に、x軸･y軸の目盛り(axis.text.x, axis.text.y)の文字の大きさを9に指定します。
ヒストグラムの描き方の詳しい説明はメダカの体長のヒストグラムを参照して下さい。図の軸の説明の文字を日本語にしたい場合は、その方法も、メダカの体長のヒストグラムを参照して下さい。
2.2. 連結関数は比例(identity)に
適用する連結関数を指定します。説明変数の値と応答変数の平均値の関係性：連結関数を用いた解析と適用する連結関数の指定の説明を読んで下さい。説明変数が範疇の場合は必ず比例(identity)を指定します。説明変数に便宜的に0, 1という数値を与え、それらを直線で結ぶ回帰式を計算するからです。この回帰式は比例の式です。
2.3. 一般化線形モデルの実行
一般化線形モデルの命令文を実行します。一般化線形モデルの命令文の実行の説明を読んで下さい。説明変数が範疇で、比較する対象が｢25°｣｢30°｣｢35°｣｢40°｣｢45°｣の5つあります。なので、1. 一般化線形モデルの命令文および2. 追加命令文を実行します。
r <- glm(formula = Flying.distance.m ~ Angle.degree, family = gaussian(link = identity), data = d)
summary(r)
# 一般化線形モデルの命令文を実行して、その結果をデータフレームrに格納します。
# 解析するデータフレーム；d
# 応答変数；飛行距離(Flying.distance.m)
# 説明変数；発射角度(Angle.degree)
# 確率分布；正規分布(gaussian)
# 連結関数；比例(identity)
# summary(r)を実行して、解析結果rの詳細を表示させます。
# 左上の赤枠内にあるのが、各発射角度における平均飛行距離を示す数値です(ただし30°以降では、平均飛行距離そのもではありません)。各発射角度の平均飛行距離は以下のようになります。
25°の平均飛行距離 = (Intercept) = 31.6125 m
30°の平均飛行距離 = Angle.degree30 + (Intercept) = 5.5375 + 31.6125 = 37.15 m
35°の平均飛行距離 = Angle.degree35 + (Intercept) = 8.9625 + 31.6125 = 40.575 m
40°の平均飛行距離 = Angle.degree40 + (Intercept) = 6.0375 + 31.6125 = 37.65 m
45°の平均飛行距離 = Angle.degree45 + (Intercept) = 3.4375 + 31.6125 = 35.05 m
連結関数が比例なので、y = ax + b (y；応答変数の平均値 x；説明変数の値 a；傾きの値 b；切片の値)という形で応答変数の平均値を示します。このとき、データの並び順の先頭に出てくるもの(この例では25°)の平均値を切片(Intercept)として扱います。つまり切片b(Intercept)の値が25°の平均飛行距離です。他のもの(30°~45°)の平均値は、x = 1のときの値、すなわちy = a + bとして表現します。このとき、25°~45°の平均飛行距離すべてを1つの回帰式に組み込むのではなく、｢25°と30°の平均飛行距離の回帰式｣｢25°と35°の平均飛行距離の回帰式｣｢25°と40°の平均飛行距離の回帰式｣｢25°と45°の平均飛行距離の回帰式｣という別々の回帰式を計算します。そして、各回帰式にx = 1を当てはめた値y = a + bを求めます。左の赤枠内の2行目以降が、その発射角度の回帰式のa(傾き)の値です。これに、25°の平均飛行距離(切片bの値)を加えた値がその発射角度における平均飛行距離です。
# 右上の赤枠内にあるのがP値です。並び順の先頭(25°)に関しては、切片(Intercept)の値(bの値)が有意に0と異なるのかどうかを示しています。他のもの(30°~45°)に関しては、切片の値との差(aの値)が有意に0と異なるのかどうかを示しています。｢e-｣とあるのは10のマイナス何乗かということです。｢2e-16｣は｢2 × 10- 16｣で、｢6.51e-08｣は｢6.51 × 10- 8｣です。
25°の平均飛行距離(Intercept)のP値；2 × 10- 16 有意に0と異なる。
25°と30°の平均飛行距離の差のP値；6.51 × 10- 8 有意に0と異なる。つまり、25°と30°の平均飛行距離に有意な差がある。平均飛行距離は、25°が31.6125mで30°が37.15mであり、30°の方が有意に長い。
25°と35°の平均飛行距離の差のP値；5.98 × 10- 13 有意に0と異なる。つまり、25°と35°の平均飛行距離に有意な差がある。平均飛行距離は、25°が31.6125mで35°が40.575mであり、35°の方が有意に長い。
25°と40°の平均飛行距離の差のP値；1.05 × 10- 8 有意に0と異なる。つまり、25°と40°の平均飛行距離に有意な差がある。平均飛行距離は、25°が31.6125mで40°が37.65mであり、40°の方が有意に長い。
25°と45°の平均飛行距離の差のP値；0.000158 有意に0と異なる。つまり、25°と45°の平均飛行距離に有意な差がある。平均飛行距離は、25°が31.6125mで45°が35.05mであり、45°の方が有意に長い。
# 下の赤枠内の数値がAICです。AICによる、確率分布と連結関数の選択を参照して下さい。この値が一番小さい、確率分布と連結関数の組合せを採用します。ガンマ分布と比例の連結関数の組合せで実行するとAICは162.4になります。正規分布と比例の連結関数の場合は158.95なので、正規分布の方が良いということです。
# 上記の解析では、30°, 35°, 40°, 45°間での平均飛行距離の差を検定できていません。そのため、2. 追加命令文を実行してこれらの間での差を検定します。
# 比較の対象が2つだけの場合は2. 追加命令文は不要で、これで解析は終わりです。たとえば、25°と30°のデータしかなく、この2つの間での差のみを検定したい場合などです。上記の解析で両者の差(30°と、切片(Intercept)である25°の差)が検定できているので、これ以上の解析は必要ありません。
pairs(emmeans(r, "Angle.degree"), adjust = "holm")
# 1. 一般化線形モデルの命令文の実行結果がデータフレームrに格納されています。
# データフレームrと説明変数Angle.degreeを指定します。Angle.degreeを""で囲みます。
# 発射角度25°, 30°, 35°, 40°, 45°間で、平均値の有意差を総当たりで検定します。切片となる発射角度(データの並び順で先頭になるもの)を自動的に入れ替え、その切片との比較をしてくれます。
# adjust="holm"を忘れずに書いて下さい。
# 実行結果です。左の赤枠内にあるのが、発射角度の各組合せでの平均飛行距離の差です。右の赤枠内にあるのが平均値の差のP値です。平均値に有意に差があるのかどうかを示しています。
# 発射角度25°と45°の差のP値が0.0009になっています。ところが1. 一般化線形モデルの命令文の実行の方では、このP値は0.000158です。25°と他の角度との差も実はP値が変わっています(< .0001と表記されているので、変わっているとわからないのですが)。変ですね。その理由を、多重比較における補正で説明しています。この説明を必ず読んで下さい。
# P値はすべて、この追加命令文によって得た値を採用して下さい。25°との比較のP値は1. 一般化線形モデルの命令文でも出ていますが、それらではなく、2. 追加命令文で得た値を採用します。
# 追加命令文
3. 多重比較における補正
同時に複数組の検定を行うことを多重比較といいます。たとえば上記の2. 追加命令文の実行では、発射角度25°, 30°, 35°, 40°, 45°間の総当たりで合計10組の検定を行っています。実は、こうした多重比較には考慮すべきことがあります。
｢下手な鉄砲も数撃ちゃ当たる｣で説明しましょう。100発中1発しか当たらない(命中率が0.01)鉄砲撃ちがいるとします。その鉄砲撃ちに1発だけ撃って貰ったところ、見事に命中しました。稀な現象が起きたので驚きます。今度は100発撃って貰いましょう。稀な現象はそうそう起きないので、実力通りに外しまくります。でも、100発も撃てば1発くらいは当たるでしょう。その、たまたま当たった1発を捉えて、｢稀な現象が起きた｣と驚くことはありません。命中率が0.01であっても、100発中1発以上当たる確率は1 - 0.99100 = 0.634ですから。｢たくさん撃てばどれかは当たる｣ことを無視して、たまたま当たった1発のことをことさら取り上げるのは無意味です。
同じことが検定にも当てはまります。ある1組だけの検定を行い、そのP値が0.01であったとします。P値が低いので稀な現象であるとして、母集団間に平均値に有意な差があると判定します。一方、100組の検定を同時に行った場合はどうでしょう。どの組合せでも本当は母集団間に差がなかったとしても、100組もデータ取りをすれば、たまたま、平均値に大きな差が出るものもあるでしょう。つまり、本当は差がなくとも、たまたまP値が0.01になる検定も出てきます。その、P値が0.01になった検定をことさら捉えて、｢0.01の確率でしか起きないことが起きた｣と判定するのは間違っています。100組の検定をすればどれかのP値が0.01になる確率を考慮していないからです。このように、複数組の検定を同時に行うと、どれかがたまたま有意になる確率が高まります。これが、多重比較において考慮すべきことです。
複数組の検定を同時に行う場合は、どれかがたまたま有意になる確率を補正する必要があります。そして、複数組の検定を同時に行ったことを考慮した上でのP値を計算します。こうして補正したP値は、同時に行った検定数に依存せず、本当にそのデータを得る確率を示したものです。
2. 追加命令文のadjust="holm"がこうした補正を行う命令文です。1. 一般化線形モデルの命令文の実行で得たP値は補正前のもので、2. 追加命令文の実行で得たP値が補正後のものです。この補正後のP値を採用して下さい。難しいので、具体的な補正方法は説明しません。論文やプレゼンでは、｢Holm法で補正した｣と書いて下さい。
4. 作図
検定結果を作図する方法を説明します。
画面3
ペットボトルロケットの飛行距離の平均と標準誤差の作図の説明に従って作図をして下さい。そしてそれをPowerPoint等に読む込んで下さい。読み込む方法は図の保存･コピーを参照して下さい。
P値が0.05以下の発射角度間では平均飛行距離に有意な差があると判定します。ここでいうP値は、比較の対象が3つ以上で複数組の検定を行った場合は補正後のもの(2. 追加命令文の実行で得たP値)、2つしかなく検定を1回しか行っていない場合は補正前のもの(1. 一般化線形モデルの命令文の実行で得たP値)です。そして以下のようにアルファベットを添えます(PowerPoint等に読み込んで書き込みます)。
平均に有意な差がないもの；同じ文字を添える
平均に有意な差があるもの；異なる文字を添える
論文では、図の説明文中に、｢同じ文字を添えたものには有意な差がない(P > 0.05)｣などと書きます。プレゼンでは口頭で説明しましょう。本解析の場合、2. 追加命令文の実行の結果、発射角度30°と40°の間には有意な差がなく、他の組合せにはすべて有意な差がありました。その場合は画面3のように、30°と40°に同じ文字を添え、他にはすべて異なる文字を添えます。慣れないと難しいかもしれません。いくつかの例を出しますので習得して下さい。｢有｣が有意差ありで、｢無｣が有意差なしの関係です。
画面4
検定結果が複雑だと、1つの平均値に複数の文字を添えて差の有無の範囲を示すことになります。この場合、結果の理解もちょっと面倒です。たとえば右下の場合、｢CとDに有意差なし｣｢DとEに有意差なし｣です。ならばCとEにも差がないと思ってしまいそうですが、｢CとEには有意差あり｣です。検定の結果を正確に理解して下さいね。