advanced_univariate.Rに対応したR Notebook

第7章の内容

データの準備

まずは第6章までと同じデータを読み込む

###For chapter 07
phyto_metadata <- readRDS("phyto_metadata.obj")
species_ryuko_data <- readRDS("phyto_ryuko_data.obj")
metadata_ecoplate <- readRDS("metadata_ecopl.obj")
summary_ecoplate <- readRDS("summary_ecopl.obj")
species_richness <- apply(species_ryuko_data > 0, 1, sum) #種数の計算
total_abundance <- apply(species_ryuko_data, 1, sum) #総個体数の計算

7.1 重回帰分析

7.1.1 説明変数が二つ以上ある場合の可視化方法

復習として二つの説明変数それぞれに対する散布図を描いてみる

#種数と温度の関係
plot(species_richness ~ phyto_metadata$temp, cex = 2) 

#種数と総個体数の関係
plot(species_richness ~ total_abundance, cex = 2)

3次元散布図について描いてみる(scatterplot3dというパッケージのインストールを済ませましょう)

library(scatterplot3d)
scatterplot3d(
  x = phyto_metadata$temp, y = total_abundance, z = species_richness,　#軸に使うベクトルの設定
  xlab = "temperature", ylab = "abundance", zlab = "richness",　#軸ラベルの設定
  angle = 30　#角度の設定
)

マウスで操作が可能なインタラクティブな3次元散布図も描いてみましょう(rglというパッケージのインストールが必要ですが，うまくいかない時は諦めて飛ばしてよいです)

library(rgl)
plot3d(
  x = phyto_metadata$temp, y = total_abundance, z = species_richness,
  xlab = "temperature", ylab = "abundance", zlab = "richness",　#軸ラベルの設定
  type = "s",　#マークのタイプ
  size = 3,　　#マークのサイズ
  col = c(2,9)[as.factor(phyto_metadata$year)]　#年ごとに色指定を変える
)
rgl.snapshot("fig7.3b.tiff")

Rスクリプトを実行したときのグラフは新しいウィンドウに表示されます．そのイメージは以下の通り．

図7.3b

7.1.2 重回帰分析のイメージ

重回帰分析に行く前に「単」回帰分析をしてみる

#温度及び総個体数による単回帰分析
model_temp <- lm(species_richness ~ phyto_metadata$temp)
model_abundance <- lm(species_richness ~ total_abundance)
#温度による回帰直線の描写
plot(species_richness ~ phyto_metadata$temp, cex = 2)
abline(model_temp) #これが回帰分析の結果から回帰直線を描くコード

#総個体数による回帰直線の描写
plot(species_richness ~ total_abundance, cex = 2)
abline(model_abundance) #これが回帰分析の結果から回帰直線を描くコード

重回帰分析をイメージしてもらうために「回帰平面」を描いてみる

#3次元散布図のコード
plot3d(
  x = phyto_metadata$temp, y = total_abundance, z = species_richness,
  xlab = "temperature", ylab = "abundance", zlab = "richness",
  type = "s",
  size = 3,
  col = c(2,9)[as.factor(phyto_metadata$year)] #年ごとに色指定を変える
)
#重回帰分析のコード
fit_model <- lm(species_richness ~ phyto_metadata$temp + total_abundance)
#重回帰分析の結果を格納したオブジェクト(fit_model)から回帰係数を抽出する
plane_coeff <- coef(fit_model)
#回帰係数を使って「回帰平面」を上描きする
planes3d(plane_coeff[2], plane_coeff[3], -1, plane_coeff[1], col = "blue", alpha = 0.5)
#rgl.snapshot("fig7.5.tiff")

Rスクリプトを実行したときのグラフは新しいウィンドウに表示されます．そのイメージは以下の通り．

図7.5 3次元散布図と回帰平面

7.1.3 重回帰分析のコーディング

#lm関数を使って重回帰分析を実行する
multi_model <- lm(species_richness ~ phyto_metadata$temp + total_abundance)
#summary関数をつかって重回帰分析の結果を表示する
summary(multi_model)

Call:
lm(formula = species_richness ~ phyto_metadata$temp + total_abundance)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.2519 -3.6671 -0.5555  1.4574 13.0203 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)
(Intercept)         8.733183   7.271159   1.201    0.253
phyto_metadata$temp 0.974593   0.548926   1.775    0.101
total_abundance     0.005174   0.011083   0.467    0.649

Residual standard error: 5.93 on 12 degrees of freedom
Multiple R-squared:  0.5009,    Adjusted R-squared:  0.4178 
F-statistic: 6.023 on 2 and 12 DF,  p-value: 0.01545

7.2 交互作用

7.2.1 交互作用ありの分散分析

交互作用を直感的に理解するためにまずは可視化してみる

参照元: https://momonoki2017.blogspot.com/2018/08/r019.html

