library(StatDA)
data(chorizon)
data(kola.background)
attach(chorizon)
X=chorizon[,"XCOO"]
Y=chorizon[,"YCOO"]
# classical
res=lm(log10(Be) ~ Al_XRF+Ca_XRF+Fe_XRF+K_XRF+Mg_XRF+Mn_XRF+Na_XRF+P_XRF+Si_XRF+Ti_XRF)
#######################################################################################
# plot with map:
# true representation of x and y axis of map for plot
xwid=diff(range(X))/12e4
ywid=diff(range(Y))/12e4
pdf("fig-16-7.pdf",width=2*xwid,height=2*ywid)
par(mfrow=c(2,2),mar=c(1.5,1.5,1.5,1.5))
### NEW MAP:
dat=res$resid
# generate plot with packground
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T)
SymbLegend(X,Y,dat,type="boxplot",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8,
leg.position="topright",leg.title="Be [mg/kg]",leg.title.cex=1.0,leg.round=2,leg.wid=4,
leg.just="right",cex.scale=0.8,xf=9e3,logscale=FALSE,accentuate=FALSE)
# Class legend
text(min(X)+diff(range(X))*5/8,max(Y),"LS residuals",cex=1.0)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.5,shiftkm=37e3,sizetext=0.8)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7838750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)
#####################################################################################
par(mar=c(4,4,2,2))
plot(res$fitted,res$resid,cex.lab=1.2,xlab="Fitted values",ylab="Residuals")
abline(h=0,col="grey",lty=2)
boxinfo=boxplot(dat,plot=F)
abline(h=boxinfo$stats[c(1,5)],lty=3,cex=1.1)
#####################################################################################
plot.lm(res,which=2,caption="",sub.caption="",id.n=2,cex.lab=1.2,pch=3,cex=0.9)
qqline(residuals(res)/sd(residuals(res)),datax=FALSE,lwd=1)
plot.lm(res,which=4,caption="",sub.caption="",id.n=4,cex.lab=1.2)
dev.off()