Last updated on 2026-02-28 11:49:24 CET.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| 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 | 190.00 | 199.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 |
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’ [39s/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 20:35:18 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) 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 { 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 } 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") 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 <- ((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")) { 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 >= 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)) 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})(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 ...
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) 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 { 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 } 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") 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 <- ((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")) { 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 >= 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)) 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})(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’ [3s/3s]
Running ‘lgm.R’ [7s/7s]
Running ‘lgmRaster.R’ [49s/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 20:16:40 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") }
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)] } if (any(class(coordinates) == "dist")) coordinates = as(as.matrix(coordinates), "dsyMatrix") if (length(grep("SpatVect", class(coordinates)))) { coordinatesOrig = 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") 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)) { } 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)))) { 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)) 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"] } else { if (is.matrix(coordinates)) { if (ncol(coordinates) == 2) { coordinates = dist(coordinates) } param = param[names(param) != "variance"] } else { if (any(names(param) == "variance")) { else { coordinates = as(coordinates, "dsyMatrix") } } if (any(class(coordinates) == "dist")) coordinates = as(as.matrix(coordinates), "dsyMatrix") estimateVariance = FALSE } } paramToEstimate = gsub("anisoAngleDegrees", "anisoAngleRadians", paramToEstimate) degToRad = function(par) { if (length(grep("anisoAngleDegrees", names(par)))) { if (!length(grep("anisoAngleRadians", names(par)))) if (length(grep("SpatVect", class(coordinates)))) { if (!nchar(crs(coordinates))) terra::crs(coordinates) = "+proj=utm +zone=1" coordinates = new("dsyMatrix", Dim = rep(length(coordinates), par["anisoAngleRadians"] = 2 * pi * par["anisoAngleDegrees"]/360 par = par[grep("anisoAngleDegrees", names(par), invert = TRUE)] } par } 2), x = as.vector(as.matrix(terra::distance(coordinates))), uplo = "L") } maxDist = max(coordinates, na.rm = TRUE) 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 } } 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))) 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 noNA = !theNA theFactors = model.frame(trend, data[noNA, ]) whichFactors = unlist(lapply(theFactors, is.factor)) theFactors = theFactors[, whichFactors, drop = FALSE] 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]] 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), 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] collapse = "|"), ")NA$", sep = ""), colnames(covariates), invert = TRUE), drop = FALSE] observations = all.vars(trend)[1] if (!any(names(data) == observations)) 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[, warning("can't find observations ", observations, "in data") observations = data[noNA, observations, drop = FALSE] } else { trend = as.matrix(trend) 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"] == 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] 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 } else { if (ncol(coordinates) == nrow(coordinates)) { coordinates = coordinates[noNA, noNA] } else { coordinates = coordinates[noNA, ] } } } } obsCov = cbind(y1 = observations, y2 = 0, y3 = 0, covariates) if (all(paramToEstimate == "variance") & param["nugget"] == 0) { estimateVariance = TRUE if (any(paramToEstimate == "variance")) { paramToEstimate = paramToEstimate[paramToEstimate != "variance"] param = param[names(param) != "variance"] } else { 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 { if (any(names(param) == "variance")) { estimateVariance = FALSE } } paramToEstimate = gsub("anisoAngleDegrees", "anisoAngleRadians", 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) { 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)] fromOptim$start["anisoRatio"] <- 1/fromOptim$start["anisoRatio"] if (fromOptim$start["anisoAngleRadians"] + pi/2 >= pi/2) { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] - } 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, pi/2 } else { fromOptim$start["anisoAngleRadians"] <- fromOptim$start["anisoAngleRadians"] + 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, 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", anisoAngleRadians = 0.01, anisoRatio = 0.01, variance = 0.01, shape = 0.01) paramDefaults[names(param)] = param parscaleDefaults[names(parscale)] = parscale lowerDefaults[names(lower)] = lower 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") names(result$betaHat) = colnames(covariates) dimnames(result$varBetaHat) = list(names(result$betaHat), names(result$betaHat)) 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) } 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), 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) c(".ml", ".reml")[reml + 1], sep = "") result$data = cbind(data.frame(observations = observations, fitted = covariates %*% result$parameters[colnames(covariates)]), data.frame(data)[noNA, ]) 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, 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"] 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 } } } else { result$data$obsBC <- obsCov[, 1] } 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)) result$data$resid = result$data$obsBC - result$data$fitted 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) } 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)) } 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) { 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)) { minFac = min(sameFac) toAdd = 0 names(toAdd) = baseParams[D, "pasted"] rparams = c(rparams[pseq < minFac], toAdd, rparams[pseq >= 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", 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)) 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})(param = dots[[1L]][[1L]], data = list(fitted = c(7.09263093277653, 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") 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$totalVarHat7.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, 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 <- ((observations^result$parameters["boxcox"]) - 1)/result$parameters["boxcox"] } } else {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, 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")) {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, 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)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, 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) 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 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, 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})(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, 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, 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, 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))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,
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) { 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, 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)) 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, 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) 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: .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))
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) {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)
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 ")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, 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 ...
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