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

副読文書
VIII-6. 一般化線形モデル入門
1種類のデータにおけるデータ分布の偏りの解析
条件の変化の影響を見る場合
誇大広告の店に行くか正直な広告の店に行くかの選択の解析

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("Advertisement.csv")

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

# 画面1-2のデータのファイル名はAdvertisement2なので、こちらを読み込む場合はd <- read.csv("Advertisement2.csv")として下さい。

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

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

d2 <- reshape(d, idvar = "State", timevar = "Shop", direction = "wide")

# データフレームdは、画面2-1の形式になっています。それを、画面2-2の形式に変換してデータフレームd2に格納します。

# idvar = "State"として、データフレームdのデータ列State(1000yen, 3000yen, 5000yenが入っている)を、データフレームd2で縦方向に並ぶように指定します。Stateを""で囲みます。

# timevar = "Shop"として、データフレームdのデータ列Shop(Hype, Honestが入っている)を、データフレームd2で横方向に並ぶように指定します。Shopを""で囲みます。

# direction = "wide"を忘れずに書いて下さい。

# データフレームの形式の変換

colnames(d2) <- c("State", "Hype", "Honest")

# 列の名称を、State, Hype, Honestに変更します。名称を""で囲みます。

# データ列の名称の変更；まとめて変更する場合

r <- glm(formula = cbind(Hype, Honest) ~ State, family = binomial(link = logit), data = d2)
summary(r)

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

# 解析するデータフレーム；d2
# 応答変数；店の選択(cbind(Hype, Honest))。データ列Hypeに、誇大広告の店の選択数が、データ列Honestに、正直な広告の店の選択数が入っています。cbindを使って、各金額の条件における、(誇大広告の店の選択数, 正直な店の選択数)という組合せからなるデータを作ります。
# 説明変数；使える金額(State)。データ行1000yen, 3000yen, 5000yenのそれぞれに、その金額を使えるときのお店の選択が入っています。
# 確率分布；二項分布(binomial)
# 連結関数；S字(logit)

# 求める選択確率は、cbindで先に並べた方のものです。cbind(Hype, Honest)としたので、Hype(誇大広告の店)を選択する確率を計算します。Honest(正直な広告の店)を選択する確率は、｢1 - Hypeを選択する確率｣となります。

# cbind；データを、横に増える(列が増える)方向に統合

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

# 左の赤枠内にあるのが、各金額における、Hype(誇大広告の店)を選択する確率を示す数値です(ただし、確率そのものではありません)。各金額における選択確率は以下のようになります。

1000円の時に誇大広告の店を選択(Hype)する確率 = e^(Intercept)/(1 + e^(Intercept)) = e^{- 0.5754}/(1 + e^{- 0.5754}) = 0.3600
3000円の時に誇大広告の店を選択(Hype)する確率 = e^{State3000yen + (Intercept)}/(1 + e^{State3000yen + (Intercept)}) = e^{0.6554 - 0.5754}/(1 + e^{0.6554 - 0.5754}) = 0.5200
5000円の時に誇大広告の店を選択(Hype)する確率 = e^{State5000yen + (Intercept)}/(1 + e^{State5000yen + (Intercept)}) = e^{1.3755 - 0.5754}/(1 + e^{1.3755 - 0.5754}) = 0.6900

連結関数がS字なので、y = e^{ax + b}/(1 + e^{ax + b}) (y；応答変数(選択確率)　x；説明変数の値　a；傾きの値　b；切片の値)という形で応答変数(選択確率)を示します。このとき、データの並び順の先頭に出てくるものの選択確率を切片(Intercept)として扱います。本解析ではデータの並び順を指定していないので、数値が小さいもの(1000円)が自動的に先頭となります。つまり、1000円の場合の選択確率は、x = 0で切片b(Intercept)の値のみを持つe^b/(1 + e^b)に当てはめたものとなります。3000円, 5000円の場合の選択確率は、x = 1のときの値、すなわちy = e^{a + b}/(1 + e^{a + b})として表現します。このとき、1000円, 3000円, 5000円すべてを1つの回帰式に組み込むのではなく、｢1000円と3000円の場合の選択確率の回帰式｣｢1000円と5000円の場合の選択確率の回帰式｣という別々の回帰式を計算します。そして、各回帰式にx = 1を当てはめた値y = e^{a + b}/(1 + e^{a + b})を求めます。左の赤枠内の2行目以降が、その金額の回帰式のa(傾き)の値です。これに、1000円の場合の値(切片bの値)を加えた値から、その金額における選択確率を得ます。

#　右の赤枠内にあるのがP値です。1番上は、切片(Intercept)の値(bの値)が有意に0と異なるのかを示しています。b = 0のときy = 1/2となるので、どちらにも偏らず1/2の確率で選択するのかどうかを見ることになります。2番目以降は、aの値(切片である1000円との差)が有意に0と異なるのかを示しています。｢e-｣とあるのは10のマイナス何乗かということです。｢4.63e-06｣は｢4.63 × 10^{- 6}｣です。

切片(Intercept)のP値；0.00575　有意に0と異なる　つまり、1000円の時に誇大広告の店を選択(Hype)する確率は0.5と有意に異なる。
3000円におけるaの値(1000円との差)のP値；0.02329　有意に0と異なる　つまり、1000円の時と3000円の時で、誇大広告の店を選択(Hype)する確率が有意に異なる。
5000円におけるaの値(1000円との差)のP値；4.63 × 10^{- 6}　有意に0と異なる　つまり、1000円の時と5000円の時で、誇大広告の店を選択(Hype)する確率が有意に異なる。