data("ToothGrowth") #Rに標準で付属するサンプルデータを読み込み
sub_TG <- subset(ToothGrowth, dose > 0.5)　#摂取量が0.5より大きいデータだけをsubset関数で抜き出す

#摂取量ｘ摂取タイプで箱ひげ図を描写する
boxplot(
  len ~ dose + supp, data = sub_TG, #公式
  outline = F,　#外れ値の扱い
  col = c("white", "white", "gray", "gray")　#列ごとの色の指定
)
#生データを箱ひげ図に重ね合わせる
stripchart(
  len ~ dose + supp, data = sub_TG, #公式
  method = "stack",
  pch = c(1,1,2,2),　#記号の指定:１と２列目(1.OJと2.OJ)に1を指定，３と４列目(1.VCと2.VC)に2を指定
  cex = 3, #マークの大きさを指定
  vertical = TRUE,
  add = TRUE
)

交互作用に特化した可視化方法もある

interaction.plot(
  x.factor = sub_TG$dose,　#主要因の指定
  response = sub_TG$len,   #応答変数(ｙ軸)の指定
  trace.factor = sub_TG$supp,　#副要因の指定
  xlab = "dose level", ylab = "growth"　#軸ラベルの指定
)

「二要因の分散分析と交互作用」についてやってみる

まずは、交互作用ナシの二要因分散分析(「応答変数 ~ 説明変数1 + 説明変数２」の形で公式を書く)

model_tooth01 <- lm(len ~ dose + supp, data = sub_TG)
anova(model_tooth01) #分散分析の結果の表示にはanova関数を使う

Analysis of Variance Table

Response: len
          Df Sum Sq Mean Sq F value   Pr(>F)    
dose       1 405.13  405.13 26.9841 7.72e-06 ***
supp       1  85.56   85.56  5.6985   0.0222 *  
Residuals 37 555.51   15.01                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

次に、交互作用アリの二要因分散分析を、「応答変数 ~ 説明変数1 + 説明変数２ + 説明変数１: 説明変数２」の形で公式を書いて実現する

model_tooth02 <- lm(len ~ dose + supp + dose:supp, data = sub_TG)
anova(model_tooth02)

Analysis of Variance Table

Response: len
          Df Sum Sq Mean Sq F value    Pr(>F)    
dose       1 405.13  405.13 31.3510 2.387e-06 ***
supp       1  85.56   85.56  6.6207   0.01435 *  
dose:supp  1  90.30   90.30  6.9878   0.01208 *  
Residuals 36 465.21   12.92                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

7.3 一般線形モデルとモデル選択

7.3.3 一般線形モデルとしてのRのlm関数の使いかた

分散分析タイプのデータを一般線形モデルで扱う

(復習)anova関数を使うと分散分析の表として結果が出力される

anova(lm(species_richness ~ as.factor(phyto_metadata$month)))

Analysis of Variance Table

Response: species_richness
                                Df Sum Sq Mean Sq F value  Pr(>F)  
as.factor(phyto_metadata$month)  3 418.83 139.611  3.5985 0.04961 *
Residuals                       11 426.77  38.797                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

一般線形モデル(General Linear　Model)の形式で出力するには、summary関数を使えばよい

####General Linear Model01####
summary(lm(species_richness ~ phyto_metadata$month))

Call:
lm(formula = species_richness ~ phyto_metadata$month)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.200  -2.167  -0.800   1.650  11.200 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)   
(Intercept)               16.500      4.404   3.746  0.00323 **
phyto_metadata$monthM04    9.300      5.211   1.785  0.10190   
phyto_metadata$monthM05   11.167      5.686   1.964  0.07531 . 
phyto_metadata$monthM10   16.700      5.211   3.205  0.00839 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.229 on 11 degrees of freedom
Multiple R-squared:  0.4953,    Adjusted R-squared:  0.3577 
F-statistic: 3.599 on 3 and 11 DF,  p-value: 0.04961

カテゴリー変数と量的変数が混じったデータを一般線形モデルで扱う

ToothGrowthデータをすべて使って、一般線形モデルで解析してみよう．

####General Linear Model02####
data("ToothGrowth")　#データの読み込み

まずは、変数のタイプ(クラス)を確認する

class(ToothGrowth$supp)

[1] "factor"

class(ToothGrowth$dose)

[1] "numeric"

interaction.plot関数を使って交互作用を可視化してみよう

interaction.plot(
  x.factor = ToothGrowth$dose,　#主要因を摂取量にする
  response = ToothGrowth$len,　 #応答変数は成長量
  trace.factor = ToothGrowth$supp,　#副要因は摂取タイプ
  xlab = "dose level", ylab = "growth"　#軸ラベルの指定
)

では、実際に一般線形モデルとして解析・結果の出力をしてみましょう

