6.2 Linear Regression Lines

# Fig. 6.3.: xy-plot with LS-line
library(StatDA)
pdf("fig-6-3.pdf",width=9,height=4.5)
par(mfrow=c(1,2),mar=c(4,4,2,2))

data(moss)
Cu=moss[,"Cu"]
Ni=moss[,"Ni"]
plot(log10(Cu),log10(Ni),xlab="Cu in Moss [mg/kg]",
ylab="Ni in Moss [mg/kg]", pch=3, cex=0.7, cex.lab=1.2,xaxt="n",
yaxt="n")
axis(1,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
axis(2,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
abline(lsfit(log10(Cu),log10(Ni)))



data(chorizon)
coun=chorizon[,"COUN"]
Sr=chorizon[coun=="RUS","Sr"]
Ca=chorizon[coun=="RUS","Ca"]

plot(log10(Sr),log10(Ca),xlab="Sr in C-horizon, Russia [mg/kg]",
ylab="Ca in C-horizon, Russia [mg/kg]", pch=3, cex=0.7, cex.lab=1.2,xaxt="n",
yaxt="n")
axis(1,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
axis(2,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
abline(lsfit(log10(Sr),log10(Ca)))

dev.off()
# Fig. 6.4.: xy-plot with robust line
library(StatDA)
pdf("fig-6-4.pdf",width=9,height=4.5)
par(mfrow=c(1,2),mar=c(4,4,2,2))

# left plot:
data(moss)
Cu=moss[,"Cu"]
Ni=moss[,"Ni"]
plot(log10(Cu),log10(Ni), pch=3, cex=0.7, cex.lab=1.2,xaxt="n", yaxt="n",
xlab="Cu in Moss [mg/kg]", ylab="Ni in Moss [mg/kg]")
axis(1,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
axis(2,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
abline(lsfit(log10(Cu),log10(Ni)),lty=2)
abline(ltsReg(log10(Cu),log10(Ni)))
legend("bottomright",legend=c("Robust line","LS line"),lty=c(1,2))

# right plot:
data(chorizon)
coun=chorizon[,"COUN"]
Sr=chorizon[coun=="RUS","Sr"]
Ca=chorizon[coun=="RUS","Ca"]
plot(log10(Sr),log10(Ca), pch=3, cex=0.7, cex.lab=1.2,xaxt="n", yaxt="n",
xlab="Sr in C-horizon, Russia [mg/kg]", ylab="Ca in C-horizon, Russia [mg/kg]")
axis(1,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
axis(2,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
abline(lsfit(log10(Sr),log10(Ca)),lty=2)
abline(ltsReg(log10(Sr),log10(Ca)))
legend("bottomright",legend=c("Robust line","LS line"),lty=c(1,2))

dev.off()