R Under development (unstable) (2025-08-19 r88650 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 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.244298024731204, -0.170258948491207, 0.0173105473447226, -0.0685886765359066, -0.149589443472457, -0.30393153180931, -0.0744849494750161, -0.0610089646585014, -0.169877422812463, -0.0250603400774691), rsph2 = 2.68494135241001, sPSD2 = c(0.188093425588121, 0.112413776792952, -0.261555808943098, -0.0901237569697608, 0.29454004199496, -0.0495893037098249, 0.11692412883546, -0.0364983638756834, 0.0556329418710303, -0.0560865939820887), rexp3 = 0.324821592199274, sPSD3 = c(0.807182513235476, 0.0670977945920792, -0.345098470157714, -0.182438480010288, -0.159981183558307, -0.409605615460391, -0.233274243198665, -0.0255766856936466, -0.109462834722641, -0.23959455345376 )) 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.93263 attr(,"fit.quality")$estimate [1] -0.24429802 -0.17025895 0.01731055 -0.06858868 -0.14958944 -0.30393153 [7] -0.07448495 -0.06100896 -0.16987742 -0.02506034 2.68494135 0.18809343 [13] 0.11241378 -0.26155581 -0.09012376 0.29454004 -0.04958930 0.11692413 [19] -0.03649836 0.05563294 -0.05608659 0.32482159 0.80718251 0.06709779 [25] -0.34509847 -0.18243848 -0.15998118 -0.40960562 -0.23327424 -0.02557669 [31] -0.10946283 -0.23959455 attr(,"fit.quality")$gradient [1] -0.25207498 -1.05596004 0.66026703 0.90129715 1.78771442 3.26354022 [7] 0.63006026 -1.55926098 -4.37817149 -0.62597682 0.09655334 0.37780685 [13] 0.07140695 -1.13281622 0.45469188 -1.38118128 0.33245011 -2.15605071 [19] 2.61659128 0.21376763 -0.66840362 0.57441820 -0.43695040 3.16570520 [25] -0.66597325 -3.04732342 2.22245924 2.12128084 -0.42278601 -0.78979618 [31] -4.28758105 1.74733349 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 12.31 3.43 15.68