8.2 Graphical Comparison of the Data Distributions of Several Data Sets

# Fig. 8.1: Density traces of some elements in 4 layers
library(StatDA)
data(moss)
data(ohorizon)
data(bhorizon)
data(chorizon)


pdf("fig-8-1.pdf",width=8,height=8)
par(mfrow=c(2,2),mar=c(4,4,2,2))

nam="Al"
M=density(log10(moss[,nam]))
O=density(log10(ohorizon[,nam]))
B=density(log10(bhorizon[,nam]))
C=density(log10(chorizon[,nam]))
plot(0,0,xlim=c(min(M$x,O$x,B$x,C$x),max(M$x,O$x,B$x,C$x)),
ylim=c(min(M$y,O$y,B$y,C$y),max(M$y,O$y,B$y,C$y)),
main="",xlab=paste(nam," [mg/kg]",sep=""),ylab="Density",
xaxt="n",cex.lab=1.2,type="n")
axis(1,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
lines(M,lty=1)
lines(O,lty=2)
lines(B,lty=4)
lines(C,lty=5)
text(1.5,1.5,"Moss",pos=4)
text(2.5,1.45,"O-horizon",pos=4)
text(4.6,0.9,"B-horizon",pos=4,srt=90)
text(3.7,1.4,"C-horizon",pos=4,srt=90)

nam="K"
M=density(log10(moss[,nam]))
O=density(log10(ohorizon[,nam]))
B=density(log10(bhorizon[,nam]))
C=density(log10(chorizon[,nam]))
plot(0,0,xlim=c(min(M$x,O$x,B$x,C$x),max(M$x,O$x,B$x,C$x)),
ylim=c(min(M$y,O$y,B$y,C$y),max(M$y,O$y,B$y,C$y)),
main="",xlab=paste(nam," [mg/kg]",sep=""),ylab="Density",
xaxt="n",cex.lab=1.2,type="n")
axis(1,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
lines(M,lty=1)
lines(O,lty=2)
lines(B,lty=4)
lines(C,lty=5)
text(3.7,4,"Moss",pos=4)
text(2.6,3.9,"O-horizon",pos=4)
text(1.8,1,"B-horizon",pos=4)
text(3.3,0.9,"C-horizon",pos=4,srt=90)


nam="Pb"
M=density(log10(moss[,nam]))
O=density(log10(ohorizon[,nam]))
B=density(log10(bhorizon[,nam]))
C=density(log10(na.omit(chorizon[,nam])))
plot(0,0,xlim=c(min(M$x,O$x,B$x,C$x),max(M$x,O$x,B$x,C$x)),
ylim=c(min(M$y,O$y,B$y,C$y),max(M$y,O$y,B$y,C$y)),
main="",xlab=paste(nam," [mg/kg]",sep=""),ylab="Density",
xaxt="n",cex.lab=1.2,type="n")
axis(1,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
lines(M,lty=1)
lines(O,lty=2)
lines(B,lty=4)
lines(C,lty=5)
text(-0.33,2.2,"Moss",pos=4)
text(1.5,1.1,"O-horizon",pos=4)
text(0.35,1.25,"B-horizon",pos=4,srt=90)
text(-0.4,0.7,"C-horizon",pos=4,srt=85)


nam="S"
M=density(log10(moss[,nam]))
O=density(log10(ohorizon[,nam]))
B=density(log10(bhorizon[,nam]))
C=density(log10(chorizon[,nam]))
plot(0,0,xlim=c(min(M$x,O$x,B$x,C$x),max(M$x,O$x,B$x,C$x)),
ylim=c(min(M$y,O$y,B$y,C$y),max(M$y,O$y,B$y,C$y)),
main="",xlab=paste(nam," [mg/kg]",sep=""),ylab="Density",
xaxt="n",cex.lab=1.2,type="n")
axis(1,at=log10(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)
lines(M,lty=1)
lines(O,lty=2)
lines(B,lty=4)
lines(C,lty=5)
text(2.3,5,"Moss",pos=4)
text(3.4,2,"O-horizon",pos=4,srt=90)
text(1.7,1.9,"B-horizon",pos=4)
text(0.8,1.7,"C-horizon",pos=4)


dev.off()


# Fig.8.2.: CP-plots of some elements in 4 layers
library(StatDA)
data(moss)
data(ohorizon)
data(bhorizon)
data(chorizon)

pdf("fig-8-2.pdf",width=8,height=8)
par(mfrow=c(2,2),mar=c(4,4,2,2))

nam="Al"
M=log10(moss[,nam])
O=log10(ohorizon[,nam])
B=log10(bhorizon[,nam])
C=log10(chorizon[,nam])
qpplot.das(M,qdist=qnorm,xlab=paste(nam," [mg/kg]",sep=""),
ylab="Probability [%]", pch=3,cex=0.7, logx=TRUE,
logfinetick=c(2,5,10),logfinelab=c(2,5,10),line=F,cex.lab=1.2,
xlim=c(min(M,O,B,C),max(M,O,B,C)))
points(sort(O),qnorm(ppoints(length(O))),pch=4,cex=0.7)
points(sort(B),qnorm(ppoints(length(B))),pch=1,cex=0.7)
points(sort(C),qnorm(ppoints(length(C))),pch=22,cex=0.7)
legend("topleft",legend=c("Moss","O-horizon","B-horizon","C-horizon"),
pch=c(3,4,1,22),bg="white")


nam="K"
M=log10(moss[,nam])
O=log10(ohorizon[,nam])
B=log10(bhorizon[,nam])
C=log10(chorizon[,nam])
qpplot.das(M,qdist=qnorm,xlab=paste(nam," [mg/kg]",sep=""),
ylab="Probability [%]", pch=3,cex=0.7, logx=TRUE,
logfinetick=c(2,5,10),logfinelab=c(2,5,10),line=F,cex.lab=1.2,
xlim=c(min(M,O,B,C),max(M,O,B,C)))
points(sort(O),qnorm(ppoints(length(O))),pch=4,cex=0.7)
points(sort(B),qnorm(ppoints(length(B))),pch=1,cex=0.7)
points(sort(C),qnorm(ppoints(length(C))),pch=22,cex=0.7)
legend("topleft",legend=c("Moss","O-horizon","B-horizon","C-horizon"),
pch=c(3,4,1,22),bg="white")



nam="Pb"
M=log10(moss[,nam])
O=log10(ohorizon[,nam])
B=log10(bhorizon[,nam])
C=log10(na.omit(chorizon[,nam]))
qpplot.das(M,qdist=qnorm,xlab=paste(nam," [mg/kg]",sep=""),
ylab="Probability [%]", pch=3,cex=0.7, logx=TRUE,
logfinetick=c(2,5,10),logfinelab=c(2,5,10),line=F,cex.lab=1.2,
xlim=c(min(M,O,B,C),max(M,O,B,C)))
points(sort(O),qnorm(ppoints(length(O))),pch=4,cex=0.7)
points(sort(B),qnorm(ppoints(length(B))),pch=1,cex=0.7)
points(sort(C),qnorm(ppoints(length(C))),pch=22,cex=0.7)
legend("topleft",legend=c("Moss","O-horizon","B-horizon","C-horizon"),
pch=c(3,4,1,22),bg="white")



nam="S"
M=log10(moss[,nam])
O=log10(ohorizon[,nam])
B=log10(bhorizon[,nam])
C=log10(chorizon[,nam])
qpplot.das(M,qdist=qnorm,xlab=paste(nam," [mg/kg]",sep=""),
ylab="Probability [%]", pch=3,cex=0.7, logx=TRUE,
logfinetick=c(2,5,10),logfinelab=c(2,5,10),line=F,cex.lab=1.2,
xlim=c(min(M,O,B,C),max(M,O,B,C)))
points(sort(O),qnorm(ppoints(length(O))),pch=4,cex=0.7)
points(sort(B),qnorm(ppoints(length(B))),pch=1,cex=0.7)
points(sort(C),qnorm(ppoints(length(C))),pch=22,cex=0.7)
legend("topleft",legend=c("Moss","O-horizon","B-horizon","C-horizon"),
pch=c(3,4,1,22),bg="white")




dev.off()


# Fig. 8.3.: boxplots of some elements in 4 layers
library(StatDA)
data(moss)
data(ohorizon)
data(bhorizon)
data(chorizon)

pdf("fig-8-3.pdf",width=8,height=8)
par(mfrow=c(2,2),mar=c(4,4,2,2))

nam="Al"
M=moss[,nam]
O=ohorizon[,nam]
B=bhorizon[,nam]
C=chorizon[,nam]
boxplotlog(C,B,O,M,horizontal=T,xlab=paste(nam," [mg/kg]",sep=""),cex=0.7,cex.axis=1.2,
log="x",cex.lab=1.2,pch=3,names=c("C-hor","B-hor","O-hor","Moss"),xaxt="n")
axis(1,at=(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)

nam="K"
M=moss[,nam]
O=ohorizon[,nam]
B=bhorizon[,nam]
C=chorizon[,nam]
boxplotlog(C,B,O,M,horizontal=T,xlab=paste(nam," [mg/kg]",sep=""),cex=0.7,cex.axis=1.2,
log="x",cex.lab=1.2,pch=3,names=c("C-hor","B-hor","O-hor","Moss"),xaxt="n")
axis(1,at=(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)



nam="Pb"
M=moss[,nam]
O=ohorizon[,nam]
B=bhorizon[,nam]
C=chorizon[,nam]
boxplotlog(C,B,O,M,horizontal=T,xlab=paste(nam," [mg/kg]",sep=""),cex=0.7,cex.axis=1.2,
log="x",cex.lab=1.2,pch=3,names=c("C-hor","B-hor","O-hor","Moss"),xaxt="n")
axis(1,at=(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)



nam="S"
M=moss[,nam]
O=ohorizon[,nam]
B=bhorizon[,nam]
C=chorizon[,nam]
boxplotlog(C,B,O,M,horizontal=T,xlab=paste(nam," [mg/kg]",sep=""),cex=0.7,cex.axis=1.2,
log="x",cex.lab=1.2,pch=3,names=c("C-hor","B-hor","O-hor","Moss"),xaxt="n")
axis(1,at=(a<-sort(c((10^(-50:50))%*%t(c(2,5,10))))),labels=a)




dev.off()