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 > data("swissRain") > swissRain = unwrap(swissRain) > swissAltitude= unwrap(swissAltitude) > swissLandType = unwrap(swissLandType) > > bob = function(x) { + thepar = x$param + pdf(tempfile("lgm", tmpdir=".", fileext=".pdf")) + plot(x$predict[["predict"]], main= + paste( + paste(names(thepar), thepar, sep="="), + collapse=", "),cex.main=0.3 + ) + dev.off() + } > > > > # specify formula name of raster layer > swissFit = lgm(data=swissRain[1:60,], formula=rain~ CHE_alt, + grid=80, covariates=swissAltitude, + shape=2, fixShape=TRUE, + boxcox=0.5, fixBoxcox=TRUE, + aniso=TRUE) > names(swissFit) [1] "predict" "optim" "betaHat" "varParam" "parameters" [6] "data" "model" "summary" > swissFit$summary[,1:4] estimate stdErr ci0.005 ci0.995 (Intercept) 5.539667532 1.5186552399 1.627870863 9.451464e+00 CHE_alt 0.001035103 0.0005719224 -0.000438072 2.508277e-03 sdNugget 0.968107658 NA 0.641881116 1.460134e+00 sdSpatial 3.151230623 NA 1.493616604 6.648463e+00 range/1000 46.171859918 NA 17.631646783 1.209099e+02 shape 2.000000000 NA NA NA anisoRatio 8.126224277 NA 2.494065801 2.647706e+01 anisoAngleRadians 0.642065353 NA 0.505367946 7.787628e-01 anisoAngleDegrees 36.787634917 NA 28.955450418 4.461982e+01 boxcox 0.500000000 NA NA NA > bob(swissFit) null device 1 > > > # specify formula using name of list element > > swissFitAgain = lgm( + data=swissRain[1:80,], + formula=rain~ elev+land, + grid=80, covariates=list(elev=swissAltitude,land=swissLandType), + shape=1, fixShape=TRUE, + boxcox=0.5, fixBoxcox=TRUE, + aniso=TRUE, verbose=FALSE) > names(swissFitAgain) [1] "predict" "optim" "betaHat" "varParam" "parameters" [6] "data" "model" "summary" > swissFitAgain$summary[,c('estimate','stdErr','Estimated')] estimate stdErr Estimated (Intercept) 5.8267433146 1.3533565416 TRUE elev 0.0006942889 0.0006197647 TRUE landCropland/natural vegetation mosaic 0.0000000000 NA FALSE landMixed forests -0.5529200263 0.3722759763 TRUE landCroplands -0.7490923590 0.4604695432 TRUE landGrasslands -0.7465372004 0.5981679226 TRUE landUrban and built-up -1.1794725318 0.6416584043 TRUE landWater bodies -1.2376473219 0.8826197492 TRUE landDeciduous broadleaf forest 0.3587127247 0.8997706864 TRUE landEvergreen needleleaf forest -1.0077078508 1.3050293552 TRUE landDeciduous needleleaf forest -2.8208590497 1.2235509390 TRUE landPermanent Wetlands -3.7422186144 2.0667075839 TRUE sdNugget 0.8932900194 NA TRUE sdSpatial 2.7000563506 NA TRUE range/1000 53.4871108877 NA TRUE shape 1.0000000000 NA FALSE anisoRatio 6.6745956225 NA TRUE anisoAngleRadians 0.6253994861 NA TRUE anisoAngleDegrees 35.8327510635 NA TRUE boxcox 0.5000000000 NA FALSE > bob(swissFitAgain) null device 1 > > swissFitAgain = lgm( + data=swissRain[1:60,], + formula="rain", + grid=80, covariates=swissAltitude, + shape=1, fixShape=TRUE, + boxcox=0.5, fixBoxcox=TRUE, + aniso=TRUE) > swissFitAgain$summary[,c('estimate','stdErr','Estimated')] estimate stdErr Estimated (Intercept) 5.725580581 1.9988235022 TRUE CHE_alt 0.001067468 0.0005765973 TRUE sdNugget 0.867983062 NA TRUE sdSpatial 3.593905716 NA TRUE range/1000 71.546074504 NA TRUE shape 1.000000000 NA FALSE anisoRatio 7.981526356 NA TRUE anisoAngleRadians 0.634302759 NA TRUE anisoAngleDegrees 36.342871044 NA TRUE boxcox 0.500000000 NA FALSE > bob(swissFitAgain) null device 1 > > > > > swissFitAgain = lgm(data=swissRain, formula="rain", + grid=80, covariates=list(elev=swissAltitude,land=swissLandType), + shape=1, fixShape=TRUE, + boxcox=0.5, fixBoxcox=TRUE, + aniso=TRUE) > swissFitAgain$summary[,c('estimate','stdErr','Estimated')] estimate stdErr Estimated (Intercept) 6.0525478421 1.1726947666 TRUE elev -0.0001028764 0.0005095701 TRUE landCropland/natural vegetation mosaic 0.0000000000 NA FALSE landMixed forests -0.5140076042 0.3557662719 TRUE landGrasslands -0.7632888385 0.5420796265 TRUE landCroplands -0.7377049020 0.4465854844 TRUE landUrban and built-up -1.3199271082 0.6450922388 TRUE landEvergreen needleleaf forest -1.4648339239 0.7622905027 TRUE landWater bodies -1.3266755496 0.8774555256 TRUE landDeciduous needleleaf forest -2.3802142458 0.9147492306 TRUE landDeciduous broadleaf forest 0.1987287296 0.9072897145 TRUE landOpen shrublands -1.2190667425 1.2732587969 TRUE landPermanent Wetlands -4.0866237520 1.7919563808 TRUE sdNugget 0.8691160151 NA TRUE sdSpatial 2.6680347252 NA TRUE range/1000 49.1183059982 NA TRUE shape 1.0000000000 NA FALSE anisoRatio 6.2702892780 NA TRUE anisoAngleRadians 0.6401670343 NA TRUE anisoAngleDegrees 36.6788692503 NA TRUE boxcox 0.5000000000 NA FALSE > bob(swissFitAgain) null device 1 > > > # land type, factor covariate > swissRes2 = lgm(rain ~ elev + factor(land), swissRain, + grid=30, + covariates=list(elev=swissAltitude,land=swissLandType), + boxcox=0.5, fixBoxcox=TRUE, + aniso=TRUE + ) > swissRes2$summary estimate stdErr ci0.005 (Intercept) 6.0525478421 1.1726947666 3.031886298 elev -0.0001028764 0.0005095701 -0.001415442 factor(land)Mixed forests -0.5140076042 0.3557662719 -1.430400793 factor(land)Grasslands -0.7632888385 0.5420796265 -2.159593425 factor(land)Croplands -0.7377049020 0.4465854844 -1.888032879 factor(land)Urban and built-up -1.3199271082 0.6450922388 -2.981574600 factor(land)Evergreen needleleaf forest -1.4648339239 0.7622905027 -3.428364139 factor(land)Water bodies -1.3266755496 0.8774555256 -3.586851205 factor(land)Deciduous needleleaf forest -2.3802142458 0.9147492306 -4.736452119 factor(land)Deciduous broadleaf forest 0.1987287296 0.9072897145 -2.138294704 factor(land)Open shrublands -1.2190667425 1.2732587969 -4.498764063 factor(land)Permanent Wetlands -4.0866237520 1.7919563808 -8.702397508 sdNugget 0.8691160151 NA 0.546757020 sdSpatial 2.6680347252 NA 1.509910320 range/1000 49.1183059982 NA 20.115531653 shape 1.0000000000 NA NA anisoRatio 6.2702892780 NA 2.563712758 anisoAngleRadians 0.6401670343 NA 0.478839868 anisoAngleDegrees 36.6788692503 NA 27.435503505 boxcox 0.5000000000 NA NA ci0.995 ci0.025 (Intercept) 9.073209386 3.754108335 elev 0.001209689 -0.001101615 factor(land)Mixed forests 0.402385584 -1.211296684 factor(land)Grasslands 0.633015748 -1.825745383 factor(land)Croplands 0.412623075 -1.612996367 factor(land)Urban and built-up 0.341720384 -2.584284663 factor(land)Evergreen needleleaf forest 0.498696291 -2.958895855 factor(land)Water bodies 0.933500106 -3.046456778 factor(land)Deciduous needleleaf forest -0.023976372 -4.173089793 factor(land)Deciduous broadleaf forest 2.535752163 -1.579526434 factor(land)Open shrublands 2.060630578 -3.714608127 factor(land)Permanent Wetlands 0.529150004 -7.598793720 sdNugget 1.381532600 0.610829519 sdSpatial 4.714458336 1.730074601 range/1000 119.937569920 24.901820267 shape NA NA anisoRatio 15.335777185 3.174956847 anisoAngleRadians 0.801494201 0.517412225 anisoAngleDegrees 45.922234996 29.645536750 boxcox NA NA ci0.975 ci0.05 (Intercept) 8.3509873495 4.1236366020 elev 0.0008958626 -0.0009410446 factor(land)Mixed forests 0.1832814756 -1.0991910469 factor(land)Grasslands 0.2991677062 -1.6549304783 factor(land)Croplands 0.1375865634 -1.4722726558 factor(land)Urban and built-up -0.0555695534 -2.3810094169 factor(land)Evergreen needleleaf forest 0.0292280072 -2.7186902221 factor(land)Water bodies 0.3931056787 -2.7699614534 factor(land)Deciduous needleleaf forest -0.5873386990 -3.8848428354 factor(land)Deciduous broadleaf forest 1.9769838936 -1.2936300481 factor(land)Open shrublands 1.2764746425 -3.3133910927 factor(land)Permanent Wetlands -0.5744537838 -7.0341297043 sdNugget 1.2366177869 0.6464630825 sdSpatial 4.1145100275 1.8548577902 range/1000 96.8848043347 27.7754603220 shape NA NA anisoRatio 12.3833266148 3.5420477207 anisoAngleRadians 0.7629218439 0.5371479507 anisoAngleDegrees 43.7122017511 30.7763105478 boxcox NA NA ci0.95 ci0.1 (Intercept) 7.9814590822 4.5496790281 elev 0.0007352918 -0.0007559167 factor(land)Mixed forests 0.0711758384 -0.9699404269 factor(land)Grasslands 0.1283528012 -1.4579918325 factor(land)Croplands -0.0031371483 -1.3100272287 factor(land)Urban and built-up -0.2588447995 -2.1466460767 factor(land)Evergreen needleleaf forest -0.2109776257 -2.4417485111 factor(land)Water bodies 0.1166103542 -2.4511800522 factor(land)Deciduous needleleaf forest -0.8755856561 -3.5525125543 factor(land)Deciduous broadleaf forest 1.6910875072 -0.9640098245 factor(land)Open shrublands 0.8752576077 -2.8508135470 factor(land)Permanent Wetlands -1.1391177998 -6.3831082572 sdNugget 1.1684544226 0.6901337903 sdSpatial 3.8377116200 2.0099352869 range/1000 86.8611341148 31.5025141130 shape NA NA anisoRatio 11.0999429508 4.0182595246 anisoAngleRadians 0.7431861180 0.5599019787 anisoAngleDegrees 42.5814279528 32.0800203181 boxcox NA NA ci0.9 pval Estimated (Intercept) 7.5554166561 2.453321e-07 TRUE elev 0.0005501639 8.400038e-01 TRUE factor(land)Mixed forests -0.0580747815 1.485168e-01 TRUE factor(land)Grasslands -0.0685858445 1.591088e-01 TRUE factor(land)Croplands -0.1653825754 9.855935e-02 TRUE factor(land)Urban and built-up -0.4932081396 4.074594e-02 TRUE factor(land)Evergreen needleleaf forest -0.4879193367 5.465338e-02 TRUE factor(land)Water bodies -0.2021710470 1.305446e-01 TRUE factor(land)Deciduous needleleaf forest -1.2079159373 9.267100e-03 TRUE factor(land)Deciduous broadleaf forest 1.3614672836 8.266223e-01 TRUE factor(land)Open shrublands 0.4126800621 3.383461e-01 TRUE factor(land)Permanent Wetlands -1.7901392468 2.257583e-02 TRUE sdNugget 1.0945162495 NA TRUE sdSpatial 3.5416111857 NA TRUE range/1000 76.5846171985 NA TRUE shape NA NA FALSE anisoRatio 9.7844669785 NA TRUE anisoAngleRadians 0.7204320900 NA TRUE anisoAngleDegrees 41.2777181825 NA TRUE boxcox NA NA FALSE > bob(swissRes2) null device 1 > > > > # simulated data (without a CRS) > # and all covariates are in 'data' object > myModel = c(intercept=0,variance=2^2,nugget=0.5^2, range=4.5,shape=2, + cov1=0.2, cov2=-0.5) > covariates = rast( + xmin=0,ymin=0,xmax=10,ymax=10, + ncols=200,nrows=200,nlyrs=2, crs='+proj=merc') > values(covariates) = c( + rep(seq(0,1,len=nrow(covariates)), ncol(covariates)), + rep(seq(0,1,len=nrow(covariates)), + rep(nrow(covariates), ncol(covariates)))) > names(covariates) = c("cov1","cov2") > > Npoints = 30 > set.seed(0) > myPoints = vect(cbind(runif(Npoints,0,10), + seq(0,10, len=Npoints)), crs = '+proj=merc') > > # don't test using the randomFields package, it's currently broken on some R builds > options(useRandomFields = FALSE) > > myPoints = RFsimulate(myModel,myPoints) > > values(myPoints) = cbind( + values(myPoints), + as.data.frame(extract(covariates, myPoints)) + ) > > myPoints$y= myModel["intercept"] + + as.matrix(values(myPoints)[,names(covariates)]) %*% + myModel[names(covariates)] + + myPoints$sim+ + rnorm(length(myPoints), 0, sqrt(myModel["nugget"])) > > fitLikfit = likfitLgm(y~cov1+cov2, myPoints, + param=c(range=1,nugget=0,shape=1)) > > > > > # run lgm without providing covariates > fitMLE = lgm( + formula=y~ cov1+cov2, + data=myPoints, + grid=10, covariates=list(), + shape=1, fixShape=TRUE) > > c(fitMLE$summary["range","estimate"], fitLikfit$summary["range","estimate"]) [1] 4.712256 4.712273 > bob(fitMLE) null device 1 > > > # now give covariates as raster brick > fitMLE = lgm( y~ cov1 + cov2, myPoints, grid=10, + covariates=covariates, + shape=1, fixShape=TRUE) > c(fitMLE$summary["range","estimate"], fitLikfit$summary["range","estimate"]) [1] 4.712256 4.712273 > bob(fitMLE) null device 1 > # now give covariates as list > fitMLE = lgm(y~ cov1 + cov2, myPoints, grid=10, + covariates=list(cov1=covariates[["cov1"]], + cov2 = covariates[["cov2"]]), + shape=1, fixShape=TRUE) > c(fitMLE$summary["range","estimate"], fitLikfit$summary["range","estimate"]) [1] 4.712256 4.712273 > bob(fitMLE) null device 1 > > # not remove covariates from data > myPoints = vect(crds(myPoints), + atts=values(myPoints)[,"y",drop=FALSE], + crs = crs(myPoints)) > > # now give covariates as raster brick > fitMLE = lgm(y~ cov1 + cov2, myPoints, grid=10, + covariates=covariates, + shape=1, fixShape=TRUE) > c(fitMLE$summary["range","estimate"], fitLikfit$summary["range","estimate"]) [1] 4.722069 4.712273 > bob(fitMLE) null device 1 > > > > proc.time() user system elapsed 17.78 0.29 18.07