# Examples for Multiple Regression
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-6.pdf",width=2*xwid,height=1*ywid)
par(mfrow=c(1,2),mar=c(1.5,1.5,1.5,1.5))
### NEW MAP:
dat=chorizon[,"Be"]
# 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=0.8,leg.round=2,leg.wid=4,
leg.just="right",cex.scale=0.75,xf=9e3,logscale=TRUE,accentuate=FALSE)
# Class legend
text(min(X)+diff(range(X))*5/8,max(Y),"Original data",cex=0.8)
# 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)
### NEW MAP:
dat=10^(res$fitted)
# 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=0.8,leg.round=2,leg.wid=4,
leg.just="right",cex.scale=0.75,xf=9e3,logscale=TRUE,accentuate=FALSE)
# Class legend
text(min(X)+diff(range(X))*5/8,max(Y),"LS prediction",cex=0.8)
# 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)
dev.off()