[Georob-commits] r12 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Jul 10 11:02:18 CEST 2013
Author: papritz
Date: 2013-07-10 11:02:14 +0200 (Wed, 10 Jul 2013)
New Revision: 12
Modified:
pkg/ChangeLog
pkg/R/georob.exported.functions.R
pkg/R/georob.private.functions.R
pkg/man/cv.georob.Rd
pkg/man/georob.Rd
Log:
computing robust intial values of parameters by minimizing sum of squared estimating equations
Modified: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog 2013-07-09 11:15:07 UTC (rev 11)
+++ pkg/ChangeLog 2013-07-10 09:02:14 UTC (rev 12)
@@ -124,3 +124,10 @@
2013-07-09 Andreas Papritz <papritz at env.ethz.ch>
* georob.private.functions.R (georob.fit): catching errors occuring when fitting anisotropic variogram models with default anisotropy parameters
+
+
+2013-07-10 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.exported.functions.R (georob): computing robust initial variogram parameter estimates by minimizing sum of squared estimating equations
+* georob.private.functions.R (georob.fit, compute.estimating.equations): computing robust initial variogram parameter estimates by minimizing sum of squared estimating equations
+
Modified: pkg/R/georob.exported.functions.R
===================================================================
--- pkg/R/georob.exported.functions.R 2013-07-09 11:15:07 UTC (rev 11)
+++ pkg/R/georob.exported.functions.R 2013-07-10 09:02:14 UTC (rev 12)
@@ -18,7 +18,7 @@
)[ 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 ),
- tuning.psi = 2, initial.param = TRUE,
+ tuning.psi = 2, initial.param = c( "minimize", "exclude", "no" ),
control = georob.control( ... ), verbose = 0,
...
)
@@ -188,7 +188,7 @@
signif( condnum, 2 ), ")\ninitial values of fixed effects coefficients are computed by 'lm'\n\n"
)
control[["initial.method"]] <- "lm"
- initial.param <- FALSE
+ initial.param <- "mininimize"
warning(
"design matrix has not full column rank (condition number of X^T X: ",
signif( condnum, 2 ), ")\ninitial values of fixed effects coefficients are computed by 'lm'"
@@ -202,7 +202,7 @@
## adjust choice for initial.method to compute regression coefficients
- if( initial.param ) control[["initial.method"]] <- "lmrob"
+ if( identical( initial.param, "exclude" ) ) control[["initial.method"]] <- "lmrob"
## compute initial guess of fixed effects parameters (betahat)
@@ -318,86 +318,162 @@
## prune data set for computing initial values of variogram and
## anisotropy parameters by reml
- if( initial.param && tuning.psi < control[["tuning.psi.nr"]] ){
-
- t.sel <- switch(
- control[["initial.method"]],
- lmrob = r.initial.fit[["rweights"]] > control[["min.rweight"]],
- stop(
- "computing initial values of covariance parameters requires 'lmrob' as ",
- "method for computing initial regression coefficients"
+ initial.param <- match.arg( initial.param )
+
+ if( tuning.psi < control[["tuning.psi.nr"]] ){
+
+ if( identical( initial.param, "exclude" ) ){
+
+ if( verbose > 0 ) cat( "\ncomputing robust initial parameter estimates ...\n" )
+
+ t.sel <- switch(
+ control[["initial.method"]],
+ lmrob = r.initial.fit[["rweights"]] > control[["min.rweight"]],
+ stop(
+ "computing initial values of covariance parameters requires 'lmrob' as ",
+ "method for computing initial regression coefficients"
+ )
)
- )
+
+ if( verbose > 0 ) cat(
+ "\ndiscarding", sum( !t.sel ), "of", length( t.sel ),
+ "observations for computing initial estimates of variogram\nand anisotropy parameter by gaussian reml\n"
+ )
+
+ ## collect.items for initial object
+
+ initial.objects <- list(
+ x = as.matrix( x[t.sel, ] ),
+ y = yy[t.sel],
+ betahat = coef( r.initial.fit ),
+ bhat = if( is.null( control[["bhat"]] ) ){
+ rep( 0., length( yy ) )[t.sel]
+ } else {
+ control[["bhat"]][t.sel]
+ },
+ initial.fit = r.initial.fit,
+ locations.objects = list(
+ locations = locations,
+ coordinates = locations.coords[t.sel, ]
+ ),
+ isotropic = aniso.missing
+ )
+
+ ## estimate model parameters with pruned data set
+
+ t.georob <- georob.fit(
+ slv = TRUE,
+ envir = envir,
+ initial.objects = initial.objects,
+ variogram.model = variogram.model, param = param, fit.param = fit.param,
+ aniso = aniso, fit.aniso = fit.aniso,
+ param.tf = control[["param.tf"]],
+ fwd.tf = control[["fwd.tf"]],
+ deriv.fwd.tf = control[["deriv.fwd.tf"]],
+ bwd.tf = control[["bwd.tf"]],
+ safe.param = control[["safe.param"]],
+ tuning.psi = control[["tuning.psi.nr"]],
+ cov.bhat = control[["cov.bhat"]], full.cov.bhat = control[["full.cov.bhat"]],
+ cov.betahat = control[["cov.betahat"]],
+ cov.bhat.betahat = control[["cov.bhat.betahat"]],
+ cov.delta.bhat = control[["cov.delta.bhat"]],
+ full.cov.delta.bhat = control[["full.cov.delta.bhat"]],
+ cov.delta.bhat.betahat = control[["cov.delta.bhat.betahat"]],
+ cov.ehat = control[["cov.ehat"]], full.cov.ehat = control[["full.cov.ehat"]],
+ cov.ehat.p.bhat = control[["cov.ehat.p.bhat"]], full.cov.ehat.p.bhat = control[["full.cov.ehat.p.bhat"]],
+ aux.cov.pred.target = control[["aux.cov.pred.target"]],
+ min.condnum = control[["min.condnum"]],
+ rankdef.x = rankdef.x,
+ psi.func = control[["psi.func"]],
+ tuning.psi.nr = tuning.psi,
+ irwls.initial = control[["irwls.initial"]],
+ irwls.maxiter = control[["irwls.maxiter"]],
+ irwls.reltol = control[["irwls.reltol"]],
+ force.gradient = control[["force.gradient"]],
+ zero.dist = control[["zero.dist"]],
+ nleqslv.method = control[["nleqslv"]][["method"]],
+ nleqslv.control = control[["nleqslv"]][["control"]],
+ optim.method = control[["optim"]][["method"]],
+ optim.lower = control[["optim"]][["lower"]],
+ optim.upper = control[["optim"]][["upper"]],
+ hessian = control[["optim"]][["hessian"]],
+ optim.control = control[["optim"]][["control"]],
+ full.output = control[["full.output"]],
+ verbose = verbose
+ )
+
+ param = t.georob[["param"]][names(fit.param)]
+ aniso = t.georob[["aniso"]][["aniso"]][names(fit.aniso)]
+
+ } else if( identical( initial.param, "minimize" ) ){
+
+ if( verbose > 0 ) cat( "\ncomputing robust initial parameter estimates ...\n" )
+
+ initial.objects <- list(
+ x = as.matrix( x ),
+ y = yy,
+ betahat = coef( r.initial.fit ),
+ bhat = if( is.null( control[["bhat"]] ) ){
+ rep( 0., length( yy ) )
+ } else {
+ control[["bhat"]]
+ },
+ initial.fit = r.initial.fit,
+ locations.objects = list(
+ locations = locations,
+ coordinates = locations.coords
+ ),
+ isotropic = aniso.missing
+ )
+
+ ## estimate model parameters by minimizing sum( gradient^2)
+
+ t.georob <- georob.fit(
+ slv = FALSE,
+ envir = envir,
+ initial.objects = initial.objects,
+ variogram.model = variogram.model, param = param, fit.param = fit.param,
+ aniso = aniso, fit.aniso = fit.aniso,
+ param.tf = control[["param.tf"]],
+ fwd.tf = control[["fwd.tf"]],
+ deriv.fwd.tf = control[["deriv.fwd.tf"]],
+ bwd.tf = control[["bwd.tf"]],
+ safe.param = control[["safe.param"]],
+ tuning.psi = tuning.psi,
+ cov.bhat = control[["cov.bhat"]], full.cov.bhat = control[["full.cov.bhat"]],
+ cov.betahat = control[["cov.betahat"]],
+ cov.bhat.betahat = control[["cov.bhat.betahat"]],
+ cov.delta.bhat = control[["cov.delta.bhat"]],
+ full.cov.delta.bhat = control[["full.cov.delta.bhat"]],
+ cov.delta.bhat.betahat = control[["cov.delta.bhat.betahat"]],
+ cov.ehat = control[["cov.ehat"]], full.cov.ehat = control[["full.cov.ehat"]],
+ cov.ehat.p.bhat = control[["cov.ehat.p.bhat"]], full.cov.ehat.p.bhat = control[["full.cov.ehat.p.bhat"]],
+ aux.cov.pred.target = control[["aux.cov.pred.target"]],
+ min.condnum = control[["min.condnum"]],
+ rankdef.x = rankdef.x,
+ psi.func = control[["psi.func"]],
+ tuning.psi.nr = control[["tuning.psi.nr"]],
+ irwls.initial = control[["irwls.initial"]],
+ irwls.maxiter = control[["irwls.maxiter"]],
+ irwls.reltol = control[["irwls.reltol"]],
+ force.gradient = control[["force.gradient"]],
+ zero.dist = control[["zero.dist"]],
+ nleqslv.method = control[["nleqslv"]][["method"]],
+ nleqslv.control = control[["nleqslv"]][["control"]],
+ optim.method = control[["optim"]][["method"]],
+ optim.lower = control[["optim"]][["lower"]],
+ optim.upper = control[["optim"]][["upper"]],
+ hessian = control[["optim"]][["hessian"]],
+ optim.control = control[["optim"]][["control"]],
+ full.output = control[["full.output"]],
+ verbose = verbose
+ )
+
+ param = t.georob[["param"]][names(fit.param)]
+ aniso = t.georob[["aniso"]][["aniso"]][names(fit.aniso)]
+
+ }
- if( verbose > 0 ) cat(
- "\ndiscarding", sum( !t.sel ), "of", length( t.sel ),
- "observations for computing initial estimates of variogram\nand anisotropy parameter by gaussian reml\n"
- )
-
- ## collect.items for initial object
-
- initial.objects <- list(
- x = as.matrix( x[t.sel, ] ),
- y = yy[t.sel],
- betahat = coef( r.initial.fit ),
- bhat = if( is.null( control[["bhat"]] ) ){
- rep( 0., length( yy ) )[t.sel]
- } else {
- control[["bhat"]][t.sel]
- },
- initial.fit = r.initial.fit,
- locations.objects = list(
- locations = locations,
- coordinates = locations.coords[t.sel, ]
- ),
- isotropic = aniso.missing
- )
-
- ## estimate model parameters with pruned data set
-
- t.georob <- georob.fit(
- envir = envir,
- initial.objects = initial.objects,
- variogram.model = variogram.model, param = param, fit.param = fit.param,
- aniso = aniso, fit.aniso = fit.aniso,
- param.tf = control[["param.tf"]],
- fwd.tf = control[["fwd.tf"]],
- deriv.fwd.tf = control[["deriv.fwd.tf"]],
- bwd.tf = control[["bwd.tf"]],
- safe.param = control[["safe.param"]],
- tuning.psi = control[["tuning.psi.nr"]],
- cov.bhat = control[["cov.bhat"]], full.cov.bhat = control[["full.cov.bhat"]],
- cov.betahat = control[["cov.betahat"]],
- cov.bhat.betahat = control[["cov.bhat.betahat"]],
- cov.delta.bhat = control[["cov.delta.bhat"]],
- full.cov.delta.bhat = control[["full.cov.delta.bhat"]],
- cov.delta.bhat.betahat = control[["cov.delta.bhat.betahat"]],
- cov.ehat = control[["cov.ehat"]], full.cov.ehat = control[["full.cov.ehat"]],
- cov.ehat.p.bhat = control[["cov.ehat.p.bhat"]], full.cov.ehat.p.bhat = control[["full.cov.ehat.p.bhat"]],
- aux.cov.pred.target = control[["aux.cov.pred.target"]],
- min.condnum = control[["min.condnum"]],
- rankdef.x = rankdef.x,
- psi.func = control[["psi.func"]],
- tuning.psi.nr = control[["tuning.psi.nr"]],
- irwls.initial = control[["irwls.initial"]],
- irwls.maxiter = control[["irwls.maxiter"]],
- irwls.reltol = control[["irwls.reltol"]],
- force.gradient = control[["force.gradient"]],
- zero.dist = control[["zero.dist"]],
- nleqslv.method = control[["nleqslv"]][["method"]],
- nleqslv.control = control[["nleqslv"]][["control"]],
- optim.method = control[["optim"]][["method"]],
- optim.lower = control[["optim"]][["lower"]],
- optim.upper = control[["optim"]][["upper"]],
- hessian = control[["optim"]][["hessian"]],
- optim.control = control[["optim"]][["control"]],
- full.output = control[["full.output"]],
- verbose = verbose
- )
-
- param = t.georob[["param"]][names(fit.param)]
- aniso = t.georob[["aniso"]][["aniso"]][names(fit.aniso)]
-
}
## collect.items for initial object
@@ -420,8 +496,11 @@
)
## estimate model parameters
-
+
+ if( verbose > 0 ) cat( "computing final parameter estimates ...\n" )
+
r.georob <- georob.fit(
+ slv = TRUE,
envir = envir,
initial.objects = initial.objects,
variogram.model = variogram.model, param = param, fit.param = fit.param,
Modified: pkg/R/georob.private.functions.R
===================================================================
--- pkg/R/georob.private.functions.R 2013-07-09 11:15:07 UTC (rev 11)
+++ pkg/R/georob.private.functions.R 2013-07-10 09:02:14 UTC (rev 12)
@@ -2348,6 +2348,7 @@
compute.estimating.equations <-
function(
adjustable.param,
+ slv,
envir,
variogram.model, fixed.param, param.name, aniso.name,
param.tf, bwd.tf, safe.param,
@@ -2589,7 +2590,19 @@
assign( "lik.item", lik.item, pos = as.environment( envir ) )
- return( eeq.emp / eeq.exp - 1. )
+ if( slv ){
+ return( eeq.emp / eeq.exp - 1. )
+ } else {
+ res <- sum( (eeq.emp / eeq.exp - 1.)^2 )
+ if( verbose > 1 ) cat(
+ " sum(EEQ^2) :",
+ format(
+ signif( res, digits = 7 ),
+ scientific = TRUE, width = 14
+ ), "\n", sep = ""
+ )
+ return( res )
+ }
} else {
@@ -3155,6 +3168,7 @@
georob.fit <-
function(
+ slv,
envir,
initial.objects,
variogram.model, param, fit.param,
@@ -3535,10 +3549,10 @@
if( fit.aniso["omega"] && aniso["f1"] == 1. ) aniso["f1"] <- aniso["f1"] - sqrt( .Machine$double.eps )
if( fit.aniso["phi"] ){
- if( aniso["f1"] == 1. ) aniso["f1"] <- aniso["f1"] - sqrt( .Machine$double.eps )
- if( aniso["f2"] == 1. ) aniso["f2"] <- aniso["f2"] - sqrt( .Machine$double.eps )
+ if( aniso["f1"] == 1. ) aniso["f1"] <- aniso["f1"] - 0.0001
+ if( aniso["f2"] == 1. ) aniso["f2"] <- aniso["f2"] - 0.0001
}
- if( fit.aniso["zeta"] && aniso["f2"] == 1. ) aniso["f2"] <- aniso["f2"] - sqrt( .Machine$double.eps )
+ if( fit.aniso["zeta"] && aniso["f2"] == 1. ) aniso["f2"] <- aniso["f2"] - 0.0001
## rearrange and check flags controlling anisotropy parameter fitting
@@ -3607,7 +3621,7 @@
if( !identical( t.exp[["message"]], "OK" ) ) stop( t.exp[["message"]] )
expectations["dpsi"] <- t.exp[["value"]]
if( verbose > 1 ) cat(
- "\nexpectation of psi'(epsilon/sigma) :",
+ "\nexpectation of psi'(epsilon/sigma) :",
signif( expectations["dpsi"] ), "\n"
)
@@ -3624,81 +3638,181 @@
if( !identical( t.exp[["message"]], "OK" ) ) stop( t.exp[["message"]] )
expectations["psi2"] <- t.exp[["value"]]
if( verbose > 1 ) cat(
- "expectation of (psi(epsilon/sigma))^2 :",
+ "expectation of (psi(epsilon/sigma))^2 :",
signif( t.exp[["value"]] ), "\n"
)
-
-
r.hessian <- NULL
if( tuning.psi < tuning.psi.nr ) {
+ ## robust REML estimation
+
if( any( c( fit.param, fit.aniso ) ) ){
- ## some variogram parameters are fitted
+ ## find roots of estimating equations
- ## find root of estimating equations
+ if( slv ){
- r.root <- nleqslv(
- x = c(
- transformed.param[ fit.param ],
- transformed.aniso[ fit.aniso ]
- ),
- fn = compute.estimating.equations,
- method = nleqslv.method,
- control = nleqslv.control,
- envir = envir,
- variogram.model = variogram.model,
- fixed.param = c(
- transformed.param[ !fit.param ],
- transformed.aniso[ !fit.aniso ]
- ),
- param.name = param.name,
- aniso.name = aniso.name,
- param.tf = param.tf,
- bwd.tf = bwd.tf,
- safe.param = safe.param,
- lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
- yy = yy, betahat = betahat, TT = TT, bhat = bhat,
- psi.function = rho.psi.etc[["psi.function"]],
- dpsi.function = rho.psi.etc[["dpsi.function"]],
- tuning.psi = tuning.psi,
- tuning.psi.nr = tuning.psi.nr,
- irwls.initial = irwls.initial,
- irwls.maxiter = irwls.maxiter,
- irwls.reltol = irwls.reltol,
- force.gradient = force.gradient,
- expectations = expectations,
- verbose = verbose
- )
+ r.root <- nleqslv(
+ x = c(
+ transformed.param[ fit.param ],
+ transformed.aniso[ fit.aniso ]
+ ),
+ fn = compute.estimating.equations,
+ method = nleqslv.method,
+ control = nleqslv.control,
+ slv = slv,
+ envir = envir,
+ variogram.model = variogram.model,
+ fixed.param = c(
+ transformed.param[ !fit.param ],
+ transformed.aniso[ !fit.aniso ]
+ ),
+ param.name = param.name,
+ aniso.name = aniso.name,
+ param.tf = param.tf,
+ bwd.tf = bwd.tf,
+ safe.param = safe.param,
+ lag.vectors = lag.vectors,
+ XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
+ yy = yy, betahat = betahat, TT = TT, bhat = bhat,
+ psi.function = rho.psi.etc[["psi.function"]],
+ dpsi.function = rho.psi.etc[["dpsi.function"]],
+ tuning.psi = tuning.psi,
+ tuning.psi.nr = tuning.psi.nr,
+ irwls.initial = irwls.initial,
+ irwls.maxiter = irwls.maxiter,
+ irwls.reltol = irwls.reltol,
+ force.gradient = force.gradient,
+ expectations = expectations,
+ verbose = verbose
+ )
+
+ # r.param <- r.root[["x"]] names( r.param ) <- names(
+ # transformed.param[ fit.param ] )
+
+ r.gradient <- r.root[["fvec"]]
+ names( r.gradient ) <- c(
+ names( transformed.param[ fit.param ] ),
+ names( transformed.aniso[ fit.aniso ] )
+ )
+
+ r.converged <- r.root[["termcd"]] == 1
+ r.convergence.code <- r.root[["termcd"]]
+
+ r.counts <- c( nfcnt = r.root[["nfcnt"]], njcnt = r.root[["njcnt"]] )
+
+ } else {
+
+ ## minimize sum of squared estimating equations
+
+ r.opt.eeq.sq <- optim(
+ par = c(
+ transformed.param[ fit.param ],
+ transformed.aniso[ fit.aniso ]
+ ),
+ fn = compute.estimating.equations,
+ # gr = gradient.negative.restricted.loglikelihood,
+ method = optim.method,
+ lower = optim.lower,
+ upper = optim.upper,
+ control = optim.control,
+ hessian = FALSE,
+ slv = slv,
+ envir = envir,
+ variogram.model = variogram.model,
+ fixed.param = c(
+ transformed.param[ !fit.param ],
+ transformed.aniso[ !fit.aniso ]
+ ),
+ param.name = param.name,
+ aniso.name = aniso.name,
+ param.tf = param.tf,
+ bwd.tf = bwd.tf,
+ safe.param = safe.param,
+ lag.vectors = lag.vectors,
+ XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
+ yy = yy, betahat = betahat, TT = TT, bhat = bhat,
+ psi.function = rho.psi.etc[["psi.function"]],
+ dpsi.function = rho.psi.etc[["dpsi.function"]],
+ tuning.psi = tuning.psi,
+ tuning.psi.nr = tuning.psi.nr,
+ irwls.initial = irwls.initial,
+ irwls.maxiter = irwls.maxiter,
+ irwls.reltol = irwls.reltol,
+ force.gradient = force.gradient,
+ expectations = expectations,
+ verbose = verbose
+ )
+
+ r.opt.neg.loglik <- r.opt.eeq.sq[["value"]]
+ r.converged <- r.opt.eeq.sq[["convergence"]] == 0
+ r.convergence.code <- r.opt.eeq.sq[["convergence"]]
+ r.counts <- r.opt.eeq.sq[["counts"]]
+
+ if( verbose > 0 ){
+ cat(
+ "\n sum(EEQ^2) :",
+ format(
+ signif( r.opt.eeq.sq[["value"]], digits = 7 ),
+ scientific = TRUE, width = 14
+ ), sep = ""
+ )
+ cat(
+ "\n convergence code :",
+ format(
+ signif( r.opt.eeq.sq[["convergence"]], digits = 0 ),
+ scientific = FALSE, width = 14
+ ), "\n\n", sep = ""
+ )
+ }
+
+ # if( hessian ) r.hessian <- r.opt.eeq.sq[["hessian"]]
+
+ r.gradient <- compute.estimating.equations(
+ adjustable.param = r.opt.eeq.sq[["par"]],
+ slv = TRUE,
+ envir = envir,
+ variogram.model = variogram.model,
+ fixed.param = c(
+ transformed.param[ !fit.param ],
+ transformed.aniso[ !fit.aniso ]
+ ),
+ param.name = param.name,
+ aniso.name = aniso.name,
+ param.tf = param.tf,
+ bwd.tf = bwd.tf,
+ safe.param = safe.param,
+ lag.vectors = lag.vectors,
+ XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
+ yy = yy, betahat = betahat, TT = TT, bhat = bhat,
+ psi.function = rho.psi.etc[["psi.function"]],
+ dpsi.function = rho.psi.etc[["dpsi.function"]],
+ tuning.psi = tuning.psi,
+ tuning.psi.nr = tuning.psi.nr,
+ irwls.initial = irwls.initial,
+ irwls.maxiter = irwls.maxiter,
+ irwls.reltol = irwls.reltol,
+ force.gradient = force.gradient,
+ expectations = expectations,
+ verbose = verbose
+ )
+
+ }
- # r.param <- r.root[["x"]]
- # names( r.param ) <- names( transformed.param[ fit.param ] )
-
- r.gradient <- r.root[["fvec"]]
- names( r.gradient ) <- c(
- names( transformed.param[ fit.param ] ),
- names( transformed.aniso[ fit.aniso ] )
- )
-
- r.converged <- r.root[["termcd"]] == 1
- r.convergence.code <- r.root[["termcd"]]
-
- r.counts <- c( nfcnt = r.root[["nfcnt"]], njcnt = r.root[["njcnt"]] )
-
} else {
## all variogram parameters are fixed
- ## compute values of estimating equations
+ ## evaluate estimating equations
r.gradient <- compute.estimating.equations(
adjustable.param = c(
transformed.param[ fit.param ],
transformed.aniso[ fit.aniso ]
),
+ slv = TRUE,
envir = envir,
variogram.model = variogram.model,
fixed.param = c(
@@ -3737,8 +3851,7 @@
if( any( fit.param ) ){
- ## some variogram parameters are fitted
- ## minimize laplace approximation of negative restricted loglikelihood
+ ## Gaussian REML estimation
r.opt.neg.restricted.loglik <- optim(
par = c(
@@ -3786,7 +3899,6 @@
if( hessian ) r.hessian <- r.opt.neg.restricted.loglik[["hessian"]]
-
r.gradient <- gradient.negative.restricted.loglikelihood(
adjustable.param = r.opt.neg.restricted.loglik[["par"]],
envir = envir,
Modified: pkg/man/cv.georob.Rd
===================================================================
--- pkg/man/cv.georob.Rd 2013-07-09 11:15:07 UTC (rev 11)
+++ pkg/man/cv.georob.Rd 2013-07-10 09:02:14 UTC (rev 12)
@@ -1,4 +1,4 @@
-% 2013-07-05 A. Papritz
+% 2013-07-10 A. Papritz
% R CMD Rdconv -t html -o bla.html cv.georob.Rd ; open bla.html; R CMD Rd2pdf --force cv.georob.Rd;
\encoding{macintosh}
@@ -105,7 +105,7 @@
\item{return.fit}{logical controlling whether information about the fit
should be returned for when re-estimating the model with the reduced data
- sets (default \code{TRUE}).}
+ sets (default \code{FALSE}).}
\item{reduced.output}{logical controlling whether the complete fitted
model objects, fitted to the reduced data sets, are returned
Modified: pkg/man/georob.Rd
===================================================================
--- pkg/man/georob.Rd 2013-07-09 11:15:07 UTC (rev 11)
+++ pkg/man/georob.Rd 2013-07-10 09:02:14 UTC (rev 12)
@@ -1,4 +1,4 @@
-% 2013-06-12 A. Papritz
+% 2013-07-10 A. Papritz
% R CMD Rdconv -t html -o bla.html georob.Rd ; open bla.html; R CMD Rd2pdf --force georob.Rd;
\encoding{macintosh}
@@ -28,7 +28,8 @@
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),
- tuning.psi = 2, initial.param = TRUE, control = georob.control(...),
+ tuning.psi = 2, initial.param = c("minimize", "exclude", "no"),
+ control = georob.control(...),
verbose = 0, ...)
}
%
@@ -144,12 +145,21 @@
\item{tuning.psi}{positive numeric. The tuning constant \eqn{c} of the
\eqn{\psi_c}-function of the robust REML algorithm.}
- \item{initial.param}{logical. If \code{TRUE} (default) robust initial
- values of variogram parameters are computed by discarding outlying
- observations based on the \dQuote{robustness weights} of the initial fit
- of the regression model by \code{\link[robustbase]{lmrob}} and fitting
- the spatial linear model by Gaussian REML to the pruned data set (see
- \emph{Details}).}
+ \item{initial.param}{character, controlling whether initial values of
+ parameters are computed for solving the estimating equations of the
+ variogram and anisotropy parameters.
+
+ If \code{initial.param = "minimize"} (defaullt) robust initial values are
+ computed by minimizing the sum of the squared robustified estimating
+ equations using \code{\link[stats]{optim}} (see \emph{Details}).
+ If \code{initial.param = "exclude"} robust initial values of parameters are
+ computed by discarding outlying observations based on the
+ \dQuote{robustness weights} of the initial fit of the regression model by
+ \code{\link[robustbase]{lmrob}} and fitting the spatial linear model by
+ Gaussian REML to the pruned data set (see \emph{Details}).
+ For \code{initial.param = "no"} no initial parameter values are computed
+ and the estimating equations are solved with the initial values passed by
+ \code{param} and \code{aniso} to \code{georob}.}
\item{control}{a list specifying parameters that control the behaviour of
\code{georob}. Use the function \code{\link{georob.control}} and see its
@@ -308,13 +318,21 @@
}
Finding the roots of the robustified estimating equations of the
- variogram parameters is more sensitive to a good choice of initial
- values than maximizing the Gaussian restricted loglikelihood with
- respect to the same parameters. To get good initial values that are
- often sufficiently close to the roots so that
- \code{\link[nleqslv]{nleqslv}} converges, one can use
- \code{initial.param = TRUE}. This has the following effects:
+ variogram and anisotropy parameters is more sensitive to a good choice
+ of initial values than maximizing the Gaussian restricted loglikelihood
+ with respect to the same parameters. Two options are implemented to
+ get good initial values that are often sufficiently close to the roots
+ so that \code{\link[nleqslv]{nleqslv}} converges:
+ Setting \code{initial.param = "minimize"} invokes
+ \code{\link[stats]{optim}} to minimize the \emph{sum of squared
+ estimating equations}. The required accuracy of the initial estimates
+ are best controlled by the argument \code{abstol} of
+ \code{\link[stats]{optim}}, e.g. by using the argument
+ \code{control = georob.control(optim = optim.control(optim.control = list(abstol = 1.e-6)))}.
+
+ Setting \code{initial.param = "exclude"} has the following effects:
+
\enumerate{
\item Initial values of the regression parameters are computed by
More information about the Georob-commits
mailing list