これから研究を始める高校生と指導教員のために　第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を終了して再び起動したときは再実行する必要があります。

# パッケージ(関数集)のRへの読み込み

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と異なるのかを示しています。x = 0でb = 0のときy = e⁰/(1 + e⁰) = 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と有意に異なる。誇大広告の店を選択する確率0.3600は0.5よりも有意に小さい。

3000円におけるaの値(1000円との差)のP値；0.02329　有意に0と異なる。つまり、1000円の時と3000円の時で、誇大広告の店を選択(Hype)する確率が有意に異なる。誇大広告の店を選択する確率は、1000円の時に0.3600で3000円の時に0.5200である。3000円の時の方が有意に高い。

5000円におけるaの値(1000円との差)のP値；4.63 × 10^{- 6}　有意に0と異なる。つまり、1000円の時と5000円の時で、誇大広告の店を選択(Hype)する確率が有意に異なる。誇大広告の店を選択する確率は、1000円の時に0.3600で5000円の時に0.6900である。5000円の時の方が有意に高い。

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

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

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

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

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

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

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

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

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

# 実行結果です。左の赤枠内にあるのが、金額の各組合せでの、選択確率(y = e^{ax + b}/(1 + e^{ax + b}))の中の係数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, 3-2のように加工します。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円を書き込みます。作図に用いるデータフレーム(この例ではデータフレームd)におけるデータの並び順(あなたが指定した並び順または数値･アルファベット順)通りに並べて下さい。並びが異なっていたらデータの対応がおかしくなります。名称を""で囲みます。

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

theme_classic(base_family = "HiraKakuPro-W3") +

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

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

# 図のデザインの指定

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