CRAN Package Check Results for Package geostatsp

Last updated on 2026-02-28 02:54:17 CET.

Flavor	Version	T_install	T_check	T_total	Status
r-devel-linux-x86_64-debian-clang	2.0.10	39.88	558.83	598.71	OK
r-devel-linux-x86_64-debian-gcc	2.0.11	24.89	368.58	393.47	OK
r-devel-linux-x86_64-fedora-clang	2.0.8	77.00	475.34	552.34	ERROR
r-devel-linux-x86_64-fedora-gcc	2.0.11	73.00	784.92	857.92	ERROR
r-devel-macos-arm64	2.0.11	9.00	188.00	197.00	ERROR
r-devel-windows-x86_64	2.0.11	40.00	430.00	470.00	OK
r-patched-linux-x86_64	2.0.8	38.37	550.70	589.07	OK
r-release-linux-x86_64	2.0.8	36.53	557.60	594.13	OK
r-release-macos-arm64	2.0.11	9.00	122.00	131.00	ERROR
r-release-macos-x86_64	2.0.11	33.00	1173.00	1206.00	OK
r-release-windows-x86_64	2.0.11	40.00	429.00	469.00	OK
r-oldrel-macos-arm64	2.0.11	8.00	148.00	156.00	NOTE
r-oldrel-macos-x86_64	2.0.11	31.00	1240.00	1271.00	NOTE
r-oldrel-windows-x86_64	2.0.11	56.00	505.00	561.00	NOTE

Check Details

Version: 2.0.8
Check: examples
Result: ERROR Running examples in ‘geostatsp-Ex.R’ failed The error most likely occurred in: > ### Name: RFsimulate > ### Title: Simulation of Random Fields > ### Aliases: RFsimulate modelRandomFields RFsimulate RFsimulate-methods > ### RFsimulate,ANY,SpatRaster-method RFsimulate,numeric,SpatRaster-method > ### RFsimulate,numeric,SpatVector-method > ### RFsimulate,RMmodel,SpatVector-method > ### RFsimulate,RMmodel,SpatRaster-method > ### RFsimulate,matrix,SpatRaster-method > ### RFsimulate,matrix,SpatVector-method RFsimulate,data.frame,ANY-method > ### Keywords: spatial > > ### ** Examples > > library('geostatsp') > > # exclude this line to use the RandomFields package > options(useRandomFields = FALSE) > > model1 <- c(var=5, range=1,shape=0.5) > > > myraster = rast(nrows=20,ncols=30,extent = ext(0,6,0,4), + crs="+proj=utm +zone=17 +datum=NAD27 +units=m +no_defs") > > set.seed(0) > > simu <- RFsimulate(model1, x=myraster, n=3) install the RandomFields package for faster simulations OMP: Warning #96: Cannot form a team with 24 threads, using 2 instead. OMP: Hint Consider unsetting KMP_DEVICE_THREAD_LIMIT (KMP_ALL_THREADS), KMP_TEAMS_THREAD_LIMIT, and OMP_THREAD_LIMIT (if any are set). Flavor: r-devel-linux-x86_64-fedora-clang

Version: 2.0.8
Check: tests
Result: ERROR Running ‘RFsimulate.R’ [90m/59m] Running ‘krige.R’ [18s/19s] Running ‘lgcp.R’ [19s/20s] Running ‘lgm.R’ [70s/59s] Running ‘lgmRaster.R’ [0m/30m] Running ‘likfitLgm.R’ Running the tests in ‘tests/RFsimulate.R’ failed. Complete output: > library("geostatsp") Loading required package: Matrix Loading required package: terra terra 1.8.93 > > model <- c(var=5, range=20,shape=0.5) > > # any old crs > theCrs = "+proj=utm +zone=17 +datum=NAD27 +units=m +no_defs" > > # don't test using the randomFields package, it's currently broken on some R builds > options(useRandomFields = FALSE) > > myraster = rast(nrows=20,ncols=20,extent = ext(100,110,100,110), + crs=theCrs) > > set.seed(0) > simu = RFsimulate(model = rbind(a=model, b=model+0.1), + x=myraster, n=4 + ) OMP: Warning #96: Cannot form a team with 24 threads, using 2 instead. OMP: Hint Consider unsetting KMP_DEVICE_THREAD_LIMIT (KMP_ALL_THREADS), KMP_TEAMS_THREAD_LIMIT, and OMP_THREAD_LIMIT (if any are set). Running the tests in ‘tests/lgmRaster.R’ failed. Complete output: > #+ setup > library('geostatsp') Loading required package: Matrix Loading required package: terra terra 1.8.93 > #' > > #' # simulated data > > # exclude this line to use the RandomFields package > options(useRandomFields = FALSE) > > Ncell = 40 > > myRaster = squareRaster(ext(0,6000,0,6000), Ncell) > > myParam=c(oneminusar=0.1, conditionalVariance=2.5^2,shape=2) > myQ = maternGmrfPrec(myRaster, param=myParam) > attributes(myQ)$info$optimalShape shape variance range cellSize 4.092496 110.524266 900.000000 150.000000 > set.seed(0) > mySim = RFsimulate(attributes(myQ)$info$optimalShape, myRaster) install the RandomFields package for faster simulations OMP: Warning #96: Cannot form a team with 24 threads, using 2 instead. OMP: Hint Consider unsetting KMP_DEVICE_THREAD_LIMIT (KMP_ALL_THREADS), KMP_TEAMS_THREAD_LIMIT, and OMP_THREAD_LIMIT (if any are set). Flavor: r-devel-linux-x86_64-fedora-clang

Version: 2.0.11
Check: examples
Result: ERROR Running examples in ‘geostatsp-Ex.R’ failed The error most likely occurred in: > ### Name: RFsimulate > ### Title: Simulation of Random Fields > ### Aliases: RFsimulate modelRandomFields RFsimulate RFsimulate-methods > ### RFsimulate,ANY,SpatRaster-method RFsimulate,numeric,SpatRaster-method > ### RFsimulate,numeric,SpatVector-method > ### RFsimulate,RMmodel,SpatVector-method > ### RFsimulate,RMmodel,SpatRaster-method > ### RFsimulate,matrix,SpatRaster-method > ### RFsimulate,matrix,SpatVector-method RFsimulate,data.frame,ANY-method > ### Keywords: spatial > > ### ** Examples > > library('geostatsp') > > # exclude this line to use the RandomFields package > options(useRandomFields = FALSE) > > model1 <- c(var=5, range=1,shape=0.5) > > > myraster = rast(nrows=20,ncols=30,extent = ext(0,6,0,4), + crs="+proj=utm +zone=17 +datum=NAD27 +units=m +no_defs") > > set.seed(0) > > simu <- RFsimulate(model1, x=myraster, n=3) install the RandomFields package for faster simulations Flavor: r-devel-linux-x86_64-fedora-gcc

