🔍 检索绿色通道
     0.ESI查询      1.JCR-IF查询      2.CiteScore-IF查询     3.SCI检索     4.中科院分区查询     5.SSCI检索     6.AHCI检索     7.ESCI检索(新增)     8.EI检索     9.EI检索备用     10.ISTP/CPCI检索备用(新增)      11.梅斯期刊智能查询     12.CSCD检索     13.CSSCI检索     14.北大核心检索      15.国自然官方查询     16.国自然梅斯查询     17.国自然科学网查询     18.论文被检索情况查询方法

Run ANCOVA correctly in R software


(本文于 2016-5-8 09:53 首发于 “科学网”)

Briefly, covariance analysis (ANCOVA) is expressed in the R language as y ~ x +/* y0, where the variable is y, the covariate is y0, the independent variable is x. Then, go ahead code test.

1
2
3
4
5
6
7
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)

> <font color=red>#Step1: ANCOVA with InterAction</font>

> mod.IA<-aov(d15N~d13C*sitef,data=eco,  contrasts=list(sitef=contr.sum))
N.B., aov() function, “contr.sum” for factor “sitef”(not numerical “site”)
> summary(mod.IA) <font color=darkgreen>???              # this summary is <font color=red>wroooong</font>, it's for type I , not type III, don't do this!</font>

> library(car)

> Anova(mod.IA,type=3)           <font color=darkgreen># this is riiiiight. Anova() function, not anova  </font>

>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.<font color=red>NIA</font>.tc)          <font color=darkgreen>#this step is riiiight,for graphics with type I is what we need.</font>

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.6545d13C+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")         <font color=darkgreen># type="n" means plot with axes,without any plots.</font>

> 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)
Aboveall,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 have correctly achieved what we need, that’s enough ^_^



<已有 次阅读>


由于本文作者水平有限,文中如有错误之处,欢迎大家批评指正!

① 本文仅代表作者个人观点,不代表任何其它立场,欢迎交流合作!

② 转载与分享请注明:本文源于为学为研网 http://meiweiping.cn

(>看完记得五星好评哦亲<)

分享到:
    友情链接
    e教学: 蓝墨云班课     传课网     慕课网     作业帮
    BBS: 小木虫     丁香园     科学网bbs     零点花园     数据狗论坛     ResearchGate
    预印本: 中国科技论文在线(教育部)     ChinaXiv(中科院)     bioRxiv(冷泉港)     arXiv(康奈尔大学)
    科研资讯: 科学网   科学之家     51science     Publons     J-STAGE     日本の研究
    信息查询: 中国研究生招生信息网     学信网(学历查询)     中国人事考试网     更多>>
    网站建设: Hexo     极简图床     sm.ms     html-unicode     UTF-8转码     颜色代码
    关于本站: 关于我们      Logo含义      广告合作      免责声明      联系我们      RSS订阅