1. トップページ
  2. Rによる群集組成の解析
  3. 群集組成と環境変数の関係：直接傾度分析

Rによる群集組成の解析

English

群集組成と環境変数の関係：直接傾度分析

直接傾度分析では、群集組成の違いを環境変数の重回帰で説明します。間接傾度分析に続けて行いますので、同じ群集データdatと環境データedatを使います。

環境変数を整理

環境変数のセットを説明変数として重回帰を行いますので、重回帰の場合の一般的な注意が必要です。名義変数や類別変数が含まれる場合も同様です。ここでは、次の３点について検討し、必要なデータの修正を行います。

  • 正規性
  • 多重共線性
  • 標準化

説明変数は正規分布をしているのが理想的ですが、実際のデータではなかなか難しいです。重回帰分析の頑健性を信じて、分布が大きく偏った変数についてのみ、変数変換で分布の偏りを修正します。hist()関数で分布図を見てみます。

hist(edat$MdG)

環境変数MdGは底質の砂の中央粒径値(mm)です。見た目で判断していますが、正規分布とは大きく外れています。そこで、底質粒度を表す際に良く用いられるファイスケールに変換します。2を底に対数変換して正負を逆転します。粒径が1ミリだとファイスケールでゼロになります。

edat$MdPh <- -log(edat$MdG, 2)
hist(edat$MdPh)

今度は、逆に偏りましたが、変換前よりは偏りがましになったのでこれで良いことにしておきます。 他にも３つの変数で偏りが見られましたので、これらは平方根変換をしておきます。

edat$sqPM <- sqrt(edat$PM)
edat$sqTN <- sqrt(edat$TN)
edat$sqTOC <- sqrt(edat$TOC)

まだまだ正規分布にはかけ離れているのですが、次に移ります。

多重共線性はvariance inflation factor (VIF)を求めて検討します。VIFが10以上だと多重共線性の問題が生じるといわれていますので、10未満のより小さなVIFをとるように変数を順次削除していきます。まず、edatのうち変数変換をした元データは使わないので、新たな変数を作っておきます。また、２変数間の相関を確認しておきます。cor()とpairs()で十分ですが、見やすさを重視してpsychパッケージのpairs.panelsで散布図を示します。

edat2 <- edat[,c(1:7,12:15)]
psych::pairs.panels(edat2)

vif <- diag(solve(cor(edat2)))
vif
##       PSU        PH      ORP3      ORP5      ORP8        IL      Chla 
##  4.638988  4.147434 12.680515 29.958645 19.160637  6.747144 22.484479 
##      MdPh      sqPM      sqTN     sqTOC 
##  6.649999  7.595803 17.796963  9.618832

明らかにVIFは10を越えています。底質の酸化還元電位（ORP3、ORP5、ORP8）の値に問題がありそうなので、ORP5を削除してみます。

vif <- diag(solve(cor(edat2[,-4])))
vif
##       PSU        PH      ORP3      ORP8        IL      Chla      MdPh 
##  3.792806  3.433859  9.190892 12.508433  6.046272 20.499752  5.710017 
##      sqPM      sqTN     sqTOC 
##  7.401397 16.304317  9.043974

順次vifを調べ、次はChla、その次はsqTOCと削除してゆきます。

vif <- diag(solve(cor(edat2[,c(-4, -7, -11)])))
vif
##      PSU       PH     ORP3     ORP8       IL     MdPh     sqPM     sqTN 
## 2.545664 1.474089 2.594832 2.617503 2.188238 1.948444 2.347553 1.865053

これでVIFは2.6175032となりました。さらにORP8を削除しても良いかもしれませんが、ひとまず変数の削除は終了としておきます。混乱しないように新たな環境変数のセットをedat2と名付け、標準化しておきます。

