これから研究を始める高校生と指導教員のために　第2版
探究活動と課題研究の進め方･論文の書き方･口頭とポスター発表の仕方
(これ研)

副読文書
VIII-3. 一般化線形モデル入門
説明変数が範疇の場合の、各条件間での平均値の差の解析
確率分布がポアソン分布の場合
製造社間での、レーズンロールに入っている干し葡萄の平均数の差の解析

setwd("/Users/sakai/Documents/書籍等原稿/これ研2版/課題研究解析")

# あなたの作業ディレクトリの位置を""で囲んで書きます。

# この命令文を実行しておきます。

# 作業ディレクトリの指定

install.packages("emmeans")
install.packages("tidyverse")

# RStudioを起動しこの命令文を実行します。emmeans, tidyverseを""で囲みます。｢作図の準備｣も参照して下さい。

library(emmeans)
library(tidyverse)

# RStudioでこの命令文を実行しておきます。emmeans, tidyverseを""で囲みません。Rstudio起動後、一度だけ実行すればよいです。Rstudioを終了して再び起動したときは再実行する必要があります。

d <- read.csv("Raisin.csv")

# csvファイル｢Raisin.csv｣を読み込んでデータフレームdに格納します。ファイル名を""で囲みます。ファイル名の拡張子｢.csv｣も忘れずに書きます。

# データフレームの名称(この例ではd)はお好みのものでよいです。

# Excelで作った解析用ファイルのRへの読み込み

# 命令文が続く場合には+で繋げ、最後の命令文の後には+を付けません。

# 必須命令文

ggplot(d, aes(x = Number)) +

# データフレームdを指定します。

# x = Numberとして、ヒストグラムを描く干し葡萄数のデータ列Numberをx軸(横軸)に指定します。xは小文字です。

# 末尾の+を忘れないで下さい。

# ヒストグラムにおける、作図に用いるデータフレームとデータの指定

geom_histogram() +

# 描く図としてヒストグラムを指定します。

# 末尾の+を忘れないで下さい。

# ヒストグラムの指定

facet_wrap(~ Maker, ncol = 1) +

# データ列Makerに入っている製造社ごとにヒストグラムを描きます。それらの図を並べて示します。Makerの前の~を忘れずに書いて下さい。

# 図の並べ方をncol(またはnrow)で指定します。ncol=1にすると縦1列に並びます。

# 末尾の+を忘れないで下さい。ただし、以降の命令文を省略する場合は+を付けず、｢facet_wrap(~ Maker, ncol = 1)｣とします。

# 複数の図の並べ方の指定

# 省略してもよい命令文

# 以降の命令文を省略するとデフォルトで自動で描きます。

labs(x = "Number of raisins", y = "Frequency") +

# x軸(横軸)の名称をNumber of raisinsに、y軸(縦軸)の名称をFrequencyにします。名称を""で囲みます。x, yは小文字です。

# x軸･y軸の名称の指定

scale_y_continuous(limits = c(0, 15)) +

# y軸(縦軸)を描く範囲を0から15に指定します。ヒストグラムの場合は、y軸の最小値を必ず0にします。

# x軸･y軸の値の最小値･最大値の指定

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に指定します。

# x軸･y軸の文字の大きさと書き方の指定

r <- glm(formula = Number ~ Maker, family = poisson(link = identity), data = d)
summary(r)

# 一般化線形モデルの命令文を実行して、その結果をデータフレームrに格納します。

# 解析するデータフレーム；d
# 応答変数；干し葡萄の数(Number)
# 説明変数；製造社(Maker)
# 確率分布；ポアソン分布(poisson)
# 連結関数；比例(identity)

# summary(r)を実行して、解析結果rの詳細を表示させます。

# 左上の赤枠内にあるのが、各社における平均干し葡萄数を示す数値です(ただしB社以降では、平均干し葡萄数そのもではありません)。各社の平均干し葡萄数は以下のようになります。

A社の平均干し葡萄数 = (Intercept) = 13.9300
B社の平均干し葡萄数 = MakerB + (Intercept) = - 0.1100 + 13.9300 = 13.82
C社の平均干し葡萄数 = MakerC + (Intercept) = 2.5100 + 13.9300 = 16.44
D社の平均干し葡萄数 = MakerD + (Intercept) = 2.0700 + 13.9300 = 16.00

連結関数が比例なので、y = ax + b (y；応答変数の平均値　x；説明変数の値　a；傾きの値　b；切片の値)という形で応答変数の平均値を示します。このとき、データの並び順の先頭に出てくるもの(この例ではA社)の平均値を切片(Intercept)として扱います。つまり切片b(Intercept)の値がA社の平均干し葡萄数です。他社(B~D社)の平均値は、x = 1のときの値、すなわちy = a + bとして表現します。このとき、A~D社の平均干し葡萄数すべてを1つの回帰式に組み込むのではなく、｢A社とB社の平均干し葡萄数の回帰式｣｢A社とC社の平均干し葡萄数の回帰式｣｢A社とD社の平均干し葡萄数の回帰式｣という別々の回帰式を計算します。そして、各回帰式にx = 1を当てはめた値y = a + bを求めます。左上の赤枠内の2行目以降が、その社の回帰式のa(傾き)の値です。これに、A社の平均干し葡萄数(切片bの値)を加えた値がその社における平均干し葡萄数です。

