tests/TestSimp.R

In compositions: Compositional Data Analysis

options(warn=1)
# ad-hod function for comparing numbers
qround <- function(x) round(x, dig=3)

library("compositions")

#js <- read.table("juraset.dat",skip=17,header=TRUE)
#js$Land <- factor( c("Wald","Weide","Wiese","Acker")[js$Land] )
#js$Rock <- factor( c("Argovian","Kimmeridgian","Sequanian","Portlandian","Quaternary")[js$Rock] )

#Rock <- js$Rock
#Land <- js$Land

#cdata <- js[,c("Cd","Cu","Pb","Co","Cr","Ni","Zn")]
#cdata <- data.matrix(js[,c("Cd","Cu","Pb","Co","Cr","Ni","Zn")])

#cd1 <- cdata[,1:3]
#cd2 <- cdata[,4:7]
data(SimulatedAmounts)
cdata <- sa.groups5
Land  <- sa.groups5.area
cd1 <- cdata[,1:3]
cd2 <- cdata[,4:5]

# Transformations
# clo 
clo(c(1,2,3))
clo(matrix(1:4,ncol=2))
data(iris)
clo(iris[1:10,1:4])
clo(0.5)
clo(matrix(0.5))
clo(matrix(0.5,nrow=5))
clo(matrix(0.5,ncol=5))

clo(iris[1:10,],c("Sepal.Length","Petal.Length"))
clo(iris[1:10,],c(2,3))

checker <- function(x,y) {
  x<-unclass(x)
  y<-unclass(y)
  if( sum(c(x-y)^2) > 1E-10 )
    stop("Wrong results in ", deparse(substitute(x)))
  x
}

qround( clr(cdata[1:10,]) )

qround( checker( clrInv(clr(cdata)) , clo(cdata) ) )
qround( ilr(cdata) )
qround( checker( ilrInv(ilr(cdata)) , clo(cdata) ) )
qround( alr(cdata) )
qround( checker( alrInv(alr(cdata)) , clo(cdata) ) )
qround( cpt(cdata) )
qround( checker( cptInv(cpt(cdata)) , clo(cdata) ) )
qround( ipt(cdata) )
qround( checker( iptInv(ipt(cdata)) , clo(cdata) ) )
qround( apt(cdata) )
qround( checker( aptInv(apt(cdata)) , clo(cdata) ) )
qround( ilt(cdata) )
qround( checker( iltInv(ilt(cdata)) , cdata ) )
qround( iit(cdata) )
qround( checker( iitInv(iit(cdata)) , cdata ) )

qround( clr(c(a=1,2,3)))
qround( ilr(c(a=1,2,3)))
qround( alr(c(a=1,2,3)))
qround( cpt(c(a=1,2,3)))
qround( ipt(c(a=1,2,3)))
qround( apt(c(a=1,2,3)))
qround( ilt(c(a=1,2,3)))
qround( iit(c(a=1,2,3)))
qround( checker( clrInv(clr(c(a=1,2,3))) , clo(c(1,2,3))) )
qround( checker( ilrInv(ilr(c(a=1,2,3))) , clo(c(1,2,3))) )
qround( checker( alrInv(alr(c(a=1,2,3))) , clo(c(1,2,3))) )
qround( checker( cptInv(cpt(c(a=1,2,3))) , clo(c(1,2,3))) )
qround( checker( iptInv(ipt(c(a=1,2,3))) , clo(c(1,2,3))) )
qround( checker( aptInv(apt(c(a=1,2,3))) , clo(c(1,2,3))) )
qround( checker( iltInv(ilt(c(a=1,2,3))) , c(1,2,3)) )
qround( checker( iitInv(iit(c(a=1,2,3))) , c(1,2,3)) )

# mean

qround( mean(acomp(cdata)) )
qround( mean(rcomp(cdata)) )
qround( mean(aplus(cdata)) )
qround( mean(rplus(cdata)) )

qround( meanCol(cdata) )
qround( meanCol(clo(cdata)) )
qround( clo(meanCol(cdata)) )

# var (Variation Matrix)
qround( var(rcomp(cdata)) )
qround( var(acomp(cdata)) )
qround( var(aplus(cdata)) )
qround( var(rplus(cdata)) )

# clr
qround( clr(mean(acomp(cdata))) )
qround( meanCol(clr(cdata)) )
qround( clrInv(meanCol(clr(cdata))) )
qround( mean(acomp(cdata)) )

# ilr
qround( ilr(mean(acomp(cdata))) )
qround( meanCol(ilr(cdata)) )
qround( ilrInv(meanCol(ilr(cdata))) )
qround( mean(acomp(cdata)) )

# alr
qround( alr(mean(acomp(cdata))) )
qround( meanCol(alr(cdata)) )
qround( alrInv(meanCol(alr(cdata))) )
qround( mean(acomp(cdata)) )

# Operations
mean(acomp(3 * (cdata - mean(acomp(cdata)))))


# barplot
barplot(acomp(cdata[1:10,]))
barplot(rcomp(cdata[1:10,]))
barplot(aplus(cdata[1:10,]))
barplot(rplus(cdata[1:10,]))