edat2 <- edat2[,c(-4, -7, -11)]
edat2 <- data.frame(scale(edat2))
colMeans(edat2)
##           PSU            PH          ORP3          ORP8            IL 
##  2.016905e-16 -2.035409e-15  9.066821e-17  5.551115e-17 -5.458597e-17 
##          MdPh          sqPM          sqTN 
##  4.070818e-17  1.850372e-17 -6.476301e-17
apply(edat2, 2, sd)
##  PSU   PH ORP3 ORP8   IL MdPh sqPM sqTN 
##    1    1    1    1    1    1    1    1

各列の平均が0、標準偏差が1となりましたので、これで説明変数の準備ができました。

dbRDAによって群集組成を環境変数で重回帰

dbRDAでは群集組成の違いを環境変数で説明します。計算上は、PCoAで求めた各群集に対応する座標をうまく再現できる環境変数の係数を求めます。得られた重回帰式によって予測が可能になるのが特長です。ただし、環境変数に対し線型な関係を示すことを想定していますので、線型でない部分は上手く扱えません。

dbRDAの計算にはveganのcapscale()とdbrda()が使えます。ここではdbrda()の例を示します。

res1.rda <- dbrda(dat~., data=edat2, distance="horn")
res1.rda
## Call: dbrda(formula = dat ~ PSU + PH + ORP3 + ORP8 + IL + MdPh +
## sqPM + sqTN, data = edat2, distance = "horn")
## 
##               Inertia Proportion Rank RealDims
## Total          4.6471     1.0000              
## Constrained    3.1141     0.6701    8        7
## Unconstrained  1.5329     0.3299    6        5
## Inertia is squared Horn distance 
## 
## Eigenvalues for constrained axes:
##  dbRDA1  dbRDA2  dbRDA3  dbRDA4  dbRDA5  dbRDA6  dbRDA7 idbRDA1 
##  1.1001  0.7293  0.5035  0.4615  0.2383  0.1006  0.0233 -0.0425 
## 
## Eigenvalues for unconstrained axes:
##    MDS1    MDS2    MDS3    MDS4    MDS5   iMDS1 
##  0.8778  0.3836  0.1772  0.1471  0.0103 -0.0631
plot(res1.rda)

Total InertiaとConstrained Inertiaの値を確認しておきます。また、何番目の軸まで固有値が1より大きいか、確認しておきます。

anova.cca(res1.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = dat ~ PSU + PH + ORP3 + ORP8 + IL + MdPh + sqPM + sqTN, data = edat2, distance = "horn")
##          Df SumOfSqs      F Pr(>F)  
## Model     8   3.1141 1.5236  0.066 .
## Residual  6   1.5329                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dbRDAの有意性を調べてみたところ、残念ながら有意とはなりませんでした。

ordistep()にて変数選択を行います。変数を増やす方向か減らす方向かのどちらかを選べますが、ここではフルモデルから変数を減らす方向で変数選択をしてみます。anova.cca()で有意性を検定しておきます。