Version: 2.0.11
Check: tests
Result: ERROR Running ‘RFsimulate.R’ [90m/46m] Running ‘krige.R’ [21s/22s] Running ‘lgcp.R’ [68s/56s] Running ‘lgm.R’ [63s/61s] Running ‘lgmRaster.R’ [0m/42m] Running ‘likfitLgm.R’ Running the tests in ‘tests/RFsimulate.R’ failed. Complete output: > library("geostatsp") Loading required package: Matrix Loading required package: terra terra 1.8.93 > > Sys.setenv( + OMP_NUM_THREADS = "2", + OPENBLAS_NUM_THREADS = "2", + MKL_NUM_THREADS = "2", + BLIS_NUM_THREADS = "2", + VECLIB_MAXIMUM_THREADS = "2", # macOS Accelerate + NUMEXPR_NUM_THREADS = "2" + ) > options(mc.cores = 2) > > model <- c(var=5, range=20,shape=0.5) > > # any old crs > theCrs = "+proj=utm +zone=17 +datum=NAD27 +units=m +no_defs" > > # don't test using the randomFields package, it's currently broken on some R builds > options(useRandomFields = FALSE) > > myraster = rast(nrows=20,ncols=20,extent = ext(100,110,100,110), + crs=theCrs) > > set.seed(0) > simu = RFsimulate(model = rbind(a=model, b=model+0.1), + x=myraster, n=4 + ) Running the tests in ‘tests/lgmRaster.R’ failed. Complete output: > #+ setup > library('geostatsp') Loading required package: Matrix Loading required package: terra terra 1.8.93 > #' > > #' # simulated data > > # exclude this line to use the RandomFields package > options(useRandomFields = FALSE) > > Ncell = 40 > > myRaster = squareRaster(ext(0,6000,0,6000), Ncell) > > myParam=c(oneminusar=0.1, conditionalVariance=2.5^2,shape=2) > myQ = maternGmrfPrec(myRaster, param=myParam) > attributes(myQ)$info$optimalShape shape variance range cellSize 4.092496 110.524266 900.000000 150.000000 > set.seed(0) > mySim = RFsimulate(attributes(myQ)$info$optimalShape, myRaster) install the RandomFields package for faster simulations Flavor: r-devel-linux-x86_64-fedora-gcc

Version: 2.0.11
Check: tests
Result: ERROR Running ‘RFsimulate.R’ [3s/3s] Running ‘krige.R’ [3s/3s] Running ‘lgcp.R’ [12s/12s] Running ‘lgm.R’ [7s/7s] Running ‘lgmRaster.R’ [38s/21s] Running ‘likfitLgm.R’ [2s/2s] Running ‘matern.R’ [2s/2s] Running ‘maternGmrfPrec.R’ [3s/3s] Running ‘profLlgm.R’ [2s/2s] Running the tests in ‘tests/profLlgm.R’ failed. Complete output: > > library('geostatsp') Loading required package: Matrix Loading required package: terra terra 1.8.93 > data('swissRain') > swissRain = unwrap(swissRain) > swissAltitude = unwrap(swissAltitude) > > Ncores = c(1,2)[1+(.Platform$OS.type=='unix')] > > > > sr2 = swissRain > sr2$elev = extract(swissAltitude, sr2) Warning message: [`[[<-`] only using the first column > swissFit = likfitLgm( + data=sr2, + formula=rain~ elev, + param=c(range=10000,shape=1,nugget=0,boxcox=0.5,anisoRatio=2,anisoAngleDegrees=45), + paramToEstimate = c("range",'anisoAngleDegrees','anisoRatio'), + reml=FALSE, + verbose=FALSE + ) > > > # calculate log-likelihood at the MLE's, but re-estimate variance > sl = loglikLgm( + swissFit$param[c('range','shape','boxcox', 'anisoRatio', 'anisoAngleRadians')], + data=sr2, + formula=rain~ elev, + reml=swissFit$model$reml) > > > # calculate log-likelihood without re-estimating variance > sigSqHat = attributes(sl)$totalVarHat > sl1 = loglikLgm( + c(attributes(sl)$param[ + c('boxcox','anisoRatio','anisoAngleRadians','shape', 'range')], + variance=sigSqHat), + data=sr2, + formula=rain~ elev, + reml=swissFit$model$reml) > > > # re=estimate the anisotropy parameters but not the range > sf2 = likfitLgm( + data=swissFit$data, + formula=swissFit$model$formula, + param= swissFit$param[c('range','nugget','shape','boxcox', 'anisoRatio', 'anisoAngleRadians')], + paramToEstimate = c('variance','anisoAngleRadians','anisoRatio'), + reml=swissFit$model$reml) > > # these should all be the same > as.numeric(sl1) [1] 644.4812 > as.numeric(sl) [1] 644.4812 > swissFit$optim$logL m2logL.ml logL.ml 644.4812 -322.2406 > sf2$optim$logL m2logL.ml logL.ml 644.4812 -322.2406 > > date() [1] "Sat Feb 28 10:56:47 2026" > x=profLlgm(swissFit, mc.cores=Ncores, + range=seq(15000, 55000 , len=12) + ) *** caught segfault *** address 0x110, cause 'invalid permissions' *** caught segfault *** address 0x110, cause 'invalid permissions' Traceback: 1: (function (formula, data, paramToEstimate = c("range", "nugget"), reml = TRUE, coordinates = data, param = NULL, upper = NULL, lower = NULL, parscale = NULL, verbose = FALSE) { param = param[!is.na(param)] coordinatesOrig = coordinates aniso = as.logical(length(grep("^aniso", c(names(param), paramToEstimate)))) | any(abs(param["anisoRatio"] - 1) > 1e-05, na.rm = TRUE) if (aniso) { if (is.matrix(coordinates)) { if (ncol(coordinates) != 2 | nrow(coordinates) != nrow(data)) stop("anisotropic model requested but coordinates appears to be a distance matrix") coordinates = vect(coordinates) } maxDist = dist(matrix(ext(coordinates), ncol = 2)) } else { if (is.matrix(coordinates)) { if (ncol(coordinates) == 2) { coordinates = dist(coordinates) } else { coordinates = as(coordinates, "dsyMatrix") } } if (any(class(coordinates) == "dist")) coordinates = as(as.matrix(coordinates), "dsyMatrix") if (length(grep("SpatVect", class(coordinates)))) { if (!nchar(crs(coordinates))) terra::crs(coordinates) = "+proj=utm +zone=1" coordinates = new("dsyMatrix", Dim = rep(length(coordinates), 2), x = as.vector(as.matrix(terra::distance(coordinates))), uplo = "L") } maxDist = max(coordinates, na.rm = TRUE) } trend = formula theFactors = NULL if (any(class(trend) == "formula")) { data = data.frame(data) theNA = apply(data[, all.vars(trend), drop = FALSE], 1, function(qq) any(is.na(qq))) noNA = !theNA theFactors = model.frame(trend, data[noNA, ]) whichFactors = unlist(lapply(theFactors, is.factor)) theFactors = theFactors[, whichFactors, drop = FALSE] theFactors = lapply(theFactors, levels) theFactors = unlist(lapply(theFactors, function(qq) qq[1])) covariates = model.matrix(trend, data[noNA, ]) covariates = covariates[, grep(paste("^(", paste(names(theFactors), collapse = "|"), ")NA$", sep = ""), colnames(covariates), invert = TRUE), drop = FALSE] observations = all.vars(trend)[1] if (!any(names(data) == observations)) warning("can't find observations ", observations, "in data") observations = data[noNA, observations, drop = FALSE] } else { trend = as.matrix(trend) theNA = is.na(data) | apply(trend, 1, function(qq) any(is.na(qq))) noNA = !theNA observations = data[noNA] covariates = trend[noNA, , drop = FALSE] } if (any(theNA)) { if (length(grep("^SpatVector", class(coordinates)))) { coordinates = coordinates[noNA] } else { if (ncol(coordinates) == nrow(coordinates)) { coordinates = coordinates[noNA, noNA] } else { coordinates = coordinates[noNA, ] } } } estimateVariance = TRUE if (any(paramToEstimate == "variance")) { paramToEstimate = paramToEstimate[paramToEstimate != "variance"] param = param[names(param) != "variance"] } else { if (any(names(param) == "variance")) { estimateVariance = FALSE } } paramToEstimate = gsub("anisoAngleDegrees", "anisoAngleRadians", paramToEstimate) degToRad = function(par) { if (length(grep("anisoAngleDegrees", names(par)))) { if (!length(grep("anisoAngleRadians", names(par)))) par["anisoAngleRadians"] = 2 * pi * par["anisoAngleDegrees"]/360 par = par[grep("anisoAngleDegrees", names(par), invert = TRUE)] } par } param = degToRad(param) lower = degToRad(lower) upper = degToRad(upper) parscale = degToRad(parscale) lowerDefaults = c(nugget = 0, range = maxDist/1000, anisoRatio = 0.01, anisoAngleRadians = -pi/2, shape = 0.1, boxcox = -1.5, variance = 0) upperDefaults = c(nugget = Inf, range = 10 * maxDist, anisoRatio = 100, anisoAngleRadians = pi/2, shape = 4, boxcox = 2.5, variance = Inf) paramDefaults = c(nugget = 0, anisoRatio = 1, anisoAngleRadians = 0, shape = 1.5, boxcox = 1, range = maxDist/10) if (any(names(paramToEstimate) == "nugget")) { paramDefaults["nugget"] = 1 } parscaleDefaults = c(range = maxDist/5, nugget = 0.05, boxcox = 0.5, anisoAngleRadians = 0.2, anisoRatio = 1, variance = 1, shape = 0.2) ndepsDefault = c(range = 0.01, nugget = 0.05, boxcox = 0.005, anisoAngleRadians = 0.01, anisoRatio = 0.01, variance = 0.01, shape = 0.01) paramDefaults[names(param)] = param parscaleDefaults[names(parscale)] = parscale lowerDefaults[names(lower)] = lower upperDefaults[names(upper)] = upper if (any(paramToEstimate == "nugget") & paramDefaults["nugget"] == lowerDefaults["nugget"]) { paramDefaults["nugget"] = min(c(0.5, upperDefaults["nugget"])) } startingParam = paramDefaults[paramToEstimate] names(startingParam) = paramToEstimate naStarting = is.na(startingParam) startingParam[naStarting] = paramDefaults[names(startingParam)[naStarting]] moreParams = paramDefaults[!names(paramDefaults) %in% paramToEstimate] allParams = c(startingParam, moreParams) allParams = fillParam(allParams) paramsForC = allParams[c("nugget", "variance", "range", "shape", "anisoRatio", "anisoAngleRadians", "boxcox")] Sparam = names(paramsForC) %in% paramToEstimate names(Sparam) = names(paramsForC) paramToEstimate = names(Sparam)[Sparam] parOptions = cbind(lower = lowerDefaults[paramToEstimate], upper = upperDefaults[paramToEstimate], parscale = parscaleDefaults[paramToEstimate], ndeps = ndepsDefault[paramToEstimate]) forO = list(scalarF = c(fnscale = -1, abstol = -1, reltol = -1, alpha = -1, beta = -1, gamma = -1, factr = 1e+06, pgtol = 0), scalarInt = c(trace = 0, maxit = 200, REPORT = 1, type = -1, lmm = 25, tmax = -1, temp = -1), pars = parOptions[, c("lower", "upper", "parscale", "ndeps"), drop = FALSE]) if (verbose) { forO$scalarInt["trace"] = 6 forO$scalarInt["REPORT"] = 200 } forO$parsInt = rep(0, nrow(forO$pars)) names(forO$parsInt) = rownames(forO$pars) forO$parsInt[forO$pars[, "lower"] != -Inf & forO$pars[, "upper"] != Inf] = 2 forO$parsInt[forO$pars[, "lower"] != -Inf & forO$pars[, "upper"] == Inf] = 1 forO$parsInt[forO$pars[, "lower"] == -Inf & forO$pars[, "upper"] != Inf] = 3 forOparsOrig = forO$pars forO$pars[!is.finite(forO$pars)] = -1 forO$pars = c(forO$pars, rep(0, ncol(covariates) + ncol(covariates)^2)) if (aniso) { xcoord = crds(coordinates)[, 1] ycoord = crds(coordinates)[, 2] } else { xcoord = as.vector(coordinates) ycoord = -99 } obsCov = cbind(y1 = observations, y2 = 0, y3 = 0, covariates) if (all(paramToEstimate == "variance") & param["nugget"] == 0) { theL = loglikLgm(param, data = observations, formula = covariates, coordinates = coordinates, reml = reml) Traceback: 1: (function (formula, data, paramToEstimate = c("range", "nugget"), reml = TRUE, coordinates = data, param = NULL, upper = NULL, lower = NULL, parscale = NULL, verbose = FALSE) { param = param[!is.na(param)] coordinatesOrig = coordinates aniso = as.logical(length(grep("^aniso", c(names(param), paramToEstimate)))) | any(abs(param["anisoRatio"] - 1) > 1e-05, na.rm = TRUE) if (aniso) { if (is.matrix(coordinates)) { fromOptim = attributes(theL) result = list(optim = list(mle = fillParam(fromOptim$param), logL = c(m2logL = as.numeric(theL), logL = -as.numeric(theL)/2), totalVarHat = fromOptim$totalVarHat, message = "numerical optimization not needed", options = NULL, detail = NULL), betaHat = fromOptim$betaHat, varBetaHat = fromOptim$varBetaHat) } if (ncol(coordinates) != 2 | nrow(coordinates) != nrow(data)) stop("anisotropic model requested but coordinates appears to be a distance matrix") coordinates = vect(coordinates) } maxDist = dist(matrix(ext(coordinates), ncol = 2)) } else { if (is.matrix(coordinates)) { if (ncol(coordinates) == 2) { coordinates = dist(coordinates) } else { coordinates = as(coordinates, "dsyMatrix") else { fromOptim = .C(C_maternLogLOpt, start = as.double(paramsForC), Sparam = as.integer(Sparam), obsCov = as.double(as.matrix(obsCov)), as.double(xcoord), as.double(ycoord), as.integer(aniso), N = as.integer(c(nrow(obsCov), 3, ncol(covariates))), Ltype = as.integer(reml + 2 * !estimateVariance), optInt = as.integer(forO$scalarInt), optF = as.double(forO$scalarF), betas = as.double(forO$pars), limType = as.integer(forO$parsInt), message = format(" ", width = 80)) names(fromOptim$start) = names(paramsForC) if (fromOptim$start["anisoRatio"] < 1) { fromOptim$start["anisoRatio"] <- 1/fromOptim$start["anisoRatio"] } } if (any(class(coordinates) == "dist")) coordinates = as(as.matrix(coordinates), "dsyMatrix") if (length(grep("SpatVect", class(coordinates)))) { if (!nchar(crs(coordinates))) terra::crs(coordinates) = "+proj=utm +zone=1" coordinates = new("dsyMatrix", Dim = rep(length(coordinates), 2), x = as.vector(as.matrix(terra::distance(coordinates))), uplo = "L") } maxDist = max(coordinates, na.rm = TRUE) } if (fromOptim$start["anisoAngleRadians"] + pi/2 >= pi/2) { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] - pi/2 } else { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] + pi/2 } } result = list(optim = list(mle = fromOptim$start, logL = c(m2logL = fromOptim$optF[1], logL = -fromOptim$optF[1]/2), totalVarHat = fromOptim$optF[2], boxcox = fromOptim$optF[3:5], determinants = fromOptim$optF[6:7], message = fromOptim$message, options = cbind(start = paramsForC[Sparam], opt = fromOptim$start[Sparam], parOptions[, c("parscale", "lower", "upper", "ndeps")]), detail = fromOptim$optInt[1:3]), betaHat = fromOptim$betas[1:ncol(covariates)], varBetaHat = new("dsyMatrix", Dim = as.integer(rep(ncol(covariates), 2)), uplo = "L", x = fromOptim$betas[seq(1 + ncol(covariates), len = ncol(covariates)^2)])) result$optim$options = cbind(result$optim$options, gradient = fromOptim$betas[seq(ncol(covariates)^2 + ncol(covariates) + 1, len = sum(Sparam))]) names(result$optim$detail) = c("fail", "fncount", "grcount") trend = formula theFactors = NULL if (any(class(trend) == "formula")) { data = data.frame(data) theNA = apply(data[, all.vars(trend), drop = FALSE], 1, function(qq) any(is.na(qq))) noNA = !theNA theFactors = model.frame(trend, data[noNA, ]) whichFactors = unlist(lapply(theFactors, is.factor)) theFactors = theFactors[, whichFactors, drop = FALSE] theFactors = lapply(theFactors, levels) theFactors = unlist(lapply(theFactors, function(qq) qq[1])) covariates = model.matrix(trend, data[noNA, ]) covariates = covariates[, grep(paste("^(", paste(names(theFactors), collapse = "|"), ")NA$", sep = ""), colnames(covariates), invert = TRUE), drop = FALSE] observations = all.vars(trend)[1] if (!any(names(data) == observations)) warning("can't find observations ", observations, names(result$betaHat) = colnames(covariates) dimnames(result$varBetaHat) = list(names(result$betaHat), names(result$betaHat)) } result$parameters = fillParam(c(result$optim$mle, result$betaHat)) if (estimateVariance) { result$parameters[c("nugget", "variance")] = result$parameters[c("nugget", "variance")] * result$optim$totalVarHat result$varBetaHat = result$varBetaHat * result$optim$totalVarHat } names(result$optim$logL) = paste(names(result$optim$logL), "in data") observations = data[noNA, observations, drop = FALSE] } else { trend = as.matrix(trend) theNA = is.na(data) | apply(trend, 1, function(qq) any(is.na(qq))) noNA = !theNA observations = data[noNA] covariates = trend[noNA, , drop = FALSE] } if (any(theNA)) { c(".ml", ".reml")[reml + 1], sep = "") result$data = cbind(data.frame(observations = observations, fitted = covariates %*% result$parameters[colnames(covariates)]), data.frame(data)[noNA, ]) if (abs(result$parameters["boxcox"] - 1) > 1e-04) { if (abs(result$parameters["boxcox"]) < 0.001) { result$data$obsBC = log(observations) } else { result$data$obsBC <- ((observations^result$parameters["boxcox"]) - 1)/result$parameters["boxcox"] } } else { result$data$obsBC <- obsCov[, 1] } result$data$resid = result$data$obsBC - result$data$fitted if (length(grep("^SpatVector", class(coordinatesOrig)))) { forDf = rep(NA, length(noNA)) forDf[noNA] = seq(1, sum(noNA)) theDf = result$data[forDf, ] result$data = vect(x = crds(coordinatesOrig), atts = theDf, crs = crs(coordinatesOrig)) } result$model = list(reml = reml, baseline = theFactors) if (any(class(trend) == "formula")) { if (length(grep("^SpatVector", class(coordinates)))) { coordinates = coordinates[noNA] } else { if (ncol(coordinates) == nrow(coordinates)) { coordinates = coordinates[noNA, noNA] } else { coordinates = coordinates[noNA, ] } } } estimateVariance = TRUE if (any(paramToEstimate == "variance")) { paramToEstimate = paramToEstimate[paramToEstimate != "variance"] param = param[names(param) != "variance"] } else { if (any(names(param) == "variance")) { estimateVariance = FALSE } result$model$formula = trend result$data[[all.vars(formula)[1]]] = result$data$observations } else { result$model$formula = names(trend) } if (length(result$model$baseline)) { rparams = result$parameters baseParams = data.frame(var = names(result$model$baseline), base = result$model$baseline, pasted = paste(names(result$model$baseline), result$model$baseline, sep = "")) for (D in which(!baseParams$pasted %in% names(rparams))) { sameFac = grep(paste("^", baseParams[D, "var"], sep = ""), names(rparams)) pseq = 1:length(rparams) if (length(sameFac)) { minFac = min(sameFac) toAdd = 0 names(toAdd) = baseParams[D, "pasted"] rparams = c(rparams[pseq < minFac], toAdd, rparams[pseq >= } paramToEstimate = gsub("anisoAngleDegrees", "anisoAngleRadians", paramToEstimate) degToRad = function(par) { if (length(grep("anisoAngleDegrees", names(par)))) { if (!length(grep("anisoAngleRadians", names(par)))) par["anisoAngleRadians"] = 2 * pi * par["anisoAngleDegrees"]/360 par = par[grep("anisoAngleDegrees", names(par), invert = TRUE)] } par } param = degToRad(param) lower = degToRad(lower) upper = degToRad(upper) parscale = degToRad(parscale) lowerDefaults = c(nugget = 0, range = maxDist/1000, anisoRatio = 0.01, anisoAngleRadians = -pi/2, shape = 0.1, boxcox = -1.5, variance = 0) upperDefaults = c(nugget = Inf, range = 10 * maxDist, anisoRatio = 100, anisoAngleRadians = pi/2, shape = 4, boxcox = 2.5, variance = Inf) paramDefaults = c(nugget = 0, anisoRatio = 1, anisoAngleRadians = 0, shape = 1.5, boxcox = 1, range = maxDist/10) minFac]) } } result$parameters = rparams } parameterTable = data.frame(estimate = result$parameters) rownames(parameterTable) = names(result$parameters) parameterTable$stdErr = NA stdErr = sqrt(Matrix::diag(result$varBetaHat)) parameterTable[rownames(result$varBetaHat), "stdErr"] = stdErr thelims = c(0.005, 0.025, 0.05, 0.1) thelims = c(rbind(thelims, 1 - thelims)) if (any(names(paramToEstimate) == "nugget")) { paramDefaults["nugget"] = 1 } parscaleDefaults = c(range = maxDist/5, nugget = 0.05, boxcox = 0.5, anisoAngleRadians = 0.2, anisoRatio = 1, variance = 1, shape = 0.2) ndepsDefault = c(range = 0.01, nugget = 0.05, boxcox = 0.005, anisoAngleRadians = 0.01, anisoRatio = 0.01, variance = 0.01, shape = 0.01) paramDefaults[names(param)] = param parscaleDefaults[names(parscale)] = parscale theQ = qnorm(thelims) toadd = outer(parameterTable$stdErr, theQ) toadd = toadd + matrix(parameterTable$estimate, ncol = length(thelims), nrow = dim(parameterTable)[1]) colnames(toadd) = paste("ci", thelims, sep = "") parameterTable = cbind(parameterTable, toadd) parameterTable[, "pval"] = pchisq(parameterTable$estimate^2/parameterTable$stdErr^2, df = 1, lower.tail = FALSE) parameterTable[, "Estimated"] = FALSE parameterTable[paramToEstimate, "Estimated"] = TRUE parameterTable[rownames(result$varBetaHat), "Estimated"] = TRUE if (estimateVariance) parameterTable["variance", "Estimated"] = TRUE rownames(parameterTable) = gsub("^variance$", "sdSpatial", rownames(parameterTable)) rownames(parameterTable) = gsub("^nugget$", "sdNugget", rownames(parameterTable)) parameterTable[c("sdSpatial", "sdNugget"), "estimate"] = sqrt(parameterTable[c("sdSpatial", lowerDefaults[names(lower)] = lower upperDefaults[names(upper)] = upper if (any(paramToEstimate == "nugget") & paramDefaults["nugget"] == lowerDefaults["nugget"]) { paramDefaults["nugget"] = min(c(0.5, upperDefaults["nugget"])) } startingParam = paramDefaults[paramToEstimate] names(startingParam) = paramToEstimate naStarting = is.na(startingParam) startingParam[naStarting] = paramDefaults[names(startingParam)[naStarting]] moreParams = paramDefaults[!names(paramDefaults) %in% paramToEstimate] allParams = c(startingParam, moreParams) allParams = fillParam(allParams) paramsForC = allParams[c("nugget", "variance", "range", "shape", "anisoRatio", "anisoAngleRadians", "boxcox")] "sdNugget"), "estimate"]) result$summary = as.data.frame(parameterTable) result$summary$Estimated = as.logical(result$summary$Estimated) result})(param = dots[[1L]][[1L]], data = list(fitted = c(7.09263093277653, 7.0632591348391, 7.03388733690166, 7.00451553896422, 6.97514374102678, 6.94577194308934, 6.9164001451519, 6.88702834721446, 6.85765654927703, 6.82828475133959, 6.79891295340215, 6.76954115546471, 6.74016935752727, 6.71079755958983, 6.68142576165239, 6.65205396371495, 6.62268216577752, 6.59331036784008, 6.56393856990264, 6.5345667719652, 6.50519497402776, 6.47582317609032, 6.44645137815288, 6.41707958021545, 6.38770778227801, 6.35833598434057, 6.32896418640313, 6.29959238846569, 6.27022059052825, 6.24084879259081, 6.21147699465337, 6.18210519671594, 6.1527333987785, 6.12336160084106, 6.09398980290362, 6.06461800496618, 6.03524620702874, Sparam = names(paramsForC) %in% paramToEstimate names(Sparam) = names(paramsForC) paramToEstimate = names(Sparam)[Sparam] parOptions = cbind(lower = lowerDefaults[paramToEstimate], upper = upperDefaults[paramToEstimate], parscale = parscaleDefaults[paramToEstimate], ndeps = ndepsDefault[paramToEstimate]) forO = list(scalarF = c(fnscale = -1, abstol = -1, reltol = -1, alpha = -1, beta = -1, gamma = -1, factr = 1e+06, pgtol = 0), scalarInt = c(trace = 0, maxit = 200, REPORT = 1, type = -1, lmm = 25, tmax = -1, temp = -1), pars = parOptions[, c("lower", "upper", "parscale", "ndeps"), drop = FALSE]) if (verbose) { forO$scalarInt["trace"] = 6 forO$scalarInt["REPORT"] = 200 } forO$parsInt = rep(0, nrow(forO$pars)) names(forO$parsInt) = rownames(forO$pars) forO$parsInt[forO$pars[, "lower"] != -Inf & forO$pars[, "upper"] != Inf] = 26.0058744090913, 5.97650261115387, 5.94713081321643, 5.91775901527899, 5.88838721734155, 5.85901541940411, 5.82964362146667, 5.80027182352923, 5.7709000255918, 5.74152822765436, 5.71215642971692, 5.68278463177948, 5.65341283384204, 5.6240410359046, 5.59466923796716, 5.56529744002973, 5.53592564209229, 5.50655384415485, 5.47718204621741, 5.44781024827997, 5.41843845034253, 5.38906665240509, 5.35969485446766, 5.33032305653022, 5.30095125859278, 5.27157946065534, 5.2422076627179, 5.21283586478046, 5.18346406684302, 5.15409226890558, 5.12472047096815, 5.09534867303071, 5.06597687509327, 5.03660507715583, 5.00723327921839, 4.97786148128095, 4.94848968334351, 4.91911788540608, 4.88974608746864, 4.8603742895312, 4.83100249159376, 4.80163069365632, 4.77225889571888, 4.74288709778144, 4.713515299844, 4.68414350190657, 4.65477170396913, 4.62539990603169, 4.59602810809425, 4.56665631015681, 4.53728451221937, 4.50791271428193, 4.4785409163445, 4.44916911840706, 4.41979732046962, 4.39042552253218, 4.36105372459474, 4.3316819266573, 4.30231012871986, 4.27293833078243, 4.24356653284499, 4.21419473490755, 4.18482293697011), ID = c(13L, 14L, 22L, 23L, 24L, 29L, 30L, 35L, 36L, 37L, 52L, 55L, 66L, 71L, forO$parsInt[forO$pars[, "lower"] != -Inf & forO$pars[, "upper"] == Inf] = 1 forO$parsInt[forO$pars[, "lower"] == -Inf & forO$pars[, "upper"] != Inf] = 3 forOparsOrig = forO$pars forO$pars[!is.finite(forO$pars)] = -1 forO$pars = c(forO$pars, rep(0, ncol(covariates) + ncol(covariates)^2)) if (aniso) { xcoord = crds(coordinates)[, 1] ycoord = crds(coordinates)[, 2] } else { xcoord = as.vector(coordinates) ycoord = -99 } obsCov = cbind(y1 = observations, y2 = 0, y3 = 0, covariates) if (all(paramToEstimate == "variance") & param["nugget"] == 0) { theL = loglikLgm(param, data = observations, formula = covariates, coordinates = coordinates, reml = reml) fromOptim = attributes(theL) result = list(optim = list(mle = fillParam(fromOptim$param), logL = c(m2logL = as.numeric(theL), logL = -as.numeric(theL)/2), totalVarHat = fromOptim$totalVarHat, message = "numerical optimization not needed", options = NULL, detail = NULL), betaHat = fromOptim$betaHat, varBetaHat = fromOptim$varBetaHat) } else {73L, 84L, 95L, 102L, 105L, 126L, 130L, 136L, 138L, 159L, 168L, 172L, 178L, 181L, 185L, 188L, 192L, 198L, 202L, 203L, 207L, 208L, 218L, 224L, 226L, 235L, 245L, 246L, 247L, 254L, 261L, 263L, 274L, 275L, 276L, 277L, 281L, 283L, 287L, 292L, 293L, 300L, 302L, 305L, 314L, 317L, 320L, 322L, 335L, 336L, 341L, 342L, 344L, 357L, 362L, 368L, 369L, 372L, 373L, 377L, 378L, 381L, 384L, 392L, 399L, 400L, 401L, 406L, 408L, 409L, 421L, 423L, 424L, 425L, 436L, 442L, 449L, 450L, 451L, 455L, 456L, 458L, 460L, 466L, 468L, 471L), rain = c(15.1, 25.5, 7.9, 19.1, 19.4, 33.4, 10.7, 29.6, 39.4, 39.4, 32.4, 10.5, 13.5, 58.5, 11.4, 33.4, 13.1, 7.8, 39.8, 14.1, 19.2, 15.1, 10.7, 14.5, 33.4, 32.7, 21.3, 33.1, 40, 32.7, 38, 9.4, 18.5, 23.9, 33, 3, 25.4, 5.3, 7.8, 7.1, 6.2, 7.1, 5.9, 6, 12.4, 15.3, 7.5, 13.7, 8.6, 12.9, 34.5, 44.1, 18.4, 12.1, 34.6, 27, 10, 4.5, 10.7, 35.9, 27.8, 7.2, 13.1, 9, 14.1, 13.1, 45.2, 1.6, 13.6, 13, 11.8, 10.9, 14.5, 25.4, 14, 15.2, 6, 28.3, 18.4, 12.7, 22, 17.8, 21.8, fromOptim = .C(C_maternLogLOpt, start = as.double(paramsForC), Sparam = as.integer(Sparam), obsCov = as.double(as.matrix(obsCov)), as.double(xcoord), as.double(ycoord), as.integer(aniso), N = as.integer(c(nrow(obsCov), 3, ncol(covariates))), Ltype = as.integer(reml + 2 * !estimateVariance), optInt = as.integer(forO$scalarInt), optF = as.double(forO$scalarF), betas = as.double(forO$pars), limType = as.integer(forO$parsInt), message = format(" ", width = 80)) names(fromOptim$start) = names(paramsForC) if (fromOptim$start["anisoRatio"] < 1) { fromOptim$start["anisoRatio"] <- 1/fromOptim$start["anisoRatio"] if (fromOptim$start["anisoAngleRadians"] + pi/2 >= 13.7, 14.4, 23, 28.2, 12.9, 6.5, 19, 17, 15.6, 13.1, 1, 9.9, 9.2, 6.7, 1.8, 2, 5.5), elev = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100)), formula = rain ~ elev, paramToEstimate = c("variance", "anisoRatio", "anisoAngleRadians" ), reml = FALSE, coordinates = c(2514336.71, 2518588.71, 2531615.71, 2533523.71, 2534125.71, 2538019.71, 2539319.71, 2545119.71, 2545607.71, 2545627.71, 2561259.71, 2562512.71, 2570067.71, 2571109.71, 2571016.71, 2579205.71, 2586092.71, 2590945.71, 2592122.71, 2608356.71, 2609399.71, 2612205.71, 2612286.71, 2620667.71, 2624915.71, 2626374.71, 2628041.71, 2629488.71, 2630082.71, 2630821.71, 2633239.71, 2635189.71, 2637581.71, 2637613.71, 2640771.71, 2641232.71, 2648275.71, 2655053.71, 2656589.71, 2661121.71, 2668821.71, 2668770.71, 2670977.71, 2676213.71, 2679163.71, 2680537.71, 2685168.71, 2686062.71, 2686015.71, 2687100.71, 2687949.71, 2688944.71, 2688673.71, 2692431.71, 2694076.71, 2696179.71, 2696965.71, 2699592.71, 2700787.71, 2702527.71, 2704185.71, 2704251.71, pi/2) { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] - pi/2 } else { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] + pi/2 } } result = list(optim = list(mle = fromOptim$start, logL = c(m2logL = fromOptim$optF[1], logL = -fromOptim$optF[1]/2), totalVarHat = fromOptim$optF[2], boxcox = fromOptim$optF[3:5], determinants = fromOptim$optF[6:7], message = fromOptim$message, options = cbind(start = paramsForC[Sparam], opt = fromOptim$start[Sparam], parOptions[, c("parscale", "lower", "upper", "ndeps")]), detail = fromOptim$optInt[1:3]), betaHat = fromOptim$betas[1:ncol(covariates)], varBetaHat = new("dsyMatrix", Dim = as.integer(rep(ncol(covariates), 2)), uplo = "L", x = fromOptim$betas[seq(1 + ncol(covariates), len = ncol(covariates)^2)])) result$optim$options = cbind(result$optim$options, gradient = fromOptim$betas[seq(ncol(covariates)^2 + ncol(covariates) + 1, len = sum(Sparam))]) names(result$optim$detail) = c("fail", "fncount", "grcount") 2709024.71, 2711362.71, 2710661.71, 2710650.71, 2713957.71, 2718384.71, 2717971.71, 2720517.71, 2720255.71, 2721055.71, 2723237.71, 2726025.71, 2725583.71, 2726727.71, 2728111.71, 2733095.71, 2736658.71, 2736759.71, 2739493.71, 2740339.71, 2742780.71, 2741684.71, 2749700.71, 2751519.71, 2751266.71, 2752986.71, 2759987.71, 2761869.71, 2767951.71, 2766435.71, 2769345.71, 2775690.71, 2777669.71, 2779806.71, 2782895.71, 2796876.71, 2800101.71, 2805720.71, 1156023.34, 1174834.34, 1177889.34, 1196758.34, 1188960.34, 1154404.34, 1182184.34, 1207655.34, 1150925.34, 1152037.34, 1172904.34, 1252956.34, 1255067.34, 1168310.34, 1106035.34, 1204901.34, 1143655.34, 1094673.34, 1206974.34, 1149002.34, 1185690.34, 1257948.34, 1269067.34, 1160040.34, 1244526.34, 1234511.34, 1268974.34, names(result$betaHat) = colnames(covariates) dimnames(result$varBetaHat) = list(names(result$betaHat), names(result$betaHat)) } result$parameters = fillParam(c(result$optim$mle, result$betaHat)) if (estimateVariance) { result$parameters[c("nugget", "variance")] = result$parameters[c("nugget", "variance")] * result$optim$totalVarHat result$varBetaHat = result$varBetaHat * result$optim$totalVarHat } names(result$optim$logL) = paste(names(result$optim$logL), c(".ml", ".reml")[reml + 1], sep = "") 1256736.34, 1218927.34, 1214477.34, 1254497.34, 1161086.34, 1195550.34, 1206669.34, 1258922.34, 1129935.34, 1248902.34, 1185517.34, 1135480.34, 1203312.34, 1154399.34, 1176638.34, 1205554.34, 1225586.34, 1243389.34, 1273417.34, 1247865.34, 1221182.34, 1230078.34, 1168925.34, 1153362.34, 1112225.34, 1292361.34, 1289049.34, 1152287.34, 1180101.34, 1282408.34, 1232389.34, 1272429.34, 1147900.34, 1132346.34, 1216858.34, 1273611.34, 1103497.34, 1259171.34, 1260283.34, 1153562.34, 1095782.34, 1278149.34, 1251488.34, 1274838.34, 1270399.34, 1211486.34, 1168149.34, 1268227.34, 1235992.34, 1246018.34, 1152668.34, 1227225.34, 1273933.34, 1187230.34, 1233949.34, 1170596.34, 1245090.34, 1261895.34, 1196309.34, 1211875.34, 1152960.34, 1146405.34, 1169794.34, 1171018.34, 1251067.34, 1137679.34, 1165608.34, 1143403.34, 1187937.34, 1185778.34, 1141601.34, 1135004.34, 1125130.34)) 2: .mapply(FUN, dots, MoreArgs) 3: result$data = cbind(data.frame(observations = observations, fitted = covariates %*% result$parameters[colnames(covariates)]), data.frame(data)[noNA, ]) if (abs(result$parameters["boxcox"] - 1) > 1e-04) { if (abs(result$parameters["boxcox"]) < 0.001) { result$data$obsBC = log(observations) } else { result$data$obsBC <- ((observations^result$parameters["boxcox"]) - 1)/result$parameters["boxcox"] } } else { result$data$obsBC <- obsCov[, 1] } result$data$resid = result$data$obsBC - result$data$fitted if (length(grep("^SpatVector", class(coordinatesOrig)))) { forDf = rep(NA, length(noNA)) forDf[noNA] = seq(1, sum(noNA))FUN(X[[i]], ...) 4: lapply(X = S, FUN = FUN, ...) 5: doTryCatch(return(expr), name, parentenv, handler) 6: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 7: tryCatchList(expr, classes, parentenv, handlers) 8: tryCatch(expr, error = function(e) { call <- conditionCall(e) if (!is.null(call)) { if (identical(call[[1L]], quote(doTryCatch))) call <- sys.call(-4L) dcall <- deparse(call, nlines = 1L) prefix <- paste("Error in", dcall, ": ") LONG <- 75L sm <- strsplit(conditionMessage(e), "\n")[[1L]] w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w") if (is.na(w)) w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], type = "b") if (w > LONG) prefix <- paste0(prefix, "\n ") theDf = result$data[forDf, ] result$data = vect(x = crds(coordinatesOrig), atts = theDf, crs = crs(coordinatesOrig)) } result$model = list(reml = reml, baseline = theFactors) if (any(class(trend) == "formula")) { result$model$formula = trend result$data[[all.vars(formula)[1]]] = result$data$observations } else { result$model$formula = names(trend) } if (length(result$model$baseline)) { rparams = result$parameters baseParams = data.frame(var = names(result$model$baseline), base = result$model$baseline, pasted = paste(names(result$model$baseline), result$model$baseline, sep = "")) for (D in which(!baseParams$pasted %in% names(rparams))) { sameFac = grep(paste("^", baseParams[D, "var"], sep = ""), names(rparams)) pseq = 1:length(rparams) if (length(sameFac)) { } else prefix <- "Error : " msg <- paste0(prefix, conditionMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && isTRUE(getOption("show.error.messages"))) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))}) 9: minFac = min(sameFac) toAdd = 0 names(toAdd) = baseParams[D, "pasted"] rparams = c(rparams[pseq < minFac], toAdd, rparams[pseq >= minFac]) } } result$parameters = rparams } parameterTable = data.frame(estimate = result$parameters) rownames(parameterTable) = names(result$parameters) parameterTable$stdErr = NA stdErr = sqrt(Matrix::diag(result$varBetaHat)) parameterTable[rownames(result$varBetaHat), "stdErr"] = stdErr thelims = c(0.005, 0.025, 0.05, 0.1)try(lapply(X = S, FUN = FUN, ...), silent = TRUE) 10: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) 11: FUN(X[[i]], ...) 12: lapply(seq_len(cores), inner.do) 13: mclapply(seq_len(n), do_one, mc.preschedule = mc.preschedule, mc.set.seed = mc.set.seed, mc.silent = mc.silent, mc.cores = mc.cores, thelims = c(rbind(thelims, 1 - thelims)) theQ = qnorm(thelims) toadd = outer(parameterTable$stdErr, theQ) toadd = toadd + matrix(parameterTable$estimate, ncol = length(thelims), nrow = dim(parameterTable)[1]) colnames(toadd) = paste("ci", thelims, sep = "") parameterTable = cbind(parameterTable, toadd) parameterTable[, "pval"] = pchisq(parameterTable$estimate^2/parameterTable$stdErr^2, df = 1, lower.tail = FALSE) parameterTable[, "Estimated"] = FALSE parameterTable[paramToEstimate, "Estimated"] = TRUE parameterTable[rownames(result$varBetaHat), "Estimated"] = TRUE if (estimateVariance) parameterTable["variance", "Estimated"] = TRUE rownames(parameterTable) = gsub("^variance$", "sdSpatial", rownames(parameterTable)) rownames(parameterTable) = gsub("^nugget$", "sdNugget", rownames(parameterTable)) parameterTable[c("sdSpatial", "sdNugget"), "estimate"] = sqrt(parameterTable[c("sdSpatial", "sdNugget"), "estimate"]) result$summary = as.data.frame(parameterTable) result$summary$Estimated = as.logical(result$summary$Estimated) result mc.cleanup = mc.cleanup, affinity.list = affinity.list) 14: parallel::mcmapply(likfitLgm, param = parList, MoreArgs = list(data = as.data.frame(fit$data), formula = fit$model$formula, paramToEstimate = reEstimate, reml = fit$model$reml, coordinates = terra::crds(fit$data)), mc.cores = mc.cores, SIMPLIFY = FALSE) 15: profLlgm(swissFit, mc.cores = Ncores, range = seq(15000, 55000, len = 12)) An irrecoverable exception occurred. R is aborting now ... })(param = dots[[1L]][[1L]], data = list(fitted = c(7.09263093277653, 7.0632591348391, 7.03388733690166, 7.00451553896422, 6.97514374102678, 6.94577194308934, 6.9164001451519, 6.88702834721446, 6.85765654927703, 6.82828475133959, 6.79891295340215, 6.76954115546471, 6.74016935752727, 6.71079755958983, 6.68142576165239, 6.65205396371495, 6.62268216577752, 6.59331036784008, 6.56393856990264, 6.5345667719652, 6.50519497402776, 6.47582317609032, 6.44645137815288, 6.41707958021545, 6.38770778227801, 6.35833598434057, 6.32896418640313, 6.29959238846569, 6.27022059052825, 6.24084879259081, 6.21147699465337, 6.18210519671594, 6.1527333987785, 6.12336160084106, 6.09398980290362, 6.06461800496618, 6.03524620702874, 6.0058744090913, 5.97650261115387, 5.94713081321643, 5.91775901527899, 5.88838721734155, 5.85901541940411, 5.82964362146667, 5.80027182352923, 5.7709000255918, 5.74152822765436, 5.71215642971692, 5.68278463177948, 5.65341283384204, 5.6240410359046, 5.59466923796716, 5.56529744002973, 5.53592564209229, 5.50655384415485, 5.47718204621741, 5.44781024827997, 5.41843845034253, 5.38906665240509, 5.35969485446766, 5.33032305653022, 5.30095125859278, 5.27157946065534, 5.2422076627179, 5.21283586478046, 5.18346406684302, 5.15409226890558, 5.12472047096815, 5.09534867303071, 5.06597687509327, 5.03660507715583, 5.00723327921839, 4.97786148128095, 4.94848968334351, 4.91911788540608, 4.88974608746864, 4.8603742895312, 4.83100249159376, 4.80163069365632, 4.77225889571888, 4.74288709778144, 4.713515299844, 4.68414350190657, 4.65477170396913, 4.62539990603169, 4.59602810809425, 4.56665631015681, 4.53728451221937, 4.50791271428193, 4.4785409163445, 4.44916911840706, 4.41979732046962, 4.39042552253218, 4.36105372459474, 4.3316819266573, 4.30231012871986, 4.27293833078243, 4.24356653284499, 4.21419473490755, 4.18482293697011), ID = c(13L, 14L, 22L, 23L, 24L, 29L, 30L, 35L, 36L, 37L, 52L, 55L, 66L, 71L, 73L, 84L, 95L, 102L, 105L, 126L, 130L, 136L, 138L, 159L, 168L, 172L, 178L, 181L, 185L, 188L, 192L, 198L, 202L, 203L, 207L, 208L, 218L, 224L, 226L, 235L, 245L, 246L, 247L, 254L, 261L, 263L, 274L, 275L, 276L, 277L, 281L, 283L, 287L, 292L, 293L, 300L, 302L, 305L, 314L, 317L, 320L, 322L, 335L, 336L, 341L, 342L, 344L, 357L, 362L, 368L, 369L, 372L, 373L, 377L, 378L, 381L, 384L, 392L, 399L, 400L, 401L, 406L, 408L, 409L, 421L, 423L, 424L, 425L, 436L, 442L, 449L, 450L, 451L, 455L, 456L, 458L, 460L, 466L, 468L, 471L), rain = c(15.1, 25.5, 7.9, 19.1, 19.4, 33.4, 10.7, 29.6, 39.4, 39.4, 32.4, 10.5, 13.5, 58.5, 11.4, 33.4, 13.1, 7.8, 39.8, 14.1, 19.2, 15.1, 10.7, 14.5, 33.4, 32.7, 21.3, 33.1, 40, 32.7, 38, 9.4, 18.5, 23.9, 33, 3, 25.4, 5.3, 7.8, 7.1, 6.2, 7.1, 5.9, 6, 12.4, 15.3, 7.5, 13.7, 8.6, 12.9, 34.5, 44.1, 18.4, 12.1, 34.6, 27, 10, 4.5, 10.7, 35.9, 27.8, 7.2, 13.1, 9, 14.1, 13.1, 45.2, 1.6, 13.6, 13, 11.8, 10.9, 14.5, 25.4, 14, 15.2, 6, 28.3, 18.4, 12.7, 22, 17.8, 21.8, 13.7, 14.4, 23, 28.2, 12.9, 6.5, 19, 17, 15.6, 13.1, 1, 9.9, 9.2, 6.7, 1.8, 2, 5.5), elev = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100)), formula = rain ~ elev, paramToEstimate = c("variance", "anisoRatio", "anisoAngleRadians" ), reml = FALSE, coordinates = c(2514336.71, 2518588.71, 2531615.71, 2533523.71, 2534125.71, 2538019.71, 2539319.71, 2545119.71, 2545607.71, 2545627.71, 2561259.71, 2562512.71, 2570067.71, 2571109.71, 2571016.71, 2579205.71, 2586092.71, 2590945.71, 2592122.71, 2608356.71, 2609399.71, 2612205.71, 2612286.71, 2620667.71, 2624915.71, 2626374.71, 2628041.71, 2629488.71, 2630082.71, 2630821.71, 2633239.71, 2635189.71, 2637581.71, 2637613.71, 2640771.71, 2641232.71, 2648275.71, 2655053.71, 2656589.71, 2661121.71, 2668821.71, 2668770.71, 2670977.71, 2676213.71, 2679163.71, 2680537.71, 2685168.71, 2686062.71, 2686015.71, 2687100.71, 2687949.71, 2688944.71, 2688673.71, 2692431.71, 2694076.71, 2696179.71, 2696965.71, 2699592.71, 2700787.71, 2702527.71, 2704185.71, 2704251.71, 2709024.71, 2711362.71, 2710661.71, 2710650.71, 2713957.71, 2718384.71, 2717971.71, 2720517.71, 2720255.71, 2721055.71, 2723237.71, 2726025.71, 2725583.71, 2726727.71, 2728111.71, 2733095.71, 2736658.71, 2736759.71, 2739493.71, 2740339.71, 2742780.71, 2741684.71, 2749700.71, 2751519.71, 2751266.71, 2752986.71, 2759987.71, 2761869.71, 2767951.71, 2766435.71, 2769345.71, 2775690.71, 2777669.71, 2779806.71, 2782895.71, 2796876.71, 2800101.71, 2805720.71, 1156023.34, 1174834.34, 1177889.34, 1196758.34, 1188960.34, 1154404.34, 1182184.34, 1207655.34, 1150925.34, 1152037.34, 1172904.34, 1252956.34, 1255067.34, 1168310.34, 1106035.34, 1204901.34, 1143655.34, 1094673.34, 1206974.34, 1149002.34, 1185690.34, 1257948.34, 1269067.34, 1160040.34, 1244526.34, 1234511.34, 1268974.34, 1256736.34, 1218927.34, 1214477.34, 1254497.34, 1161086.34, 1195550.34, 1206669.34, 1258922.34, 1129935.34, 1248902.34, 1185517.34, 1135480.34, 1203312.34, 1154399.34, 1176638.34, 1205554.34, 1225586.34, 1243389.34, 1273417.34, 1247865.34, 1221182.34, 1230078.34, 1168925.34, 1153362.34, 1112225.34, 1292361.34, 1289049.34, 1152287.34, 1180101.34, 1282408.34, 1232389.34, 1272429.34, 1147900.34, 1132346.34, 1216858.34, 1273611.34, 1103497.34, 1259171.34, 1260283.34, 1153562.34, 1095782.34, 1278149.34, 1251488.34, 1274838.34, 1270399.34, 1211486.34, 1168149.34, 1268227.34, 1235992.34, 1246018.34, 1152668.34, 1227225.34, 1273933.34, 1187230.34, 1233949.34, 1170596.34, 1245090.34, 1261895.34, 1196309.34, 1211875.34, 1152960.34, 1146405.34, 1169794.34, 1171018.34, 1251067.34, 1137679.34, 1165608.34, 1143403.34, 1187937.34, 1185778.34, 1141601.34, 1135004.34, 1125130.34)) 2: .mapply(FUN, dots, MoreArgs) 3: FUN(X[[i]], ...) 4: lapply(X = S, FUN = FUN, ...) 5: doTryCatch(return(expr), name, parentenv, handler) 6: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 7: tryCatchList(expr, classes, parentenv, handlers) 8: tryCatch(expr, error = function(e) { call <- conditionCall(e) if (!is.null(call)) { if (identical(call[[1L]], quote(doTryCatch))) call <- sys.call(-4L) dcall <- deparse(call, nlines = 1L) prefix <- paste("Error in", dcall, ": ") LONG <- 75L sm <- strsplit(conditionMessage(e), "\n")[[1L]] w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w") if (is.na(w)) w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], type = "b") if (w > LONG) prefix <- paste0(prefix, "\n ") } else prefix <- "Error : " msg <- paste0(prefix, conditionMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && isTRUE(getOption("show.error.messages"))) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))}) 9: try(lapply(X = S, FUN = FUN, ...), silent = TRUE) 10: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) 11: FUN(X[[i]], ...) 12: lapply(seq_len(cores), inner.do) 13: mclapply(seq_len(n), do_one, mc.preschedule = mc.preschedule, mc.set.seed = mc.set.seed, mc.silent = mc.silent, mc.cores = mc.cores, mc.cleanup = mc.cleanup, affinity.list = affinity.list) 14: parallel::mcmapply(likfitLgm, param = parList, MoreArgs = list(data = as.data.frame(fit$data), formula = fit$model$formula, paramToEstimate = reEstimate, reml = fit$model$reml, coordinates = terra::crds(fit$data)), mc.cores = mc.cores, SIMPLIFY = FALSE) 15: profLlgm(swissFit, mc.cores = Ncores, range = seq(15000, 55000, len = 12)) An irrecoverable exception occurred. R is aborting now ... Error in resL[grep("^m2logL", rownames(resL)), ] : incorrect number of dimensions Calls: profLlgm In addition: Warning message: In mclapply(seq_len(n), do_one, mc.preschedule = mc.preschedule, : scheduled cores 1, 2 did not deliver results, all values of the jobs will be affected Execution halted Flavor: r-devel-macos-arm64

Version: 2.0.11
Check: tests
Result: ERROR Running ‘RFsimulate.R’ [3s/3s] Running ‘krige.R’ [3s/3s] Running ‘lgcp.R’ [2s/2s] Running ‘lgm.R’ [7s/7s] Running ‘lgmRaster.R’ [48s/22s] Running ‘likfitLgm.R’ [2s/2s] Running ‘matern.R’ [2s/2s] Running ‘maternGmrfPrec.R’ [3s/3s] Running ‘profLlgm.R’ [2s/2s] Running the tests in ‘tests/profLlgm.R’ failed. Complete output: > > library('geostatsp') Loading required package: Matrix Loading required package: terra terra 1.8.93 > data('swissRain') > swissRain = unwrap(swissRain) > swissAltitude = unwrap(swissAltitude) > > Ncores = c(1,2)[1+(.Platform$OS.type=='unix')] > > > > sr2 = swissRain > sr2$elev = extract(swissAltitude, sr2) Warning message: [`[[<-`] only using the first column > swissFit = likfitLgm( + data=sr2, + formula=rain~ elev, + param=c(range=10000,shape=1,nugget=0,boxcox=0.5,anisoRatio=2,anisoAngleDegrees=45), + paramToEstimate = c("range",'anisoAngleDegrees','anisoRatio'), + reml=FALSE, + verbose=FALSE + ) > > > # calculate log-likelihood at the MLE's, but re-estimate variance > sl = loglikLgm( + swissFit$param[c('range','shape','boxcox', 'anisoRatio', 'anisoAngleRadians')], + data=sr2, + formula=rain~ elev, + reml=swissFit$model$reml) > > > # calculate log-likelihood without re-estimating variance > sigSqHat = attributes(sl)$totalVarHat > sl1 = loglikLgm( + c(attributes(sl)$param[ + c('boxcox','anisoRatio','anisoAngleRadians','shape', 'range')], + variance=sigSqHat), + data=sr2, + formula=rain~ elev, + reml=swissFit$model$reml) > > > # re=estimate the anisotropy parameters but not the range > sf2 = likfitLgm( + data=swissFit$data, + formula=swissFit$model$formula, + param= swissFit$param[c('range','nugget','shape','boxcox', 'anisoRatio', 'anisoAngleRadians')], + paramToEstimate = c('variance','anisoAngleRadians','anisoRatio'), + reml=swissFit$model$reml) > > # these should all be the same > as.numeric(sl1) [1] 644.4812 > as.numeric(sl) [1] 644.4812 > swissFit$optim$logL m2logL.ml logL.ml 644.4812 -322.2406 > sf2$optim$logL m2logL.ml logL.ml 644.4812 -322.2406 > > date() [1] "Sat Feb 28 10:35:13 2026" > x=profLlgm(swissFit, mc.cores=Ncores, + range=seq(15000, 55000 , len=12) + ) *** caught segfault *** address 0x110, cause 'invalid permissions' *** caught segfault *** address 0x110, cause 'invalid permissions' Traceback: 1: (function (formula, data, paramToEstimate = c("range", "nugget"), reml = TRUE, coordinates = data, param = NULL, upper = NULL, lower = NULL, parscale = NULL, verbose = FALSE) { param = param[!is.na(param)] coordinatesOrig = coordinates aniso = as.logical(length(grep("^aniso", c(names(param), paramToEstimate)))) | any(abs(param["anisoRatio"] - 1) > 1e-05, na.rm = TRUE) if (aniso) { if (is.matrix(coordinates)) { if (ncol(coordinates) != 2 | nrow(coordinates) != nrow(data)) stop("anisotropic model requested but coordinates appears to be a distance matrix") coordinates = vect(coordinates) } maxDist = dist(matrix(ext(coordinates), ncol = 2)) } else { if (is.matrix(coordinates)) { if (ncol(coordinates) == 2) { coordinates = dist(coordinates) } else { coordinates = as(coordinates, "dsyMatrix") } } if (any(class(coordinates) == "dist")) coordinates = as(as.matrix(coordinates), "dsyMatrix") if (length(grep("SpatVect", class(coordinates)))) { if (!nchar(crs(coordinates))) terra::crs(coordinates) = "+proj=utm +zone=1" coordinates = new("dsyMatrix", Dim = rep(length(coordinates), 2), x = as.vector(as.matrix(terra::distance(coordinates))), uplo = "L") } maxDist = max(coordinates, na.rm = TRUE) } trend = formula theFactors = NULL if (any(class(trend) == "formula")) { data = data.frame(data) theNA = apply(data[, all.vars(trend), drop = FALSE], 1, function(qq) any(is.na(qq))) noNA = !theNA Traceback: 1: (function (formula, data, paramToEstimate = c("range", "nugget"), reml = TRUE, coordinates = data, param = NULL, upper = NULL, lower = NULL, parscale = NULL, verbose = FALSE) { param = param[!is.na(param)] coordinatesOrig = coordinates aniso = as.logical(length(grep("^aniso", c(names(param), paramToEstimate)))) | any(abs(param["anisoRatio"] - 1) > 1e-05, na.rm = TRUE) if (aniso) { if (is.matrix(coordinates)) { if (ncol(coordinates) != 2 | nrow(coordinates) != theFactors = model.frame(trend, data[noNA, ]) whichFactors = unlist(lapply(theFactors, is.factor)) theFactors = theFactors[, whichFactors, drop = FALSE] theFactors = lapply(theFactors, levels) theFactors = unlist(lapply(theFactors, function(qq) qq[1])) covariates = model.matrix(trend, data[noNA, ]) covariates = covariates[, grep(paste("^(", paste(names(theFactors), collapse = "|"), ")NA$", sep = ""), colnames(covariates), invert = TRUE), drop = FALSE] observations = all.vars(trend)[1] if (!any(names(data) == observations)) nrow(data)) stop("anisotropic model requested but coordinates appears to be a distance matrix") coordinates = vect(coordinates) } maxDist = dist(matrix(ext(coordinates), ncol = 2)) } else { if (is.matrix(coordinates)) { if (ncol(coordinates) == 2) { coordinates = dist(coordinates) } else { coordinates = as(coordinates, "dsyMatrix") } warning("can't find observations ", observations, "in data") observations = data[noNA, observations, drop = FALSE] } else { trend = as.matrix(trend) theNA = is.na(data) | apply(trend, 1, function(qq) any(is.na(qq))) noNA = !theNA observations = data[noNA] covariates = trend[noNA, , drop = FALSE] } if (any(theNA)) { if (length(grep("^SpatVector", class(coordinates)))) { coordinates = coordinates[noNA] } if (any(class(coordinates) == "dist")) coordinates = as(as.matrix(coordinates), "dsyMatrix") if (length(grep("SpatVect", class(coordinates)))) { if (!nchar(crs(coordinates))) terra::crs(coordinates) = "+proj=utm +zone=1" coordinates = new("dsyMatrix", Dim = rep(length(coordinates), 2), x = as.vector(as.matrix(terra::distance(coordinates))), uplo = "L") } maxDist = max(coordinates, na.rm = TRUE) } else { if (ncol(coordinates) == nrow(coordinates)) { coordinates = coordinates[noNA, noNA] } else { coordinates = coordinates[noNA, ] } } } estimateVariance = TRUE if (any(paramToEstimate == "variance")) { paramToEstimate = paramToEstimate[paramToEstimate != "variance"] param = param[names(param) != "variance"] } else { } trend = formula theFactors = NULL if (any(class(trend) == "formula")) { data = data.frame(data) theNA = apply(data[, all.vars(trend), drop = FALSE], 1, function(qq) any(is.na(qq))) noNA = !theNA if (any(names(param) == "variance")) { estimateVariance = FALSE } } paramToEstimate = gsub("anisoAngleDegrees", "anisoAngleRadians", paramToEstimate) degToRad = function(par) { if (length(grep("anisoAngleDegrees", names(par)))) { if (!length(grep("anisoAngleRadians", names(par)))) par["anisoAngleRadians"] = 2 * pi * par["anisoAngleDegrees"]/360 par = par[grep("anisoAngleDegrees", names(par), invert = TRUE)] } par } param = degToRad(param) lower = degToRad(lower) upper = degToRad(upper) parscale = degToRad(parscale) lowerDefaults = c(nugget = 0, range = maxDist/1000, anisoRatio = 0.01, anisoAngleRadians = -pi/2, shape = 0.1, boxcox = -1.5, variance = 0) upperDefaults = c(nugget = Inf, range = 10 * maxDist, anisoRatio = 100, theFactors = model.frame(trend, data[noNA, ]) whichFactors = unlist(lapply(theFactors, is.factor)) theFactors = theFactors[, whichFactors, drop = FALSE] theFactors = lapply(theFactors, levels) theFactors = unlist(lapply(theFactors, function(qq) qq[1])) covariates = model.matrix(trend, data[noNA, ]) covariates = covariates[, grep(paste("^(", paste(names(theFactors), collapse = "|"), ")NA$", sep = ""), colnames(covariates), invert = TRUE), drop = FALSE] observations = all.vars(trend)[1] if (!any(names(data) == observations)) warning("can't find observations ", observations, "in data") observations = data[noNA, observations, drop = FALSE] } else { trend = as.matrix(trend) anisoAngleRadians = pi/2, shape = 4, boxcox = 2.5, variance = Inf) paramDefaults = c(nugget = 0, anisoRatio = 1, anisoAngleRadians = 0, shape = 1.5, boxcox = 1, range = maxDist/10) if (any(names(paramToEstimate) == "nugget")) { paramDefaults["nugget"] = 1 } parscaleDefaults = c(range = maxDist/5, nugget = 0.05, boxcox = 0.5, anisoAngleRadians = 0.2, anisoRatio = 1, variance = 1, shape = 0.2) ndepsDefault = c(range = 0.01, nugget = 0.05, boxcox = 0.005, theNA = is.na(data) | apply(trend, 1, function(qq) any(is.na(qq))) noNA = !theNA observations = data[noNA] covariates = trend[noNA, , drop = FALSE] } if (any(theNA)) { if (length(grep("^SpatVector", class(coordinates)))) { anisoAngleRadians = 0.01, anisoRatio = 0.01, variance = 0.01, coordinates = coordinates[noNA] shape = 0.01) } else { if (ncol(coordinates) == nrow(coordinates)) { coordinates = coordinates[noNA, noNA] } else { coordinates = coordinates[noNA, ] } } } estimateVariance = TRUE if (any(paramToEstimate == "variance")) { paramDefaults[names(param)] = param paramToEstimate = paramToEstimate[paramToEstimate != "variance"] param = param[names(param) != "variance"] } else { if (any(names(param) == "variance")) { estimateVariance = FALSE } } paramToEstimate = gsub("anisoAngleDegrees", "anisoAngleRadians", parscaleDefaults[names(parscale)] = parscale lowerDefaults[names(lower)] = lower upperDefaults[names(upper)] = upper if (any(paramToEstimate == "nugget") & paramDefaults["nugget"] == lowerDefaults["nugget"]) { paramDefaults["nugget"] = min(c(0.5, upperDefaults["nugget"])) } startingParam = paramDefaults[paramToEstimate] paramToEstimate) degToRad = function(par) { if (length(grep("anisoAngleDegrees", names(par)))) { if (!length(grep("anisoAngleRadians", names(par)))) par["anisoAngleRadians"] = 2 * pi * par["anisoAngleDegrees"]/360 par = par[grep("anisoAngleDegrees", names(par), invert = TRUE)] } par } names(startingParam) = paramToEstimate naStarting = is.na(startingParam) startingParam[naStarting] = paramDefaults[names(startingParam)[naStarting]] moreParams = paramDefaults[!names(paramDefaults) %in% paramToEstimate] allParams = c(startingParam, moreParams) allParams = fillParam(allParams) paramsForC = allParams[c("nugget", "variance", "range", "shape", "anisoRatio", "anisoAngleRadians", "boxcox")] Sparam = names(paramsForC) %in% paramToEstimate names(Sparam) = names(paramsForC) paramToEstimate = names(Sparam)[Sparam] parOptions = cbind(lower = lowerDefaults[paramToEstimate], upper = upperDefaults[paramToEstimate], parscale = parscaleDefaults[paramToEstimate], ndeps = ndepsDefault[paramToEstimate]) forO = list(scalarF = c(fnscale = -1, abstol = -1, reltol = -1, alpha = -1, beta = -1, gamma = -1, factr = 1e+06, pgtol = 0), scalarInt = c(trace = 0, maxit = 200, REPORT = 1, type = -1, param = degToRad(param) lower = degToRad(lower) upper = degToRad(upper) parscale = degToRad(parscale) lowerDefaults = c(nugget = 0, range = maxDist/1000, anisoRatio = 0.01, anisoAngleRadians = -pi/2, shape = 0.1, boxcox = -1.5, variance = 0) upperDefaults = c(nugget = Inf, range = 10 * maxDist, anisoRatio = 100, anisoAngleRadians = pi/2, shape = 4, boxcox = 2.5, variance = Inf) paramDefaults = c(nugget = 0, anisoRatio = 1, anisoAngleRadians = 0, shape = 1.5, boxcox = 1, range = maxDist/10) if (any(names(paramToEstimate) == "nugget")) { paramDefaults["nugget"] = 1 } parscaleDefaults = c(range = maxDist/5, nugget = 0.05, boxcox = 0.5, lmm = 25, tmax = -1, temp = -1), pars = parOptions[, c("lower", "upper", "parscale", "ndeps"), drop = FALSE]) if (verbose) { forO$scalarInt["trace"] = 6 forO$scalarInt["REPORT"] = 200 } forO$parsInt = rep(0, nrow(forO$pars)) names(forO$parsInt) = rownames(forO$pars) forO$parsInt[forO$pars[, "lower"] != -Inf & forO$pars[, "upper"] != Inf] = 2 forO$parsInt[forO$pars[, "lower"] != -Inf & forO$pars[, "upper"] == Inf] = 1 anisoAngleRadians = 0.2, anisoRatio = 1, variance = 1, shape = 0.2) ndepsDefault = c(range = 0.01, nugget = 0.05, boxcox = 0.005, anisoAngleRadians = 0.01, anisoRatio = 0.01, variance = 0.01, shape = 0.01) paramDefaults[names(param)] = param parscaleDefaults[names(parscale)] = parscale lowerDefaults[names(lower)] = lower upperDefaults[names(upper)] = upper if (any(paramToEstimate == "nugget") & paramDefaults["nugget"] == lowerDefaults["nugget"]) { forO$parsInt[forO$pars[, "lower"] == -Inf & forO$pars[, "upper"] != Inf] = 3 forOparsOrig = forO$pars forO$pars[!is.finite(forO$pars)] = -1 forO$pars = c(forO$pars, rep(0, ncol(covariates) + ncol(covariates)^2)) if (aniso) { xcoord = crds(coordinates)[, 1] ycoord = crds(coordinates)[, 2] } else { xcoord = as.vector(coordinates) paramDefaults["nugget"] = min(c(0.5, upperDefaults["nugget"])) } startingParam = paramDefaults[paramToEstimate] names(startingParam) = paramToEstimate naStarting = is.na(startingParam) startingParam[naStarting] = paramDefaults[names(startingParam)[naStarting]] moreParams = paramDefaults[!names(paramDefaults) %in% paramToEstimate] allParams = c(startingParam, moreParams) allParams = fillParam(allParams) paramsForC = allParams[c("nugget", "variance", "range", "shape", "anisoRatio", "anisoAngleRadians", "boxcox")] Sparam = names(paramsForC) %in% paramToEstimate names(Sparam) = names(paramsForC) paramToEstimate = names(Sparam)[Sparam] ycoord = -99 } obsCov = cbind(y1 = observations, y2 = 0, y3 = 0, covariates) if (all(paramToEstimate == "variance") & param["nugget"] == 0) { theL = loglikLgm(param, data = observations, formula = covariates, coordinates = coordinates, reml = reml) fromOptim = attributes(theL) result = list(optim = list(mle = fillParam(fromOptim$param), logL = c(m2logL = as.numeric(theL), logL = -as.numeric(theL)/2), totalVarHat = fromOptim$totalVarHat, message = "numerical optimization not needed", options = NULL, detail = NULL), betaHat = fromOptim$betaHat, parOptions = cbind(lower = lowerDefaults[paramToEstimate], upper = upperDefaults[paramToEstimate], parscale = parscaleDefaults[paramToEstimate], ndeps = ndepsDefault[paramToEstimate]) forO = list(scalarF = c(fnscale = -1, abstol = -1, reltol = -1, alpha = -1, beta = -1, gamma = -1, factr = 1e+06, pgtol = 0), scalarInt = c(trace = 0, maxit = 200, REPORT = 1, type = -1, lmm = 25, tmax = -1, temp = -1), pars = parOptions[, c("lower", "upper", "parscale", "ndeps"), drop = FALSE]) if (verbose) { forO$scalarInt["trace"] = 6 forO$scalarInt["REPORT"] = 200 } forO$parsInt = rep(0, nrow(forO$pars)) names(forO$parsInt) = rownames(forO$pars) forO$parsInt[forO$pars[, "lower"] != -Inf & forO$pars[, "upper"] != varBetaHat = fromOptim$varBetaHat) } else { fromOptim = .C(C_maternLogLOpt, start = as.double(paramsForC), Sparam = as.integer(Sparam), obsCov = as.double(as.matrix(obsCov)), as.double(xcoord), as.double(ycoord), as.integer(aniso), N = as.integer(c(nrow(obsCov), 3, ncol(covariates))), Ltype = as.integer(reml + 2 * !estimateVariance), optInt = as.integer(forO$scalarInt), optF = as.double(forO$scalarF), betas = as.double(forO$pars), limType = as.integer(forO$parsInt), message = format(" ", width = 80)) names(fromOptim$start) = names(paramsForC) if (fromOptim$start["anisoRatio"] < 1) { fromOptim$start["anisoRatio"] <- 1/fromOptim$start["anisoRatio"] if (fromOptim$start["anisoAngleRadians"] + pi/2 >= pi/2) { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] - pi/2 Inf] = 2 forO$parsInt[forO$pars[, "lower"] != -Inf & forO$pars[, "upper"] == Inf] = 1 forO$parsInt[forO$pars[, "lower"] == -Inf & forO$pars[, "upper"] != Inf] = 3 forOparsOrig = forO$pars forO$pars[!is.finite(forO$pars)] = -1 forO$pars = c(forO$pars, rep(0, ncol(covariates) + ncol(covariates)^2)) if (aniso) { xcoord = crds(coordinates)[, 1] ycoord = crds(coordinates)[, 2] } else { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] + pi/2 } } result = list(optim = list(mle = fromOptim$start, logL = c(m2logL = fromOptim$optF[1], logL = -fromOptim$optF[1]/2), totalVarHat = fromOptim$optF[2], boxcox = fromOptim$optF[3:5], determinants = fromOptim$optF[6:7], message = fromOptim$message, options = cbind(start = paramsForC[Sparam], opt = fromOptim$start[Sparam], parOptions[, c("parscale", "lower", "upper", "ndeps")]), detail = fromOptim$optInt[1:3]), betaHat = fromOptim$betas[1:ncol(covariates)], varBetaHat = new("dsyMatrix", Dim = as.integer(rep(ncol(covariates), 2)), uplo = "L", x = fromOptim$betas[seq(1 + ncol(covariates), len = ncol(covariates)^2)])) result$optim$options = cbind(result$optim$options, gradient = fromOptim$betas[seq(ncol(covariates)^2 + } else { xcoord = as.vector(coordinates) ycoord = -99 } obsCov = cbind(y1 = observations, y2 = 0, y3 = 0, covariates) if (all(paramToEstimate == "variance") & param["nugget"] == 0) { theL = loglikLgm(param, data = observations, formula = covariates, coordinates = coordinates, reml = reml) fromOptim = attributes(theL) result = list(optim = list(mle = fillParam(fromOptim$param), logL = c(m2logL = as.numeric(theL), logL = -as.numeric(theL)/2), totalVarHat = fromOptim$totalVarHat, message = "numerical optimization not needed", options = NULL, detail = NULL), betaHat = fromOptim$betaHat, varBetaHat = fromOptim$varBetaHat) ncol(covariates) + 1, len = sum(Sparam))]) names(result$optim$detail) = c("fail", "fncount", "grcount") names(result$betaHat) = colnames(covariates) dimnames(result$varBetaHat) = list(names(result$betaHat), names(result$betaHat)) } result$parameters = fillParam(c(result$optim$mle, result$betaHat)) } else { fromOptim = .C(C_maternLogLOpt, start = as.double(paramsForC), Sparam = as.integer(Sparam), obsCov = as.double(as.matrix(obsCov)), as.double(xcoord), as.double(ycoord), as.integer(aniso), N = as.integer(c(nrow(obsCov), 3, ncol(covariates))), Ltype = as.integer(reml + 2 * !estimateVariance), optInt = as.integer(forO$scalarInt), optF = as.double(forO$scalarF), betas = as.double(forO$pars), limType = as.integer(forO$parsInt), message = format(" ", width = 80)) names(fromOptim$start) = names(paramsForC) if (fromOptim$start["anisoRatio"] < 1) { fromOptim$start["anisoRatio"] <- 1/fromOptim$start["anisoRatio"] if (fromOptim$start["anisoAngleRadians"] + pi/2 >= pi/2) { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] - pi/2 } if (estimateVariance) { result$parameters[c("nugget", "variance")] = result$parameters[c("nugget", "variance")] * result$optim$totalVarHat result$varBetaHat = result$varBetaHat * result$optim$totalVarHat } names(result$optim$logL) = paste(names(result$optim$logL), c(".ml", ".reml")[reml + 1], sep = "") else { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] + pi/2 } } result = list(optim = list(mle = fromOptim$start, logL = c(m2logL = fromOptim$optF[1], logL = -fromOptim$optF[1]/2), totalVarHat = fromOptim$optF[2], boxcox = fromOptim$optF[3:5], determinants = fromOptim$optF[6:7], message = fromOptim$message, options = cbind(start = paramsForC[Sparam], opt = fromOptim$start[Sparam], parOptions[, c("parscale", "lower", "upper", "ndeps")]), detail = fromOptim$optInt[1:3]), betaHat = fromOptim$betas[1:ncol(covariates)], varBetaHat = new("dsyMatrix", Dim = as.integer(rep(ncol(covariates), 2)), uplo = "L", x = fromOptim$betas[seq(1 + ncol(covariates), len = ncol(covariates)^2)])) result$optim$options = cbind(result$optim$options, gradient = fromOptim$betas[seq(ncol(covariates)^2 + ncol(covariates) + 1, len = sum(Sparam))]) names(result$optim$detail) = c("fail", "fncount", "grcount") result$data = cbind(data.frame(observations = observations, fitted = covariates %*% result$parameters[colnames(covariates)]), data.frame(data)[noNA, ]) if (abs(result$parameters["boxcox"] - 1) > 1e-04) { if (abs(result$parameters["boxcox"]) < 0.001) { result$data$obsBC = log(observations) } else { result$data$obsBC <- ((observations^result$parameters["boxcox"]) - 1)/result$parameters["boxcox"] } names(result$betaHat) = colnames(covariates) dimnames(result$varBetaHat) = list(names(result$betaHat), names(result$betaHat)) } result$parameters = fillParam(c(result$optim$mle, result$betaHat)) if (estimateVariance) { result$parameters[c("nugget", "variance")] = result$parameters[c("nugget", "variance")] * result$optim$totalVarHat result$varBetaHat = result$varBetaHat * result$optim$totalVarHat } names(result$optim$logL) = paste(names(result$optim$logL), c(".ml", ".reml")[reml + 1], sep = "") result$data = cbind(data.frame(observations = observations, fitted = covariates %*% result$parameters[colnames(covariates)]), data.frame(data)[noNA, ]) if (abs(result$parameters["boxcox"] - 1) > 1e-04) { if (abs(result$parameters["boxcox"]) < 0.001) { result$data$obsBC = log(observations) } else { result$data$obsBC <- obsCov[, 1] } result$data$resid = result$data$obsBC - result$data$fitted if (length(grep("^SpatVector", class(coordinatesOrig)))) { forDf = rep(NA, length(noNA)) forDf[noNA] = seq(1, sum(noNA)) theDf = result$data[forDf, ] result$data = vect(x = crds(coordinatesOrig), atts = theDf, crs = crs(coordinatesOrig)) } result$model = list(reml = reml, baseline = theFactors) if (any(class(trend) == "formula")) { result$model$formula = trend result$data[[all.vars(formula)[1]]] = result$data$observations } else { result$model$formula = names(trend) } } else { result$data$obsBC <- ((observations^result$parameters["boxcox"]) - 1)/result$parameters["boxcox"] } } else { result$data$obsBC <- obsCov[, 1] } result$data$resid = result$data$obsBC - result$data$fitted if (length(grep("^SpatVector", class(coordinatesOrig)))) { if (length(result$model$baseline)) { rparams = result$parameters baseParams = data.frame(var = names(result$model$baseline), base = result$model$baseline, pasted = paste(names(result$model$baseline), result$model$baseline, sep = "")) for (D in which(!baseParams$pasted %in% names(rparams))) { sameFac = grep(paste("^", baseParams[D, "var"], sep = ""), names(rparams)) pseq = 1:length(rparams) if (length(sameFac)) { minFac = min(sameFac) toAdd = 0 names(toAdd) = baseParams[D, "pasted"] rparams = c(rparams[pseq < minFac], toAdd, rparams[pseq >= minFac]) forDf = rep(NA, length(noNA)) forDf[noNA] = seq(1, sum(noNA)) theDf = result$data[forDf, ] result$data = vect(x = crds(coordinatesOrig), atts = theDf, crs = crs(coordinatesOrig)) } result$model = list(reml = reml, baseline = theFactors) if (any(class(trend) == "formula")) { result$model$formula = trend result$data[[all.vars(formula)[1]]] = result$data$observations } else { result$model$formula = names(trend) } } if (length(result$model$baseline)) { rparams = result$parameters baseParams = data.frame(var = names(result$model$baseline), base = result$model$baseline, pasted = paste(names(result$model$baseline), } result$parameters = rparams } parameterTable = data.frame(estimate = result$parameters) result$model$baseline, sep = "")) for (D in which(!baseParams$pasted %in% names(rparams))) { sameFac = grep(paste("^", baseParams[D, "var"], sep = ""), names(rparams)) pseq = 1:length(rparams) if (length(sameFac)) { minFac = min(sameFac) toAdd = 0 names(toAdd) = baseParams[D, "pasted"] rparams = c(rparams[pseq < minFac], toAdd, rparams[pseq >= minFac]) rownames(parameterTable) = names(result$parameters) parameterTable$stdErr = NA stdErr = sqrt(Matrix::diag(result$varBetaHat)) parameterTable[rownames(result$varBetaHat), "stdErr"] = stdErr thelims = c(0.005, 0.025, 0.05, 0.1) thelims = c(rbind(thelims, 1 - thelims)) } } result$parameters = rparams } parameterTable = data.frame(estimate = result$parameters) rownames(parameterTable) = names(result$parameters) parameterTable$stdErr = NA stdErr = sqrt(Matrix::diag(result$varBetaHat)) theQ = qnorm(thelims) toadd = outer(parameterTable$stdErr, theQ) toadd = toadd + matrix(parameterTable$estimate, ncol = length(thelims), nrow = dim(parameterTable)[1]) colnames(toadd) = paste("ci", thelims, sep = "") parameterTable = cbind(parameterTable, toadd) parameterTable[, "pval"] = pchisq(parameterTable$estimate^2/parameterTable$stdErr^2, df = 1, lower.tail = FALSE) parameterTable[rownames(result$varBetaHat), "stdErr"] = stdErr thelims = c(0.005, 0.025, 0.05, 0.1) thelims = c(rbind(thelims, 1 - thelims)) theQ = qnorm(thelims) parameterTable[, "Estimated"] = FALSE parameterTable[paramToEstimate, "Estimated"] = TRUE parameterTable[rownames(result$varBetaHat), "Estimated"] = TRUE toadd = outer(parameterTable$stdErr, theQ) toadd = toadd + matrix(parameterTable$estimate, ncol = length(thelims), nrow = dim(parameterTable)[1]) colnames(toadd) = paste("ci", thelims, sep = "") if (estimateVariance) parameterTable["variance", "Estimated"] = TRUE rownames(parameterTable) = gsub("^variance$", "sdSpatial", rownames(parameterTable)) rownames(parameterTable) = gsub("^nugget$", "sdNugget", rownames(parameterTable)) parameterTable = cbind(parameterTable, toadd) parameterTable[, "pval"] = pchisq(parameterTable$estimate^2/parameterTable$stdErr^2, df = 1, lower.tail = FALSE) parameterTable[, "Estimated"] = FALSE parameterTable[paramToEstimate, "Estimated"] = TRUE parameterTable[rownames(result$varBetaHat), "Estimated"] = TRUE if (estimateVariance) parameterTable["variance", "Estimated"] = TRUE rownames(parameterTable) = gsub("^variance$", "sdSpatial", rownames(parameterTable)) rownames(parameterTable) = gsub("^nugget$", "sdNugget", rownames(parameterTable)) parameterTable[c("sdSpatial", "sdNugget"), "estimate"] = sqrt(parameterTable[c("sdSpatial", "sdNugget"), "estimate"]) result$summary = as.data.frame(parameterTable) result$summary$Estimated = as.logical(result$summary$Estimated) parameterTable[c("sdSpatial", "sdNugget"), "estimate"] = sqrt(parameterTable[c("sdSpatial", "sdNugget"), "estimate"]) result$summary = as.data.frame(parameterTable) result$summary$Estimated = as.logical(result$summary$Estimated) result})(param = dots[[1L]][[1L]], data = list(fitted = c(7.09263093277653, 7.0632591348391, 7.03388733690166, 7.00451553896422, 6.97514374102678, 6.94577194308934, 6.9164001451519, 6.88702834721446, 6.85765654927703, 6.82828475133959, 6.79891295340215, 6.76954115546471, 6.74016935752727, 6.71079755958983, 6.68142576165239, 6.65205396371495, 6.62268216577752, 6.59331036784008, 6.56393856990264, 6.5345667719652, 6.50519497402776, 6.47582317609032, 6.44645137815288, 6.41707958021545, 6.38770778227801, result})(param = dots[[1L]][[1L]], data = list(fitted = c(7.09263093277653, 7.0632591348391, 7.03388733690166, 7.00451553896422, 6.97514374102678, 6.94577194308934, 6.9164001451519, 6.88702834721446, 6.85765654927703, 6.82828475133959, 6.79891295340215, 6.76954115546471, 6.74016935752727, 6.71079755958983, 6.68142576165239, 6.65205396371495, 6.62268216577752, 6.59331036784008, 6.56393856990264, 6.5345667719652, 6.50519497402776, 6.47582317609032, 6.44645137815288, 6.41707958021545, 6.38770778227801, 6.35833598434057, 6.32896418640313, 6.29959238846569, 6.27022059052825, 6.24084879259081, 6.21147699465337, 6.18210519671594, 6.1527333987785, 6.12336160084106, 6.09398980290362, 6.06461800496618, 6.03524620702874, 6.0058744090913, 5.97650261115387, 5.94713081321643, 5.91775901527899, 5.88838721734155, 5.85901541940411, 5.82964362146667, 5.80027182352923, 6.35833598434057, 6.32896418640313, 6.29959238846569, 6.27022059052825, 6.24084879259081, 6.21147699465337, 6.18210519671594, 6.1527333987785, 6.12336160084106, 6.09398980290362, 6.06461800496618, 6.03524620702874, 6.0058744090913, 5.97650261115387, 5.94713081321643, 5.91775901527899, 5.88838721734155, 5.85901541940411, 5.82964362146667, 5.80027182352923, 5.7709000255918, 5.74152822765436, 5.71215642971692, 5.68278463177948, 5.7709000255918, 5.74152822765436, 5.71215642971692, 5.68278463177948, 5.65341283384204, 5.6240410359046, 5.59466923796716, 5.56529744002973, 5.53592564209229, 5.50655384415485, 5.47718204621741, 5.44781024827997, 5.41843845034253, 5.38906665240509, 5.35969485446766, 5.33032305653022, 5.30095125859278, 5.27157946065534, 5.2422076627179, 5.21283586478046, 5.18346406684302, 5.15409226890558, 5.12472047096815, 5.09534867303071, 5.06597687509327, 5.03660507715583, 5.00723327921839, 4.97786148128095, 5.65341283384204, 5.6240410359046, 5.59466923796716, 5.56529744002973, 5.53592564209229, 5.50655384415485, 5.47718204621741, 5.44781024827997, 5.41843845034253, 5.38906665240509, 5.35969485446766, 5.33032305653022, 5.30095125859278, 5.27157946065534, 5.2422076627179, 5.21283586478046, 5.18346406684302, 5.15409226890558, 5.12472047096815, 5.09534867303071, 5.06597687509327, 5.03660507715583, 5.00723327921839, 4.97786148128095, 4.94848968334351, 4.91911788540608, 4.88974608746864, 4.8603742895312, 4.83100249159376, 4.80163069365632, 4.77225889571888, 4.74288709778144, 4.713515299844, 4.68414350190657, 4.65477170396913, 4.62539990603169, 4.59602810809425, 4.56665631015681, 4.53728451221937, 4.50791271428193, 4.4785409163445, 4.44916911840706, 4.41979732046962, 4.39042552253218, 4.36105372459474, 4.3316819266573, 4.30231012871986, 4.27293833078243, 4.24356653284499, 4.21419473490755, 4.18482293697011), ID = c(13L, 14L, 22L, 23L, 24L, 29L, 30L, 35L, 36L, 37L, 52L, 55L, 66L, 71L, 4.94848968334351, 4.91911788540608, 4.88974608746864, 4.8603742895312, 4.83100249159376, 4.80163069365632, 4.77225889571888, 4.74288709778144, 4.713515299844, 4.68414350190657, 4.65477170396913, 4.62539990603169, 4.59602810809425, 4.56665631015681, 4.53728451221937, 4.50791271428193, 4.4785409163445, 4.44916911840706, 4.41979732046962, 4.39042552253218, 4.36105372459474, 4.3316819266573, 4.30231012871986, 4.27293833078243, 4.24356653284499, 4.21419473490755, 4.18482293697011), ID = c(13L, 14L, 22L, 23L, 24L, 29L, 30L, 35L, 36L, 37L, 52L, 55L, 66L, 71L, 73L, 84L, 95L, 102L, 105L, 126L, 130L, 136L, 138L, 159L, 168L, 172L, 178L, 181L, 185L, 188L, 192L, 198L, 202L, 203L, 207L, 208L, 218L, 224L, 226L, 235L, 245L, 246L, 247L, 254L, 261L, 263L, 274L, 73L, 84L, 95L, 102L, 105L, 126L, 130L, 136L, 138L, 159L, 168L, 172L, 178L, 181L, 185L, 188L, 192L, 198L, 202L, 203L, 207L, 208L, 218L, 224L, 226L, 235L, 245L, 246L, 247L, 254L, 261L, 263L, 274L, 275L, 276L, 277L, 281L, 283L, 287L, 292L, 293L, 300L, 302L, 305L, 314L, 317L, 320L, 322L, 335L, 336L, 341L, 342L, 344L, 357L, 362L, 368L, 369L, 372L, 373L, 377L, 378L, 381L, 384L, 392L, 399L, 400L, 401L, 406L, 408L, 409L, 421L, 423L, 424L, 425L, 436L, 442L, 449L, 450L, 451L, 455L, 456L, 458L, 460L, 466L, 468L, 471L), rain = c(15.1, 25.5, 7.9, 19.1, 19.4, 33.4, 10.7, 29.6, 39.4, 39.4, 32.4, 10.5, 13.5, 58.5, 11.4, 33.4, 13.1, 7.8, 39.8, 14.1, 19.2, 15.1, 10.7, 14.5, 33.4, 32.7, 21.3, 33.1, 40, 32.7, 38, 9.4, 18.5, 23.9, 33, 3, 25.4, 5.3, 7.8, 7.1, 6.2, 7.1, 5.9, 6, 12.4, 15.3, 7.5, 13.7, 8.6, 12.9, 34.5, 44.1, 18.4, 12.1, 34.6, 27, 10, 4.5, 10.7, 35.9, 27.8, 7.2, 13.1, 9, 14.1, 13.1, 45.2, 1.6, 13.6, 13, 11.8, 10.9, 14.5, 25.4, 14, 15.2, 6, 28.3, 18.4, 12.7, 22, 17.8, 21.8, 13.7, 14.4, 23, 28.2, 12.9, 6.5, 19, 17, 15.6, 13.1, 1, 9.9, 9.2, 6.7, 1.8, 2, 5.5), elev = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 275L, 276L, 277L, 281L, 283L, 287L, 292L, 293L, 300L, 302L, 305L, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100)), formula = rain ~ elev, 314L, 317L, 320L, 322L, 335L, 336L, 341L, 342L, 344L, 357L, 362L, 368L, 369L, 372L, 373L, 377L, 378L, 381L, 384L, 392L, 399L, 400L, 401L, 406L, 408L, 409L, 421L, 423L, 424L, 425L, 436L, 442L, 449L, 450L, 451L, 455L, 456L, 458L, 460L, 466L, 468L, 471L), rain = c(15.1, 25.5, 7.9, 19.1, 19.4, 33.4, 10.7, 29.6, 39.4, 39.4, 32.4, 10.5, 13.5, 58.5, 11.4, 33.4, 13.1, 7.8, 39.8, 14.1, 19.2, 15.1, 10.7, 14.5, 33.4, 32.7, 21.3, 33.1, 40, 32.7, 38, 9.4, 18.5, 23.9, paramToEstimate = c("variance", "anisoRatio", "anisoAngleRadians" ), reml = FALSE, coordinates = c(2514336.71, 2518588.71, 2531615.71, 2533523.71, 2534125.71, 2538019.71, 2539319.71, 2545119.71, 2545607.71, 2545627.71, 2561259.71, 2562512.71, 2570067.71, 2571109.71, 2571016.71, 2579205.71, 2586092.71, 2590945.71, 2592122.71, 2608356.71, 2609399.71, 2612205.71, 2612286.71, 2620667.71, 2624915.71, 2626374.71, 2628041.71, 2629488.71, 2630082.71, 2630821.71, 2633239.71, 2635189.71, 33, 3, 25.4, 5.3, 7.8, 7.1, 6.2, 7.1, 5.9, 6, 12.4, 15.3, 7.5, 13.7, 8.6, 12.9, 34.5, 44.1, 18.4, 12.1, 34.6, 27, 10, 4.5, 10.7, 35.9, 27.8, 7.2, 13.1, 9, 14.1, 13.1, 45.2, 1.6, 13.6, 13, 11.8, 10.9, 14.5, 25.4, 14, 15.2, 6, 28.3, 18.4, 12.7, 22, 17.8, 21.8, 13.7, 14.4, 23, 28.2, 12.9, 6.5, 19, 17, 15.6, 13.1, 1, 9.9, 9.2, 6.7, 1.8, 2, 5.5), elev = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100)), formula = rain ~ elev, 2637581.71, 2637613.71, 2640771.71, 2641232.71, 2648275.71, 2655053.71, 2656589.71, 2661121.71, 2668821.71, 2668770.71, 2670977.71, 2676213.71, 2679163.71, 2680537.71, 2685168.71, 2686062.71, 2686015.71, 2687100.71, 2687949.71, 2688944.71, 2688673.71, 2692431.71, 2694076.71, 2696179.71, 2696965.71, 2699592.71, 2700787.71, 2702527.71, 2704185.71, 2704251.71, 2709024.71, 2711362.71, 2710661.71, 2710650.71, 2713957.71, 2718384.71, 2717971.71, 2720517.71, 2720255.71, 2721055.71, 2723237.71, 2726025.71, 2725583.71, 2726727.71, 2728111.71, 2733095.71, 2736658.71, 2736759.71, 2739493.71, 2740339.71, 2742780.71, 2741684.71, 2749700.71, 2751519.71, 2751266.71, paramToEstimate = c("variance", "anisoRatio", "anisoAngleRadians" ), reml = FALSE, coordinates = c(2514336.71, 2518588.71, 2531615.71, 2533523.71, 2534125.71, 2538019.71, 2539319.71, 2545119.71, 2545607.71, 2545627.71, 2561259.71, 2562512.71, 2570067.71, 2571109.71, 2571016.71, 2579205.71, 2586092.71, 2590945.71, 2592122.71, 2608356.71, 2609399.71, 2612205.71, 2612286.71, 2620667.71, 2624915.71, 2626374.71, 2628041.71, 2629488.71, 2630082.71, 2630821.71, 2633239.71, 2635189.71, 2637581.71, 2637613.71, 2640771.71, 2641232.71, 2648275.71, 2655053.71, 2656589.71, 2661121.71, 2668821.71, 2668770.71, 2752986.71, 2759987.71, 2761869.71, 2767951.71, 2766435.71, 2769345.71, 2775690.71, 2777669.71, 2779806.71, 2782895.71, 2796876.71, 2800101.71, 2805720.71, 1156023.34, 1174834.34, 1177889.34, 1196758.34, 1188960.34, 1154404.34, 1182184.34, 1207655.34, 1150925.34, 1152037.34, 1172904.34, 1252956.34, 1255067.34, 1168310.34, 1106035.34, 1204901.34, 1143655.34, 1094673.34, 1206974.34, 1149002.34, 1185690.34, 1257948.34, 1269067.34, 1160040.34, 1244526.34, 1234511.34, 1268974.34, 1256736.34, 1218927.34, 1214477.34, 1254497.34, 1161086.34, 1195550.34, 1206669.34, 1258922.34, 1129935.34, 1248902.34, 1185517.34, 1135480.34, 1203312.34, 1154399.34, 1176638.34, 2670977.71, 2676213.71, 2679163.71, 2680537.71, 2685168.71, 2686062.71, 2686015.71, 2687100.71, 2687949.71, 2688944.71, 2688673.71, 2692431.71, 2694076.71, 2696179.71, 2696965.71, 2699592.71, 2700787.71, 2702527.71, 2704185.71, 2704251.71, 2709024.71, 2711362.71, 2710661.71, 2710650.71, 2713957.71, 2718384.71, 2717971.71, 2720517.71, 2720255.71, 2721055.71, 2723237.71, 2726025.71, 2725583.71, 2726727.71, 2728111.71, 2733095.71, 2736658.71, 2736759.71, 2739493.71, 2740339.71, 2742780.71, 2741684.71, 2749700.71, 2751519.71, 2751266.71, 2752986.71, 2759987.71, 2761869.71, 2767951.71, 2766435.71, 1205554.34, 1225586.34, 1243389.34, 1273417.34, 1247865.34, 1221182.34, 1230078.34, 1168925.34, 1153362.34, 1112225.34, 1292361.34, 1289049.34, 1152287.34, 1180101.34, 1282408.34, 1232389.34, 1272429.34, 1147900.34, 1132346.34, 1216858.34, 1273611.34, 1103497.34, 1259171.34, 1260283.34, 1153562.34, 1095782.34, 1278149.34, 1251488.34, 1274838.34, 1270399.34, 1211486.34, 1168149.34, 1268227.34, 1235992.34, 1246018.34, 1152668.34, 1227225.34, 1273933.34, 1187230.34, 1233949.34, 1170596.34, 1245090.34, 1261895.34, 1196309.34, 1211875.34, 1152960.34, 1146405.34, 1169794.34, 1171018.34, 1251067.34, 1137679.34, 1165608.34, 1143403.34, 1187937.34, 1185778.34, 2769345.71, 2775690.71, 2777669.71, 2779806.71, 2782895.71, 2796876.71, 2800101.71, 2805720.71, 1156023.34, 1174834.34, 1177889.34, 1196758.34, 1188960.34, 1154404.34, 1182184.34, 1207655.34, 1150925.34, 1152037.34, 1172904.34, 1252956.34, 1255067.34, 1168310.34, 1106035.34, 1204901.34, 1143655.34, 1094673.34, 1206974.34, 1149002.34, 1185690.34, 1257948.34, 1269067.34, 1160040.34, 1244526.34, 1234511.34, 1268974.34, 1256736.34, 1218927.34, 1214477.34, 1254497.34, 1161086.34, 1195550.34, 1206669.34, 1258922.34, 1129935.34, 1248902.34, 1185517.34, 1135480.34, 1203312.34, 1154399.34, 1176638.34, 1205554.34, 1225586.34, 1243389.34, 1273417.34, 1247865.34, 1221182.34, 1230078.34, 1168925.34, 1153362.34, 1112225.34, 1292361.34, 1289049.34, 1152287.34, 1180101.34, 1282408.34, 1232389.34, 1272429.34, 1147900.34, 1132346.34, 1216858.34, 1141601.34, 1135004.34, 1125130.34)) 2: .mapply(FUN, dots, MoreArgs) 3: FUN(X[[i]], ...) 4: lapply(X = S, FUN = FUN, ...) 5: doTryCatch(return(expr), name, parentenv, handler) 6: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 1273611.34, 1103497.34, 1259171.34, 1260283.34, 1153562.34, 1095782.34, 1278149.34, 1251488.34, 1274838.34, 1270399.34, 1211486.34, 1168149.34, 1268227.34, 1235992.34, 1246018.34, 1152668.34, 1227225.34, 1273933.34, 1187230.34, 1233949.34, 1170596.34, 1245090.34, 1261895.34, 1196309.34, 1211875.34, 1152960.34, 1146405.34, 1169794.34, 1171018.34, 1251067.34, 1137679.34, 1165608.34, 1143403.34, 1187937.34, 1185778.34, 1141601.34, 1135004.34, 1125130.34)) 2: .mapply(FUN, dots, MoreArgs) 7: tryCatchList(expr, classes, parentenv, handlers) 8: tryCatch(expr, error = function(e) { call <- conditionCall(e) if (!is.null(call)) { if (identical(call[[1L]], quote(doTryCatch))) call <- sys.call(-4L) dcall <- deparse(call, nlines = 1L) prefix <- paste("Error in", dcall, ": ") LONG <- 75L sm <- strsplit(conditionMessage(e), "\n")[[1L]] w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w") if (is.na(w)) w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], type = "b") if (w > LONG) prefix <- paste0(prefix, "\n ") } else prefix <- "Error : " msg <- paste0(prefix, conditionMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && isTRUE(getOption("show.error.messages"))) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))}) 9: try(lapply(X = S, FUN = FUN, ...), silent = TRUE) 10: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) 3: FUN(X[[i]], ...) 4: lapply(X = S, FUN = FUN, ...) 5: doTryCatch(return(expr), name, parentenv, handler)11: FUN(X[[i]], ...) 12: lapply(seq_len(cores), inner.do) 13: mclapply(seq_len(n), do_one, mc.preschedule = mc.preschedule, mc.set.seed = mc.set.seed, mc.silent = mc.silent, mc.cores = mc.cores, mc.cleanup = mc.cleanup, affinity.list = affinity.list) 14: parallel::mcmapply(likfitLgm, param = parList, MoreArgs = list(data = as.data.frame(fit$data), formula = fit$model$formula, paramToEstimate = reEstimate, reml = fit$model$reml, coordinates = terra::crds(fit$data)), mc.cores = mc.cores, SIMPLIFY = FALSE) 15: profLlgm(swissFit, mc.cores = Ncores, range = seq(15000, 55000, len = 12)) An irrecoverable exception occurred. R is aborting now ... 6: tryCatchOne(expr, names, parentenv, handlers[[1L]]) 7: tryCatchList(expr, classes, parentenv, handlers) 8: tryCatch(expr, error = function(e) { call <- conditionCall(e) if (!is.null(call)) { if (identical(call[[1L]], quote(doTryCatch))) call <- sys.call(-4L) dcall <- deparse(call, nlines = 1L) prefix <- paste("Error in", dcall, ": ") LONG <- 75L sm <- strsplit(conditionMessage(e), "\n")[[1L]] w <- 14L + nchar(dcall, type = "w") + nchar(sm[1L], type = "w") if (is.na(w)) w <- 14L + nchar(dcall, type = "b") + nchar(sm[1L], type = "b") if (w > LONG) prefix <- paste0(prefix, "\n ") } else prefix <- "Error : " msg <- paste0(prefix, conditionMessage(e), "\n") .Internal(seterrmessage(msg[1L])) if (!silent && isTRUE(getOption("show.error.messages"))) { cat(msg, file = outFile) .Internal(printDeferredWarnings()) } invisible(structure(msg, class = "try-error", condition = e))}) 9: try(lapply(X = S, FUN = FUN, ...), silent = TRUE) 10: sendMaster(try(lapply(X = S, FUN = FUN, ...), silent = TRUE)) 11: FUN(X[[i]], ...) 12: lapply(seq_len(cores), inner.do) 13: mclapply(seq_len(n), do_one, mc.preschedule = mc.preschedule, mc.set.seed = mc.set.seed, mc.silent = mc.silent, mc.cores = mc.cores, mc.cleanup = mc.cleanup, affinity.list = affinity.list) 14: parallel::mcmapply(likfitLgm, param = parList, MoreArgs = list(data = as.data.frame(fit$data), formula = fit$model$formula, paramToEstimate = reEstimate, reml = fit$model$reml, coordinates = terra::crds(fit$data)), mc.cores = mc.cores, SIMPLIFY = FALSE) 15: profLlgm(swissFit, mc.cores = Ncores, range = seq(15000, 55000, len = 12)) An irrecoverable exception occurred. R is aborting now ... Error in resL[grep("^m2logL", rownames(resL)), ] : incorrect number of dimensions Calls: profLlgm In addition: Warning message: In mclapply(seq_len(n), do_one, mc.preschedule = mc.preschedule, : scheduled cores 1, 2 did not deliver results, all values of the jobs will be affected Execution halted Flavor: r-release-macos-arm64

Version: 2.0.11
Check: package dependencies
Result: NOTE Package suggested but not available for checking: ‘RandomFields’ Package which this enhances but not available for checking: ‘INLA’ Flavors: r-oldrel-macos-arm64, r-oldrel-macos-x86_64

Version: 2.0.11
Check: package dependencies
Result: NOTE Package suggested but not available for checking: 'RandomFields' Flavor: r-oldrel-windows-x86_64