barplot(mean(acomp(cdata)))
barplot(mean(rcomp(cdata)))
barplot(mean(aplus(cdata)))
barplot(mean(rplus(cdata)))

# piechart
pie(mean(acomp(cdata)))
pie(mean(rcomp(cdata)))
pie(mean(aplus(cdata)))



# Triangular Diagrams
plot(acomp(cdata[,1:3]))
qround( mean(acomp(cdata[,1:3])) )

# In this file we suppress unnecessary warnings targeting end-users
# Human: if you want to see them, replace suppressWarnings() by I()
suppressWarnings( plot(acomp(cdata),margin="rcomp",pca=TRUE) )
plot(acomp(cdata),margin="acomp",pca=TRUE)
plot(acomp(cdata),margin="Cd",pca=TRUE) # bug corrected
qround( mean(acomp(cdata)) )

boxplot(acomp(cdata))              # boxplotscale
suppressWarnings( boxplot(acomp(cdata),Land,notch=TRUE) )# notch
#boxplot(acomp(cdata),Rock)

boxplot(acomp(cdata),log=FALSE)
boxplot(acomp(cdata),Land,log=FALSE)
#boxplot(acomp(cdata),Rock,log=FALSE)

boxplot(acomp(cdata),log=FALSE,ylim=c(0,5))
boxplot(acomp(cdata),Land,log=FALSE,ylim=c(0,5))
#boxplot(acomp(cdata),Rock,log=FALSE,ylim=c(0,5))

qqnorm(acomp(cdata),alpha=100)
qqnorm(acomp(cdata),alpha=0.05)
qqnorm(acomp(cdata[,-3]),alpha=0.05)
qqnorm(acomp(cdata[,-3]),alpha=0.05)
 
#boxplot(acomp(cdata[,-1]),js$Cd)
plot(Land,data.matrix(cdata)%*% rep(1,ncol(cdata)))


# rcomp.plots

boxplot(rcomp(cdata))
boxplot(rcomp(cdata),Land)
#boxplot(rcomp)(cdata,Rock)

qqnorm(rcomp(cdata))
qqnorm(rcomp(cdata),alpha=0.05)
qqnorm(rcomp(cdata[,1:3]),alpha=0.05)
plot(acomp(cdata[,1:3]))
ellipses(mean(acomp(cdata[,1:3])), var(acomp(cdata,1:3)),col="red",r=2)

ellipses(mean(rcomp(cdata[,1:3])), var(rcomp(cdata[,1:3])),col="blue",r=2)


plot(rplus(cdata[,1:2]))
ellipses(rplus(mean(rplus(cdata[,1:2]))), var(rplus(cdata[,1:2])),col="blue",r=2)
ellipses(aplus(mean(aplus(cdata[,1:2]))), var(aplus(cdata[,1:2])),col="red",r=2)

plot(aplus(cdata[,1:2]))
ellipses(aplus(mean(aplus(cdata[,1:2]))), var(aplus(cdata[,1:2])),col="red",r=2)
ellipses(rplus(mean(rplus(cdata[,1:2]))), var(rplus(cdata[,1:2])),col="blue",r=2)


straight(acomp(c(1,1,1)),c(2,1,3))


boxplot(rcomp(cdata[,-1]),cdata[,"Cd"])

# biplot, princomp

biplot(princomp(cdata))
biplot(princomp(acomp(cdata)))
biplot(princomp(rcomp(cdata)))
biplot(princomp(aplus(cdata)),choice=c(2,3))
biplot(princomp(rplus(cdata)),choice=c(2,3))

summary(princomp(cdata))
summary(princomp(acomp(cdata)))
summary(princomp(rcomp(cdata)))
summary(princomp(aplus(cdata)))
summary(princomp(rplus(cdata)))


# The following lines produce loadings, which are known to be defined
#   only up to a sign; to ensure meaningful comparison, they are
#   converted to absolute values
# Human: if you want to see the actual numbers, replace abs() by I()
abs( loadings(princomp(cdata)) )
abs( loadings(princomp(acomp(cdata))) )
abs( loadings(princomp(rcomp(cdata))) )
abs( loadings(princomp(aplus(cdata))) )
abs( loadings(princomp(rplus(cdata))) )


#names
qround( meanCol(cdata[,1:3]) )
qround( mean(acomp(cdata[,1:3])) )
oneOrDataset(c(a=1,b=2,c=3))

# covariance
qround( cov(acomp(cd1),acomp(cd2)) )
qround( cov(rcomp(cd1),rcomp(cd2)) )
qround( cov(aplus(cd1),aplus(cd2)) )
qround( cov(rplus(cd1),rplus(cd2)) )

# tmp <- princov(acomp(cd1,cd2))
# The following line produce loadings and scores, which are known to be defined
#   only up to a sign; to ensure meaningful comparison, they are
#   commented
# Human: if you want to see the actual numbers, uncomment next lines
#tmp

tmp <- princomp(acomp(cdata))
# tmp
plot(tmp)
biplot(tmp)
biplot(tmp,ch=2:3)

#pplot(acomp(cdata))