smatr: (Standardised) Major Axis Estimation and Testing Routines

• CRAN: http://cran.r-project.org/web/packages/smatr/index.html

> library(smatr)

Loading required package: plyr
-------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
-------------------------------------------------------------------------

Attaching package: 'plyr'

The following object is masked from 'package:maps':

    ozone

The following objects are masked from 'package:Hmisc':

    is.discrete, summarize

The following objects are masked from 'package:dplyr':

    arrange, count, desc, failwith, id, mutate, rename, summarise,
    summarize

> library(dplyr)

バージョン: 3.4.3

leaflife （データセット）

> data("leaflife")
> str(leaflife)

'data.frame':    67 obs. of  5 variables:
 $ site  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
 $ rain  : Factor w/ 2 levels "high","low": 1 1 1 1 1 1 1 1 1 1 ...
 $ soilp : Factor w/ 2 levels "high","low": 1 1 1 1 1 1 1 1 1 1 ...
 $ longev: num  1.115 0.516 0.972 0.672 1.095 ...
 $ lma   : num  125.5 82.3 71 94.7 119.7 ...

> # site... 種の形質データを調査した地点コード
> # rain... 各siteの年降雨量の区分（high or low）
> # solip... 各siteにおける土壌窒素濃度の水準（high or low）
> # longev... 葉寿命（単位は年）
> # lma... 葉の面積あたりの重量(m2/kg)
> summary(leaflife)

 site     rain     soilp        longev            lma        
 1:17   high:34   high:40   Min.   :0.5162   Min.   : 56.94  
 2:17   low :33   low :27   1st Qu.:1.0100   1st Qu.:116.05  
 3:23                       Median :1.4591   Median :153.60  
 4:10                       Mean   :1.6676   Mean   :162.12  
                            3rd Qu.:2.1318   3rd Qu.:202.79  
                            Max.   :4.7304   Max.   :333.89

sma:

II型回帰モデル。複数ある形式にも対応。

leaflife %>% filter(soilp == "low", rain == "low") %>%

> leaflife %$% ma(longev ~ lma, log = 'xy', data = .) # MA

Call: sma(formula = ..1, data = ..3, log = "xy", method = "MA") 

Fit using Major Axis 

These variables were log-transformed before fitting: xy 

Confidence intervals (CI) are at 95%

------------------------------------------------------------
Coefficients:
            elevation    slope
estimate    -3.085214 1.492616
lower limit -3.968020 1.146777
upper limit -2.202407 2.001084

H0 : variables uncorrelated
R-squared : 0.4544809 
P-value : 0.00000000040171

回帰が有意なものであるか（slope = 1）を検証

> leaflife %$% ma(longev ~ lma, log = 'xy', slope.test = 1, data = .) %T>% 
+   summary(.) -> res

Call: sma(formula = ..1, data = ..4, log = "xy", method = "MA", slope.test = 1) 

Fit using Major Axis 

These variables were log-transformed before fitting: xy 

Confidence intervals (CI) are at 95%

------------------------------------------------------------
Coefficients:
            elevation    slope
estimate    -3.085214 1.492616
lower limit -3.968020 1.146777
upper limit -2.202407 2.001084

H0 : variables uncorrelated
R-squared : 0.4544809 
P-value : 0.00000000040171 

------------------------------------------------------------
H0 : slope not different from 1 
Test statistic : r= 0.3515 with 65 degrees of freedom under H0
P-value : 0.0035393

0.0035393 であるので１より有意（P < 0.05）に異なる傾き slopeであるという結果となった

残差プロットを描く

> plot(res, type = "residual")

（なにも描画されない？？）

> # 土壌窒素濃度が低いサイトのデータを抽出
> leaf.low.soilp <- filter(leaflife, soilp == "low")

Error in filter(leaflife, soilp == "low"): object 'soilp' not found

> # nrow(leaf.low.soilp) # 27

降水量の高いサイト、低いサイトそれぞれでSMAによる当てはめ、傾きの違いを比較

> leaf.low.soilp %>% sma(formula = "longev ~ lma * rain", 
+                        log = "xy", data = .) %>% {
+   print(.)
+   plot(.)
+ }

年降雨量の多さ（high and low）が葉寿命に効いていないという帰無仮説H0に対して、P = 0.12395が得られた。帰無仮説を棄却することはできないということになる。