★手計算とRで学ぶ統計学: 線形モデリング (回帰分析)

まずは，実験や調査の前にどれだけのデータを集めればよいのかを把握しておきましょう。

以下のコードをR上で実行すると

# (初めてパッケージを使用するときのみ行う) パッケージをインストールする  
install.packages("pwr")  
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む 
library(pwr) 
# 単回帰モデル
# サンプルサイズ設計 (厳密には残差の自由度の推定)
pwr_est_simple <- pwr.f2.test(u = 1, # パラメータ数
           v = NULL, # 残差の自由度。サンプルサイズによって決定される値なので設定しない
           f2 = 0.35, # 大きな効果量
           sig.level = 0.05, # 有意水準
           power = 0.8 # 検定力
　　　　　　)
# 残差の自由度の推定結果の表示
pwr_est_simple

以下のような結果が出力されます。

     Multiple regression power calculation 

              u = 1
              v = 22.50313
             f2 = 0.35
      sig.level = 0.05
          power = 0.8

計算結果は無事に出力されましたが，このコードだけでサンプルサイズを計算することはできません (サンプルサイズ設計に必要な値が計算されるだけ)。そのため，以下のコードを実行して必要なサンプルサイズを求めます。

# サンプルサイズ設計のための下準備
v_simple <- pwr_est_simple$v # 変数pwr_est_simpleのv列の値を変数v_simpleに代入
u_simple <- pwr_est_simple$u # 変数pwr_est_simpleのu列の値を変数u_simpleに代入

# サンプルサイズの決定
n_simple <- v_simple + u_simple + 1
# サンプルサイズを表示
n_simple

上記のコードを実行することで，コンソール上に[1]　24.50313という値が表示されます。 そのため，今回の分析では最低でも25個のデータがあればよいとわかりました。

それでは，実際にデータをRに読み込み，データの確認を行ってみましょう。

まず，以下のコードでデータを読み込み，列名でデータを指定できるようにします。

# 使用するデータの設定。作業ディレクトリとして指定したフォルダに入れておかないとエラーを吐きます
data_simple_reg <- read.csv("simple_regression.csv")
# 列名でデータを指定できるようにする
attach(data_simple_reg)

データを読み込んだら，必ずデータの中身を確認しましょう。

# データ内容の確認
head(data_simple_reg)

データの頭から6行を確認すると，今回行う解析に用いるデータを正しく読み込めたことがわかります。

> head(data_simple_reg)
   ID speed dist 
1  1     4    2 
2  2     4   10 
3  3     7    4 
4  4     7   22 
5  5     8   16 
6  6     9   10

データを読み込んだあとは， ヒストグラムを表示して分布の形状を確認します。

# ヒストグラムの表示
hist(data_simple_reg$speed)
hist(data_simple_reg$dist)

ヒストグラムを見る限りでは，いずれの変数にも2つ以上のピークがないことがわかります。

実際にRで制動距離を目的変数，停止直前の速度を説明変数とした単回帰モデルを作成し，2つの変数の関係を検討しましょう。

以下のコードをR上で実行すると

# 単回帰モデルの作成
# (目的変数 ~ 説明変数, データ) と指定
simple_reg_model <- lm(dist ~ speed, data = data_simple_reg)

# 回帰モデルを要約した結果の表示　 
summary(simple_reg_model)

たちまちに結果が得られます。

Call: 
lm(formula = dist ~ speed, data = data_simple_reg) 
Residuals: 
    Min      1Q  Median      3Q     Max 
-19.333  -9.177  -1.153   3.112  43.956 

Coefficients: 
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -10.0031    10.2663  -0.974  0.34001   
speed         3.2891     0.8932   3.683  0.00123 ** 
--- 
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 
# Estimate: 推定値 (傾き，切片の値)
# Std.Error: 傾きや切片の標準誤差

Residual standard error: 13.66 on 23 degrees of freedom # 残差の標準誤差
Multiple R-squared:  0.3709 # R2乗値,    Adjusted R-squared:  0.3436  # 修正済みR2乗値
F-statistic: 13.56 on 1 and 23 DF,  p-value: 0.001233

ブレーキをかける直前の速度と制動距離の関係を表す回帰モデル式は\(y = - 10.0 +3.29x\) と求められました。
また，モデル全体の調整済み決定係数は\(\mathrm{adjusted}\; R^2= 0.34 \), モデル全体に意味があるかどうかを検討する分散分析の結果は\(F(1, 23) = 13.56,\;p < .001\) とわかります。
以上の結果から，5%有意水準のもとでモデル全体が有意であることがわかりました。

