R Under development (unstable) (2024-02-19 r85946 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. > library('geostatsp') Loading required package: Matrix Loading required package: terra terra 1.7.71 > > # exclude this line to use the RandomFields package > options(useRandomFields = FALSE) > > if(interactive() | Sys.info()['user'] =='patrick') { + n=300 + } else { + n=100 + } > > set.seed(0) > mydat = vect(cbind(seq(0,1,len=n), runif(n)), + atts=data.frame(cov1 = rnorm(n), cov2 = rpois(n, 0.5)) + ) > > > > > trueParamAniso = c( + variance=2^2, range=0.2, shape=2, + nugget=1^2,anisoRatio=4,anisoAngleDegrees=10) > > trueParamIso = c( + variance=2^2, range=0.2, shape=2, + nugget=1^2) > > > mydat$U = geostatsp::RFsimulate(trueParamAniso,mydat)$sim > mydat$Uiso = geostatsp::RFsimulate(trueParamIso,mydat)$sim > > mydat$Y = -3 + 0.5*mydat$cov1 + 0.2*mydat$cov2 + + mydat$U + rnorm(length(mydat), 0, sd=sqrt(trueParamAniso["nugget"])) > > mydat$Yiso = -3 + 0.5*mydat$cov1 + 0.2*mydat$cov2 + + mydat$Uiso + rnorm(length(mydat), 0, sd=sqrt(trueParamIso["nugget"])) > > mydat$Ybc = (mydat$Y*0.5+1)^2 > > mydat$YbcIso = (mydat$Yiso*0.5+1)^2 > > print(range(mydat$Ybc, na.rm=TRUE)) [1] 0.0002490923 8.1045536632 > print(range(mydat$YbcIso, na.rm=TRUE)) [1] 8.595556e-06 8.800440e+00 > > date() [1] "Tue Feb 20 17:38:49 2024" > myres = likfitLgm( + formula=Ybc ~ cov1 + cov2, + data=mydat, + coordinates=mydat, + param=c(range=0.01,nugget=2,shape=2, + anisoAngleDegrees=20, anisoRatio=2, + boxcox=1), + paramToEstimate = c("range","nugget", + "anisoRatio","anisoAngleDegrees", + "boxcox","shape"), + reml=TRUE + ) > date() [1] "Tue Feb 20 17:38:50 2024" > myres$opt$logL m2logL.reml logL.reml 193.65745 -96.82873 > myres$par range shape variance nugget 0.41123688 4.00000000 1.35971060 1.55529409 anisoRatio anisoAngleRadians boxcox (Intercept) 2.92424026 0.02122102 0.19819599 -0.59791543 cov1 cov2 anisoAngleDegrees -0.21847661 0.13996281 1.21587462 > > date() [1] "Tue Feb 20 17:38:50 2024" > myresIso = likfitLgm( + formula=YbcIso ~ cov1 + cov2, + data=mydat, + param=c(shape=2,boxcox=0.4), + paramToEstimate = c( + "range","nugget", + "boxcox","shape"), + reml=TRUE + ) > date() [1] "Tue Feb 20 17:38:51 2024" > myresIso$opt$logL m2logL.reml logL.reml 160.22317 -80.11158 > myresIso$par range shape variance nugget 0.14605288 4.00000000 1.54069991 1.16246249 anisoRatio anisoAngleRadians boxcox (Intercept) 1.00000000 0.00000000 0.22859290 -0.51411657 cov1 cov2 anisoAngleDegrees 0.01856578 -0.20523131 0.00000000 > > myres$summary[,grep("^ci", colnames(myres$summary),invert=TRUE)] estimate stdErr pval Estimated range 0.41123688 NA NA TRUE shape 4.00000000 NA NA TRUE sdSpatial 1.16606629 NA NA TRUE sdNugget 1.24711430 NA NA TRUE anisoRatio 2.92424026 NA NA TRUE anisoAngleRadians 0.02122102 NA NA TRUE boxcox 0.19819599 NA NA TRUE (Intercept) -0.59791543 0.6358818 0.3470671 TRUE cov1 -0.21847661 0.1524077 0.1517146 TRUE cov2 0.13996281 0.1632619 0.3912846 TRUE anisoAngleDegrees 1.21587462 NA NA FALSE > > loglikLgm(formula=Ybc ~ cov1 + cov2, + data=mydat, + param=myres$param + ) minusTwoLogRestrictedLik 193.6575 attr(,"param") range shape variance nugget 0.41123688 4.00000000 1.35971060 1.55529409 anisoRatio anisoAngleRadians boxcox (Intercept) 2.92424026 0.02122102 0.19819599 -0.59791543 cov1 cov2 anisoAngleDegrees -0.21847661 0.13996281 1.21587462 attr(,"totalVarHat") [1] 1 attr(,"betaHat") (Intercept) cov1 cov2 -0.5979154 -0.2184766 0.1399628 attr(,"varBetaHat") 3 x 3 Matrix of class "dsyMatrix" (Intercept) cov1 cov2 (Intercept) 0.404345664 0.002847572 -0.013569374 cov1 0.002847572 0.023228096 -0.003517577 cov2 -0.013569374 -0.003517577 0.026654452 attr(,"reml") [1] TRUE attr(,"jacobian") boxcox -155.0984 attr(,"Ltype") [1] 3 attr(,"Lorig") [1] 348.7559 attr(,"aniso") [1] TRUE attr(,"determinants") [1] 32.57578 4.16511 > > # only estimate variance > myres = likfitLgm( + formula=Ybc ~ cov1 + cov2, + data=mydat, + coordinates=mydat, + param=c(range=0.01,nugget=2,shape=2, + anisoAngleDegrees=20, anisoRatio=2, + boxcox=1), + paramToEstimate = c("variance"), + reml=TRUE + ) > > > pdf("ligfitLgm.pdf") > par(mfrow=c(1,2)) > > myraster = rast(nrows=30,ncols=30,xmin=0,xmax=1,ymin=0,ymax=1) > covEst = matern(myraster, y=c(0.5, 0.5), par=myres$param) > covTrue = matern(myraster, y=c(0.5, 0.5), par=trueParamAniso) > > plot(covEst, main="estimate") > plot(covTrue, main="true") > > dev.off() null device 1 > > if(interactive() | Sys.info()['user'] =='patrick') { + + + + library("geostatsp") + data("swissRain") + swissAltitude = unwrap(swissAltitude) + + sr2 = unwrap(swissRain) + sr2$elev = extract(swissAltitude, sr2) + swissFitAgain = likfitLgm(data=sr2, + formula=rain~ elev, + param=c(range=1000,shape=1,nugget=0.1,boxcox=0.5), + paramToEstimate = c("range","nugget") + ) + swissFitAgain$par + } > > > proc.time() user system elapsed 6.29 0.15 6.45