res0.rda <- dbrda(dat~1, data=edat2, distance="horn")
res2.rda <- ordistep(res1.rda, scope=formula(res0.rda), direction="both")
## 
## Start: dat ~ PSU + PH + ORP3 + ORP8 + IL + MdPh + sqPM + sqTN 
## 
##        Df    AIC      F Pr(>F)  
## - ORP3  1 22.216 0.3468  0.905  
## - ORP8  1 22.365 0.4105  0.865  
## - IL    1 24.001 1.1492  0.370  
## - sqPM  1 24.180 1.2346  0.285  
## - sqTN  1 25.044 1.6637  0.140  
## - MdPh  1 25.209 1.7486  0.080 .
## - PSU   1 25.427 1.8619  0.070 .
## - PH    1 28.043 3.3598  0.020 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PSU + PH + ORP8 + IL + MdPh + sqPM + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - ORP8  1 20.833 0.2943  0.930   
## - IL    1 22.617 1.2153  0.305   
## - sqPM  1 22.752 1.2898  0.250   
## - MdPh  1 23.915 1.9581  0.085 . 
## - sqTN  1 23.591 1.7665  0.075 . 
## - PSU   1 23.934 1.9695  0.060 . 
## - PH    1 27.119 4.0906  0.010 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PSU + PH + IL + MdPh + sqPM + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - sqPM  1 21.302 1.4309  0.175   
## - IL    1 21.874 1.7980  0.115   
## - sqTN  1 22.046 1.9108  0.090 . 
## - PSU   1 22.390 2.1408  0.065 . 
## - MdPh  1 22.555 2.2531  0.050 * 
## - PH    1 25.491 4.4692  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PSU + PH + IL + MdPh + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - PSU   1 21.790 1.6242  0.150   
## - IL    1 21.885 1.6916  0.150   
## - sqTN  1 21.640 1.5180  0.145   
## - MdPh  1 22.804 2.3672  0.045 * 
## - PH    1 25.445 4.5555  0.010 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PH + IL + MdPh + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - IL    1 22.052 1.6276   0.14   
## - sqTN  1 22.478 1.9623   0.07 . 
## - MdPh  1 22.825 2.2420   0.03 * 
## - PH    1 24.876 4.0363   0.01 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PH + MdPh + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - MdPh  1 22.653 2.0825  0.075 . 
## - sqTN  1 23.014 2.4014  0.015 * 
## - PH    1 24.547 3.8428  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res2.rda
## Call: dbrda(formula = dat ~ PH + MdPh + sqTN, data = edat2,
## distance = "horn")
## 
##               Inertia Proportion Rank RealDims
## Total          4.6471     1.0000              
## Constrained    1.9129     0.4116    3        3
## Unconstrained  2.7341     0.5884   11        8
## Inertia is squared Horn distance 
## 
## Eigenvalues for constrained axes:
## dbRDA1 dbRDA2 dbRDA3 
## 0.9743 0.6191 0.3195 
## 
## Eigenvalues for unconstrained axes:
##    MDS1    MDS2    MDS3    MDS4    MDS5    MDS6    MDS7    MDS8   iMDS1 
##  1.0463  0.6487  0.4820  0.3541  0.2177  0.1143  0.0509  0.0016 -0.0159 
##   iMDS2   iMDS3 
## -0.0579 -0.1077
plot(res2.rda)

anova.cca(res2.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = dat ~ PH + MdPh + sqTN, data = edat2, distance = "horn")
##          Df SumOfSqs      F Pr(>F)    
## Model     3   1.9129 2.5654  0.001 ***
## Residual 11   2.7341                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(res2.rda, by="axis")
## Permutation test for dbrda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = dat ~ PH + MdPh + sqTN, data = edat2, distance = "horn")
##          Df SumOfSqs      F Pr(>F)  
## dbRDA1    1  0.97431 3.9199  0.014 *
## dbRDA2    1  0.61912 2.4909  0.106  
## dbRDA3    1  0.31950 1.2854  0.283  
## Residual 11  2.73412                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(res2.rda, by="terms")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = dat ~ PH + MdPh + sqTN, data = edat2, distance = "horn")
##          Df SumOfSqs      F Pr(>F)   
## PH        1  0.85141 3.4254  0.010 **
## MdPh      1  0.46465 1.8694  0.080 . 
## sqTN      1  0.59687 2.4014  0.018 * 
## Residual 11  2.73412                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

変数を減らした結果、説明できた部分(Constrained Proportion)は減りましたが、dbRDAは有意となりました。 係数はcoef()で表示できます。

coef(res2.rda)
##           dbRDA1      dbRDA2     dbRDA3
## PH    0.27375601 -0.02673265  0.0433661
## MdPh -0.03684029  0.22352641  0.1695031
## sqTN  0.10816218  0.22639967 -0.1415342

参考文献

  • Borcard D, Gille F, Legendre P (2011) Numerical Ecology with R. Springer
  • Legendre P, Legendre L (2012) Numerical Ecology. Elsevier
  • Magurran AE, McGill BJ (2011) Biological Diversity. Oxford Univ.

Rによる群集組成の解析に戻る