単回帰分析の効果量\(f^2\)は，決定係数\(R^2\)から求めることができます \( \left( f^2 = \frac{R^2}{1-R^2} \right) \)。
ここでは，決定係数\(R^2\)を単回帰モデルの要約結果から抜き出して効果量を求めます。

# 効果量の計算 
# R2乗値を変数r_squaredに代入
r_squared_simple <- summary(simple_reg_model)["r.squared"]$r.squared 
# 効果量の計算
simple_reg_effect <- r_squared_simple / (1-r_squared_simple)
# 効果量の表示
simple_reg_effect

上記のコードを実行すると，以下のように計算結果が表示されます。

> simple_reg_effect
[1] 0.5896113

サンプルサイズ設計の段階では \(f^2 = 0.35\) (大きな効果を検出できることを想定) と設定しましたが，実際にはそれを上回る効果があったことがわかります。

どのような仮説のもと単回帰分析を行ったか，その結果はどうだったか，を論文やレポートに書く際は以下のようにします。

(※) この書き方は心理学をはじめとする社会科学分野で一般的なものです。実際にレポートや論文を書く際には，担当教員や指導教員の指示，それぞれの分野の慣例に従ってください。  

まずは，実験や調査の前にどれだけのデータを集めればよいのかを把握しておきましょう。

以下のコードをR上で実行すると

# サンプルサイズ設計 
# (初めてパッケージを使用するときのみ行う) パッケージをインストールする  
install.packages("pwr")  
# (RやRStudioを新しく起動したら毎回行う) パッケージを読み込む  
library(pwr) 

# サンプルサイズ設計 (厳密には残差の自由度の推定) 
pwr_est_multi <- pwr.f2.test(u = 3, # パラメータ数 
                       v = NULL, # 残差の自由度。サンプルサイズによって決定される値なので設定しない 
                       f2 = 0.15, # 中程度の効果量 
                       sig.level = 0.05, # 有意水準 
                       power = 0.8 ) # 検定力 
# 残差の自由度の推定結果の表示
pwr_est_multi

以下のような結果が出力されます。

> pwr_est_multi    
 Multiple regression power calculation
               u = 3
               v = 72.70583
              f2 = 0.15
       sig.level = 0.05
           power = 0.8

単回帰分析と同様，現段階ではサンプルサイズ設計に必要な値だけがわかっている状態です。そのため，以下のコードを実行して必要なサンプルサイズを求めます。

# サンプルサイズ設計のための下準備
v_multi <- pwr_est_multi$v # 変数pwr_est_multiのv列の値を変数v_multiに代入
u_multi <- pwr_est_multi$u # 変数pwr_est_multiのu列の値を変数u_multiに代入

# サンプルサイズの決定
n_multi <- v_multi + u_multi + 1
# サンプルサイズを表示
n_multi

上記のコードを実行することで，コンソール上に[1]　76.70583という値が表示されます。 そのため，今回の分析では最低でも77個のデータがあればよいとわかりました。
通常，企業等の市場調査ではキリの良い数字が使われるため，今回は100人のデータをモデル作成に使用します。

それでは，実際にデータをRに読み込み，データの確認を行ってみましょう。

まず，作業ディレクトリを設定します。次に，以下のコードでデータを読み込み，列名でデータを指定できるようにします。

# 作業ディレクトリの設定
setwd("csvファイルが入っているフォルダ")
# 使用するデータの設定
data_multiple_reg <- read.csv("multiple_regression.csv")
# 列名でデータを指定できるようにする
attach(data_multiple_reg)

データを読み込んだら，必ずデータの中身を確認しましょう。

# データ内容の確認
head(data_multiple_reg)

データの頭から6行を確認すると，今回行う解析に用いるデータを正しく読み込めたことがわかります。

> head(data_multiple_reg)
   ID sex age tend value
 1  1   m  55    2   671
 2  2   f  72    7 21741
 3  3   m  30    2  4986
 4  4   f  26    4  9417
 5  5   f  35    6 18438
 6  6   m  45    4 15317

実際にRで一人当たりのビールの年間消費額を目的変数，年齢，性別，週あたりの飲酒日数を説明変数とした重回帰モデルを作成し，変数の関係を検討しましょう。

以下のコードをR上で実行すると

