統計解析ノート（R）

soc-ja

統計解析ノート（R）

1. Rのインストール
  以下のページを参考にインストールする．
  Rのインストール方法（Windows）

1. dplyrパッケージのインストール
  以下のページを参考にインストールする
  dplyrパッケージのインストール方法（Windows – R）

1. 棒グラフの書き方
  例えば2群の値を棒グラフで示すとき、有意差があるかどうかを直感的に分かるようにするために、エラーバーを示すようにする。エラーバーは、標準誤差 (SE) を示すことが一般的であるが、それぞれの群の平均値の推定の95%信頼区間を示すこともある。
  標準誤差は、標準偏差を√nで割ることで得られる。平均値の推定の95%信頼区間は、
  t.test(a)
  t.test(b)
  のように、それぞれの群のベクトルデータをt.test関数に入力すれば得られる。
  多重比較等を行っていない場合は、95%信頼区間が重なっていなければ、有意差があるとみなして良い。標準誤差の場合は、重なっていなくても、統計的有意差が得られない場合があるため、注意が必要である。
  参考：https://www.igaku-shoin.co.jp/paper/archive/y2012/PA02967_03
2. 図のeps化
  dev.copy2eps関数を使うと良い．
  a = rnorm(50, 50, 15)
  hist(a)
  dev.copy2eps(file = “figure.eps”, family = “Japan1”)
3. 2変数（定量データ）の相関分析
  以下の手順に従うこと
  １．ヒストグラムにてそれぞれデータの分布を目視で確認
  ２．散布図でそれぞれのデータの平面上での分布を目視で確認
  ３．データが正規分布に従っているかどうかの確認（コルモゴロフ–スミルノフ検定・シャピロ・ウィルク検定など）
  （正規分布に従っている場合（パラメトリック手法を適用）
  ４A－１：相関係数（Pearsonの積率相関係数）を算出
  ４A－２：無相関検定を実施
  （正規分布に従っていない場合（ノンパラメトリック手法を適用）
  ４B－１：順位相関係数（Spearmanの順位相関係数など）を算出
  ４B－２：無相関検定（t分布を利用）を実施

1. ヒストグラムの表示
  定量データのヒストグラムの確認する．
  > hist(x)
  ※xはベクトルデータ

1. 散布図の表示
  2変数の散布図の確認する．
  > plot(x, y)
  ※x, yはベクトルデータ

