R Under development (unstable) (2024-01-24 r85824 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # source("acheck.R") > > # try(library(compositions,lib.loc="compositions.Rcheck")) > # options(error=dump.frames) > options(warn=1) > library(compositions) Welcome to compositions, a package for compositional data analysis. Find an intro with "? compositions" Attaching package: 'compositions' The following objects are masked from 'package:stats': anova, cor, cov, dist, var The following object is masked from 'package:graphics': segments The following objects are masked from 'package:base': %*%, norm, scale, scale.default > data(juraset) > > X <- with(juraset,cbind(X,Y)) > comp <- acomp(juraset,c("Cd","Cu","Pb","Co","Cr")) > > > lrv <- logratioVariogram(comp,X,maxdist=1,nbins=10) > > plot(lrv) > > > > > fff <- CompLinModCoReg(~nugget()+sph(0.5)+exp(0.7),comp) > #fff <- CompLinModCoReg(~nugget()+sph(0.5),comp) > fit <- fit.lmc(lrv,fff,iterlim=100,print.level=0) > fit function (h, sPSD1 = c(-0.245450591724301, -0.168154709789692, 0.01860280635602, -0.0736817867554658, -0.149655578832826, -0.304594956147885, -0.0757401634322841, -0.0589150432944755, -0.171222228453802, -0.0302944166192619), rsph2 = 2.68460626319117, sPSD2 = c(0.184843199892676, 0.126729855384512, -0.257181209661489, -0.096557226789793, 0.297312881307108, -0.0522204235745241, 0.110863884787551, -0.0351320064170399, 0.0589031354367264, -0.0631964050281496), rexp3 = 0.325097298627938, sPSD3 = c(0.807249701794204, 0.066498057570053, -0.345552850541565, -0.183821683117109, -0.159973437372879, -0.406283794937714, -0.233153261339591, -0.0246234469526221, -0.108237899753436, -0.239059435120841)) vgram.nugget(h = h) %o% parametricPosdefClrMat(sPSD1) + vgram.sph(h = h, range = rsph2) %o% parametricPosdefClrMat(sPSD2) + vgram.exp(h = h, range = rexp3) %o% parametricPosdefClrMat(sPSD3) attr(,"fit.quality") attr(,"fit.quality")$minimum [1] 59.9583 attr(,"fit.quality")$estimate [1] -0.24545059 -0.16815471 0.01860281 -0.07368179 -0.14965558 -0.30459496 [7] -0.07574016 -0.05891504 -0.17122223 -0.03029442 2.68460626 0.18484320 [13] 0.12672986 -0.25718121 -0.09655723 0.29731288 -0.05222042 0.11086388 [19] -0.03513201 0.05890314 -0.06319641 0.32509730 0.80724970 0.06649806 [25] -0.34555285 -0.18382168 -0.15997344 -0.40628379 -0.23315326 -0.02462345 [31] -0.10823790 -0.23905944 attr(,"fit.quality")$gradient [1] -2.24950302 -0.17323085 -0.78284982 -1.09509236 1.74475646 3.19613633 [7] -1.54259745 -1.35325590 -5.13244142 -1.29018736 -0.03007521 0.54477530 [13] 0.22471518 -0.96319219 -0.42388764 0.04004963 -0.16811607 -1.27233890 [19] 0.92200550 0.77790909 -1.06857870 0.02421064 2.86511550 -0.92123032 [25] -1.05432882 -0.79887797 -0.15013512 1.66693021 -0.42764547 -0.80618602 [31] -2.46352216 -0.58604568 attr(,"fit.quality")$code [1] 4 attr(,"fit.quality")$iterations [1] 100 attr(,"class") [1] "CompLinModCoReg" "function" > fff(1:3) , , 1 1 2 3 4 5 [1,] 2.388922 -0.5972304 -0.5972304 -0.5972304 -0.5972304 [2,] 2.399847 -0.5999616 -0.5999616 -0.5999616 -0.5999616 [3,] 2.399998 -0.5999995 -0.5999995 -0.5999995 -0.5999995 , , 2 1 2 3 4 5 [1,] -0.5972304 2.388922 -0.5972304 -0.5972304 -0.5972304 [2,] -0.5999616 2.399847 -0.5999616 -0.5999616 -0.5999616 [3,] -0.5999995 2.399998 -0.5999995 -0.5999995 -0.5999995 , , 3 1 2 3 4 5 [1,] -0.5972304 -0.5972304 2.388922 -0.5972304 -0.5972304 [2,] -0.5999616 -0.5999616 2.399847 -0.5999616 -0.5999616 [3,] -0.5999995 -0.5999995 2.399998 -0.5999995 -0.5999995 , , 4 1 2 3 4 5 [1,] -0.5972304 -0.5972304 -0.5972304 2.388922 -0.5972304 [2,] -0.5999616 -0.5999616 -0.5999616 2.399847 -0.5999616 [3,] -0.5999995 -0.5999995 -0.5999995 2.399998 -0.5999995 , , 5 1 2 3 4 5 [1,] -0.5972304 -0.5972304 -0.5972304 -0.5972304 2.388922 [2,] -0.5999616 -0.5999616 -0.5999616 -0.5999616 2.399847 [3,] -0.5999995 -0.5999995 -0.5999995 -0.5999995 2.399998 > plot(lrv,lrvg=vgram2lrvgram(fit)) > > x <- (0:60/60)*6 > y <- (0:40/40)*6 > Xnew <- cbind(rep(x,length(y)),rep(y,each=length(x))) > > > # Grid example > erg <- compOKriging(comp,X,Xnew,fit,err=FALSE) > image(x,y,structure(unclass(balance(erg$Z,~Cd/Cu)),dim=c(length(x),length(y))), asp=1) > par(mar=c(0,0,1,0)) > suppressWarnings( + pairwisePlot(erg$Z, + panel=function(a,b,xlab,ylab){ + image(x,y, asp=1, structure(log(a/b),dim=c(length(x),length(y))),main=paste("log(",xlab,"/",ylab,")",sep="") + ) + points(X,pch=".") + } + ) + ) > > #erg <- compOKriging(comp,X,Xnew,fit$vg,err=TRUE) > > > # Checking > ergR <- compOKriging(comp,X,X+1E-7,fit,err=FALSE) > plot(ergR$Z-comp) > ergR <- compOKriging(comp,X,X,fit,err=FALSE) > plot(ergR$Z-comp) > > pairwisePlot(ilr(comp),ilr(ergR$Z)) Warning in .S3methods(generic.function, class, envir, all.names = all.names, : generic function 'panel' dispatches methods for generic 'plot' > > ergR <- compOKriging(comp,X,X[rev(1:31),],fit,err=FALSE) > pairwisePlot(ilr(comp)[rev(1:31),],ilr(ergR$Z)) Warning in .S3methods(generic.function, class, envir, all.names = all.names, : generic function 'panel' dispatches methods for generic 'plot' > > plot(ergR$Z) > plot(comp) > # > > # Small example > X <- rbind(c(0,1),c(1,0),c(1,1),c(0,0)) > Z <- acomp(rbind(c(2,1,1),c(1,2,1),c(1,1,2),c(1,1,1))) > Xnew <- rbind(c(0,1),c(0.5,0.5)) > > vg <- CompLinModCoReg(~nugget()+sph(3),Z) > vg function (h, sPSD1 = c(1, 0, 1), rsph2 = 3, sPSD2 = c(1, 0, 1 )) vgram.nugget(h = h) %o% parametricPosdefClrMat(sPSD1) + vgram.sph(h = h, range = rsph2) %o% parametricPosdefClrMat(sPSD2) attr(,"class") [1] "CompLinModCoReg" "function" > (ergR <- compOKriging(Z,X,X,vg,err=FALSE)) $X [,1] [,2] [1,] 0 1 [2,] 1 0 [3,] 1 1 [4,] 0 0 $Z [,1] [,2] [,3] [1,] "0.5000000" "0.2500000" "0.2500000" [2,] "0.2500000" "0.5000000" "0.2500000" [3,] "0.2500000" "0.2500000" "0.5000000" [4,] "0.3333333" "0.3333333" "0.3333333" attr(,"class") [1] "acomp" $err NULL > > ## Explicit vgmodel; check the following functions prefixed with compositions::: > #vgmodel <- function(h,nugget=0.2*parameterPosDefMat(diag(5)),sill=0.6*parameterPosDefMat(diag(5)),range1=1) { > # (h>0) %o% parametricPosdefMat(nugget)+ vgram.sph(h,range=range1)%o%parametricPosdefMat(sill) > #} > # plot(lrv,lrvg=vgram2lrvgram(vgmodel)) > > #fit <- vgmFit(lrv,vgmodel) > #fit <- vgmFit(lrv,vgmodel,mode="ls") > #fit > #vgmGof(emp=lrv,vg=vgmodel,mode="ls") > #vgmGof(emp=lrv,vg=vgmodel,mode="log") > #lrv$h > #lrv$vg > #vgmGof(emp=lrv,vg=fit$vg,mode="ls") > #vgmGof(emp=lrv,vg=fit$vg,mode="log") > > #plot(lrv,lrvg=vgram2lrvgram(fit$vg)) > > > > > proc.time() user system elapsed 7.68 1.35 9.03