# 重回帰モデルの作成 
multiple_reg_model <- lm(value ~ sex + age + tend, data = data_multiple_reg)　# (目的変数 ~ 説明変数, データ) と指定
# 計算結果の要約の表示 
summary(multiple_reg_model)

たちまちに結果が得られます。

> summary(multiple_reg_model) 
Call: lm(formula = value ~ sex + age + tend, data = data_multiple_reg) 
Residuals: 
     Min       1Q   Median       3Q      Max 
-11508.0  -2944.2   -263.7   2860.1  10295.6 
Coefficients:
             Estimate Std. Error t value Pr(>|t|)     
(Intercept)  2518.62    1842.54   1.367  0.17484     
sexm        -3131.32    1044.70  -2.997  0.00347 ** 
age            21.89      30.58   0.716  0.47594     
tend         1790.59     227.69   7.864  5.5e-12 *** 
--- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 5143 on 96 degrees of freedom 
Multiple R-squared:  0.4324,    Adjusted R-squared:  0.4147 
F-statistic: 24.38 on 3 and 96 DF,  p-value: 8.21e-12

一人当たりのビールの年間消費額と年齢，性別，週あたりの飲酒日数の関係を表す回帰モデル式は\(y = 2518.62 -3131.32x_{sexm}+21.89x_{age} + 1790.59x_{tend}\) と求められました。
しかし，切片 (Intercept) とageはそれぞれの係数がモデルに与える影響力についての\(t\)検定において5%有意水準のもとで有意でない\(p > .05 \) ため，これらの変数が回帰モデルに影響を与えているかどうかについての判断は保留します。
また，モデル全体の調整済み決定係数は\(\mathrm{adjusted}\; R^2=0.41 \), モデル全体に意味があるかどうかを検討する分散分析の結果は\(F(3, 96) = 24.38,\;p < .001\) とわかります。
以上の結果から，5%有意水準のもとでモデル全体が有意であることがわかりました。

今回は生データをそのまま使って重回帰分析を行いました。しかし，単位が異なる場合は偏回帰係数を直接比較してモデルへの影響を検討することが難しいです。
そのため，生データを標準化 (データを平均0，分散1に変換する) して回帰分析を行うことで，それぞれの係数の単位がなくなり値どうしを直接比較することができるようになります。

Rを用いてデータの標準化と重回帰分析を行うには以下のコードを実行します (データの読み込み等は省略)。

# 標準化するためのデータオブジェクトを作成 
data_multiple_reg_scaled <- data_multiple_reg 
# 年齢，頻度，消費額について標準化 
data_multiple_reg_scaled$age <- scale(data_multiple_reg_scaled$age) 
data_multiple_reg_scaled$tend <- scale(data_multiple_reg_scaled$tend) 
data_multiple_reg_scaled$value <- scale(data_multiple_reg_scaled$value) 

# 重回帰モデルの作成 
multiple_reg_model_scaled <- lm(value ~ sex + age + tend, data = data_multiple_reg_scaled)　# (目的変数 ~ 説明変数, データ) と指定 
summary(multiple_reg_model_scaled)

上記のコードを実行すると，\(y = 0.265 + -0.466x_{sexm} + 0.05x_{age} + 0.60x_{tend}\) というモデル式であることがわかり，係数同士を直接比較して議論できるようになりました ( \(x_{sexm} = 1\) (つまり男性) であることはビールの年間消費額に-0.466の影響を与える，など)。
論文やレポートには，標準化しない値で重回帰分析を行った時の係数と，標準化した時の係数の両方を記載します。

単回帰分析と同様，決定係数\(R^2\)をモデルの要約結果から抜き出して効果量\(f^2\)を求めます。

# 効果量の計算 
# R2乗値を変数r_squaredに代入 
r_squared_multiple <- summary(multiple_reg_model)["r.squared"]$r.squared 
# 効果量の計算 
multiple_reg_effect <- r_squared_multiple / (1-r_squared_multiple)
# 効果量の表示
multiple_reg_effect

上記のコードを実行すると，以下のように計算結果が表示されます。

> multiple_reg_effect
[1] 0.7617872

サンプルサイズ設計の段階では \(f^2 = 0.15\)  (中程度の効果を検出できることを想定) と設定しましたが，実際にはそれを上回る効果があったことがわかります。

どのような仮説のもと重回帰分析を行ったか，その結果はどうだったか，を論文やレポートに書く際は以下のようにします。

(※) 数表を作成するためのJavaScriptライブラリの都合上表題が掲載できませんでした…。実際のレポートでは，表の上部に必ずタイトルを記載しましょう。