#　右上の赤枠内にあるのがP値です。並び順の先頭(A社)に関しては、切片(Intercept)の値(bの値)が有意に0と異なるのかどうかを示しています。他のもの(B, C, D社)に関しては、切片の値との差(aの値)が有意に0と異なるのかどうかを示しています。｢e-｣とあるのは10のマイナス何乗かということです。｢2e-16｣は｢2 × 10^{- 16}｣で、｢5.25e-06｣は｢5.25 × 10^{- 6}｣です。

A社の平均干し葡萄数(Intercept)のP値；2 × 10^{- 16}　　有意に0と異なる
A社とB社の平均干し葡萄数の差のP値；0.834593　　有意に0と異ならない　つまり、A社とB社の平均干し葡萄数に差があるとは言えない
A社とC社の平均干し葡萄数の差のP値；5.25 × 10^{- 6}　　有意に0と異なる　つまり、A社とC社の平均干し葡萄数に有意な差がある
A社とD社の平均干し葡萄数の差のP値；0.000155　　有意に0と異なる　つまり、A社とD社の平均干し葡萄数に有意な差がある

#　下の赤枠内の数値がAICです。AICによる、確率分布と連結関数の選択を参照して下さい。確率分布を、正規分布やガンマ分布にするとAICが2175よりも大きくなります。ポアソン分布が確かに合っているということです。

#　上記の解析では、B, C, D社間での平均干し葡萄数の差を検定できていません。そのため、2. 追加命令文を実行してこれらの間での差を検定します。

#　比較の対象が2つだけの場合は2. 追加命令文は不要で、これで解析は終わりです。たとえば、A社とB社のデータしかなく、この2つの間での差のみを検定したい場合などです。上記の解析で両者の差(B社と、切片(Intercept)であるA社の差)が検定できているので、これ以上の解析は必要ありません。

# 一般化線形モデルの命令文

pairs(emmeans(r,"Maker"), adjust="holm")

# 1. 一般化線形モデルの命令文の実行結果がデータフレームrに格納されています。

# データフレームrと説明変数Makerを指定します。Makerを""で囲みます。

#　A, B, C, D社間で、平均値の有意差を総当たりで検定します。切片となる製造社(データの並び順で先頭になるもの)を自動的に入れ替え、その切片との比較をしてくれます。

#　adjust="holm"を忘れずに書いて下さい。

# 実行結果です。左の赤枠内にあるのが、製造社の各組合せでの平均干し葡萄数の差です。右の赤枠内にあるのが平均値の差のP値です。平均値に有意に差があるのかどうかを示しています。

# A社とB社の差のP値が0.8796になっています。ところが1. 一般化線形モデルの命令文の実行の方では、このP値は0.834593です。A社とC, D社との差もP値が変わっています。変ですね。その理由を、多重比較における補正で説明しています。この説明を必ず読んで下さい。

# 追加命令文

# 命令文が続く場合には+で繋げ、最後の命令文の後には+を付けません。

# 必須命令文

ggplot(d, aes(x = Maker, y = Number)) +

# データフレームdを指定します。

# x = Makerとして、製造社名(Maker)が入っているデータ列Makerをx軸(横軸)に指定します。xは小文字です。

# y = Numberとして、干し葡萄の数(Number)が入っているデータ列Numberをy軸(縦軸)に指定します。yは小文字です。

# データ列Makerに入っている製造社ごとに、データ列Numberに入っている干し葡萄の数のデータを作図します。

# 末尾の+を忘れないで下さい。

# 箱ひげ図における、作図に用いるデータフレームとデータの指定

geom_boxplot() +

# 描く図として箱ひげ図を指定します。

# 末尾の+を忘れないで下さい。ただし、以降の命令文を省略する場合は+を付けず、｢geom_boxplot()｣とします。

# 箱ひげ図の指定

# 省略してもよい命令文

# 以降の命令文を省略するとデフォルトで自動で描きます。

labs(x = "Maker", y = "Number") +

# x軸(横軸)の名称をMakerに、y軸(縦軸)の名称をNumberにします。名称を""で囲みます。x, yは小文字です。

# x軸･y軸の名称の指定

scale_y_continuous(limits = c(0, 30)) +

# y軸(縦軸)を描く範囲を0から30に指定します。

# x軸･y軸の値の最小値･最大値の指定

theme_classic() +

# 図の背景色を白にします。

# 図のデザインの指定

theme(
axis.title.x = element_text(size = 18),
axis.text.x = element_text(size = 18),
axis.title.y = element_text(size = 18),
axis.text.y = element_text(size = 9)
)

# ｘ軸(横軸)の目盛り(axis.text.x)には製造社名(A, B, C, D)が入っています。製造社名とx軸･y軸(縦軸)の名称(axis.title.y)の文字の大きさを18に、y軸の目盛り(axis.text.y)の文字の大きさを9に指定します。

# x軸･y軸の文字の大きさと書き方の指定