getwd() library(ade4) read.table("http://www.bondy.ird.fr/~blouin/data/datanitrogenearthworm.txt", header=TRUE, na.string=".")->ver names(ver) azo<-as.ordered(ver$nitrogen) levels(azo) dec<-as.factor(ver$chuniodrilus) levels(dec) levels(dec)<-c("abs", "pre") dec com<-as.factor(ver$millsonia) levels(com)<-c("abs", "pre") com tot<-ver$total bas<-ver$shoots brs<-ver$roots fau<-ver$fauna lmtot<-lm(tot~azo*com*dec) lmtot anova(lmtot) lmtot<-lm(tot~azo+com) lmtot anova(lmtot) lmbas<-lm(bas~azo*com*dec) lmbas anova(lmbas) lmbas<-lm(bas~azo+com) lmbas anova(lmbas) lmbrs<-lm(brs~azo*com*dec) lmbrs anova(lmbrs) lmbrs<-lm(brs~azo+com) lmbrs<-lm(brs~com) lmbrs anova(lmbrs) azoq<-log(ver$nitrogen+1) lmlogtot<-lm(tot~azoq+I(azoq^2)+I(azoq^3)+com) anova(lmlogtot) coefficients(lmlogtot) lmlogtot<-lm(tot~azoq+I(azoq^2)+com) coefficients(lmlogtot) lmlogbas<-lm(bas~azoq+I(azoq^2)+I(azoq^3)+com) anova(lmlogbas) lmlogbas<-lm(bas~azoq+I(azoq^2)+com) anova(lmlogbas) coefficients(lmlogbas) lmlogbrs<-lm(brs~azoq+I(azoq^2)+I(azoq^3)+com) anova(lmlogbrs) lmlogbrs<-lm(brs~com) anova(lmlogbrs) par(mfrow=c(3,1), mar=c(2,5,1,1)) x0<-seq(0,8,le=100) plot(azoq, tot, pch=c(1,3,2,4)[fau], xlab="", ylab="Total dry biomass (g)", axes=F, cex=1.3, cex.lab=1.5 ) axis(2, 3:15, las=1, cex.axis=1.2) coefficients(lmlogtot) lines(x0, 6.3697747-0.7193702*x0+0.1562327*x0^2,lwd=2, lty=1) lines(x0, 6.3697747-0.7193702*x0+0.1562327*x0^2+1.4283333,lwd=2, lty=2) plot(azoq, bas, pch=c(1,3,2,4)[fau], las=1, xlab="", ylab="Shoot dry biomass (g)", axes=F, cex=1.3, cex.lab=1.5 ) axis(2, 0:14, las=1, cex.axis=1.2) coefficients(lmlogbas) lines(x0,4.7990329-0.6798388*x0+0.1485498*x0^2,lwd=2, lty=1) lines(x0,4.7990329-0.6798388*x0+0.1485498*x0^2+0.8733333,lwd=2, lty=2) plot(azoq, brs, pch=c(1,3,2,4)[fau], las=1, xlab="", ylab="Root dry biomass (g)", axes=F, cex=1.3, cex.lab=1.5 ) axis(2, 0:4, las=1, cex.axis=1.2) coefficients(lmlogbrs) lines(x0, 1.5922633+0*x0, lwd=2, lty=1) lines(x0, 1.5922633+0*x0+0.5529433,lwd=2, lty=2)