#　AICも出ているのですが、確率分布と連結関数の組合せは二項分布とS字の一択です。だから気にする必要はありません。

#　3000円と5000円間での選択確率の差は検定できていません。そのため、2. 追加命令文を実行してこれらの間での差を検定します。

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

# 上記の解析では、選択する確率が1/2かどうかの検定は1000円についてしかできていません。3000円, 5000円についても検定したい場合は、データフレーム中での並び順を変えて、3000円, 5000円を先頭(切片になる)にしてもう一度解析します(こちらを参照)。

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

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

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

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

#　1000円, 3000円, 5000円間で、誇大広告の店(Hype)を選択する確率の有意差を総当たりで検定します。切片となる金額(データの並び順で先頭になるもの)を自動的に入れ替え、その切片との比較をしてくれます。

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

# 実行結果です。左の赤枠内にあるのが、金額の各組合せでの、選択確率の中の係数aの値です。右の赤枠内にあるのが、aの値が有意に0と異なるかのP値です。選択確率に有意に差があるのかどうかを示しています。

# 1000円と3000円の差のP値が0.0291になっています。ところが1. 一般化線形モデルの命令文の実行の方では、このP値は0.002329です。1000円と5000円との差も実はP値が変わっています(< .0001と表記されているので、変わっているとわからないのですが)。変ですね。その理由を、多重比較における補正で説明しています。この説明を必ず読んで下さい。

# 追加命令文

d2$State <- factor(d2$State, levels=c("3000yen", "1000yen", "5000yen"))

# データフレームd2に入っているデータ列Stateにおけるデータの並び順を、3000円, 1000円, 5000円にします。金額は文字情報なので、金額を""で囲みます。

# # d2$Stateの中身を上書きするためにd2$Stateに再格納します。

# データの並び順の指定

r <- glm(formula = cbind(Hype, Honest) ~ State, family = binomial(link = logit), data = d2)
summary(r)

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

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

#　切片(Intercept)の値(bの値)のp値は0.6892なので、誇大広告の店(Hype)を選ぶ確率が有意に1/2と異なるとはいえません。

#　d2$State <- factor(d2$State, levels=c("5000yen", "1000yen", "3000yen"))として5000円を切片にして同様の解析を行うと、p値 = 0.000215となり、有意に1/2と異なるという結果になります。

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

# 必須命令文

ggplot(d, aes(x = State, y = Number, fill = Shop)) +

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

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

# y = Numberとして、誇大広告の店(Hype)と正直な広告の店(Honest)を選択した数が入っているデータ列Numberをy軸(縦軸)に指定します。yは小文字です。

# fill = Shopとして、描き分けに用いるデータ列Shopを指定します。

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

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

geom_bar(stat = "identity", position = "dodge") +

# 描く図として棒グラフを指定します。

# stat = "identity"を忘れずに書いて下さい。

# position = "dodge"と指定して、誇大広告の店(Hype)と正直な広告の店(Honest)を選んだ数の棒グラフを横に並べます。

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

# 棒グラフの指定

# 省略してもよい命令文

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

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

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

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

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

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

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

theme_classic() +

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

# 図のデザインの指定

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

# ｘ軸(横軸)の目盛り(axis.text.x)には金額(1000yenなど)が入っています。金額とy軸(縦軸)の名称(axis.title.y)の文字の大きさを18に、y軸の目盛り(axis.text.y)の文字の大きさを9に指定します。ｘ軸(横軸)の名称(axis.title.x)は書かないので、その文字の大きさも指定しません。

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

#　こうして描いた図をPowerPoint等に読み込み、画面3-1のように加工します。PowerPoint等への図の読み込み方は、図の保存･コピーの仕方を参照して下さい。

　画面3-2のように軸の説明を日本語にする場合は、2つの命令文(下記の1つ目と3つ目)を書き替え、1つの命令文(下記の2つ目)を書き足します。他の命令文はそのままでよいです。

labs(x = "", y = "選択数") +

# x軸(横軸)の名称は書かず、y軸(縦軸)の名称を選択数にします。名称を""で囲みます。x, yは小文字です。

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

scale_x_discrete(labels = c("1000円", "3000円", "5000円")) +

#　書き足す命令文です。labs(x = "", y = "選択数") +の次に書いて下さい。

# x軸(横軸)の各データ群名として1000円･3000円･5000円を書き込みます。データの並び順(あなたが指定した並び順または数値･アルファベット順)通りに並べて下さい。並びが異なっていたらデータの対応がおかしくなります。名称を""で囲みます。

# x軸にデータ群名を書き込む

theme_classic(base_family = "HiraKakuPro-W3") +

# 日本語のフォントとしてHiraKakuPro-W3を指定します。フォントの名称を""で囲みます。

#　HiraKakuPro-W3の他にも、HiraKakuPro-W6, Meiryo, MS Gothic, Osakaなど色々なフォントがあります。お好みのフォントを設定して下さい。

# 図のデザインの指定

#　こうして描いた図をPowerPoint等に読み込み、画面3-2のように加工します。PowerPoint等への図の読み込み方は、図の保存･コピーの仕方を参照して下さい。