### R code from vignette source 'Chapter04_Descriptive.prnw' ################################################### ### code chunk number 1: descriptive1 ################################################### library(compositions) GeoChemSed=read.csv("geochemsed.csv",header=TRUE,skip=1) x = acomp(GeoChemSed, parts=c("Cr","Zn","Pb")) ################################################### ### code chunk number 2: Chapter04_Descriptive.prnw:63-65 ################################################### x = acomp(x) mean(x) ################################################### ### code chunk number 3: Chapter04_Descriptive.prnw:70-71 ################################################### mean(x-mean(x)) ################################################### ### code chunk number 4: Chapter04_Descriptive.prnw:96-97 ################################################### mvar(x) ################################################### ### code chunk number 5: Chapter04_Descriptive.prnw:108-109 ################################################### msd(x) ################################################### ### code chunk number 6: Chapter04_Descriptive.prnw:135-136 ################################################### variation(x) ################################################### ### code chunk number 7: variation ################################################### par(mfcol=c(4,4),mar=c(2,2,0,0),oma=c(4,4,0.5,0.5)) lr <- function(x) log(x/(1-x)) xx <- seq(0.001,0.999,by=0.001) slist <- sqrt(c(0.1,0.5,1,2)) mlist <- c(0.5,0.25,0.125,0.0625) for(m in mlist){ for(s in slist) { plot(xx,dnorm(lr(xx),lr(m),sd=s)/xx/(1-xx),type="l") abline(v=m) x0=qnorm(0.025,mean=lr(m), sd=s) x1=qnorm(0.975,mean=lr(m), sd=s) x0 = exp(x0)/(1+exp(x0)) x1 = exp(x1)/(1+exp(x1)) segments(x0,0,x1,0, lwd=5, col="gray")# red? text(x=c(x0,x1),y=c(0,0),label=c("(",")")) } } mtext(side=c(1,2), at =c(1,1)/2, text =c("mean","variation"), line=2.5, out=TRUE) mtext(side=2, at = c(0:3)/4+1/8, text = rev(slist)^2, line=1,outer=TRUE) mtext(side=1, at = c(0:3)/4+1/8, text = mlist, line=1, outer = TRUE) ################################################### ### code chunk number 8: variation2 ################################################### par(mfrow=c(3,3),mar=c(2,2,2,0.5)) xx <- seq(0.00,10,by=0.01) for(s in sqrt(c(0.01,0.02,0.05,0.1,0.2,0.5,1,5,10))) { yy = dlnorm(xx,sd=s) plot(xx,yy,type="l",ylab="",xlab="") title(paste("variation=",s^2)) abline(v=1) x0=qlnorm(0.025,sd=s) x1=qlnorm(0.975,sd=s) segments(x0,max(yy),x1,max(yy), lwd=5, col="gray")#red? text(x=c(x0,x1),y=rep(max(yy),2),label=c("(",")")) } ################################################### ### code chunk number 9: Chapter04_Descriptive.prnw:219-220 ################################################### summary(x) ################################################### ### code chunk number 10: logratio-boxplots ################################################### opar <- par(mar=c(3,3,1,1)) boxplot(x) par(opar) ################################################### ### code chunk number 11: centering1 ################################################### # get the statistics: mn = mean(x) mvr = mvar(x) # center the data: dat1 = scale(x,center=TRUE,scale=TRUE) dat2 = (x-mn)/sqrt(mvr) ################################################### ### code chunk number 12: centering2 ################################################### dat3 = scale(x,center=TRUE,scale=FALSE) dat4 = x-mn ################################################### ### code chunk number 13: centered ################################################### par(mfrow=c(1,2),mar=c(3,3,0,0)) plot(x, pch=".",main="non-centered") plot(dat4, pch=".",main="centered") ################################################### ### code chunk number 14: Chapter04_Descriptive.prnw:374-382 ################################################### # get the mean: mn = mean(x) # compute variance and radii: vr = var(x) df1 = ncol(x)-1 df2 = nrow(x)-ncol(x)+1 rconf = sqrt( qf(p=0.95, df1, df2)*df1/df2 ) rprob = sqrt( qchisq(p=0.95, df=df1) ) ################################################### ### code chunk number 15: elliptic-regions ################################################### # plot the data: opar = par(mar=c(3,3,1,1)) plot(x, cex=0.25) plot(mn, pch=19, cex=0.5, add=TRUE) # plot the ellipses: ellipses(mean=mn, var=vr, r=rconf, col="red") ellipses(mean=mn, var=vr, r=rprob, lty=2) par(opar) ################################################### ### code chunk number 16: descriptive1 ################################################### library("compositions") data("Hydrochem") subcomp = acomp( Hydrochem[,c(7:10,14,17,18)] ) # average in mass proportions mean(subcomp) ################################################### ### code chunk number 17: descriptive2 ################################################### # weight of a mol of each species: mw = c(22.99, 39.098, 34.305, 40.078, 35.453, 96.063, 61.017) mw = acomp(mw) # recast mass proportions to molar proportions subcomp = subcomp-mw mean(subcomp) ################################################### ### code chunk number 18: descriptive3 ################################################### # ionic charges of each species: charges = c(1,1,2,2,-1,-2,-1) unclass(mean(subcomp)) %*% charges ################################################### ### code chunk number 19: descriptive4 ################################################### variation(subcomp) ################################################### ### code chunk number 20: Chapter04_Descriptive.prnw:453-454 ################################################### aux<-variation(subcomp) ################################################### ### code chunk number 21: matrixplot ################################################### cations = acomp( subcomp[,1:4] ) plot(cations, cex=0.5) r = sqrt( qchisq(p=0.95,df=2) ) mn = mean(cations) vr = var(cations) ellipses(mean=mn, var=vr, r=r, lwd=3 ,col="gray")#red? ################################################### ### code chunk number 22: matrixHarker ################################################### opar <- par(mar=c(3,3,1,1)) plot(as.data.frame(cations)) par(opar) ################################################### ### code chunk number 23: matrixplotMg ################################################### plot(cations, cex=0.5, center=TRUE, margin="Mg") mn = mean(cations) ellipses(mean=mn, var=vr, r=r, lwd=3,col="gray")#red? ################################################### ### code chunk number 24: matrixplotrcomp ################################################### plot(cations, cex=0.5, margin="rcomp") ################################################### ### code chunk number 25: balance-example1 (eval = FALSE) ################################################### ## balance(X=x, expr=~CaO/(Na2O*P2O5)) ################################################### ### code chunk number 26: balance-example2 ################################################### x = acomp(GeoChemSed,4:13) balanceBase(X=x, expr=~CaO/(Na2O*P2O5)) ################################################### ### code chunk number 27: codadendrogram ################################################### V = balanceBase(x, expr=~((Al2O3/SiO2)/((Na2O/(CaO/P2O5))/K2O)) / (TiO2/(MnO/(Fe2O3t/MgO))) ) #V = gsi.buildilrBase(gsi.merge2signary(hc[["merge"]]) ) CoDaDendrogram(x, V=V, lwd=2, type="lines", range=c(-5,5)) CoDaDendrogram(x, V=V, box.space=2, col="white",add=TRUE)