1. 離散データの散布図の表示
  離散データの場合，異なるデータであっても，同じ座標に点が表示されてしまい，相関が読み取れなくなる．そこで，jitter関数により，各データに人工的に測定誤差を追加して可視化する．以下は，ランダムに発生させたデータに対する利用例．7段階のリッカート尺度で回答されたデータを想定．
  a = ceiling(runif(100)*7)
  b = ceiling(runif(100)*7)
  plot(a,b)　　←　jitterなし
  plot(jitter(a),jitter(b)　　←　デフォルトのjitter．パラメータの設定なしで利用
  plot(jitter(a, factor=3, amount=0),jitter(b, factor=3, amount=0))　　←　パラメータfactorの値を設定し，より連続データに近い散布図に整形
  ※factorの値を大きくすると，測定誤差を大きくすることができる

1. コルモゴロフ–スミルノフ検定
  分布の正規性，または2標本の代表値や分布の違いを検定する方法．
  (1) 分布の正規性の検定の場合
  データが想定した確率分布に従っているかどうかを検定する．正確には，帰無仮説「データの分布が仮定した確率分布に一致していること」を棄却できるかを検定する．正規分布の場合には，得られたデータが正規分布でないことは保証できる．しかし，正規分布であることは保証できない．もし，パラメトリックな手法の適用の判断材料とするのであれば，「正規分布でないとは言えないため，パラメトリックな手法を適用する．」という説明になる．以下のように実行する．
  > ks.test(x, “pnorm”, mean=mean(x), sd=sd(x))
  ※xはベクトルデータ
  P値が0.05（有意水準5%）を上回れば，正規分布とみなすことができる．（正規分布だと仮定した場合，このようなデータ群が得られることは十分にあり得る）
  P値が0.05（有意水準5%）を下回れば，正規分布とは言えない．（正規分布だと仮定した場合，このようなデータ群が得られることはほとんどない）
  (2) 2標本の代表値の分布の違いを検定する場合
  帰無仮説「2つの標本の確率分布が一致していること」を棄却できるかを検定する．すなわち，「2つの標本の確率分布が一致しているとは言えないかどうか」を確かめる．
  > ks.test(x, y)
  ※x,yはベクトルデータ
  参考記事：http://darusmart.hatenablog.com/entry/2015/02/17/062141

1. シャピロ・ウィルク検定 (Shapiro-Wilk test)
  分布の正規性を検定する方法．コルモゴロフ・スミルノフ検定等と同様に，得られたデータが正規分布に従うものか否かを調べる検定法であるが，シャピロ・ウィルク検定の方がやや慎重に検定する傾向がある（正規性がないと判断されやすい）．そのため，大規模データを扱うことの多い今日においては，コルモゴロフ・スミルノフ検定を使うほうが実用的と言える．シャピロウィルク検定はデータ数が100以下の場合に用いられることが多い．
  コルモゴロフ・スミルノフ検定等と同様に，帰無仮説「データの分布が正規分布分布に一致していること」を棄却できるかを検定する．得られたデータが正規分布でないことは保証できるが，正規分布であることは保証できない．以下のように実行する．
  > shapiro.test(x=data)
  ※dataはベクトルデータ
  P値が0.05（有意水準5%）を上回れば，正規分布とみなすことができる．（正規分布だと仮定した場合，このようなデータ群が得られることは十分にあり得る）
  P値が0.05（有意水準5%）を下回れば，正規分布とは言えない．（正規分布だと仮定した場合，このようなデータ群が得られることはほとんどない）
  参考記事：https://data-science.gr.jp/implementation/ist_r_shapiro_wilk_test.html

1. F検定（分散の等質性の検定・等分散性の検定）
  t検定や分散分析を行う前に，分散の等質性を必ず確認しておく必要がある．2群のデータの間で，分散が等しいとみなして良いかどうかを検定する．
  以下のデータは，分散分析の実施を検討するため，3つの群（A, B, C）からなるデータ（性格を表す心理指標を想定）を例とする．
  > psyA
  [1] 15 9 18 14 18
  > psyB
  [1] 10 7 3 5 7
  > psyC
  [1] 10 6 11 7 12
  
  > var.test(psyA,psyB)
  
  F test to compare two variances
  
  data: psyA and psyB
  F = 2.0147, num df = 4, denom df = 4, p-value = 0.5142
  alternative hypothesis: true ratio of variances is not equal to 1
  95 percent confidence interval:
  0.2097662 19.3503029
  sample estimates:
  ratio of variances
  2.014706
  
  > var.test(psyA,psyC)
  
  F test to compare two variances
  
  data: psyA and psyC
  F = 2.0448, num df = 4, denom df = 4, p-value = 0.5055
  alternative hypothesis: true ratio of variances is not equal to 1
  95 percent confidence interval:
  0.2128971 19.6391133
  sample estimates:
  ratio of variances
  2.044776
  
  > var.test(psyB,psyC)
  
  F test to compare two variances
  
  data: psyB and psyC
  F = 1.0149, num df = 4, denom df = 4, p-value = 0.9889
  alternative hypothesis: true ratio of variances is not equal to 1
  95 percent confidence interval:
  0.1056715 9.7478811
  sample estimates:
  ratio of variances
  1.014925
  
  この検定は，帰無仮説が「分散が等質である」になる．上記の例では，いずれの比較も，p値が 0.05を上回ったので，帰無仮説が採択され，分散は等質であると言えることになる．p値が0.05を下回った場合は，対立仮説「分散は等質とは言えない」が採択される．
  参考記事：https://data-science.gr.jp/implementation/ist_r_f_test.html
  参考記事：http://www.math.kobe-u.ac.jp/HOME/higuchi/13.pdf

1. 相関係数（パラメトリック：Pearsonの積率相関係数）
  量的変数同士の関係を確かめるのに使用．
  > cor(x, y)
  ※x,yはベクトルデータ
  ※欠損値(NA)のある場合は，use = “complete.obs”を引数として追加すると，データが完全な行のみを用いて計算してくれる
  ※相関係数が意味を持つのは，データ数が少なくとも10個以上ある場合

1. すべての特徴間の相関係数（パラメトリック：Pearsonの積率相関係数）
  複数の質問項目から構成される心理尺度において，質問項目間の相関を調べることがある．その際，すべての質問項目間の相関を算出する必要があるが，列に各質問（特徴）を割り当て，行が回答レコードに割り当てれば，以下のように実行することで，すべての列間の相関係数を一発で出してくれる．
  > D = matrix(c(1,2,3,4,5,6,7,8,9,4,3,2), nrow=3, ncol=4)
  > D
  [,1] [,2] [,3] [,4]
  [1,] 1 4 7 4
  [2,] 2 5 8 3
  [3,] 3 6 9 2
  > cor(D)
  [,1] [,2] [,3] [,4]
  [1,] 1 1 1 -1
  [2,] 1 1 1 -1
  [3,] 1 1 1 -1
  [4,] -1 -1 -1 1

1. 無相関検定（パラメトリック）
  量的変数同士の相関の有無を統計的に検定．
  > cor.test(x, y)
  ※x,yはベクトルデータ
  p-valueが，0.05未満であれば，「相関がないとは言えない（相関係数は0ではない，相関がある）」と判断する（有意水準5%時）
  データ数が多い場合，「有意（相関はゼロではない）」と判断されやすいので，注意が必要である．相関係数は低いが，無相関検定では有意と判定された場合，「○○と○○の間には，正の相関があるが，その関連性は弱いと言える．」や「○○と○○の間には，ごくわずかに正の相関がある．」というような結論にすること．
  無相関検定を行うのに必要なデータ数（標本の大きさ）は，母相関係数をrと想定している場合，n = 4 + 8 / |r| で計算できる．例えば，r が0.4の場合は，24個のデータが必要である．

1. Spearmanの順位相関係数（ノンパラメトリック）
  正規性が保証されない二つの変数（または順序尺度）の相関を得る方法（ノンパラメトリック手法）．ピアソンの積率相関係数の特別な（相関係数を計算する前にデータを順位に変換した）場合に当たる．
  > cor(x, y, method=”spearman”)
  ※x,yはベクトルデータ
  ※同じデータでも，ピアソン相関（Pearsonの積率相関係数）に比べると，低い相関係数が出る．スピアマンの順位相関係数は，ピアソン相関よりも，より慎重に（保守的に）相関係数を算出してくれる．そのため，査読でもケチをつけられにくい．
  参考記事：http://cse.naro.affrc.go.jp/takezawa/r-tips/r/67.html
  参考記事：http://www.tamagaki.com/math/Statistics609.html

1. Spearmanの無相関検定（ノンパラメトリック）
  正規性が保証されない二つの変数（または順序尺度）の相関の有無を検定する．スピアマンの順位相関係数$\rho$が0でないかどうかを検定する．
  > cor.test ( x, y, “spearman”)
  ※x,yはベクトルデータ
  ※同じ値のデータがあると，「タイのため正確な p 値を計算することができません」という警告がでるが，無視して良い．cor.test(x, y, “spearman”, exact=FALSE) とすれば，警告を消すことができる（算出されるp値は変わらない）．
  ※p-valueが，0.05未満であれば，「相関がないとは言えない（相関係数は0ではない，相関がある）」と判断する（有意水準5%時）
  データ数が多い場合，「有意（相関はゼロではない）」と判断されやすいので，注意が必要である．相関係数は低いが，無相関検定では有意と判定された場合，「○○と○○の間には，正の相関があるが，その関連性は弱いと言える．」や「○○と○○の間には，ごくわずかに正の相関がある．」というような結論にすること．
  参考記事：https://www.mk-mode.com/blog/2020/02/17/fortran95-rcc-spearman/
  参考記事：http://rstudio-pubs-static.s3.amazonaws.com/327506_fa678800ada44ab19dbb6c5ec35df6c0.html
  参考記事：https://www.sbj.or.jp/wp-content/uploads/file/sbj/9508/9508_tokei_kaiseki.pdf

1. Kendallの順位相関係数
  正規性が保証されない二つの変数（または順序尺度）の相関を得る方法（ノンパラメトリック手法）．順位関係の一致（または不一致）する組の数より算出する．
  > cor(x, y, method=”kendall”)
  ※x,yはベクトルデータ
  ※同じ値のデータがあると，「タイのため正確な p 値を計算することができません」という警告がでるが，無視して良い．
  参考記事：http://cse.naro.affrc.go.jp/takezawa/r-tips/r/67.html

1. 相関係数の判定（パラメトリック，ノンパラメトリックの両方）
  0～0.2 （ほとんど）相関なし
  0.2～0.4 弱い相関あり
  0.4～0.7 中程度の相関あり
  0.7～1.0 強い相関あり
  と言われている（厳密な定義は難しいため，公式な判定基準は存在しない）．

1. 2変数（定性データ）の相関分析
  質的変数（カテゴリカル変数）同士の関係を分析する場合は，以下の手順に従う．
  1. クロス集計表を算出
  2. クラメールの連関係数を算出
  3. カイ二乗検定を実施
  4. 残差分析を実施（変数に3水準以上のものがある場合のみ）
  ※カイ二乗検定を行わずに、直接に残差分析を行っても良い。カイ二乗検定は表全体での検定。残差分析は個別の行（セル）の検定。カイ二乗検定で有意差が出なくても、残差分析で有意差が出ることは、たまにある。

1. クロス集計表の算出
  table(x,y)
  カテゴリカル変数x, yにデータが入っているものとする．それぞれの同一要素は，同一オブジェクトに対応．
  具体的にはデータフレーム中の二つの列をそれぞれ別々の変数(x,y)で取得すればよい．

1. 連関の程度（クラメールの連関係数・Cramer’s V）
  質的変数（カテゴリカル変数）同士の関係を確かめるのに使用．
  1. vcdパッケージをインストール
  > install.packages(“vcd”, dependencies = TRUE)
  2．vcdパッケージを有効化
  > library(vcd)
  3. 関数assocstatsを使用
  > assocstats(クロス集計表の行列）
  （例）
  > D = matrix(c(2,2,1,4,2,10,11,3,15,16,9,2,21,3,2,29),nrow=4,ncol=4)
  > assocstats(D)
  X^2 df P(> X^2)
  Likelihood Ratio 66.24 9 8.2727e-11
  Pearson 58.69 9 2.3954e-09
  
  Phi-Coefficient : NA
  Contingency Coeff.: 0.555
  Cramer’s V : 0.385

1. カイ二乗検定・フィッシャーの直接確率検定
  量的変数同士の連関の有無を統計的に検定．
  以下の方針に従う．
  ・標本の大きさが50以上で、セルの実測度数で5未満のものはない→カイ二乗検定
  ・標本の大きさが50以上で、セルの実測度数で5未満のものがある→フィッシャーの直接確率検定
  ・標本の大きさが50未満→カイ二乗検定（イェーツの連続性補正）
  
  （カイ二乗検定）
  > chisq.test(D, correct=FALSE)
  
  （カイ二乗検定（イェーツの連続性補正））
  > chisq.test(D, correct=TRUE)
  
  （フィッシャーの直接確率検定）
  > fisher.test(D)
  
  ※Dはクロス集計表の入った行列
  
  （フィッシャーの直接確率検定（詳細））
  大きなクロス集計表を対象にフィッシャーの直接確率検定を行う場合は，メモリ不足になることがある．
  > fisher.test(D, workspace=10000000)
  のように，大きなworkspace値を設定するか，
  > fisher.test(D, simulate.p.value=TRUE)
  として近似を行うかすること．
  
  ※カテゴリ値のパターンが多くなりすぎないように注意すること．度数の少ないカテゴリ値があれば，他のカテゴリ値と結合できないかどうか検討する．

1. 残差分析
  クロス集計表において，どのセルが特に期待度数からのずれ（残差）が大きかったのかを統計的に検証する．変数のうち1つでも3水準以上を取りうるものがある場合には，実行する必要がある．なぜなら，カイ二乗検定のみでは，どの水準に連関が見られ，どの水準に連関が見られないかを明らかにできないからである．残差とは，「実測度数−期待度数」であり，残差分析を行うことで期待度数と実測度数のずれが特に大きかったセルを発見することが出来る．Rでは，以下のように実行する．
  > result = chisq.test(D)
  > result$stdres
  [,1] [,2]
  [1,] 4.1833804 -4.1833804
  [2,] 0.2741955 -0.2741955
  [3,] -1.5688628 1.5688628
  [4,] -0.9970697 0.9970697
  > pnorm(abs(result$stdres), lower.tail = FALSE) * 2
  [,1] [,2]
  [1,] 2.872062e-05 2.872062e-05
  [2,] 7.839344e-01 7.839344e-01
  [3,] 1.166799e-01 1.166799e-01
  [4,] 3.187307e-01 3.187307e-01
  ※Dは，クロス集計表の行列とする．
  ※result$stdresには，調整済み残差が行列の形式で入っている．
  ※この結果（最後の結果：p値を示したもの）では，1行目のカテゴリのみ，列の頻度に有意差（p値＜0.05）があることが検証できた．
  ※p値は出さなくても，調整済み残差を見るだけでも，判定できる．調整済み残差は，標準正規分布に従う．そのため，調整済み残差が1.96より大きければ（有意水準0.05の両側検定），そのセルは統計的に有意差がある（期待度数と異なる）とみて良い．
  参考記事（残差分析）：https://to-kei.net/hypothesis-testing/chi2-test-residual-analysis
  参考文献（Rでの残差分析）：http://langstat.hatenablog.com/entry/20150211/1423634853

1. 2群の代表値の検定
  以下の方法に従う．
  １．ヒストグラムにてそれぞれの群の分布を目視で確認
  ２．それぞれの群の平均値と標準偏差を求める（どちらの群の方が値が高いか，どれぐらいばらつきがあるのか，おおむね理解する）
  ３．データが正規分布に従っているかどうかの確認（コルモゴロフ–スミルノフ検定・シャピロ・ウィルク検定など）
  （正規分布に従っている場合（パラメトリック手法を適用）
  ４Ａ．等分散かどうかの検定
  （等分散であれば）
  ４Ａａ：t検定
  （等分散でなければ）
  ４Ａｂ：ウェルチ(Welch)の検定
  （正規分布に従っていない場合（ノンパラメトリック手法を適用）
  ４B：ウィルコクソンの順位和検定（マンホイットニーのU検定）
  ※データに対応がある場合（例：同じ被験者が条件Aも条件Bも実施している場合）は，対応のある検定を行う．

1. t検定
  正規分布に従うデータ（等分散が仮定できる）について，2つの群の代表値（平均値）に違いがあるかどうかを統計的に検定する．
  実行例は以下のとおりである．
  > psyA
  [1] 15 9 18 14 18
  > psyB
  [1] 10 7 3 5 7
  > t.test(psyA,psyB,var.equal=TRUE)
  
  Two Sample t-test
  
  data: psyA and psyB
  t = 4.1485, df = 8, p-value = 0.003216
  alternative hypothesis: true difference in means is not equal to 0
  95 percent confidence interval:
  3.730698 13.069302
  sample estimates:
  
  ※デフォルトは両側検定．片側検定を行いたい場合は，パラメータとして alternative=”greater” （または，alternative=”less”）を追加．
  ※デフォルトは対応無し．データに対応がある場合は，パラメータとして paired=TRUE を追加．
  
  論文に結果を載せるときには、「群Aと群Bの間に有意差が見られた（t(9)=2.618, p=0.028）のように、検定統計量tの値とp値を書く。

1. ウェルチの検定
  正規分布に従うデータ（等分散が仮定できない）について，2つの群の代表値（平均値）に違いがあるかどうかを統計的に検定する．
  実行例は以下のとおりである．
  > psyA
  [1] 15 9 18 14 18
  > psyB
  [1] 10 7 3 5 7
  > t.test(psyA,psyB)
  
  Welch Two Sample t-test
  
  data: psyA and psyB
  t = 4.1485, df = 7.1859, p-value = 0.004064
  alternative hypothesis: true difference in means is not equal to 0
  95 percent confidence interval:
  3.636998 13.163002
  sample estimates:
  mean of x mean of y
  14.8 6.4
  
  ※デフォルトは両側検定．片側検定を行いたい場合は，パラメータとして alternative=”greater” （または，alternative=”less”）を追加．
  ※デフォルトは対応無し．データに対応がある場合は，パラメータとして paired=TRUE を追加．

1. ウィルコクソンの順位和検定（マンホイットニーのU検定）
  正規分布を仮定できない場合（正規分布に従わないデータ）について，2つの群の代表値（中央値）に違いがあるかどうかを統計的に検定する．
  実行例は以下のとおりである．
  （パッケージのインストール）
  > install.packages(“exactRankTests”, repos=”http://cran.ism.ac.jp/”)
  > library(exactRankTests)
  （実行）
  > psyA
  [1] 15 9 18 14 18
  > psyB
  [1] 10 7 3 5 7
  > wilcox.exact(psyA,psyB,paired=F)
  
  Exact Wilcoxon rank sum test
  
  data: psyA and psyB
  W = 24, p-value = 0.01587
  alternative hypothesis: true mu is not equal to 0
  ※デフォルトは両側検定．片側検定を行いたい場合は，パラメータとして alternative=”greater” （または，alternative=”less”）を追加．
  ※デフォルトは対応無し．データに対応がある場合は，パラメータとして paired=TRUE を追加．
  ※検定統計量はW（Mann-WhitneyのU値に相当）．自由度は一般には指定されない．
  
  効果量の求め方
  result = wilcox.test(dat$Score~dat$Class, correct=FALSE)
  # correct=TRUEとしても正確なp値は得られないので使用しないようにする
  result
  
  # z値の計算
  pval = result$p.value
  z = qnorm(1-(pval/2))
  z
  
  # 効果量rの計算（r = z値/sqrt(サンプル数)で求める）
  r = z/sqrt(length(dat$Score)) #lengthはデータ数を確認する関数
  r
  
  検定統計量は W として出力される。これは，前述の U1 U2 と比べると， W は大きいほうの U2 なのである。この点に注意する必要がある。
  
  参考：https://mizumot.com/handbook/?page_id=422
  参考：https://corvus-window.com/excel_mann-whitney-test/
  参考：https://biolab.sakura.ne.jp/mann-whitny-u-test-statistic.html （統計量W）

1. 分散分析
  分散分析が必要なのは，
  1. 目的変数（定量データ）と説明変数（カテゴリカルデータ）ともに1つずつなのであるが，説明変数のカテゴリ（水準）が3つ以上の場合．（一元配置分散分析．1要因の分散分析）
  2. 目的変数（定量データ）が1つであるが，説明変数（カテゴリカルデータ）が2つの場合．2つの説明変数のカテゴリ（水準）は2つであっても，3つ以上であっても良い．（二元配置分散分析．2要因の分散分析）
  3. 目的変数（定量データ）が1つであるが，説明変数（カテゴリカルデータ）が3つ(以上)の場合．2つの説明変数のカテゴリ（水準）は2つであっても，3つ以上であっても良い．（三元配置分散分析（多元配置分散分析）．3要因(以上)の分散分析）
  のいずれか．
  ※基本的には，目的変数（定量データ）は1つで，調べたい条件（カテゴリカル変数）がいくつかある場合に使う．
  ※三元配置分散分析（多元配置分散分析）については，以下のURLを参照
  https://www.hulinks.co.jp/support/sigmaplot/v14/usersguide/p190700.html
  ※分散分析とt検定の利用場面の違いについては，以下のURLを参照
  https://www.recruit-ms.co.jp/issue/column/0000000645/?theme=hrdepartment

1. 1要因2水準の分散分析
  1要因2水準の分散分析（すなわち2水準の一元配置分散分析）と対応のないt検定（2群の平均値の検定）は，同じことを行っている（すなわち検定結果は一致する）．t^2=Fとなる．
  参考記事：https://bellcurve.jp/statistics/blog/20333.html

1. 分散分析の考え方と結果の表記方法
  分散分析の結果の表記は
  (F (1, 88)=2.03, p<.05) のように書き， F（群間の自由度, 群内の自由度）=2.03, p<.05 と読む．群間の自由度＝群の数-1 群内の自由度＝(群1のサンプル数-1)＋ (群2のサンプル数-1)＋ (群3のサンプル数-1)＋・・・＋(群nのサンプル数-1) である．（以下，参考）
  
  分散分析の記述について－F( )内の数字の意味
   エクセルを使った分散分析のやり方
   エクセルを使った分散分析のやり方
   分散分析表の見方を理解して仮説を検証
2. 3群以上の代表値の検定
  以下の手順に従うこと
  １．箱ひげ図（ボックスプロット）にてそれぞれデータの分布を目視で確認（平均も算出）
  ２．データが正規分布に従っているかどうかの確認（コルモゴロフ–スミルノフ検定など）
  ３．データの分散の等質性を確認（F検定など）
  （正規分布に従っている場合（パラメトリック手法を適用）
  ４A－１：一元配置分散分析（one-way ANOVA）を実施
  ４A－２：多重比較（パラメトリック）を実施
  （正規分布に従っていない場合（ノンパラメトリック手法を適用）
  ４B－１：クラスカル・ウォリス検定（Kruskal-Wallis test）を実施
  ４B－２：多重比較（ノンパラメトリック）を実施
  
  参考記事：https://www.jstage.jst.go.jp/article/kagakutoseibutsu/51/7/51_483/_pdf
3. 一元配置分散分析（対応なし）（one-way ANOVA）
  3群以上において平均値の差を検定する方法である．対象のデータが正規分布に従う場合のみ適用できる．例えば，3つの群（A, B, C）からなるデータ（性格を表す心理指標を想定）があるとする．対応がない場合なので，同じユーザ（個体）が3つの群のすべてを試行しているのではなく，一人のユーザ（個体）はA, B, Cのいずれかを試行している．
  > psy
  [1] 10 15 6 9 10 7 3 18 14 18 11 5 7 12 7
  > group
  [1] C A C A B B B A A A C B C C B
  Levels: A B C
  > oneway.test(psy~group, var.equal=TRUE)
  
  One-way analysis of means
  
  data: psy and group
  F = 10.088, num df = 2, denom df = 12, p-value = 0.002691
  のように検定する．p値が0.05（有意水準0.05）を下回っているので，有意差ありとなる．
  var.equalは，分散の等質性が仮定できる場合はTRUE，そうでない場合はFALSEとなる．
  
  以下の方法でも実行できる．
  > summary(aov(psy~group))
  Df Sum Sq Mean Sq F value Pr(>F)
  group 2 182.9 91.47 10.09 0.00269 **
  Residuals 12 108.8 9.07
  —
  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  p値が，oneway.testを用いた時と同じになっていることがわかる．
  
  分散分析の結果の解釈は，群の種類という変数によって，対象データの平均値に差が発生するということだけである．群に対するカテゴリ値，すなわちどの群とどの群（例えばA群とB群）の間で，平均値に差があるかということまでは分からない．それを明らかにするのが多重比較になる．
  
  参考記事：http://english-writing-theaters.com/6_statistics.html
4. クラスカル・ウォリス検定（対応なし）（Kruskal–Wallis test）
  3群以上において平均値の差を検定するノンパラメトリックな方法である．一元配置分散分析のノンパラメトリック版となる．例えば，3つの群（A, B, C）からなるデータ（性格を表す心理指標を想定）があるとする．
  > psy = c(10,15,6,9,10,7,3,18,14,18,11,5,7,12,7)
  > grp = c(‘C’,’A’,’C’,’A’,’B’,’B’,’B’,’A’,’A’,’A’,’C’,’B’,’C’,’C’,’B’)
  > group = factor(grp)
  > psy
  [1] 10 15 6 9 10 7 3 18 14 18 11 5 7 12 7
  > group
  [1] C A C A B B B A A A C B C C B
  Levels: A B C
  > kruskal.test(psy ~ groupr)
  
  Kruskal-Wallis rank sum test
  
  data: psy by groupr
  Kruskal-Wallis chi-squared = 7.9805, df = 2, p-value = 0.0185
  
  のように検定する．p値が0.05（有意水準0.05）を下回っているので，有意差ありとなる．
  一元配置分散分析と同様，クラスカル・ウォリス検定の結果の解釈は，群の種類という変数によって，対象データの代表値に差が発生するということだけである．群に対するカテゴリ値，すなわちどの群とどの群（例えばA群とB群）の間で，代表値に差があるかということまでは分からない．それを明らかにするのが多重比較になる．
  
  なお，二元配置分散分析以上のノンパラメトリックな方法は存在しない．
5. 一元配置分散分析（対応あり）（one-way ANOVA）
  3群以上において平均値の差を検定する方法である．対象のデータが正規分布に従う場合のみ適用できる．一人のユーザ（個体）がすべての群の実験を試行している場合を対象とする．要因1＝群, 要因2=ユーザとなっているので，繰り返しのない二元配置分散分析とも呼ばれる．
  例えば，3つの群（A, B, C）からなる生理指標データ（センサーで測定した値）があるとする．具体的には，薬A, B, Cを飲んだ後の生理指標データを測定する．対応がある場合なので，一人のユーザが，すべての群の実験に参加する．ユーザは5人いる（uA, uB, uC, uD, uE）．
  
  以下のように検定する．
  > a = c(59, 78, 86, 97, 51)
  > b = c(56, 90, 65, 98, 48)
  > c = c(87, 99, 111, 130, 100)
  > all = c(a,b,c)
  > tp = c(rep(“A”,5), rep(“B”,5), rep(“C”,5))
  > type = factor(tp)
  > us = c(rep(c(“uA”, “uB”, “uC”, “uD”, “uE”), 3))
  > user = factor(us)
  > all
  [1] 59 78 86 97 51 56 90 65 98 48 87 99 111 130 100
  > type
  [1] A A A A A B B B B B C C C C C
  Levels: A B C
  > user
  [1] uA uB uC uD uE uA uB uC uD uE uA uB uC uD uE
  Levels: uA uB uC uD uE
  > summary(aov(all~type+user))
  Df Sum Sq Mean Sq F value Pr(>F)
  type 2 3562 1781.1 19.94 0.000779 ***
  user 4 3653 913.2 10.22 0.003114 **
  Residuals 8 715 89.3
  —
  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  typeのp値が0.05（有意水準0.05）を下回っているので，薬に対して有意差ありとなる．userのp値が0.05（有意水準0.05）を下回っているので，人に対しても有意差があるが，ここでは薬についてのみ注目しているので，無視することにする．
  
  参考文献：http://english-writing-theaters.com/7_statistics.html
6. 二元配置分散分析（Two-way analysis of variance）
  目的変数に対する2つの因子が与える影響に関する分析を行う方法である．それぞれの因子の影響だけでなく，それら因子の交互作用もみることができる．対象のデータが正規分布に従う場合のみ適用できる．二元配置分散分析にはいくつかの種類があるが，ここでは最も基本的な繰り返しのある二元配置分散分析（対応なし）を紹介する．（繰り返し数が等しい二元配置分散分析）
  例えば，3つの群（A, B, C）からなる生理指標データ（センサーで測定した値）があるとする．具体的には，薬A, B, Cを飲んだ後の生理指標データを測定する．なお，薬は朝に飲む場合（Mo）と夕方に飲む場合（Ev）の2通りがあるとする．
  すなわち，因子I（A, B, Cの3水準）と，因子II（Mo, Evの2水準）があるとする．この組み合わせの場合に対して，センサーでの測定値Yがある．
  
  実行方法は以下のようになる．
  
  > a1 = c(28,33,54,39,58,43,51,57,63,75,88,72)
  > a2 = c(34,52,19,37,29,44,55,28,33,43,32,21)
  > Y = c(a1,a2)
  > Time = factor(c(rep(“Mo”,12),rep(“Ev”,12)))
  > Medcine =factor(rep(c(rep(“A”,4), rep(“B”,4), rep(“C”,4)),2))
  > Y
  [1] 28 33 54 39 58 43 51 57 63 75 18 72 34 52 19 37 29 44 55 28 33 43 32 21
  > Time
  [1] Mo Mo Mo Mo Mo Mo Mo Mo Mo Mo Mo Mo Ev Ev Ev Ev Ev Ev Ev Ev Ev Ev Ev Ev
  Levels: Ev Mo
  > Medcine
  [1] A A A A B B B B C C C C A A A A B B B B C C C C
  Levels: A B C
  
  補足すると，データのうちの前半が朝に服薬した場合で，こう後半が夕方に服薬した場合．データの前半部分のうち，さらに前半の1/3が薬A，中盤の1/3が薬B，後半の1/3が薬Cである．また，データの後半部分のうち，さらに前半の1/3が薬A，中盤の1/3が薬B，後半の1/3が薬Cである．因子Iと因子IIの各水準の組み合わせに対して（例えば薬Aで朝に服薬），4回の繰り返しがある（4つのデータがある）．
  
  交互作用を考慮しない分散分析
  > summary(aov(Y~Medcine+Time))
  Df Sum Sq Mean Sq F value Pr(>F)
  Medcine 2 1074 536.8 2.828 0.08287 .
  Time 1 2281 2281.5 12.020 0.00243 **
  Residuals 20 3796 189.8
  —
  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  となり，Timeのみ有意差あり．朝に服薬するか，夕方に服薬するかで，測定値は異なる．
  
  交互作用まで見る場合．
  > summary(aov(Y~Medcine*Time))
  Df Sum Sq Mean Sq F value Pr(>F)
  Medcine 2 1074 536.8 4.518 0.025698 *
  Time 1 2282 2281.5 19.204 0.000359 ***
  Medcine:Time 2 1658 828.9 6.977 0.005712 **
  Residuals 18 2138 118.8
  —
  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  となり，Medcine, Timeとも有意差あり．交互作用も有意差あり．薬の種類によっても測定値は異なり，朝に服薬するか，夕方に服薬するかで，測定値は異なる．
  
  なお，各群の平均値は以下のようになる．
  朝に服薬したデータの平均値
  > mean(Y[c(1:12)])
  [1] 55.08333
  
  夕方に服薬したデータの平均値
  > mean(Y[c(13:24)])
  [1] 35.58333
  
  薬Aのデータの平均値
  > mean(c(Y[c(1:4)],Y[c(13:16)]))
  [1] 37
  
  薬Bのデータの平均値
  > mean(c(Y[c(5:8)],Y[c(17:20)]))
  [1] 45.625
  
  薬Cのデータの平均値
  > mean(c(Y[c(9:12)],Y[c(21:24)]))
  [1] 53.375
  
  最後に，
  > interaction.plot(Medcine, Time, Y)
  とすると，測定値Yの平均が，MedcineとTimeと水準ごとに計算した際，どのようなグラフになるかを表示してくれる．交互作用がなければ，すべての線が平行に近くなるが，上記のデータのように交互作用があると，線が大きくかけ離れた（あるいはクロスした）グラフになる．
  
  なお，以下のようにデータフレームのまま分析することもできる．
  > DF = data.frame(Medcine=factor(rep(c(rep(“A”,4), rep(“B”,4), rep(“C”,4)),2)), Time=factor(c(rep(“Mo”,12),rep(“Ev”,12))), Y = c(a1,a2))
  
  > DF
  Medcine Time Y
  1 A Mo 28
  2 A Mo 33
  3 A Mo 54
  4 A Mo 39
  5 B Mo 58
  ・・・
  20 B Ev 28
  21 C Ev 33
  22 C Ev 43
  23 C Ev 32
  24 C Ev 21
  
  > summary(aov(Y~Medcine+Time,data=DF))
  Df Sum Sq Mean Sq F value Pr(>F)
  Medcine 2 1074 536.8 2.828 0.08287 .
  Time 1 2281 2281.5 12.020 0.00243 **
  Residuals 20 3796 189.8
  —
  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  > summary(aov(Y~Medcine*Time,data=DF))
  Df Sum Sq Mean Sq F value Pr(>F)
  Medcine 2 1074 536.8 4.518 0.025698 *
  Time 1 2282 2281.5 19.204 0.000359 ***
  Medcine:Time 2 1658 828.9 6.977 0.005712 **
  Residuals 18 2138 118.8
  —
  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘
  
  参考文献1：https://www1.doshisha.ac.jp/~mjin/R/Chap_13/13.html
  参考文献2：https://www1.doshisha.ac.jp/~mjin/R/Chap_13/13.html
  参考文献3：https://www.recruit-ms.co.jp/issue/column/0000000667/?theme=hrdepartment
7. 多重比較（multiple comparison）：ボンフェローニの調整によるt検定（t test with Bonferroni correction）
  3群（正規性が保証されている）以上のデータの平均値の比較において，一元配置分散分析により有意差が確認できた場合に，どの群とどの群の間に平均値に差があるかを検定する方法である．実施方法としては，単純に2群の平均値の差をt検定により検定すればよいだけである．ただし，ボンフェローニの調整では，有意水準を必要な検定回数で割って低めに調整する．
  例えば，有意水準が0.05の時を考える．群の数が3の場合，群間の比較はA群-B群，B群-C群，C群-A群の3つになる．この場合，3回の比較が必要になり，有意水準0.05を3で割った値（0.017）を，新たな有意水準と設定する．群の数が4の場合，群間の比較は，A群-B群，A群-C群，A群-D群，B群-C群，B群-D群，C群-D群の6つになる．この場合，6回の比較が必要になり，有意水準0.05を6で割った値（0.017）を，新たな有意水準と設定する．
  各群間で普通にt検定を行い，上記の新たな有意水準で判定すればよい．
  実行例は以下のとおりである．
  > psyA
  [1] 15 9 18 14 18
  > psyB
  [1] 10 7 3 5 7
  > psyC
  [1] 10 6 11 7 12
  
  > t.test(psyA,psyB,var.equal=TRUE)
  
  Two Sample t-test
  
  data: psyA and psyB
  t = 4.1485, df = 8, p-value = 0.003216
  alternative hypothesis: true difference in means is not equal to 0
  95 percent confidence interval:
  3.730698 13.069302
  sample estimates:
  mean of x mean of y
  14.8 6.4
  
  > t.test(psyA,psyC,var.equal=TRUE)
  
  Two Sample t-test
  
  data: psyA and psyC
  t = 2.7724, df = 8, p-value = 0.02421
  alternative hypothesis: true difference in means is not equal to 0
  95 percent confidence interval:
  0.9421 10.2579
  sample estimates:
  mean of x mean of y
  14.8 9.2
  
  > t.test(psyB,psyC,var.equal=TRUE)
  
  Two Sample t-test
  
  data: psyB and psyC
  t = -1.704, df = 8, p-value = 0.1268
  alternative hypothesis: true difference in means is not equal to 0
  95 percent confidence interval:
  -6.5891514 0.9891514
  sample estimates:
  mean of x mean of y
  6.4 9.2
  
  > 0.003216*3
  [1] 0.009648
  > 0.02421*3
  [1] 0.07263
  > 0.1268*3
  [1] 0.3804
  
  上記のように算出したp値を3倍にしてこれを0.05と比較する．この結果，A群-B群間で有意差があることが分かる．
  
  ボンフェローニの調整は，群の数が増えると有意水準が飛躍的に小さくなり，かなり厳しい検定となる．群の数が5以上になると，その厳しさが顕著になってくるので，その場合はチューキー・クレーマー検定（パラメトリック法のみ）を用いる方が良い．
  
  （参考）
  https://www.youtube.com/watch?v=HLzS5wPqWR0
8. 多重比較（multiple comparison）：チューキー・クレーマー検定 (Tukey-Kramer test)
  3群（正規性が保証されている）以上のデータの平均値の比較において，一元配置分散分析により有意差が確認できた場合に，どの群とどの群の間に平均値に差があるかを検定する方法である．群の数が多い場合にも適用できる．
  > TukeyHSD(aov(rev~group))
  Tukey multiple comparisons of means
  95% family-wise confidence level
  
  Fit: aov(formula = rev ~ group)
  
  $group
  diff lwr upr p adj
  B-A -8.4 -13.480629 -3.3193714 0.0022607
  C-A -5.6 -10.680629 -0.5193714 0.0308278
  C-B 2.8 -2.280629 7.8806286 0.3386342
  この結果から，A群-B群間，A群-C群間で有意差（p値＜0.05）があることが分かる．
  
  参考記事：https://data-science.gr.jp/implementation/ist_r_tukey_kramer_test.html
  参考記事：http://turtlewalk.hatenablog.jp/entry/2016/01/03/002050
  参考記事：http://turtlewalk.hatenablog.jp/entry/2016/01/04/020246
9. 多重比較（multiple comparison）：ホルム検定 (Holm test)
  3群（正規性が保証されている）以上のデータの平均値の比較において，一元配置分散分析により有意差が確認できた場合に，どの群とどの群の間に平均値に差があるかを検定する方法である．ボンフェローニの方法は，組み合わせの数（定数）で有意水準（危険率）の割り算をして調整したが，これでは調整方法としては厳しくなりがちである．
  ホルムの方法では，すべての組み合わせでp値を計算する．最も低いp値の組み合わせに対しては，有意水準（危険率）に対して，組み合わせ数cで割り算する（ここまではボンフェローニの方法と同じである）．次に小さなp値を持つ組み合わせに対しては，有意水準（危険率）に対して，組み合わせ数(c-1)で割り算する．すなわち，やや緩めに調整することになる．さらに次に小さなp値を持つ組み合わせに対しては，有意水準（危険率）に対して，組み合わせ数(c-2)で割り算する．このように，割る数を1ずつ減らしていく．有意な結果が得られなくなった時点で，検定を終了する．

1. 多重比較（multiple comparison）：スティールドゥワス (Steel-Dwass) の方法
  上記の方法では、いずれも群の数が多くなると、有意差が出にくくなってしまう。5群以上の場合はスティールドゥワスの方法が最も良いと思われる。
  スティールドゥワスの方法参考ページ
2. 分散分析（クラスカルウォリス検定）と多重比較
  分散分析（クラスカルウォリス検定）と多重比較の結果が異なることは，よくある．そもそも多重比較における補正方法にはいくつかの方法があり，それにより結果が変わってくるので，分散分析の結果と多重比較の結果が異なることは起こりうることである．検定で最も知りたいことは，何の群が特に値が高いかであると思われるので，分散分析を行わなくても，多重比較だけで目的を達成することはできる．分散分析の結果と多重比較の結果が異なるときは，多重比較の結果だけを載せるという手もありかと思われる（ただし，論文全体を通じて一貫している必要がある）．データ数なども考慮して，検定結果の信頼性を判断して，結論付けることが重要である．
  参考：https://personal.hs.hirosaki-u.ac.jp/pteiki/research/stat/multi.pdf (p.17)

1. フリードマン検定 (Friedmanの検定)
  ノンパラメトリックな対応のある分散分析．3変数以上の中央値を比較し差があるかないかを検定する．一人の被験者が，条件1, 2, 3, ….とすべての条件を実施した場合に，条件間で目的変数の値が異なるかどうかを検定する．（実際には，条件間に違いがあるかどうかまでは分からず，条件と目的変数に関係があるかどうかのみ検定する）
  # 7人が3つの条件すべてを実施．ベクトルのn番目の値はすべて，同じ人のデータである．
  > vx=c(3.7, 3.1, 2.5, 2.6, 3.2, 3.2, 2.5)
  > vy=c(3.8, 2.7, 4.0, 2.4, 2.2, 3.4, 5.5)
  > vz=c(2.8, 3.4, 4.5, 2.2, 3.0, 3.1, 3.4)
  > D=matrix(c(vx,vy,vz), ncol=3, nrow=7)
  > friedman.test(D)
  
  Friedman rank sum test
  
  data: D
  Friedman chi-squared = 0.28571, df = 2, p-value = 0.8669

1. 外れ値の削除
  外れ値を削除するには，以下の3つの方法がある．
  (1) ヒューリスティックスを用いる方法
  （箱ひげ図や標準偏差（2σや3σ）を用いる方法）
  (2) 検定を行う方法
  （Smirnov-Grubbs検定（正規分布ベース））
  (3) 計算論的アルゴリズムを用いる方法
  （階層的クラスタリング）
  以下，(1)の詳細を説明する．
  ■箱ひげ図を用いる方法
  「第三四分位数+1.5×IQR」より上にあるものを外れ値，また「第一四分位数-1.5×IQR」より下にあるものを外れ値とみなす．
  ここで，IQRは四分位範囲「75パーセンタイル（第三四分位数）－25パーセンタイル（第一四分位数）」を表す．基本的には，データが正規分布に従うことを前提にしている．
  ■標準偏差（2σや3σ）を用いる方法
  単純に取得したデータから算出した標本の標準偏差σを用いて，+2σ（または3σ）より上にあるものを外れ値，また-2σ（または3σ）より下にあるものを外れ値とみなす．有意水準の判定で5%や1%という値が使われることがあるが，これに近い考え方と言える．分布の十分外側にあるという理由（例えば，+2σの外側には2.3%，-2σの外側には2.3%しかデータが存在しない）から除去しているに過ぎない．基本的には，データが正規分布に従うことを前提にしている．
  ※機械的な方法ではないが，ヒストグラムを書き，視覚的にデータが中心のクラスタから外れていることを確認して，手動で削除することも考えられる．その外れ値が物理的に（常識的に）実現可能かどうかを考慮しつつ決定できる点が優れているが，査読者を納得させにくい欠点もある．

1. 箱ひげ図を用いて外れ値を削除するプログラム
  > DFH = read.table(“height.csv”, sep=”,”, header=TRUE)
  > data = DFH$height
  > data
  [1] 176.2202 177.4435 173.8930 165.1857 166.7457 168.0908 179.7985 173.0341
  ・・・
  [97] 169.5849 175.7688 166.1920 172.0525
  > source(“extremeData.r”)
  > extremeData(data, mode=TRUE)
  [1] 176.2202 177.4435 173.8930 165.1857 166.7457 168.0908 179.7985 173.0341
  [9] 171.7688 167.0790 175.0525 167.7561 176.0641 168.7115 177.3558 173.0537
  ・・・
  [89] 178.6754 168.0634 169.5849 175.7688 166.1920 172.0525 0.0000 0.0000
  [97] 0.0000 0.0000 0.0000 0.0000
  > extremeData(data, mode=FALSE)
  [1] 192.43 192.55 154.33 148.55 152.12 203.44 0.00 0.00 0.00 0.00
  ・・・
  [91] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
  
  プログラムとデータのファイル：extremeData.zip （ただし暗号化圧縮ファイル．土方ゼミ学生のみ利用可能）

1. 程度を問う質問の選択肢（リッカート尺度ほか）
  以下を基準に，聞きたい内容に合わせてカスタマイズすること（絶対基準ではない）．対立極がある場合は，”中立”を設けるのが一般的．対立極がない場合（片側単一極）の場合は，”中立”を設けるか否か，選択肢を奇数にするか偶数にするかは，研究者の考え方に依存し，決まったルールは存在しない．
  
  ＜”中立”がある回答の場合＞
  リッカート尺度（Likert scale）に従う．
  
  質問例：「自分は気晴らしに，コーヒーを飲むことが多い」
  (7段階)
  全く当てはまらない－当てはまらない－どちらかと言えば当てはまらない－どちらとも言えない－どちらかと言えば当てはまる－当てはまる－非常によく当てはまる
  Strongly disagree – Disagree – Somewhat disagree – Neither agree nor disagree – Somewhat agree – Agree – Strongly agree
  Disagree strongly – Disagree moderately – Disagree a little – Neither agree nor disagree – Agree a little – Agree moderately – Agree strongly
  (5段階)
  当てはまらない－やや当てはまらない－どちらとも言えない－やや当てはまる－当てはまる
  Disagree – Somewhat disagree – Neither agree nor disagree – Somewhat agree – Agree
  
  質問例：「政策に意見をするのであれば，選挙に行くべきだ」
  (7段階)
  全く同意できない－同意できない－どちらかと言えば同意できない－どちらともいえない－どちらかと言えば同意できる－同意できる－非常に同意できる
  Strongly disagree – Disagree – Somewhat disagree – Neither agree nor disagree – Somewhat agree – Agree – Strongly agree
  Disagree strongly – Disagree moderately – Disagree a little – Neither agree nor disagree – Agree a little – Agree moderately – Agree strongly
  (5段階)
  同意できない－やや同意できない－どちらともいえない－やや同意できる－同意できる
  Disagree – Somewhat disagree – Neither agree nor disagree – Somewhat agree – Agree
  
  質問例：「あなたは自分のことを，ひかえめで，おとなしいと思う」
  (7段階)
  全く違うと思う－おおよそ違うと思う－少し違うと思う－どちらともいえない－少しそう思う－まあまあそう思う－強くそう思う
  Strongly disagree – Disagree – Somewhat disagree – Neither agree nor disagree – Somewhat agree – Agree – Strongly agree
  Disagree strongly – Disagree moderately – Disagree a little – Neither agree nor disagree – Agree a little – Agree moderately – Agree strongly
  (5段階)
  違うと思う－やや違うと思う－どちらともいえない－ややそう思う－そう思う
  Disagree – Somewhat disagree – Neither agree nor disagree – Somewhat agree – Agree
  
  質問例：「あなたはこの映画を」
  (7段階)
  非常に嫌い－嫌い－どちらかと言えば嫌いーどちらとも言えない－どちらかと言えば好き－好き－非常に好き
  全く好きではないー好きではない－どちらかと言えば好きではないーどちらとも言えない－どちらかと言えば好き－好き－非常に好き
  Definitely deslike – Dislike – Somewhat dislike – Neither like nor dislike – Somewhat like – Like – Extremely like
  (5段階)
  嫌い－やや嫌い－どちらとも言えない－やや好き－好き
  Dislike – Somewhat dislike – Neither like nor dislike – Somewhat like – Like
  
  質問例：「あなたにとって，この政策は，」
  (7段階)
  まったく重要ではないー重要ではないーどちらかと言えば重要ではない－どちらともいえないーどちらかと言えば重要である－重要であるー非常に重要である
  Definitely unimportant – Umimportant – Somewhat unimportant – Neither important nor unimportant – Somewhat important – Important – Definitely important
  (5段階)
  重要ではない－やや重要ではない－どちらともいえないーやや重要である－重要である
  Umimportant – Somewhat unimportant – Neither important nor unimportant – Somewhat important – Important
  
  質問例：「あなたは，このサービスに，」
  (7段階)
  まったく満足していないー満足していないーどちらかと言えば満足していない－どちらともいえないーどちらかと言えば満足である－満足であるー非常に満足である
  Completely dissatisfied – Mostly dissatisfied – Somewhat dissatisfied – Neither satisfied nor dissatisfied – Somewhat satisfied – Mostly satisfied – Completely satisfied
  (5段階)
  満足していない－やや満足していない－どちらともいえないーやや満足である－満足である
  Dissatisfied – Somewhat dissatisfied – Neither satisfied nor dissatisfied – Somewhat satisfied – Satisfied
  
  ＜”中立”がない回答の場合＞
  質問例：「あなたはYouTubeのお薦め機能を」
  (4段階)
  まったく利用しない－ほとんど利用しない－少し利用する－よく利用する
  Never use it at all – Hardly use it at all – Use it a little – Use it very often
  (5段階)
  まったく利用しない－ほとんど利用しない－少し利用する－（かなり）利用する－非常によく利用する　（「かなり」はない方が良いが，どうしても英語に対応させる必要があれば，入れても良い）
  Never use it at all – Hardly use it at all – Use it a little – Use it often – Use it very often
  (6段階)
  まったく利用しない－ほとんど利用しない－少し利用する－利用する－かなり利用する－非常によく利用する
  （該当する英語選択肢はない）
  
  質問例：「あなたは買い物をするときに，クレジットカードをどの程度の割合で使いますか？」
  (4段階)
  まったく利用しない－ほとんど利用しない－ほぼ毎回利用する－毎回利用する
  Never – Almost never – Almost every time – Every time
  (5段階)
  まったく利用しない－ほとんど利用しない－時々利用する－ほぼ毎回利用する－毎回利用する
  Never – Almost never – Sometimes – Almost every time – Every time
  
  質問例：「最近1か月の間に，緊張やストレスを感じたことがどれくらいありましたか？」
  (5段階)
  全くない – ほとんどない – 時々ある – かなりある – 非常にある
  
  1つの質問で尋ねた場合は，いずれも順序尺度として扱う（ノンパラメトリック）．
  多くの質問（個人的な経験では4つ程度以上）で尋ねた場合は，比例尺度として扱っても良い場合がある（パラメトリック）
  私見ではあるが，英語の場合は，Agreeのように単独の動詞で回答させるよりも，Moderately agreeのように副詞を伴った選択肢を用意する方が良い気がする．
  
  参考記事：https://www.nhk.or.jp/bunken/summary/yoron/method/pdf/040901.pdf
  https://insight.rakuten.co.jp/knowledge/researchcolumn/vol3.html

1. 天井効果と床効果
  前記項目のように，リッカート尺度等でアンケートを取った時に，回答者からの回答値が段階の高い方，または低い方に極端に偏った場合をそれぞれ天井効果と床効果と呼ぶ．例えば1-7の7段階でアンケートを取った時に，ユーザからの回答が6や7に偏っている場合（天井効果），1や2に偏っている場合（床効果）と呼ぶ．統計的には，平均＋標準偏差が7を超えている場合や，平均－標準偏差が1を下回っている場合は，それぞれ天井効果と床効果を引き起こしているといえる．
  直感からわかるように，回答にばらつきがないことから，質問としてほとんど情報量のない内容になっている．このような現象が起きる理由としては，被験者の詳細な差を検出したいのにもかかわらず，5段階や7段階でアンケートを取った時に，そのような詳細な差を表現する回答選択肢になっていないことが挙げられる（例：「1. 自分はサルより知能が劣っていると感じる」，「2. 自分の知能はサル並みであると感じる」，「3. 自分の知能はサルより少し高い程度であると感じる」，「5. 自分の知能は人間より少し劣っている程度であると感じる」，「5. 自分の知能は人間並みであると感じる」で尋ねたら，ほとんどの人が5と答えるであろう）．
  心理尺度は，個人の差が出るように計測されるべきである．そのため，例えば7段階でアンケートを取ったのであれば，1と回答する人もいれば，7と回答する人もいなければ，それだけの粒度を持たせた意味がない．

1. 検者間信頼性評価
  検者間信頼性を評価するには，カッパ係数とICC(2,1)がある．以下の基準で使い分ける．
  
  間隔尺度，比例尺度，段階の多い順序尺度：ICC(2,1)
  名義尺度，段階の少ない順序尺度：カッパ係数
  
  参考：
  級内相関係数 (ICC)について
2. Cohenのカッパ係数
  複数人（2名のみ）が付与したラベルの妥当性の検証（ラベルは2以上なら何でもOK)
  評価者間の一致度 (inter-rater reliability)を測る指標の一つ
  （例）Aさん，Bさんが，画像に対して，ネコが写っているかどうかを判定（isCat = True or False）（2値判定）
  （例）Aさん，Bさんが，画像に対して，ネコかイヌかその他が写っているかどうかを判定（typeAnimal = cat or dog or other）（3値判定）
  以下，実行例（2値判定）
  > library(irr)
  > rater01=c(1,2,1,2,1,1,2,2,1,2,1,1,1,2,1,2)
  > rater02=c(2,2,2,1,1,1,2,1,1,1,2,1,2,2,2,2)
  > D2=cbind(rater01, rater02)
  > kappa2(D2)
  
  以下，実行例（3値判定）
  > DF = read.table(“rater2.csv”,header=TRUE,sep=”,”)
  > D = cbind(DF$u1,DF$u2)
  > kappa2(D)
  
  データソース：https://soc-research.org/ja/wp-content/uploads/2019/07/rater2.zip
  
  一致の判断（Landis and Koch (1977)）
  0.0〜0.2: わずかに一致(slight agreement)
  0.21〜0.40 まずまずの一致(fair agreement)
  0.41〜0.60 中等度の一致(moderate agreement)
  0.61〜0.80 かなりの一致(substantial agreement)
  0.81〜1.0 ほぼ完全、完全一致(almost perfect or perfect agreement)
  
  Cohen, J. (1960). A coefficient of agreement for nominal scales. Educational and Psychological Measurement, 20, 37-46.
  Landis JR, Koch GG. The measurement of observer agreement for categorical data. Biometrics. 1977; 33(1):159-174.
  参考記事：http://ides.hatenablog.com/entry/20170417/1492409363

1. 1. Fleissのカッパ係数
    複数人（3名以上）が付与したラベルの妥当性の検証
    評価者間の一致度 (inter-rater reliability)を測る指標の一つ
    （例）Aさん，Bさん，Cさんが，画像に対して，ネコが移っているかどうかを判定（isCat = True or False）（2値判定）
    以下，実行例（2値判定の場合）
    > library(irr)
    > rater01=c(1,2,1,2,1,1,2,2,1,2,1,1,1,2,1,2)
    > rater02=c(2,2,2,1,1,1,2,1,1,1,2,1,2,2,2,2)
    > rater03=c(2,2,2,1,1,1,2,1,1,1,2,2,2,2,2,1)
    > D3=cbind(rater01, rater02, rater03)
    > kappam.fleiss(D3)
    
    3値以上の判定も可
    
    Fleiss, J.L.: Measuring nominal scale agreement among many raters, Psychological Bulletin, Vol. 76, pp. 378-382, 1971.
    
    参考記事：https://www.datanovia.com/en/lessons/fleiss-kappa-in-r-for-multiple-categorical-variables/
  2. Krippendorff’s alpha (クリッペンドルフのα係数)
    Krippendorff’s alpha は，複数人の評価者が付与したラベルや評価値の一致度（一貫性） (inter-rater reliability)を測る指標の一つである．
    ・評価者が3人以上いても使える
    ・データが名義尺度・順序尺度・間隔尺度・比例尺度のどれでも使える
    ・欠損値があっても使える
    と，汎用性の高い手法である．
    
    install.packages(“irr”)
    library(irr)
    R1 = c(0,2,2,1,3,0,4,4,6)
    R2 = c(1,4,1,2,3,1,5,3,5)
    R3 = c(0,2,0,0,4,0,4,4,6)
    R4 = c(1,2,2,1,3,0,4,5,4)
    R5 = c(0,3,2,1,2,0,4,4,6)
    data = rbind(R1,R2,R3,R4,R5)
    kripp.alpha(data, method = “interval”)
    Krippendorff’s alpha
    
    Subjects = 9
    Raters = 5
    alpha = 0.851
    
    (参考文献)
    https://zenn.dev/hellorusk/articles/a75b374bb77c87 (ただしPython)
  3. ICC(2,1)
    検者間信頼性を測る方法．評価尺度は，間隔尺度，比例尺度，段階の多い順序尺度に限る．
    
    実行方法
    R1 = c(0,2,2,1,3,0,4,4,6)
    R2 = c(1,4,1,2,3,1,5,3,5)
    R3 = c(0,2,0,0,4,0,4,4,6)
    R4 = c(1,2,2,1,3,0,4,5,4)
    R5 = c(0,3,2,1,2,0,4,4,6)
    data = cbind(R1,R2,R3,R4,R5)
    rownames(data) = c(1:9)
    data
    library(irr)
    icc(data, model=”twoway”, type=”agreement”)
    
    参考記事：https://toukeier.hatenablog.com/entry/how-to-determine-sample-size-for-inter-rater-reliability-icc-case2/
  4. コンジョイント分析
    商品やサービスのどの特徴が消費者の購買行動につながるのかを、商品の事例を挙げて、それに対して被験者に評価付けを行わせることで、分析する方法である。購買行動と書いたが、SNSでのシェア、いいね！するか、妬みを抱くかなど、あらゆる行動に対して適用することができる。
    
    後日、詳細を更新する
    
    https://www.macromill.com/service/data_analysis/d012.html
    https://istat.co.jp/ta_commentary/conjoint_analysis
  5. 実験計画法（ラテン方格）
    作物が大きくなるかどうか，人が物を購入するか否か，それに影響する要因/特徴（因子）を調べる時には，いくつかの因子を候補に挙げ，その要因の特徴量/カテゴリ値（水準）を決め，実験により影響する因子を明らかにする．しかし，因子と水準の数が多くなると，すべての組み合わせで実験を行おうとすると，その組み合わせ数は膨大なものになる．例えば，3因子3水準で実験を行うとすると，その組み合わせ数は3×3×3=27パターンになる．このすべてのパターンに対して実験を行うとコストがかかるため，これを削減する方法が必要となる．その方法が，実験計画法である．
    理論的妥当性は，下記の参考ページや専門書に任せるが，機械的に削減する方法としてラテン方格法がある．このラテン方格法は，すでに用意された直交表を用いて，実験するパターンを削減するが，総当たりの場合と変わらぬ効果を期待できるものである．直交表は，2水準の場合と3水準の場合に分けられる．
    2水準の場合は，
    ・直交表:2水準 3因子
    ・直交表:2水準 7因子
    ・直交表:2水準 15因子
    ・直交表:2水準 31因子
    などを使う．それぞれ，因子の数を3, 7, 15まで試行することができる．
    3水準の場合は，
    ・L9直交表:3水準 4因子
    ・直交表:3水準 13因子
    などを使う．それぞれ，因子の数を4, 13まで試行することができる．
    参考：http://www.placeon.jp/blog/words/words_monozukuri/experimental_design/
    
    ただし，実験計画法は，これに代わる方法論があるわけではないが，使える場面は，非常に限定される．すなわち，すべての因子の水準が同一でないと，基本的には使えないのである．ただし，以下のようにいくつかの例外はある．この例外で賄えるのであれば，実験計画法は利用可能である．
    (1) 2水準の直交表において，1つだけ4水準の因子を持つ直交表を作る場合．
    2水準の直交表において，3列を統合して4水準を割り付けることで実現できる．これを，「多水準作成法」と呼ぶ．ただし，2つ以上4水準の因子を割り当てることはできない（と思われる）．
    (2) 2水準の直交表において，1つだけ3水準の因子を持つ直交表を作る場合．
    3列を使って一旦4水準をつくり、それに対して(3)のダミー法で使い3水準の因子を作る．
    (3) 3水準の直交表において，いくつか2水準の因子を作る場合．
    3水準の直交表において，一つの該当の因子を3水準のうちの3水準目のカテゴリ値に対して，1水準目のカテゴリ値か，2水準目のカテゴリ値を割り当てることで実現できる．「ダミー法」と呼ばれる．
    参考：https://www.monodukuri.com/gihou/article/75
    主要な直交表の種類は，以下を参照
    参考：https://www.monodukuri.com/gihou/article/74
    
    注) おそらく3種類の水準が混ざった因子群に対して，実験計画法を実行することはできないと思われる．
  6. クロンバックのα係数
    いくつかの質問により一つの心理指標を計算したいときに，その質問の妥当性（信頼性）を算出するための統計的指標である．例えば，鬱（うつ）の程度を一つの指標で取得するのに，
    Q1: 私の心はいつも晴れない
    Q2: 胸が張り裂けるほど苦しいときがある
    Q3: 死にたいと思うことがある
    Q4: 人と会いたくないことがある
    のような質問をリッカート尺度（例えば5段階）で尋ねておき，その総和で鬱（うつ）の程度を取得する．
    クロンバックのα係数を計算するライブラリはたくさんあるが，ltmパッケージが最も使い方が簡単である．
    > DF = read.csv(“depression_score.csv”)
    > library(ltm) # ltmパッケージの読み込み
    > cronbach.alpha(DF[,-1]) # [,-1]は1列目（ユーザ名）を含まない指示
    Cronbach’s alpha for the ‘DF[, -1]’ data-set
    Items: 4
    Sample units: 20
    alpha: 0.891
    
    ちゃんと確立された心理尺度を利用した場合は，クロンバックのα係数は0.9前後の値になるはずである．自分たちで作成した心理尺度や評価指標であれば，0.8以上あれば概ね信頼できるといえる（もちろん意味的に保障できるものではない）．0.7を下回ると，質問項目を設計し直すか，そもそも計測可能でない概念を扱おうとしている可能性があるので，検討し直した方が良い．
    
    データソース：https://soc-research.org/ja/wp-content/uploads/2020/10/depression_score.zip
  7. 単回帰分析と重回帰分析
    一つの説明変数と目的変数の間の線形の関係を推定したい場合は単回帰分析を行う．複数の説明変数と目的変数の間の線形の関係を推定したい場合は重回帰分析を行う．重回帰分析では，複数の説明変数を同時に考慮した際に，どの説明変数が有意に目的変数に相関するのかを明らかにすることができ，またP値や（標準化後の説明変数で分析した場合の）標準化偏回帰係数の大きさより，どれほどの大きさで相関しているのかを明らかにすることができる．なお，目的変数は量的データである．
    
    本説明では，哺乳類の寿命を目的変数として，行動上の特徴を説明変数とした，以下のファイルを基に説明する．
    mammals.csv [ダウンロード]
    Source：http://www.statsci.org/data/general/sleep.txt を修正したもの
    
    このファイルは，Species：哺乳類の種類，BodyWt：体重，BrainWt：脳の大きさ，NonDreaming：夢を見ていない睡眠時間，Dreaming：夢を見ている睡眠時間，TotalSleep：1日の睡眠時間，LifeSpan：寿命，Gestation：妊娠期間(日)，Predation：捕食されやすさ(5段階)，Exposure：就寝中の襲われやすさ(5段階)，Danger：他の動物からの危険の程度(5段階)，を表す．
    
    ＜データフレーム読み込み＞
    > DF = read.csv(“mammals.csv”)
    
    ・単回帰分析
    ここでは，LifeSpanを目的変数，TotalSleepを説明変数とする．
    ＜データフレームから別変数に保存して分析＞
    > LifeSpan = DF$LifeSpan
    > TotalSleep = DF$TotalSleep
    > plot(TotalSleep, LifeSpan)
    > m = lm(目的変数 ~ 説明変数)
    > m = lm(LifeSpan ~ TotalSleep)
    > summary(m)
    
    Call:
    lm(formula = LifeSpan ~ TotalSleep)
    
    Residuals:
    Min 1Q Median 3Q Max
    -21.706 -10.609 -2.382 6.310 76.440
    
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 36.3163 5.9993 6.053 1.69e-07 ***
    TotalSleep -1.5946 0.5214 -3.058 0.00354 **
    —
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 17.46 on 51 degrees of freedom
    (8 observations deleted due to missingness)
    Multiple R-squared: 0.155, Adjusted R-squared: 0.1384
    F-statistic: 9.353 on 1 and 51 DF, p-value: 0.003541
    
    上記結果において，Adjusted R-squaredは，自由度調整済み決定係数（R*2）のことで，この回帰分析の正確さ（モデルの良さ）を表す重要な指標である．
    
    # plot(TotalSleep, LifeSpan)に対して回帰直線を引く
    > abline(m)
    
    ＜データフレームに対して分析＞
    > m = lm(目的変数名 ~ 説明変数名, data=データフレーム名)
    ※dataでデータフレームを指定する
    > m = lm(LifeSpan ~ TotalSleep, data=DF)
    > summary(m)
    
    Call:
    lm(formula = LifeSpan ~ TotalSleep, data = DF)
    
    Residuals:
    Min 1Q Median 3Q Max
    -21.706 -10.609 -2.382 6.310 76.440
    
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 36.3163 5.9993 6.053 1.69e-07 ***
    TotalSleep -1.5946 0.5214 -3.058 0.00354 **
    —
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 17.46 on 51 degrees of freedom
    (8 observations deleted due to missingness)
    Multiple R-squared: 0.155, Adjusted R-squared: 0.1384
    F-statistic: 9.353 on 1 and 51 DF, p-value: 0.003541
    
    ・重回帰分析
    ここでは，LifeSpanを目的変数，TotalSleepとBodyWtを説明変数とする．
    ＜データフレームから別変数に保存して分析＞
    > BodyWt = DF$BodyWt
    > plot(BodyWt,LifeSpan)
    > m = lm(目的変数 ~ 説明変数1 + 説明変数2 + ・・・)
    > m = lm(LifeSpan ~ TotalSleep + BodyWt)
    > summary(m)
    
    Call:
    lm(formula = LifeSpan ~ TotalSleep + BodyWt)
    
    Residuals:
    Min 1Q Median 3Q Max
    -17.497 -9.849 -3.587 5.966 78.244
    
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 29.458288 6.012660 4.899 1.05e-05 ***
    TotalSleep -1.117298 0.509073 -2.195 0.03285 *
    BodyWt 0.019934 0.006583 3.028 0.00389 **
    —
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 16.21 on 50 degrees of freedom
    (8 observations deleted due to missingness)
    Multiple R-squared: 0.2859, Adjusted R-squared: 0.2574
    F-statistic: 10.01 on 2 and 50 DF, p-value: 0.0002206
    
    ＜データフレームに対して分析＞
    > m = lm(目的変数 ~ 説明変数1 + 説明変数2 + ・・・, data=データフレーム名)
    ※dataでデータフレームを指定する
    m = lm(LifeSpan ~ BodyWt + TotalSleep, data=DF)
    > summary(m)
    
    Call:
    lm(formula = LifeSpan ~ BodyWt + TotalSleep, data = DF)
    
    Residuals:
    Min 1Q Median 3Q Max
    -17.497 -9.849 -3.587 5.966 78.244
    
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 29.458288 6.012660 4.899 1.05e-05 ***
    BodyWt 0.019934 0.006583 3.028 0.00389 **
    TotalSleep -1.117298 0.509073 -2.195 0.03285 *
    —
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 16.21 on 50 degrees of freedom
    (8 observations deleted due to missingness)
    Multiple R-squared: 0.2859, Adjusted R-squared: 0.2574
    F-statistic: 10.01 on 2 and 50 DF, p-value: 0.0002206
    
    ＜説明変数を標準化して分析＞
    説明変数を標準化（基準化）しておくと，どの説明変数が目的変数と強く相関しているかを標準化偏回帰係数（Estimate）から判断しやすい．
    標準化（基準化）してzスコアを求めるには，scale関数を使うと便利．
    
    > zBodyWt = scale(BodyWt)
    > zTotalSleep = scale(TotalSleep)
    > m = lm(LifeSpan ~ zBodyWt + zTotalSleep)
    > summary(m)
    
    Call:
    lm(formula = LifeSpan ~ zBodyWt + zTotalSleep)
    
    Residuals:
    Min 1Q Median 3Q Max
    -17.497 -9.849 -3.587 5.966 78.244
    
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 19.401 2.227 8.711 1.35e-11 ***
    zBodyWt 6.792 2.243 3.028 0.00389 **
    zTotalSleep -5.077 2.313 -2.195 0.03285 *
    —
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 16.21 on 50 degrees of freedom
    (8 observations deleted due to missingness)
    Multiple R-squared: 0.2859, Adjusted R-squared: 0.2574
    F-statistic: 10.01 on 2 and 50 DF, p-value: 0.0002206
    
    ※上の標準化前と標準化後を比べると，P値は変化はないが，偏回帰係数（Estimate）だけ変化していることが分かる．
    
    ＜重要な補足事項＞
    ・mammals.csvには，いくつかのセルで欠損値（NA）がある．lm()関数は，分析に関与する変数（目的変数と説明変数）のどれかが NA のものは分析に使われない．
    ・scale()関数は，変数ごとに NA を自動的に除いて，平均値=0，標準偏差=1 になるように標準化する．
    ・説明変数と目的変数の値を標準化して，厳密に分析するには（切片項を完全に0にするには），元のデータフレームからNAがある行をあらかじめ削除しておき，それに対してscale()関数で標準化を行い，重回帰分析（または単回帰分析）を行う
    
    なお，重回帰分析を行うときには，標準化せずに行った場合と，標準化して行った場合の2つの方法の結果を載せた方が良い．具体的には，以下の3つを載せておくと安心である．（変数の効果の大きさ（絶対値）も分かるうえに，変数間の効果の相対的な違いも理解できるようになる）
    ・重回帰分析におけるB、SE B、β
    B（回帰係数、非標準化回帰係数）
    独立変数が従属変数に与える影響の大きさを示す。具体的には、独立変数が1単位増加したときに従属変数がどれだけ変化するかを表す。これは、非標準化された値であり、変数のスケールに依存しているといえる。
    
    SE B（標準誤差）
    回帰係数Bの標準誤差を表す。推定された回帰係数がどれだけ不確実であるか、つまり推定にどれだけ誤差が含まれているかを示すものである。SE Bが小さいほど、Bの推定値が信頼できるものになる。
    
    β（標準化回帰係数）
    標準化された回帰係数のこと。独立変数と従属変数が標準化された値（平均が0、標準偏差が1）に基づいて計算されるため、異なる単位やスケールを持つ変数同士の影響の大きさを比較することが可能になる。
    
    重回帰分析の結果を載せるとき，
    ＜各目的変数に対して＞
    標準化偏回帰係数，偏回帰係数（場合による：特に測定尺度の絶対値が物理量の場合に必要），標準誤差，ｔ値、p値
    ＜分析結果全体に対して＞
    F値（書き方は，F(説明変数の数, 自由度 (=データ数-(説明変数の数+切片)），F値の有意確率，修正済みR2乗値
    を出しておけば完璧．
    F値：回帰の分散の大きさを統計的に判断した結果を示す数値
    
    参考：https://ai-trend.jp/programming/r-beginner/r-3/
    参考：https://stats.biopapyrus.jp/glm/lm-r.html
偏回帰係数や独立変数に対するp値の意味（何の検定をしているか）は、以下のHPが参考になる。
参考：https://best-biostatistics.com/summary/henkaiki.html
単回帰分析・重回帰分析の説明で用いたmammals.csvを用いる．

＜データフレーム読み込み＞
> DF = read.csv(“mammals.csv”)

LifeSpanを目的変数，TotalSleepとBodyWtを説明変数とする．
ただし，LifeSpanは，以下のように30年以上生きる場合を1に，それ未満の場合を0にする．
> bLifeSpan = ifelse(LifeSpan>=30, 1, 0)

なお，以下のように論理値型に変換しても，ロジスティック回帰分析を行ことができる．
> bLifeSpan = ifelse(LifeSpan>=30, TRUE, FALSE)

＜データフレームから別変数に保存して分析＞
> BodyWt = DF$BodyWt
> plot(BodyWt,bLifeSpan)
> res = glm(目的変数 ~ 説明変数1 + 説明変数2 + ・・・, family=binomial)
> res = glm(bLifeSpan ~ BodyWt + TotalSleep, family=binomial)
警告メッセージ:
glm.fit: 数値的に 0 か 1 である確率が生じました
> summary(res)

Call:
glm(formula = bLifeSpan ~ BodyWt + TotalSleep, family = binomial)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8522 -0.4121 -0.2703 -0.1298 2.1760

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.355446 1.311101 -0.271 0.7863
BodyWt 0.019789 0.007692 2.573 0.0101 *
TotalSleep -0.229420 0.140046 -1.638 0.1014
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 54.133 on 52 degrees of freedom
Residual deviance: 28.736 on 50 degrees of freedom
(8 observations deleted due to missingness)
AIC: 34.736

Number of Fisher Scoring iterations: 7

＜説明変数を標準化して分析した場合＞
> zBodyWt = scale(BodyWt)
> zTotalSleep = scale(TotalSleep)

> res = glm(bLifeSpan ~ zBodyWt + zTotalSleep, family=binomial)
警告メッセージ:
glm.fit: 数値的に 0 か 1 である確率が生じました
> summary(res)

Call:
glm(formula = bLifeSpan ~ zBodyWt + zTotalSleep, family = binomial)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8522 -0.4121 -0.2703 -0.1298 2.1760

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9612 0.6886 -1.396 0.1627
zBodyWt 6.7423 2.6205 2.573 0.0101 *
zTotalSleep -1.0426 0.6364 -1.638 0.1014
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 54.133 on 52 degrees of freedom
Residual deviance: 28.736 on 50 degrees of freedom
(8 observations deleted due to missingness)
AIC: 34.736

Number of Fisher Scoring iterations: 7

※P値は変わらず，標準化偏回帰係数のみ値が変化していることが分かる．

＜データフレームに対して分析＞
単回帰分析と重回帰分析と同様に，data=データフレーム名のオプションを追記するだけで良い．
1. 1. ロジスティック回帰分析における特徴選択
    ロジスティック回帰分析を機械学習アルゴリズムのようにTrue of Falseの識別モデルとして扱いたい場合，最も識別性能の高い特徴量（説明変数）の組み合わせを知りたくなる．
    step()関数を使うと，すべての特徴量を使った分析（フルモデル）から1個ずつ変数を外すことを試み，外すと一番AICが良く（小さく）なるものを実際に外す。外すことによってAICが改善しなくなるまで続ける。AICとは，赤池情報量規準のことで，この値が小さいほど，良い識別モデルであると考えられる．
    
    本説明では，哺乳類の寿命を目的変数として，行動上の特徴を説明変数とした，以下のファイルを基に説明する．
    mammals2.csv [ダウンロード]
    Source：http://www.statsci.org/data/general/sleep.txt を修正したもの
    これは，単回帰分析・重回帰分析で用いた，mammals.csv から，BodyWt, BrainWt, TotalSleep, Gestation, LifeSpanの列を取り出し，NAの存在する行を取り除いたものである．
    
    ＜データフレーム読み込み＞
    > DF = read.csv(“mammals2.csv”)
    
    ＜すべての特徴で分析＞
    > res = glm(LifeSpan ~ BodyWt + BrainWt + TotalSleep + Gestation, data = DF)
    > summary(res)
    
    Call:
    glm(formula = LifeSpan ~ BodyWt + BrainWt + TotalSleep + Gestation,
    data = DF)
    
    Deviance Residuals:
    Min 1Q Median 3Q Max
    -18.677 -5.860 -2.544 3.434 39.023
    
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 9.098397 6.681148 1.362 0.18004
    BodyWt -0.079823 0.014514 -5.500 1.71e-06 ***
    BrainWt 0.050294 0.008246 6.099 2.23e-07 ***
    TotalSleep -0.123847 0.442548 -0.280 0.78087
    Gestation 0.068792 0.022216 3.096 0.00337 **
    —
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    (Dispersion parameter for gaussian family taken to be 120.5849)
    
    Null deviance: 17769.2 on 49 degrees of freedom
    Residual deviance: 5426.3 on 45 degrees of freedom
    AIC: 388.24
    
    Number of Fisher Scoring iterations: 2
    
    この時，AICは388.24である．
    
    step()関数で，最適な特徴を削除してもAICが向上しないところを探す．
    > res2 = step(res)
    Start: AIC=388.24
    LifeSpan ~ BodyWt + BrainWt + TotalSleep + Gestation
    
    Df Deviance AIC
    – TotalSleep 1 5435.8 386.33
    5426.3 388.24
    – Gestation 1 6582.5 395.90
    – BodyWt 1 9073.7 411.95
    – BrainWt 1 9911.9 416.37
    
    Step: AIC=386.33
    LifeSpan ~ BodyWt + BrainWt + Gestation
    
    Df Deviance AIC
    5435.8 386.33
    – Gestation 1 7377.4 399.60
    – BodyWt 1 9075.0 409.96
    – BrainWt 1 9935.5 414.49
    
    この結果から，TotalSleepを削除すれば，AICが386.33まで減少することが分かる．
    以下のように，実際にTotalSleepを削除して，もう一度ロジスティック回帰分析を行ってみる．
    
    > res3 = glm(LifeSpan ~ BodyWt + BrainWt + Gestation, data = DF)
    > summary(res3)
    
    Call:
    glm(formula = LifeSpan ~ BodyWt + BrainWt + Gestation, data = DF)
    
    Deviance Residuals:
    Min 1Q Median 3Q Max
    -18.423 -6.090 -2.164 3.397 39.596
    
    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 7.365246 2.481324 2.968 0.004742 **
    BodyWt -0.079693 0.014360 -5.549 1.36e-06 ***
    BrainWt 0.050022 0.008106 6.171 1.61e-07 ***
    Gestation 0.072418 0.017865 4.053 0.000193 ***
    —
    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    (Dispersion parameter for gaussian family taken to be 118.1688)
    
    Null deviance: 17769.2 on 49 degrees of freedom
    Residual deviance: 5435.8 on 46 degrees of freedom
    AIC: 386.33
    
    Number of Fisher Scoring iterations: 2
    
    確かに，AICが386.33まで減少し，すべての説明変数において，0.1%の有意水準で有意差が確認できる．
    
    便利な関数ではあるが，科学的な分析を行うときは（目的変数に対して，何の説明変数が支配的かを確かめるときには），step()関数ばかりに頼ることなく，多重共線性にも注意して，分析前に特徴量を絞り込んだ方が良い．
    また，検証用データがあり，精度・再現率・F値などの評価指標がある場合は，すべての組み合わせに対して（グリッドサーチ），性能を評価して変数選択を行う方が良い．
    
    参考情報：https://oku.edu.mie-u.ac.jp/~okumura/stat/140921.html
  2. 多重共線性の確認
    重回帰分析やロジスティック回帰分析を行うときは，説明変数間で高い相関のある組み合わせがないかどうか，事前に確認しておくこと．以下のようにすべての説明変数間の組み合わせにおけるピアソンの積率相関係数を求め，相関係数の高いものがないかどうかを確認する．高い相関係数のある組み合わせがある場合は，その説明変数の意味を確認し，どちらかを削除すること．
    
    > Atr = matrix(numeric(4*61),nrow=61,ncol=4)
    > Atr[,1] = DF$BodyWt
    > Atr[,2] = DF$BrainWt
    > Atr[,3] = DF$TotalSleep
    > Atr[,4] = DF$Gestation
    > source(“calcCorAll.r”)
    > calcCorAll(Atr)
    [,1] [,2] [,3] [,4]
    [1,] 1 0.955218 -0.3102933 0.6800555
    [2,] 0 1.000000 -0.3118036 0.6884701
    [3,] 0 0.000000 1.0000000 -0.6194160
    [4,] 0 0.000000 0.0000000 1.0000000
    
    calcCorAll.r [ダウンロード]
    ※上記ファイルのダウンロードと利用は，土方ゼミの学生に限ります．
  3. 離散値のヒストグラムの書き方
    bに全データ（例えば，-3～3）のデータが入っているとして，
    > b
    [1] 3 0 0 2 0 1 -1 -1　・・・
    c = table(b)
    barplot(c, main=”Overall trust”, xlab=”Answer”, ylab=”Frequency”, ylim=c(0,200), axis.lty=1)
    box(lty=1)
    
    参考情報：https://data-science.gr.jp/implementation/ida_r_barplot.html
</li=””>

1. ステップワイズ法と階層的重回帰分析
  （参考資料：階層的重回帰分析）
  https://cogpsy.educ.kyoto-u.ac.jp/personal/Kusumi/datasem17/lee1.pdf
2. 刺激提示前後（実験前後）の測定値の変化量の検定
  2種類の刺激（実験条件）を提示したときに，刺激提示前に測定した値と，刺激提示後に測定した値の変化量が，刺激の種類により違いがあるかどうかの検定方法には，以下の2つがある．
  1. 実験前後の値の差，すなわち「実験後の値－実験前の値」を求めておき，それらのt検定を実施する．
  2. しっかりランダム化が行われている前提で，実験後の値についてt検定を実施する．
  どちらの方法が良いかは，実験に依存する．
  たとえば，買ってきた苗を2種類の肥料で育て，買ってきたときの苗の高さと，3か月後の植物の高さから，成長に差があるかどうかを見るときには，「3か月後の植物の高さ－買ってきたときの苗の高さ」を計算し，t検定するのは問題ないと思われる．肥料の違いが成長の差になって出てくる可能性が高い．
  しかし，ある商品の好き嫌い（-3～3のカテゴリ値）を，広告を見る前と見た後で違いがあるかどうかを見るときには，「広告見た後の嗜好－広告見る前の嗜好」を計算し，t検定するのは問題ないと思うが，経験上有意差が出にくい．理由はいくつかあるが，広告を見る前の嗜好は，個人差が大きい（分散が大きい）ことが想定されるため．
  個人的には，分析手法や高度なアルゴリズムというテクニカルなことに頼る前に，しっかり実験計画を立てて，適切にランダム化を行っていることの方が重要だと思う．
  
  あと，うまくいくケースは限られるが，実験前と後の比較は，対応のあるデータになるので，各刺激（実験条件）において，実験後の値と実験前の値に差があるかどうかを，対応のあるt検定を行って有意差を求める．もし，一方の刺激（実験条件）で有意差があり，もう一方の刺激（実験条件）で有意差がなければ，刺激（実験条件）により変化量に違いがあることを検証することができる．
  
  実験前後の値の検定については以下を参照
  http://www.f.kpu-m.ac.jp/c/kouza/joho/ijou/taiouari.html
3. 共変量
  統計やデータ分析において、ある結果に影響を与える可能性がある変数のことを指す。主に、他の要因（説明変数や独立変数）とは区別されるが、分析対象に影響を与える潜在的な変数を指す。共変量をモデルに取り入れることで、分析結果の精度を上げたり、バイアスを取り除くことができる。例えば，査読者に年齢の影響って大きいんじゃないの？みたいな質問があったときにも，その影響の有無や，それが他の変数に与える影響を取り除いた効果を示すことができる．

1. 単純傾斜分析
  調整変数が、その平均よりも１標準偏差高い場合と，１標準偏差低い場合に，独立変数と従属変数の関係（すなわち傾き）が，統計的に有意に異なるかどうかを検定する方法．
  
  重回帰式として、y = a + b1x + b2w + b3xw （yは従属変数、aは切片、xは独立変数、wは調整変数、b1～b3が偏回帰係数）を考えたとき，この式をxについて変形すると，y = (a + b2w) + (b1 + b3w) x となる。つまり、yとxの関係について、切片が(a + b2w)で、傾きが(b1 + b3w)であると解釈することができる。上記の式でいうと、分析の結果b3が統計的に有意である場合、すなわちwの係数がゼロではないので交互作用がある（調整効果が存在する）ことが示唆されるので、単純傾斜分析で求めたいのは、wが特定の値（例えば平均＋１標準偏差）のときに、傾き(b1 + b3w)がゼロではないか、すなわち統計的に有意かどうかを検証することである。
  
  Rでは，以下のように実行する
  
  sampledata
  
  のサンプルデータを例にする．
  
  dat = read.csv(“sampledata.csv”)
  library(pequod)
  model1 = lmres(sat~talk*per, centered =c(“talk”, “per”), data =dat)
  summary(mode1)
  
  sat ~ talk*perというモデル式は，talkとperの主効果と，talkとperの交互作用を投入する，という意味です。つまり，sat~talk+per+talk:perと同じ意味です。
  次に，centered=では，中心化する（平均を0にする）変数を指定します。
  以下が実行結果
  
  ===
  Models
  R R^2 Adj. R^2 F df1 df2 p.value
  Model 0.467 0.218 0.211 27.576 3.000 296 9.4e-16 ***
  —
  Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  
  Residuals
  Min. 1st Qu. Median Mean 3rd Qu. Max.
  -2.3160 -0.4539 -0.0764 0.0000 0.6840 2.2810
  
  Coefficients
  Estimate StdErr t.value beta p.value
  (Intercept) 3.41647 0.05116 66.77491 <2e-16 *** talk 0.26572 0.05192 5.11783 0.2649 <2e-16 *** per 0.14054 0.02941 4.77880 0.2482 <2e-16 *** talk.XX.per 0.13022 0.03035 4.29012 0.2234 2e-05 *** — Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Collinearity VIF Tolerance talk 1.0148 0.9854 per 1.0216 0.9789 talk.XX.per 1.0270 0.9737 === talk xx perの項が有意であるので，交互作用項が有意であることがわかった．ここで単純主効果を検討する。単純主効果は，以下のコードで検討する。 model2 <- simpleSlope(model1, pred=”talk”, mod1 = “per”) summary(model2) simpleSlopeの二つ目のＳが大文字なのに注意。predは予測変数，つまり独立変数を指定します。つまり，talkがそれになります。mod1は，調整変数の一つ目，ということで，perを指定します。結果は以下のようになります。 === 　** Estimated points of sat ** Low talk (-1 SD) High talk (+1 SD) Low per (-1 SD) 3.1329 3.2064 High per (+1 SD) 3.1731 4.1534 ** Simple Slopes analysis ( df= 296 ) ** simple slope standard error t-value p.value Low per (-1 SD) 0.0370 0.0779 0.48 0.63 High per (+1 SD) 0.4944 0.0708 6.99 <2e-16 *** — Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ** Bauer & Curran 95% CI ** lower CI upper CI per -4.1656 -1.078 === perが高群（+1SD)の場合に，talkの効果は有意ですが，低群(-1SD)の場合，非有意であることがわかります。次に，単純効果をグラフ表示したい場合には，PlotSlope関数を使います。 PlotSlope(model2) 参考ページ：https://norimune.net/1856
1. 1. ロジスティック回帰分析のR^2値
    疑似R^2や，Cox and Snell R^2，Nagelkerke R^2などがある．