model_tooth03 <- lm(len ~ dose + supp + dose:supp, data = ToothGrowth) #解析の結果は一旦、オブジェクトmodel_tooth03に格納する
summary(model_tooth03) #結果の出力

Call:
lm(formula = len ~ dose + supp + dose:supp, data = ToothGrowth)

Residuals:
    Min      1Q  Median      3Q     Max 
-8.2264 -2.8462  0.0504  2.2893  7.9386 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   11.550      1.581   7.304 1.09e-09 ***
dose           7.811      1.195   6.534 2.03e-08 ***
suppVC        -8.255      2.236  -3.691 0.000507 ***
dose:suppVC    3.904      1.691   2.309 0.024631 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.083 on 56 degrees of freedom
Multiple R-squared:  0.7296,    Adjusted R-squared:  0.7151 
F-statistic: 50.36 on 3 and 56 DF,  p-value: 6.521e-16

二群比較用のデータを一般線形モデルで扱う

####General Linear Model03####
summary(lm(species_richness ~ as.factor(phyto_metadata$year))) #一般線形モデルの出力を直接summary関数の入力にしている

Call:
lm(formula = species_richness ~ as.factor(phyto_metadata$year))

Residuals:
    Min      1Q  Median      3Q     Max 
-12.429  -4.964  -0.500   4.536  14.571 

Coefficients:
                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                          28.429      3.023   9.403 3.66e-07 ***
as.factor(phyto_metadata$year)2019   -1.929      4.140  -0.466    0.649    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.999 on 13 degrees of freedom
Multiple R-squared:  0.01642,   Adjusted R-squared:  -0.05924 
F-statistic: 0.217 on 1 and 13 DF,  p-value: 0.649

7.3.4 モデル選択

ますは可視化してみる

年の間の種数の違い、温度による種数の違い、総個体数による種数の違いをそれぞれ可視化する

#2018 vs 2019
#箱ひげ図
boxplot(
  species_richness ~ phyto_metadata$year, #モデル式
  outline = F,　#外れ値の扱い
  col = c("white", "gray")　#色の指定
)
#散布図
stripchart(
  species_richness ~ phyto_metadata$year,　#モデル式
  method = "stack",　#同一値の描画設定
  pch = c(1,2),　#描画マークの種類指定
  cex = 3,　　　#描画マークの大きさ指定
  vertical = TRUE,　#横向きか縦向きかの指定
  add = TRUE　#箱ひげ図に重ね合わせるかどうかの指定
)

#回帰直線を描くための準備としての単回帰分析
model_temp <- lm(species_richness ~ phyto_metadata$temp)
model_abundance <- lm(species_richness ~ total_abundance)
#単回帰分析の結果を使って、散布図と回帰直線を重ね合わせる
plot(species_richness ~ phyto_metadata$temp, cex = 2) #散布図
abline(model_temp) #回帰直線

#単回帰分析の結果を使って、散布図と回帰直線を重ね合わせる
plot(species_richness ~ total_abundance, cex = 2) #散布図
abline(model_abundance) #回帰直線

step関数を使ったモデル選択

#ステップ関数でのモデル選択の結果をオブジェクトselected_modelに格納(選択している途上の結果もコンソール上に出力される)
selected_model <- step(lm(species_richness ~ phyto_metadata$temp + total_abundance + phyto_metadata$year))

Start:  AIC=52.99
species_richness ~ phyto_metadata$temp + total_abundance + phyto_metadata$year

                      Df Sum of Sq    RSS    AIC
- total_abundance      1    28.372 329.45 52.340
<none>                             301.08 52.990
- phyto_metadata$year  1   120.932 422.01 56.055
- phyto_metadata$temp  1   168.925 470.00 57.670

Step:  AIC=52.34
species_richness ~ phyto_metadata$temp + phyto_metadata$year

                      Df Sum of Sq    RSS    AIC
<none>                             329.45 52.340
- phyto_metadata$year  1    100.22 429.67 54.325
- phyto_metadata$temp  1    502.27 831.71 64.232

選択されたベストモデルのまとめを出力する

→温度と年が入ったモデルがベストであることがわかる

summary(selected_model)

Call:
lm(formula = species_richness ~ phyto_metadata$temp + phyto_metadata$year)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2495 -3.4168 -0.1866  1.8900 11.7505 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)   
(Intercept)         -1.286e+04  6.735e+03  -1.910  0.08037 . 
phyto_metadata$temp  1.588e+00  3.712e-01   4.277  0.00107 **
phyto_metadata$year  6.371e+00  3.334e+00   1.911  0.08023 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.24 on 12 degrees of freedom
Multiple R-squared:  0.6104,    Adjusted R-squared:  0.5455 
F-statistic:   9.4 on 2 and 12 DF,  p-value: 0.003497