匠人府分享 http://blog.sciencenet.cn/u/meiweipingg

博文

如何作 " 协方差分析 (ANCOVA)" 之 R (或 RStudio)应用篇

已有 15286 次阅读 2016-5-8 09:53 |个人分类:R语言|系统分类:科研笔记| and, 协方差分析, RStudio, ANCOVA, 协变量

简单应用来说,协方差分析 (ANCOVA) 就是在R语言中表达为y~x+/*y0,其中因变量y,协变量y0,自变量x。然后,开始写代码咯

> data.of.Wada2013.practice <- read.csv("~data of Wada2013 practice.csv")

> View(data.of.Wada2013.practice)

> eco<-data.of.Wada2013.practice

> eco

     site  d13C  d15N

1     1 -28.40  4.05

2     1 -26.90  6.70

3     1 -25.90  7.30

4     1 -25.85  7.85

5     1 -25.60  9.90

6     1 -26.35  10.30

7     1 -26.25  10.70

8     1 -26.55  12.20

9     1 -26.60  12.65

10      1 -26.50  13.25

11      1 -26.40  13.70

12      1 -24.90  14.10

13      2 -24.00   9.02

14      2 -23.50  11.80

15      2 -22.71 15.78

16      2 -22.10 15.10

17      2 -20.50 14.00

18      2 -19.75 17.30

19      2 -20.35 17.95

20      3 -23.80  3.95

21      3 -23.50  4.80

22      3 -23.80  5.20

23      3 -23.70  6.10

24      3 -22.80  6.95

25      3 -21.95  7.10


> attach(eco)
>sitef<-factor(eco$site)

> #Step1: ANCOVA with InterAction

> mod.IA<-aov(d15N~d13C*sitef,data=eco,  contrasts=list(sitef=contr.sum))

                   #aov() function, "contr.sum" for factor "sitef"(not  numerical "site")

> summary(mod.IA) ???              # this summary is wroooong, it's for type I , not type III, don't do this!

> library(car)

> Anova(mod.IA,type=3)           # this is riiiiight. Anova() function, not anova  

>Anova Table (Type III tests)

Response: d15N

                 Sum Sq       Df           F value      Pr(>F)

(Intercept) 70.098        1            13.0456   0.001858 **

d13C           43.563       1             8.1073    0.010302 *

sitef            3.341         2           0.3109      0.736451  

d13C:sitef   1.771        2             0.1648     0.849262

Residuals   102.093      19                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


# covariate d13C is significantly differ with d15N ( p = 0.001858),

 and interation of  d13C:sitef is not significantly differ ( p = 0.849262 ),

------->

then we can proceed ANCOVA next.


> # Step2: then ANCOVA NonInterAction, for equal slope, type III

> mod.NIA<-aov(d15N~d13C+sitef,data=eco,contrasts = list(sitef=contr.sum))

> Anova(mod.NIA, type=3)

Anova Table (Type III tests)

Response: d15N

                    Sum Sq     Df      F value      Pr(>F)  

(Intercept)   117.10       1        23.677     8.232e-05 ***

d13C              74.32       1         15.027     0.0008727 ***

sitef               258.75     2         26.158     1.989e-06 ***

Residuals      103.86     21                      

---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> # p value of 8.232e-05 *** means intercept is significantly not equal to 0.

> summary.lm(mod.NIA)

Call:

aov(formula = d15N ~ d13C + sitef, data = eco, contrasts = list(sitef = contr.sum))

Residuals:

   Min        1Q          Median       3Q          Max

-3.6695    -1.5658      0.1179     1.1474    3.5577

Coefficients:

                    Estimate        Std. Error     t value       Pr(>|t|)

(Intercept)   49.51511         0.1760       4.866       8.23e-05 ***

d13C            1.6545              0.4268        3.876      0.000873 ***

sitef1           4.3049             1.2335        3.490       0.002183 **

sitef2            1.0468             1.0774        0.972       0.342273

---Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.224 on 21 degrees of freedom

Multiple R-squared:  0.7556,    Adjusted R-squared:  0.7207

F-statistic: 21.65 on 3 and 21 DF,  p-value: 1.252e-06

#Based on above, we could get the equation with equal slope and common intercept of 3 sites which suitable with Wada_2013 as: d15N=1.6545[±0.4268]*d13C+49.51511[±0.1760 ]+[ecosystem specific sites constant] ( p = 0.000873).


> # Step3: then  ANCOVA NonInterAction with treatment contrasts, or type I.

> mod.NIA.tc<-aov(d15N~d13C+sitef,data=eco,contrasts = list(sitef=contr.treatment))

> summary.lm(mod.NIA.tc)          #this step is riiiight,for graphics with type I is what we need.

Call:

aov(formula = d15N ~ d13C + sitef, data = eco, contrasts = list(sitef = contr.treatment))

Residuals:

  Min          1Q          Median      3Q           Max

-3.6695     -1.5658     0.1179      1.1474    3.5577

Coefficients:

                         Estimate     Std. Error      t value       Pr(>|t|)    

(Intercept)         53.8200         11.2645       4.778      0.000101 ***

d13C                  1.6545           0.4268         3.876      0.000873 ***

sitef2                  -3.2581          2.1947        -1.485     0.152529    

sitef3                 -9.6567           1.7256        -5.596     1.49e-05 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.224 on 21 degrees of freedom

Multiple R-squared:  0.7556,    Adjusted R-squared:  0.7207

F-statistic: 21.65 on 3 and 21 DF,  p-value: 1.252e-06

#Based on above, we could get the equation for  3 sites respectively with equal slope and specific intercepts which suitable with Wada_2013 as:

site1: d15N=1.6545*d13C+53.8200 ;

site2: d15N=1.6545*d13C+53.8200-3.2581 ;

site1: d15N=1.6545*d13C+53.8200-9.6567 .

> # Step4: Graphics

> dat1<-subset(eco,sitef == "1")                  

> dat2<-subset(eco,sitef == "2")

> dat3<-subset(eco,sitef=="3")

> reg1<-lm(d15N~d13C,data=dat1)

> reg2<-lm(d15N~d13C,data=dat2)

> reg3<-lm(d15N~d13C,data=dat3)

> plot(d15N~d13C,type="n")         # type="n" means plot with axes,without any plots.

> points(dat1$d13C,dat1$d15N,pch=1)

> points(dat2$d13C,dat2$d15N,pch=2)

> points(dat3$d13C,dat3$d15N,pch=3)

> abline(reg1,lty=1)

> abline(reg2,lty=2)

> abline(reg3,lty=3)

> legend("topleft",c("L.Baikal","L.Biwa","Mongolian grassland"),lty=c(1,2,3),pch=c(1,2,3))

> #Above procedures output a graph with unequal slopes which we also need.

> summod.NIA.tc<-summary.lm(mod.NIA.tc)        #type I

> coeffs.NIA.tc<-coef(summod.NIA.tc)

> coeffs.NIA.tc

                       Estimate           Std. Error         t value            Pr(>|t|)

(Intercept)      53.820034         11.2645326    4.777831       1.014464e-04

d13C              1.654461           0.4268016      3.876416        8.726695e-04

sitef2             -3.258098          2.1947260     -1.484512       1.525294e-01

sitef3             -9.656707          1.7255811      -5.596206      1.488283e-05

> I1<-coeffs.NIA.tc[1,1]

> I2<-coeffs.NIA.tc[3,1]+I1

> I3<-coeffs.NIA.tc[4,1]+I1

> solpeAll<-coeffs.NIA.tc[2,1]

> plot(d15N~d13C,type="n")

> points(dat1$d13C,dat1$d15N,pch=1)

> points(dat2$d13C,dat2$d15N,pch=2)

> points(dat3$d13C,dat3$d15N,pch=3)

> abline(I1,slopeAll,lty=1)

> abline(I2,slopeAll,lty=2)

> abline(I3,slopeAll,lty=3)

#From above,we could clearly understand why we can not use summary.lm(mod.NIA) for graphics.

> legend("topleft",c("L.Baikal","L.Biwa","Mongolian grassland"),lty=c(1,2,3),pch=c(1,2,3))


OMG, so long so tired. Anyway, we correctly get what we need, that's enough ^_^



http://blog.sciencenet.cn/blog-651374-975758.html

上一篇:免费正版科研软件,重磅推荐
下一篇:Steel.Dwass test in R: Kruskal-Wallis post hoc test

1 梅志平

该博文允许注册用户评论 请点击登录 评论 (5 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-9-16 18:01

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部