[Georob-commits] r18 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 15 18:20:15 CEST 2014
Author: papritz
Date: 2014-05-15 18:20:15 +0200 (Thu, 15 May 2014)
New Revision: 18
Added:
pkg/man/georobModelBuilding.Rd
Modified:
pkg/ChangeLog
pkg/DESCRIPTION
pkg/NEWS
pkg/R/georob.S3methods.R
pkg/R/georob.cv.R
pkg/R/georob.exported.functions.R
pkg/R/georob.predict.R
pkg/R/georob.private.functions.R
pkg/R/variogram.R
pkg/man/fit.variogram.model.Rd
pkg/man/georob-package.Rd
pkg/man/georob.Rd
pkg/man/georob.control.Rd
pkg/man/lgnpp.Rd
pkg/man/param.names.Rd
pkg/man/plot.georob.Rd
pkg/man/predict.georob.Rd
pkg/man/sample.variogram.Rd
Log:
version 0-1.3
Modified: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog 2014-03-05 13:13:43 UTC (rev 17)
+++ pkg/ChangeLog 2014-05-15 16:20:15 UTC (rev 18)
@@ -188,5 +188,37 @@
* georob.private.functions.R (compute.semivariance, gar): changes for RandomFields version 3
* georob.predict.R (predict.georobcd R. ): changes for RandomFields version 3
+* georob.S3methods.R (add1.georob, drop1.georob, step.georob): changes in functions for stepwise selection of fixed-effects terms
+* georob.cv.R (cv.georob): catching attempt to re-estimate variogram parameters when all parameters in object are fixed
+2014-03-12 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.S3methods.R (add1.georob, deviance.georob, drop1.georob, extractAIC.georob, logLik.georob, step.georob): changes in functions for parallelized stepwise selection of fixed-effects terms
+
+
+2014-04-23 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.predict.R (predict.georob): correcting error when newdata contains data locations
+
+
+2014-04-23 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.S3methods.R (add1.georob, drop1.georob, step.georob): changes in functions for reducing memory demand
+
+
+2014-04-23 Andreas Papritz <papritz at env.ethz.ch>
+
+* variogram.R (plot.georob,lines.georob): changes for plotting covariances and correlations
+
+
+2014-05-15 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.cv.R (cv.georob): changes for version 3 of RandomFields
+* georob.exported.functions.R (georob, georob.control, param.transf, param.names, param.bounds): changes for version 3 of RandomFields
+* georob.private.functions.R (gcr, prepare.likelihood.calculations, dcorr.dparam, compute.semivariance): changes for version 3 of RandomFields
+* georob.S3methods.R (waldtest.georob, add1.georob, drop1.georob): changes for version 3 of RandomFields
+* variogram.R (fit.variogram.model): changes for version 3 of RandomFields
+
+
+
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2014-03-05 13:13:43 UTC (rev 17)
+++ pkg/DESCRIPTION 2014-05-15 16:20:15 UTC (rev 18)
@@ -1,8 +1,8 @@
Package: georob
Type: Package
Title: Robust Geostatistical Analysis of Spatial Data
-Version: 0.1-2
-Date: 2013-09-06
+Version: 0.1-4
+Date: 2014-05-15
Authors at R: c(
person( "Andreas", "Papritz", role = c( "cre", "aut" ),
email = "andreas.papritz at env.ethz.ch" ),
Modified: pkg/NEWS
===================================================================
--- pkg/NEWS 2014-03-05 13:13:43 UTC (rev 17)
+++ pkg/NEWS 2014-05-15 16:20:15 UTC (rev 18)
@@ -1,3 +1,4 @@
2012-12-14 georob Version 0.1-0 released!
2013-09-06 georob Version 0.1-1 released to CRAN
+2014-05-15 georob Version 0.1-3 released to CRAN
Modified: pkg/R/georob.S3methods.R
===================================================================
--- pkg/R/georob.S3methods.R 2014-03-05 13:13:43 UTC (rev 17)
+++ pkg/R/georob.S3methods.R 2014-05-15 16:20:15 UTC (rev 18)
@@ -831,6 +831,7 @@
## 2012-02-08 AP change for anisotropic variograms
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-02-19 AP correcting option for verbose output
+ ## 2014-05-15 AP changes for version 3 of RandomFields
test <- match.arg( test )
@@ -846,8 +847,8 @@
aniso = object[["aniso"]][["aniso"]],
fit.param = c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
- a = FALSE, alpha = FALSE, beta = FALSE, delta = FALSE,
- gamma = FALSE, lambda = FALSE, n = FALSE, nu = FALSE
+ alpha = FALSE, beta = FALSE, delta = FALSE,
+ gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[names( object[["param"]] )],
fit.aniso = c(
f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE
@@ -879,7 +880,7 @@
val <- if( REML ){
val <- object[["loglik"]]
} else {
- D <- deviance( object )
+ D <- deviance( object, ... )
-0.5 * (
D + attr( D, "log.det.covmat" ) + length( object[["residuals"]] ) * log( 2 * pi )
)
@@ -900,7 +901,7 @@
deviance.georob <-
function(
- object, ...
+ object, warn = TRUE, ...
)
{
## deviance method for class georob
@@ -910,6 +911,7 @@
## 2013-05-31 AP revised expansion of covariance matrices
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2014-02-27 AP computing 'pseudo' deviance for robustly fitted models
+ ## 2014-03-12 AP option for suppressing warnings
## redefine na.action component of object
@@ -920,7 +922,7 @@
## determine factor to compute heteroscedastic nugget
if( object[["tuning.psi"]] < georob.control()[["tuning.psi.nr"]] ){
- warning(
+ if( warn ) warning(
"deviance for robustly fitted model approximated by deviance of\n",
" equivalent model with heteroscedastic nugget"
)
@@ -967,7 +969,7 @@
edf <- sum( !is.na( fit[["coefficients"]] ) ) +
sum( fit[["initial.objects"]][["fit.param"]] ) +
sum( fit[["initial.objects"]][["fit.aniso"]] )
- loglik <- logLik( fit )
+ loglik <- logLik( fit, ... )
c(edf, -2 * loglik + k * edf)
}
@@ -986,84 +988,177 @@
## ##############################################################################
add1.georob <- function( object, scope, scale = 0, test=c("none", "Chisq" ),
- k = 2, trace = FALSE, fixed = TRUE, use.fitted.param = TRUE, verbose = 0, ... )
+ k = 2, trace = FALSE, data = NULL, fixed = TRUE, use.fitted.param = TRUE, verbose = 0,
+ ncores = 1, ... )
{
## add1 method for class georob based on add1.default{stats}
- ## 2014-03-01 AP
+ ## 2014-03-12 AP
+ ## 2014-04-24 AP function returns only variogram parameters of best fitted object
+ ## 2014-05-15 AP changes for version 3 of RandomFields
- ## check scope
+ ## auxiliary function for fitting model and extracting aic
- if( missing(scope ) || is.null( scope ) ) stop("no terms in scope")
- if( !is.character( scope ) ) scope <- add.scope( object, update.formula(object, scope) )
- if( !length(scope)) stop( "no terms in scope for adding to object" )
-
- ## initialize result
-
- ns <- length( scope )
- ans <- matrix(
- nrow = ns + 1L, ncol = 2L, dimnames = list(c("<none>", scope), c("df", "AIC"))
- )
- ans[1L, ] <- extractAIC( object, scale, k = k, ... )
-
- ## loop over all components of scope
-
- n0 <- nobs( object, use.fallback = TRUE )
- env <- environment( formula(object) )
-
- all.converged <- TRUE
-
- for( i in seq(ns) ) {
- tt <- scope[i]
+ f.aux <- function(
+ tt, scale, k, trace,
+ fixed, use.fitted.param, object, data, param, aniso, fit.param, fit.aniso,
+ verbose, ...
+ ){
+ converged <- TRUE
if(trace > 1) {
- cat("trying +", tt, "\n", sep = "")
+ cat("\ntrying +", tt, "\n", sep = "")
utils::flush.console()
}
- ## was:
- ## nfit <- update( object, as.formula(paste("~ . +", tt)), evaluate = FALSE )
- ## nfit <- eval( nfit, envir=env ) # was eval.parent(nfit)
if( fixed ){
- nfit <- update(
- object, as.formula(paste("~ . +", tt)),
- param = object[["param"]],
- aniso = object[["aniso"]][["aniso"]],
- fit.param = c(
- variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
- a = FALSE, alpha = FALSE, beta = FALSE, delta = FALSE,
- gamma = FALSE, lambda = FALSE, n = FALSE, nu = FALSE
- )[names( object[["param"]] )],
- fit.aniso = c(
- f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE
- ),
- verbose = verbose
- )
- } else {
- if( use.fitted.param ){
+ if( is.null( data ) ){
nfit <- update(
object, as.formula(paste("~ . +", tt)),
- param = object[["param"]], aniso = object[["aniso"]][["aniso"]],
+ param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
initial.param = "no", verbose = verbose
)
} else {
- nfit <- update( object, as.formula(paste("~ . +", tt)), verbose = verbose )
+ nfit <- update(
+ object, as.formula(paste("~ . +", tt)), data = data,
+ param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
+ initial.param = "no", verbose = verbose
+ )
}
+ } else {
+ if( use.fitted.param ){
+ if( is.null( data ) ){
+ nfit <- update(
+ object, as.formula(paste("~ . +", tt)),
+ param = param, aniso = aniso, initial.param = "no", verbose = verbose
+ )
+ } else {
+ nfit <- update(
+ object, as.formula(paste("~ . +", tt)), data = data,
+ param = param, aniso = aniso, initial.param = "no", verbose = verbose
+ )
+ }
+ }
+ else {
+ if( is.null( data ) ){
+ nfit <- update( object, as.formula(paste("~ . +", tt)),
+ verbose = verbose
+ )
+ } else {
+ nfit <- update( object, as.formula(paste("~ . +", tt)), data = data,
+ verbose = verbose
+ )
+ }
+ }
if( !nfit$converged ){
- all.converged <- FALSE
+ converged <- FALSE
if( verbose > 0 ) cat(
- "lack of convergence when fitting model", as.character(formula(nfit)),
+ "lack of convergence when fitting model", paste("~ . +", tt),
"\nconvergence code:", nfit$convergence.code, "\n"
)
}
}
- ans[i+1L, ] <- extractAIC( nfit, scale, k = k, ... )
nnew <- nobs( nfit, use.fallback = TRUE )
if( all(is.finite(c(n0, nnew))) && nnew != n0 ) stop(
"number of rows in use has changed: remove missing values?"
)
+ df.aic <- c( extractAIC( nfit, scale, k = k, ... ), as.numeric(converged))
+ names( df.aic ) <- c("df", "AIC", "converged")
+ attr( df.aic, "nfit" ) <- list(
+ param = nfit[["param"]],
+ aniso = nfit[["aniso"]][["aniso"]],
+ initial.objects = nfit[["initial.objects"]],
+ call = nfit[["call"]]
+ )
+ df.aic
}
- if( !all.converged ) warning( "lack of convergence when fitting models" )
+ ## check scope
+
+ if( missing(scope ) || is.null( scope ) ) stop("no terms in scope")
+ if( !is.character( scope ) ) scope <- add.scope( object, update.formula(object, scope) )
+ if( !length(scope)) stop( "no terms in scope for adding to object" )
+
+ ## initialize result
+
+ ns <- length( scope )
+ ans <- matrix(
+ nrow = ns + 1L, ncol = 2L, dimnames = list(c("<none>", scope), c("df", "AIC"))
+ )
+ ans[1L, ] <- extractAIC( object, scale, k = k, ... )
+
+ ## loop over all components of scope
+
+ n0 <- nobs( object, use.fallback = TRUE )
+ env <- environment( formula(object) )
+
+ param <- object[["param"]]
+ aniso <- object[["aniso"]][["aniso"]]
+ fit.param <- c(
+ variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
+ alpha = FALSE, beta = FALSE, delta = FALSE,
+ gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
+ )[names( object[["param"]] )]
+ fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
+
+ ## prepare for parallel execution
+
+ if( ncores > 1 ){
+ ncores <- min( ncores, ns, detectCores() )
+ trace <- 0
+ verbose <- 0
+ }
+
+ if( .Platform[["OS.type"]] == "windows" && ncores > 1 ){
+
+ if( is.null( data ) ) stop(
+ "argument 'data' required for parallel execution on windows OS"
+ )
+
+ ## create a SNOW cluster on windows OS
+
+ cl <- makePSOCKcluster( ncores, outfile = "")
+
+ ## export required items to workers
+
+ junk <- clusterEvalQ( cl, require( georob, quietly = TRUE ) )
+
+ result <- parLapply(
+ cl,
+ scope, f.aux,
+ scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
+ object = object, data = data, param = param, aniso = aniso,
+ fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose, ...
+ )
+
+ stopCluster(cl)
+
+ } else {
+
+ ## fork child processes on non-windows OS
+
+ result <- mclapply(
+ scope, f.aux,
+ scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
+ object = object, data = data, param = param, aniso = aniso,
+ fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose,
+ mc.cores = ncores, mc.allow.recursive = FALSE, ...
+ )
+
+ }
+
+ fits <- list(
+ list(
+ param = object[["param"]], aniso = object[["aniso"]][["aniso"]],
+ initial.objects = object[["initial.objects"]],
+ call = object[["call"]]
+ )
+ )
+ fits[2:nrow(ans)] <- lapply( result, function(x) attr( x, "nfit" ) )
+ result <- t( simplify2array( result ) )
+ ans[2:NROW(ans), ] <- result[, 1:2]
+ if(!all( as.logical(result[, "converged"]) ) ) warning(
+ "lack of convergence when fitting models"
+ )
## store results
@@ -1091,21 +1186,99 @@
)
class(aod) <- c("anova", "data.frame")
attr( aod, "heading") <- head
+ attr( aod, "nfit" ) <- fits[[which.min( aod[, "AIC"])]]
+
aod
}
-
## ##############################################################################
drop1.georob <- function( object, scope, scale = 0, test=c( "none", "Chisq" ),
- k = 2, trace = FALSE, fixed = TRUE, use.fitted.param = TRUE, verbose = 0, ... )
+ k = 2, trace = FALSE, data = NULL, fixed = TRUE, use.fitted.param = TRUE, verbose = 0,
+ ncores = 1, ... )
{
## drop1 method for class georob based on drop1.default{stats}
- ## 2014-02-27 AP
+ ## 2014-03-12 AP
+ ## 2014-04-24 AP function returns only variogram parameters of best fitted object
+ ## 2014-05-15 AP changes for version 3 of RandomFields
+ ## auxiliary function for fitting model and extracting aic
+
+ f.aux <- function(
+ tt, scale, k, trace,
+ fixed, use.fitted.param, object, data, param, aniso, fit.param, fit.aniso,
+ verbose, ...
+ ){
+ converged <- TRUE
+ if(trace > 1) {
+ cat("\ntrying -", tt, "\n", sep = "")
+ utils::flush.console()
+ }
+ if( fixed ){
+ if( is.null( data ) ){
+ nfit <- update(
+ object, as.formula(paste("~ . -", tt)),
+ param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
+ initial.param = "no", verbose = verbose
+ )
+ } else {
+ nfit <- update(
+ object, as.formula(paste("~ . -", tt)), data = data,
+ param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
+ initial.param = "no", verbose = verbose
+ )
+ }
+ } else {
+ if( use.fitted.param ){
+ if( is.null( data ) ){
+ nfit <- update(
+ object, as.formula(paste("~ . -", tt)),
+ param = param, aniso = aniso, initial.param = "no", verbose = verbose
+ )
+ } else {
+ nfit <- update(
+ object, as.formula(paste("~ . -", tt)), data = data,
+ param = param, aniso = aniso, initial.param = "no", verbose = verbose
+ )
+ }
+ }
+ else {
+ if( is.null( data ) ){
+ nfit <- update( object, as.formula(paste("~ . -", tt)),
+ verbose = verbose
+ )
+ } else {
+ nfit <- update( object, as.formula(paste("~ . -", tt)), data = data,
+ verbose = verbose
+ )
+ }
+ }
+ if( !nfit$converged ){
+ converged <- FALSE
+ if( verbose > 0 ) cat(
+ "lack of convergence when fitting model", paste("~ . -", tt),
+ "\nconvergence code:", nfit$convergence.code, "\n"
+ )
+ }
+ }
+ nnew <- nobs( nfit, use.fallback = TRUE )
+ if( all(is.finite(c(n0, nnew))) && nnew != n0 ) stop(
+ "number of rows in use has changed: remove missing values?"
+ )
+ df.aic <- c( extractAIC( nfit, scale, k = k, ... ), as.numeric(converged))
+ names( df.aic ) <- c("df", "AIC", "converged")
+ attr( df.aic, "nfit" ) <- list(
+ param = nfit[["param"]],
+ aniso = nfit[["aniso"]][["aniso"]],
+ initial.objects = nfit[["initial.objects"]],
+ call = nfit[["call"]]
+ )
+ df.aic
+ }
+
## check scope
tl <- attr( terms(object), "term.labels" )
@@ -1131,59 +1304,74 @@
n0 <- nobs( object, use.fallback = TRUE )
env <- environment( formula(object) )
- all.converged <- TRUE
+ param <- object[["param"]]
+ aniso <- object[["aniso"]][["aniso"]]
+ fit.param <- c(
+ variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
+ alpha = FALSE, beta = FALSE, delta = FALSE,
+ gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
+ )[names( object[["param"]] )]
+ fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
+
+ ## prepare for parallel execution
- for( i in seq(ns) ) {
- tt <- scope[i]
- if(trace > 1) {
- cat("trying -", tt, "\n", sep = "")
- utils::flush.console()
- }
- ## was:
- ## nfit <- update( object, as.formula(paste("~ . -", tt)), evaluate = FALSE )
- ## nfit <- eval( nfit, envir=env ) # was eval.parent(nfit)
- if( fixed ){
- nfit <- update(
- object, as.formula(paste("~ . -", tt)),
- param = object[["param"]],
- aniso = object[["aniso"]][["aniso"]],
- fit.param = c(
- variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
- a = FALSE, alpha = FALSE, beta = FALSE, delta = FALSE,
- gamma = FALSE, lambda = FALSE, n = FALSE, nu = FALSE
- )[names( object[["param"]] )],
- fit.aniso = c(
- f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE
- ),
- verbose = verbose
- )
- } else {
- if( use.fitted.param ){
- nfit <- update(
- object, as.formula(paste("~ . -", tt)),
- param = object[["param"]], aniso = object[["aniso"]][["aniso"]],
- initial.param = "no", verbose = verbose
- )
- }
- else {
- nfit <- update( object, as.formula(paste("~ . -", tt)), verbose = verbose )
- }
- if( !nfit$converged ){
- all.converged <- FALSE
- if( verbose > 0 ) cat(
- "lack of convergence when fitting model", as.character(formula(nfit)),
- "\nconvergence code:", nfit$convergence.code, "\n"
- )
- }
- }
- ans[i+1, ] <- extractAIC( nfit, scale, k = k, ... )
- nnew <- nobs( nfit, use.fallback = TRUE )
- if( all(is.finite(c(n0, nnew))) && nnew != n0 ) stop(
- "number of rows in use has changed: remove missing values?"
+ if( ncores > 1 ){
+ ncores <- min( ncores, ns, detectCores() )
+ trace <- 0
+ verbose <- 0
+ }
+
+ if( .Platform[["OS.type"]] == "windows" && ncores > 1 ){
+
+ if( is.null( data ) ) stop(
+ "argument 'data' required for parallel execution on windows OS"
)
+
+ ## create a SNOW cluster on windows OS
+
+ cl <- makePSOCKcluster( ncores, outfile = "")
+
+ ## export required items to workers
+
+ junk <- clusterEvalQ( cl, require( georob, quietly = TRUE ) )
+
+ result <- parLapply(
+ cl,
+ scope, f.aux,
+ scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
+ object = object, data = data, param = param, aniso = aniso,
+ fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose, ...
+ )
+
+ stopCluster(cl)
+
+ } else {
+
+ ## fork child processes on non-windows OS
+
+ result <- mclapply(
+ scope, f.aux,
+ scale = scale, k= k, trace = trace, fixed = fixed, use.fitted.param = use.fitted.param,
+ object = object, data = data, param = param, aniso = aniso,
+ fit.param = fit.param, fit.aniso = fit.aniso, verbose = verbose,
+ mc.cores = ncores, mc.allow.recursive = FALSE, ...
+ )
+
}
- if( !all.converged ) warning( "lack of convergence when fitting models" )
+ fits <- list(
+ list(
+ param = object[["param"]], aniso = object[["aniso"]][["aniso"]],
+ initial.objects = object[["initial.objects"]],
+ call = object[["call"]]
+ )
+ )
+ fits[2:nrow(ans)] <- lapply( result, function(x) attr( x, "nfit" ) )
+ result <- t( simplify2array( result ) )
+ ans[2:NROW(ans), ] <- result[, 1:2]
+ if(!all( as.logical(result[, "converged"]) ) ) warning(
+ "lack of convergence when fitting models"
+ )
## store results
@@ -1208,7 +1396,9 @@
head <- c("Single term deletions", "\nModel:", deparse(formula(object)),
if(scale > 0) paste("\nscale: ", format(scale), "\n"))
class(aod) <- c("anova", "data.frame")
- attr(aod, "heading") <- head
+ attr( aod, "heading") <- head
+ attr( aod, "nfit" ) <- fits[[which.min( aod[, "AIC"])]]
+
aod
}
@@ -1227,20 +1417,21 @@
step.georob <- function( object, scope, scale = 0,
direction = c( "both", "backward", "forward" ), trace = 1, keep = NULL, steps = 1000,
- k = 2, fixed = TRUE, use.fitted.param =TRUE, verbose = 0, ... )
+ k = 2, data = NULL, fixed = TRUE, use.fitted.param =TRUE, verbose = 0,
+ ncores = 1, ... )
{
## step method for class georob
- ## 2014-01-03 AP
+ ## 2014-03-12 AP
## code of step{stats} complemented by argument fixed to control whether
## variogram parameters should be kept fixed
mydeviance <- function(x, ...)
{
- dev <- deviance(x)
- if(!is.null(dev)) dev else extractAIC(x, k=0)[2L]
+ dev <- deviance(x, ...)
+ if(!is.null(dev)) dev else extractAIC(x, k=0, ...)[2L]
}
cut.string <- function(string)
@@ -1329,7 +1520,7 @@
}
## FIXME think about df.residual() here
- models[[nm]] <- list(deviance = mydeviance(fit), df.resid = n - edf,
+ models[[nm]] <- list(deviance = mydeviance(fit, ...), df.resid = n - edf,
change = "", AIC = bAIC)
if(!is.null(keep)) keep.list[[nm]] <- keep(fit, bAIC)
usingCp <- FALSE
@@ -1341,8 +1532,12 @@
aod <- NULL
change <- NULL
if(backward && length(scope$drop)) {
- aod <- drop1.georob(fit, scope$drop, scale = scale,
- trace = trace, k = k, fixed = fixed, verbose = verbose, ...)
+ aod <- drop1.georob(
+ fit, scope$drop, scale = scale,
+ k = k, trace = trace, data = data, fixed = fixed,
+ use.fitted.param = use.fitted.param, verbose = verbose,
+ ncores = ncores, ...
+ )
rn <- row.names(aod)
row.names(aod) <- c(rn[1L], paste("-", rn[-1L], sep=" "))
## drop zero df terms first: one at time since they
@@ -1354,8 +1549,12 @@
}
if(is.null(change)) {
if(forward && length(scope$add)) {
- aodf <- add1.georob(fit, scope$add, scale = scale,
- trace = trace, k = k, fixed = fixed, verbose = verbose, ...)
+ aodf <- add1.georob(
+ fit, scope$add, scale = scale,
+ k = k, trace = trace, data = data, fixed = fixed,
+ use.fitted.param = use.fitted.param, verbose = verbose,
+ ncores = ncores, ...
+ )
rn <- row.names(aodf)
row.names(aodf) <- c(rn[1L], paste("+", rn[-1L], sep=" "))
aod <-
@@ -1380,32 +1579,36 @@
## was:
## fit <- update(fit, paste("~ .", change), evaluate = FALSE)
## fit <- eval.parent(fit)
- if( fixed ){
+
+ nfit <- switch(
+ substr( change, 1, 1 ),
+ "-" = attr( aod, "nfit" ),
+ "+" = attr( aodf, "nfit" )
+ )
+ param <- nfit[["param"]]
+ aniso <- nfit[["aniso"]]
+ fit.param <- c(
+ variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
+ alpha = FALSE, beta = FALSE, delta = FALSE,
+ gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
+ )[names( fit[["param"]] )]
+ fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
+ if( is.null( data ) ){
fit <- update(
fit, paste("~ .", change),
- param = fit[["param"]],
- aniso = fit[["aniso"]][["aniso"]],
- fit.param = c(
- variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
- a = FALSE, alpha = FALSE, beta = FALSE, delta = FALSE,
- gamma = FALSE, lambda = FALSE, n = FALSE, nu = FALSE
- )[names( fit[["param"]] )],
- fit.aniso = c(
- f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE
- ),
- verbose = verbose
+ param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
+ initial.param = "no", verbose = verbose
)
} else {
- if( use.fitted.param ){
- fit <- update(
- fit, paste("~ .", change),
- param = fit[["param"]], aniso = fit[["aniso"]][["aniso"]],
- verbose = verbose
- )
- } else {
- fit <- update( fit, paste("~ .", change), verbose = verbose )
- }
+ fit <- update(
+ fit, paste("~ .", change), data = data,
+ param = param, aniso = aniso, fit.param = fit.param, fit.aniso = fit.aniso,
+ initial.param = "no", verbose = verbose
+ )
}
+ fit[["call"]] <- nfit[["call"]]
+ fit[["initial.objects"]] <- nfit[["initial.objects"]]
+
nnew <- nobs(fit, use.fallback = TRUE)
if(all(is.finite(c(n, nnew))) && nnew != n)
stop("number of rows in use has changed: remove missing values?")
@@ -1423,11 +1626,10 @@
nm <- nm + 1
## FIXME: think about using df.residual() here.
models[[nm]] <-
- list(deviance = mydeviance(fit), df.resid = n - edf,
+ list(deviance = mydeviance(fit, ...), df.resid = n - edf,
change = change, AIC = bAIC)
if(!is.null(keep)) keep.list[[nm]] <- keep(fit, bAIC)
}
if(!is.null(keep)) fit$keep <- re.arrange(keep.list[seq(nm)])
step.results(models = models[seq(nm)], fit, object, usingCp)
}
-
Modified: pkg/R/georob.cv.R
===================================================================
--- pkg/R/georob.cv.R 2014-03-05 13:13:43 UTC (rev 17)
+++ pkg/R/georob.cv.R 2014-05-15 16:20:15 UTC (rev 18)
@@ -61,6 +61,7 @@
## 2013-07-02 AP passing initial values of aniso and fit.aniso to georob via update
## 2013-07-05 AP return "variogram.model" as part of fit componnent
## 2014-02-21 AP catching problem when factor are very unbalanced
+ ## 2014-05-15 AP changes for version 3 of RandomFields
## auxiliary function that fits the model and computes the predictions of
## a cross-validation set
@@ -77,8 +78,8 @@
if( !re.estimate ){
fit.param <- c(
variance = FALSE, snugget = FALSE, nugget = FALSE, scale = FALSE,
- a = FALSE, alpha = FALSE, beta = FALSE, delta = FALSE,
- gamma = FALSE, lambda = FALSE, n = FALSE, nu = FALSE,
+ alpha = FALSE, beta = FALSE, delta = FALSE,
+ gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE,
f1 = FALSE, f2 =FALSE, omega = FALSE, phi = FALSE, zeta = FALSE
)[names( param )]
fit.aniso <- c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE )
Modified: pkg/R/georob.exported.functions.R
===================================================================
--- pkg/R/georob.exported.functions.R 2014-03-05 13:13:43 UTC (rev 17)
+++ pkg/R/georob.exported.functions.R 2014-05-15 16:20:15 UTC (rev 18)
@@ -4,17 +4,17 @@
model = TRUE, x = FALSE, y = FALSE,
contrasts = NULL, offset,
locations,
- variogram.model = c( "exponential", "bessel", "cauchy", "cauchytbm",
- "circular", "cubic", "dagum", "dampedcosine", "DeWijsian", "fractalB",
- "gauss", "genB", "gencauchy", "gengneiting", "gneiting", "lgd1",
- "matern", "penta", "power", "qexponential", "spherical", "stable",
- "wave", "whittle"
+ variogram.model = c( "RMexp", "RMbessel", "RMcauchy",
+ "RMcircular", "RMcubic", "RMdagum", "RMdampedcos", "RMdewijsian", "RMfbm",
+ "RMgauss", "RMgenfbm", "RMgencauchy", "RMgengneiting", "RMgneiting", "RMlgd",
+ "RMmatern", "RMpenta", "RMaskey", "RMqexp", "RMspheric", "RMstable",
+ "RMwave", "RMwhittle"
),
param,
fit.param = c(
variance = TRUE, snugget = FALSE, nugget = TRUE, scale = TRUE,
- a = FALSE, alpha = FALSE, beta = FALSE, delta = FALSE,
- gamma = FALSE, lambda = FALSE, n = FALSE, nu = FALSE
+ alpha = FALSE, beta = FALSE, delta = FALSE,
+ gamma = FALSE, kappa = FALSE, lambda = FALSE, mu = FALSE, nu = FALSE
)[ names(param) ],
aniso = c( f1 = 1., f2 = 1., omega = 90., phi = 90., zeta = 0. ),
fit.aniso = c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE ),
@@ -82,6 +82,7 @@
## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
## 2013-09-06 AP exclusive use of nleqslv for solving estimating equations
## 2014-02-18 AP correcting error when fitting models with offset
+ ## 2014-05-15 AP changes for version 3 of RandomFields
## check whether input is complete
@@ -179,6 +180,7 @@
)
initial.param <- match.arg( initial.param )
+ if( tuning.psi < control[["tuning.psi.nr"]] ) initial.param <- "no"
## check whether design matrix has full column rank
@@ -642,6 +644,7 @@
## 2013-06-12 AP changes in stored items of Valpha object
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
+ ## 2014-05-15 AP changes for version 3 of RandomFields
if(
!( all( param.tf %in% names( fwd.tf ) ) &&
@@ -673,7 +676,7 @@
cov.ehat.p.bhat = cov.ehat.p.bhat, full.cov.ehat.p.bhat = full.cov.ehat.p.bhat,
aux.cov.pred.target = aux.cov.pred.target,
min.condnum = min.condnum,
- irf.models = c( "DeWijsian", "fractalB", "genB" ),
+ irf.models = c( "RMdewijsian", "RMfbm", "RMgenfbm" ),
rq = rq, lmrob = lmrob, nleqslv = nleqslv,
## bbsolve = bbsolve,
optim = optim,
@@ -686,8 +689,9 @@
param.transf <-
function(
variance = "log", snugget = "log", nugget = "log", scale = "log",
- a = "identity", alpha = "identity", beta = "identity", delta = "identity",
- gamma = "identity", lambda = "identity", n = "identity", nu = "identity",
+ alpha = "identity", beta = "log", delta = "identity",
+ gamma = "identity", kappa = "identity", lambda = "identity",
+ mu = "log", nu = "log",
f1 = "log", f2 ="log", omega = "rad", phi = "rad", zeta = "rad"
)
{
@@ -696,11 +700,12 @@
## parameters
## 2013-07-02 A. Papritz
+ ## 2014-05-15 AP changes for version 3 of RandomFields
c(
variance = variance, snugget = snugget, nugget = nugget, scale = scale,
- a = a, alpha = alpha, beta = beta, delta = delta, gamma = gamma,
- lambda = lambda, n = n, nu = nu,
+ alpha = alpha, beta = beta, delta = delta, gamma = gamma,
+ kappa = kappa, lambda = lambda, mu = mu, nu = nu,
f1 = f1, f2 = f2, omega = omega, phi = phi, zeta = zeta
)
@@ -976,33 +981,33 @@
## models (cf. Variogram{RandomFields})
## 2012-01-24 A. Papritz
+ ## 2014-05-15 AP changes for version 3 of RandomFields
switch(
model,
- "bessel" = "nu",
- "cauchy" = "beta",
- "cauchytbm" = c( "alpha", "beta", "gamma" ),
- "circular" = NULL,
- "cubic" = NULL,
- "dagum" = c( "beta", "gamma" ),
- "dampedcosine" = "lambda",
- "DeWijsian" = "alpha",
- "exponential" = NULL,
- "fractalB" = "alpha",
- "gauss" = NULL,
- "genB" = c( "alpha", "delta" ),
- "gencauchy" = c( "alpha", "beta" ),
- "gengneiting" = c( "n", "alpha" ),
- "gneiting" = NULL,
- "lgd1" = c( "alpha", "beta" ),
- "matern" = "nu",
- "penta" = NULL,
- "power" = "a",
- "qexponential" = "alpha",
- "spherical" = NULL,
- "stable" = "alpha",
- "wave" = NULL,
- "whittle" = "nu",
+ "RMbessel" = "nu",
+ "RMcauchy" = "gamma",
+ "RMcircular" = NULL,
+ "RMcubic" = NULL,
+ "RMdagum" = c( "beta", "gamma" ),
+ "RMdampedcos" = "lambda",
+ "RMdewijsian" = "alpha",
+ "RMexp" = NULL,
+ "RMfbm" = "alpha",
+ "RMgauss" = NULL,
+ "RMgenfbm" = c( "alpha", "delta" ),
+ "RMgencauchy" = c( "alpha", "beta" ),
+ "RMgengneiting" = c( "kappa", "mu" ),
+ "RMgneiting" = NULL,
+ "RMlgd" = c( "alpha", "beta" ),
+ "RMmatern" = "nu",
+ "RMpenta" = NULL,
+ "RMaskey" = "alpha",
+ "RMqexp" = "alpha",
+ "RMspheric" = NULL,
+ "RMstable" = "alpha",
+ "RMwave" = NULL,
+ "RMwhittle" = "nu",
stop( model, " variogram not implemented" )
)
}
@@ -1010,46 +1015,46 @@
## ##############################################################################
param.bounds <-
- function( model, d, param )
+function( model, d )
{
## function returns range of parameters for which variogram models are
## valid (cf. Variogram{RandomFields})
## 2012-03-30 A. Papritz
+ ## 2014-05-15 AP changes for version 3 of RandomFields
switch(
model,
- "bessel" = list( nu = c( 0.5 * (d - 2.), Inf ) ),
- "cauchy" = list( beta = c( 1.e-18, Inf ) ),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/georob -r 18
More information about the Georob-commits
mailing list