[Georob-commits] r2 - in pkg: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 26 08:58:24 CEST 2013
Author: papritz
Date: 2013-04-26 08:58:23 +0200 (Fri, 26 Apr 2013)
New Revision: 2
Added:
pkg/ChangeLog
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/NEWS
pkg/R/
pkg/R/georob.S3methods.R
pkg/R/georob.exported.functions.R
pkg/R/georob.lgnpp.R
pkg/R/georob.predict.R
pkg/R/georob.private.functions.R
pkg/R/georob.xvalid.R
pkg/R/variogram.R
pkg/demo/
pkg/demo/00Index
pkg/demo/georob_example.R
pkg/man/
pkg/man/S3methods.georob.Rd
pkg/man/compress.Rd
pkg/man/cv.Rd
pkg/man/cv.georob.Rd
pkg/man/fit.variogram.model.Rd
pkg/man/georob-package.Rd
pkg/man/georob.Rd
pkg/man/georob.control.Rd
pkg/man/georobObject.Rd
pkg/man/internal.functions.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
pkg/man/validate.predictions.Rd
Log:
Initial import
Added: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog (rev 0)
+++ pkg/ChangeLog 2013-04-26 06:58:23 UTC (rev 2)
@@ -0,0 +1,45 @@
+2012-12-18 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.xvalid.R (print.cv.georob, print.summary.cv.georob): return invisible(x)
+* georob.S3methods.R (print.georob, print.summary.georob): return invisible(x)
+* variogram.R (print.summary.sample.variogram, print.fitted.variogram)
+(print.summary.fitted.variogram): return invisible(x)
+
+
+2012-12-18 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.S3methods.R (model.frame.georob, model.matrix.georob)
+(nobs.georob): new respective methods for class "georob"
+
+
+2012-12-18 Andreas Papritz <papritz at env.ethz.ch>
+
+* variogram.R (plot.sample.variogram, plot.georob): correction of error in
+processing col and pch arguments
+
+
+2012-12-22 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.S3methods.R (logLik.georob): computation of unrestricted loglikelihood
+(deviance.georob): new function
+
+
+2013-01-19 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.predict.R (predict.georob): computation of lag vectors
+
+
+2013-01-20 Andreas Papritz <papritz at env.ethz.ch>
+
+* NAMESPACE (--): changed imports for constrainedKriging and spatialCovariance
+* DESCRIPTION (--): changed imports for constrainedKriging and spatialCovariance
+* internal.functions.Rd (--): new help file for unexported functions
+
+
+2013-04-23 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.exported.functions.R (georob, georob.control): new names for robustness weights
+* georob.predict.R (predict.georob): new names for robustness weights
+* georob.private.functions.R (compute.covariances, update.betahat.bhat, estimate.betahat.bhat, compute.estimating equations, georob.fit): new names for robustness weights
+* georob.S3methods.R (ranef.georob, rstandard.georob, summary.georob, print.summary.georob): new names for robustness weights
+* georob.xvalid.R (cv.georob): changes for parallelization on windows os
Added: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION (rev 0)
+++ pkg/DESCRIPTION 2013-04-26 06:58:23 UTC (rev 2)
@@ -0,0 +1,21 @@
+Package: georob
+Type: Package
+Title: Robust Geostatistical Analysis of Spatial Data
+Version: 0.1-0
+Date: 2012-12-14
+Authors at R: c(
+ person( "Andreas", "Papritz", role = c( "cre", "aut" ),
+ email = "andreas.papritz at env.ethz.ch" ),
+ person( "Cornelia", "Schwierz", role = "ctb" ))
+Depends: R(>= 2.14.0), lmtest, nlme,
+ robustbase, sp(>= 0.9-60)
+Imports: constrainedKriging(>= 0.1-9), nleqslv, parallel, quantreg,
+ RandomFields(>= 2.0.55), spatialCovariance(>= 0.6-4)
+Suggests: geoR
+Description: The georob package provides functions for fitting linear models
+ with spatially correlated errors by robust and Gaussian Restricted
+ Maximum Likelihood and for computing robust and customary point
+ and block kriging predictions, along with utility functions for
+ cross-validation and for unbiased back-transformation of kriging
+ predictions of log-transformed data.
+License: GPL
Added: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE (rev 0)
+++ pkg/NAMESPACE 2013-04-26 06:58:23 UTC (rev 2)
@@ -0,0 +1,126 @@
+import( parallel, stats )
+
+importFrom( constrainedKriging, covmodel, f.point.block.cov, K, preCKrige )
+importFrom( lmtest, waldtest, waldtest.default )
+importFrom( nlme, fixef, fixed.effects, ranef, random.effects )
+importFrom( nleqslv, nleqslv )
+importFrom( quantreg, rq.fit )
+importFrom( RandomFields, Variogram )
+
+# exported functions
+
+export(
+ bwd.transf, # ok
+ compress, # ok
+ cv, # ok
+ dfwd.transf, # ok
+ expand, # ok
+ fwd.transf, # ok
+ param.bounds, # ok
+ param.names, # ok
+ fit.variogram.model, # ok
+ georob.control, # ok
+ georob, # ok
+ K, # ok re-export of K{constrainedKriging}
+ lgnpp, # ok
+ nleqslv.control, # ok
+ optim.control, # ok
+ param.transf, # ok
+ rq.control, # ok
+ sample.variogram, # ok
+ validate.predictions # ok
+)
+
+# documented but unexported functions
+#
+# deviance.georob, # ok
+# fixed.effects.georob, # ok
+# fixef.georob, # ok
+# lines.fitted.variogram, # ok
+# lines.georob, # ok
+# logLik.georob, # ok
+# model.frame.georob, # ok
+# model.matrix.georob, # ok
+# nobs.georob, # ok
+# plot.georob # ok
+# plot.sample.variogram, # ok
+# plot.cv.georob, # ok
+# predict.georob, # ok
+# print.fitted.variogram, # ok
+# print.georob, # ok
+# print.cv.georob, # ok
+# print.summary.cv.georob, # ok
+# print.summary.fitted.variogram, # ok
+# print.summary.georob, # ok
+# print.summary.sample.variogram, # ok
+# random.effects.georob, # ok
+# ranef.georob, # ok
+# resid.georob, # ok
+# residuals.georob, # ok
+# rstandard.georob, # ok
+# rstudent.georob, # ok
+# rstudent.cv.georob, # ok
+# summary.fitted.variogram, # ok
+# summary.georob, # ok
+# summary.sample.variogram, # ok
+# summary.cv.georob, # ok
+# variogram.georob, # ok
+# vcov.georob, # ok
+# waldtest.georob, # ok
+# cv.georob # ok
+# cv.likGRF # ok
+# cv.variomodel # ok
+
+# non-documented internal functions
+##
+## compute.covariances,
+## compute.estimating.equations,
+## compute.semivariance,
+## compute.U,
+## dcorr.dparam,
+## estimate.betahat.bhat,
+## gcr,
+## georob.fit,
+## getCall.georob,
+## gradient.negative.restricted.loglikelihood,
+## negative.restr.loglikelihood,
+## prepare.likelihood.calculations,
+## update.betahat.bhat
+
+S3method( deviance, georob )
+S3method( fixed.effects, georob )
+S3method( fixef, georob )
+S3method( lines, fitted.variogram )
+S3method( lines, fitted.variogram )
+S3method( lines, georob )
+S3method( logLik, georob )
+S3method( model.frame, georob )
+S3method( model.matrix, georob )
+S3method( nobs, georob )
+S3method( plot, georob )
+S3method( plot, sample.variogram )
+S3method( plot, cv.georob )
+S3method( predict, georob )
+S3method( print, fitted.variogram )
+S3method( print, georob )
+S3method( print, summary.fitted.variogram )
+S3method( print, summary.georob )
+S3method( print, summary.sample.variogram )
+S3method( print, summary.cv.georob )
+S3method( print, cv.georob )
+S3method( random.effects, georob )
+S3method( ranef, georob )
+S3method( resid, georob )
+S3method( residuals, georob )
+S3method( rstandard, georob )
+S3method( rstudent, georob )
+S3method( summary, fitted.variogram )
+S3method( summary, georob )
+S3method( summary, sample.variogram )
+S3method( summary, cv.georob )
+S3method( vcov, georob )
+S3method( waldtest, georob )
+S3method( cv, georob )
+# S3method( cv, likGRF )
+# S3method( cv, variomodel )
+
Added: pkg/NEWS
===================================================================
--- pkg/NEWS (rev 0)
+++ pkg/NEWS 2013-04-26 06:58:23 UTC (rev 2)
@@ -0,0 +1 @@
+2012-12-14 georob Version 0.1-0 released!
Added: pkg/R/georob.S3methods.R
===================================================================
--- pkg/R/georob.S3methods.R (rev 0)
+++ pkg/R/georob.S3methods.R 2013-04-26 06:58:23 UTC (rev 2)
@@ -0,0 +1,901 @@
+#####################################
+# #
+# Methoden fuer Klasse "georob" #
+# #
+#####################################
+
+# model.frame.georob
+# model.matrix.georob
+# nobs.georob
+# print.georob
+# ranef.georob
+# residuals.georob
+# resid.georob
+# rstandard.georob
+# rstudent.georob
+# summary.georob,
+# print.summary.georob
+# vcov.georob
+# waldtest.georob
+
+# cv.georob (in separatem File)
+# rstudent.cv.georob
+# summary.cv.georob
+# print.summary.cv.georob
+
+# 2011-08-11 A. Papritz
+
+# ToDos:
+# - ...$xy durch ...[["xy"]] ersetzen
+
+## ##############################################################################
+
+model.frame.georob <-
+ function(
+ formula, ...
+ )
+{
+ ## model.frame method for class georob
+
+ ## 2012-12-19 A. Papritz
+
+ class( formula ) <- "lm"
+ model.frame( formula, ... )
+}
+
+## ##############################################################################
+
+model.matrix.georob <-
+ function(
+ object, ...
+ )
+{
+ ## model.matrix method for class georob
+
+ ## 2012-12-19 A. Papritz
+
+ class( object ) <- "lm"
+ model.matrix( object, ... )
+}
+
+## ##############################################################################
+
+nobs.georob <-
+ function(
+ object, ...
+ )
+{
+ ## nobs method for class georob
+
+ ## 2012-12-19 A. Papritz
+
+ class( object ) <- "lm"
+ nobs( object, ... )
+
+}
+
+## ##############################################################################
+
+print.georob <-
+ function(
+ x, digits = max(3, getOption("digits") - 3), ...
+ )
+{
+
+ ## Print method for class "georob".
+
+ ## Arguments:
+
+ ## x an object generated by f.georob.initial.guess
+ ## digits number of digits
+
+ ## 2011-08-13 A. Papritz
+ ## 2012-02-07 AP change for anisotropic variograms
+ ## 2012-12-18 AP invisible(x)
+
+ ## code borrowed from print.lmrob for printing fixed effects coeffcients
+
+ cat("\nTuning constant: ", x[["tuning.psi"]], "\n" )
+
+ cat("\nFixed effects coefficients:\n")
+
+ print(
+ format( coef( x), digits = digits ), print.gap = 2,
+ quote = FALSE
+ )
+
+ ## print variogram parameters
+
+ cat("\n")
+ cat( "Variogram: ", x$variogram.model, "\n" )
+ param <- x$param
+ names( param ) <- ifelse(
+ x$initial.objects$fit.param,
+ names( param ),
+ paste( names( param ), "(fixed)", sep = "" )
+ )
+ print(
+ format( param, digits = digits ), print.gap = 2,
+ quote = FALSE
+ )
+
+ ## print anisotropy parameters
+
+ if( !x$aniso$isotropic ){
+
+ cat("\n")
+ cat( "Anisotropy parameters: ", "\n" )
+ aniso <- x$aniso$aniso * c( rep(1, 2), rep( 180/pi, 3 ) )
+ names( aniso ) <- ifelse(
+ x$initial.objects$fit.aniso,
+ names( aniso ),
+ paste( names( aniso ), "(fixed)", sep = "" )
+ )
+ print(
+ format( aniso, digits = digits ), print.gap = 2,
+ quote = FALSE
+ )
+
+ }
+
+ invisible( x )
+
+}
+
+## ##############################################################################
+
+ranef.georob <- random.effects.georob <-
+ function(
+ object,
+ standard = FALSE,
+ ...
+ )
+{
+
+ ## Function extracts the random effects (bhat) from georob object
+ ## (similar to ranef.lme{nlme})
+
+ ## Arguments:
+
+ ## object fitted georob object
+ ## standard an optional logical value indicating whether the estimated random effects
+ ## should be "standardized" (i.e. divided by the estimated standard error.
+ ## Defaults to FALSE.
+ ## ... further arguments passed to method (currently not used)
+
+ ## 2011-10-13 A. Papritz
+ ## 2011-12-14 AP modified for replicated observations
+ ## 2012-01-05 AP modified for compress storage of matrices
+ ## 2012-10-18 AP changes for new definition of eta
+ ## 2012-11-26 AP method for random.effects
+ ## 2013-04-23 AP new names for robustness weights
+
+ object$Valpha.objects <- expand( object$Valpha.objects )
+ object$cov <- expand( object$cov )
+
+ bhat <- object$bhat
+
+ if( standard ){
+
+ if( is.null( object$cov$cov.bhat ) ){
+
+ ## compute standard errors of residuals
+
+ if( is.null( object$Valpha.objects$Valpha.inverse ) ||
+ is.null( object$Valpha.objects$Valpha.ilcf )
+ ) stop(
+ "'Valpha.objects' incomplete or missing in georob object;\n",
+ "'Valpha.objects' must include components 'Valpha.inverse' and 'Valpha.ilcf'"
+ )
+ if( is.null( object$expectations ) ) stop(
+ "'expectations' missing in georob object;\n",
+ "use 'full.output = TRUE' when fitting the object"
+ )
+
+
+ if( is.null( object$Valpha.objects$Valpha.ucf ) ){
+
+ ## compute upper cholesky factor of correlation matrix Valpha
+ ## which is needed to compute cov( bhat )
+
+ object$Valpha.objects$Valpha.ucf <- t( solve( object$Valpha.objects$Valpha.ilcf ) )
+
+ }
+
+ X <- model.matrix(
+ terms( object ),
+ model.frame( object )
+ )[!duplicated( object$Tmat ), , drop = FALSE]
+
+ r.cov <- compute.covariances(
+ Valpha.objects = object$Valpha.objects,
+ rweights = object$rweights,
+ XX = X, TT = object$Tmat, names.yy = rownames( model.frame( object ) ),
+ nugget = object$param["nugget"],
+ eta = sum( object$param[c( "variance", "snugget")] ) / object$param["nugget"],
+ expectations = object$expectations,
+ cov.bhat = TRUE, full.cov.bhat = FALSE,
+ cov.betahat = FALSE,
+ cov.bhat.betahat = FALSE,
+ cov.delta.bhat = FALSE, full.cov.delta.bhat = FALSE,
+ cov.delta.bhat.betahat = FALSE,
+ cov.ehat = FALSE, full.cov.ehat = FALSE,
+ cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
+ aux.cov.pred.target = FALSE,
+ verbose = 0
+ )
+
+ if( r.cov$error ) stop(
+ "an error occurred when computing the variances of the random effects"
+ )
+
+ se <- sqrt( r.cov$cov.bhat )
+
+ } else {
+
+ ## extract standard errors of residuals from georob object
+
+ if( is.matrix( object$cov$cov.bhat ) ){
+ se <- sqrt( diag( object$cov$cov.bhat ) )
+ } else {
+ se <- sqrt( object$cov$cov.bhat )
+ }
+
+ }
+
+ bhat <- bhat / se
+
+ }
+
+ return( bhat )
+
+}
+
+## ##############################################################################
+
+fixef.georob <- fixed.effects.georob <-
+ function(
+ object,
+ ...
+ )
+{
+
+ ## Function extracts residuals from georob object (based on residuals.lm {stats})
+
+ ## Arguments:
+
+ ## object fitted georob object
+ ## type character, type of resdiuals to computed
+ ## ... further arguments passed to methods
+
+ ## 2012-11-26 A. Papritz
+
+ object$coef
+
+}
+
+
+
+## ##############################################################################
+
+residuals.georob <- resid.georob <-
+ function(
+ object,
+ type = c("working", "response", "deviance", "pearson", "partial" ),
+ level = 1,
+ ...
+ )
+{
+
+ ## Function extracts residuals from georob object (based on residuals.lm {stats})
+
+ ## Arguments:
+
+ ## object fitted georob object
+ ## type character, type of resdiuals to computed
+ ## level integer scalar to select whether ehat (level == 1) or
+ ## ehat + bhat (level == 0) is returned,
+ ## only effective for type %in% c( "working", "response", "partial" )
+ ## ... further arguments passed to methods
+
+ ## 2011-10-13 A. Papritz
+ ## 2011-12-14 AP modified for replicated observations
+
+ type <- match.arg( type )
+
+ if( !level %in% 1:0 ) stop( "wrong level: must be either 1 or 0" )
+
+ r <- object$residuals
+ res <- switch(
+ type,
+ working = ,
+ response = r,
+ deviance = ,
+ pearson = if( is.null(object$weights) ) r else r * sqrt(object$weights),
+ partial = r
+ )
+
+ if( level == 0 && any( type %in% c( "working", "response", "partial" ) ) ){
+ res <- res + ranef( object, standard = FALSE )[object$Tmat]
+ }
+
+ res <- naresid(object$na.action, res)
+ if( type == "partial" )
+ res <- res + predict( object, type = "terms" )$fit
+ res
+}
+
+
+## ##############################################################################
+
+rstandard.georob <-
+ function( model, level = 1, ... )
+{
+
+ ## Function extracts standardized residuals from georob object
+
+ ## Arguments:
+
+ ## model fitted georob object
+ ## level integer scalar to select whether ehat (level == 1) or
+ ## ehat + bhat (level == 0) is returned,
+
+ ## ... further arguments (currently not used)
+
+ ## 2011-10-13 A. Papritz
+ ## 2011-12-14 AP modified for replicated observations
+ ## 2012-01-05 AP modified for compress storage of matrices
+ ## 2012-10-18 AP changes for new definition of eta
+ ## 2013-04-23 AP new names for robustness weights
+
+ object <- model
+ object$Valpha.objects <- expand( object$Valpha.objects )
+ object$cov <- expand( object$cov )
+
+ if( !level %in% 1:0 ) stop( "wrong level: must be either 1 or 0" )
+
+ if(
+ ( is.null( object$cov$cov.ehat ) & level == 1 ) ||
+ ( is.null( object$cov$cov.ehat.p.bhat ) & level == 0 )
+ ){
+
+ ## compute standard errors of residuals
+
+ if( is.null( object$Valpha.objects$Valpha.inverse ) ||
+ is.null( object$Valpha.objects$Valpha.ilcf )
+ ) stop(
+ "'Valpha.objects' incomplete or missing in georob object;\n",
+ "'Valpha.objects' must include components 'Valpha.inverse' and 'Valpha.ilcf'"
+ )
+ if( is.null( object$expectations ) ) stop(
+ "'expectations' missing in georob object;\n",
+ "use 'full.output = TRUE' when fitting the model"
+ )
+
+ X <- model.matrix(
+ terms( object),
+ model.frame( object )
+ )[!duplicated( object$Tmat ), , drop = FALSE]
+
+ if( is.null( object$Valpha.objects$Valpha.ucf ) ){
+ object$Valpha.objects$Valpha.ucf <- t( solve( object$Valpha.objects$Valpha.ilcf ) )
+ }
+
+ r.cov <- compute.covariances(
+ Valpha.objects = object$Valpha.objects,
+ rweights = object$rweights,
+ XX = X, TT = object$Tmat, names.yy = rownames( model.frame( object ) ),
+ nugget = object$param["nugget"],
+ eta = sum( object$param[c( "variance", "snugget")] ) / object$param["nugget"],
+ expectations = object$expectations,
+ cov.bhat = FALSE, full.cov.bhat = FALSE,
+ cov.betahat = FALSE,
+ cov.bhat.betahat = FALSE,
+ cov.delta.bhat = FALSE, full.cov.delta.bhat = FALSE,
+ cov.delta.bhat.betahat = FALSE,
+ cov.ehat = level == 1, full.cov.ehat = FALSE,
+ cov.ehat.p.bhat = level == 0, full.cov.ehat.p.bhat = FALSE,
+ aux.cov.pred.target = FALSE,
+ verbose = 0
+ )
+
+ if( r.cov$error ) stop(
+ "an error occurred when computing the variance of the residuals"
+ )
+
+ if( level == 1 ){
+ object$cov$cov.ehat <-r.cov$cov.ehat
+ } else {
+ object$cov$cov.ehat.p.bhat <-r.cov$cov.ehat.p.bhat
+ }
+
+ }
+
+ ## extract standard errors of residuals from georob object
+
+ if( level == 1 ){
+ se <- object$cov$cov.ehat
+ } else {
+ se <- object$cov$cov.ehat.p.bhat
+ }
+ if( is.matrix( se ) ){
+ se <- sqrt( diag( se ) )
+ } else {
+ se <- sqrt( se )
+ }
+
+ ## compute standardized residuals
+
+ residuals( model, level = level ) / se
+
+}
+
+## ##############################################################################
+
+rstudent.georob <-
+ function( model, ... )
+{
+
+ ## Function computes studentized residuals for fitted georob object
+
+ ## Arguments:
+
+ ## model fitted georob object
+ ## data data frame that was used to fit model
+ ## ... further arguments passed to cv.georob
+
+ ## 2011-12-22 A. Papritz
+
+ if( !identical( class( model )[1], "georob" ) ) stop(
+ "model is not of class 'georob'"
+ )
+
+ r.cv <- cv( model, ... )
+
+ rstudent( model = r.cv )
+
+}
+
+## ##############################################################################
+
+summary.georob <-
+ function (
+ object, correlation = FALSE,
+ signif = 0.95,
+ ...
+ )
+{
+
+ ## ToDos:
+
+ ## - Terms Objekt einfuegen
+ ## - ausgewaehlte Angaben zur Fitmethode ausgeben
+ ## - Wald-Test des Modells y ~ 1
+
+ ## 2012-01-05 A. Papritz
+ ## 2012-01-05 AP modified for compress storage of matrices
+ ## 2012-01-31 AP corretion of error for computing CI for variogram parameters
+ ## 2012-02-07 AP change for anisotropic variograms
+ ## 2012-03-29 AP warning messages inserted
+ ## 2012-05-23 ap correction of error for computing correlation matrix of variogram parameters
+ ## 2012-11-04 AP handling compressed cov.betahat
+ ## 2012-11-27 AP changes in parameter back-transformation
+ ## 2013-04-23 AP new names for robustness weights
+
+ object$cov <- expand( object$cov )
+
+ ans <- object[c(
+ "call", "residuals", "bhat", "rweights", "converged", "convergence.code",
+ "iter", "loglik", "variogram.model", "gradient",
+ "tuning.psi", "df.residual"
+ )]
+ ans <- c( ans, object$initial.objects["fit.param"] )
+
+ if( !object$aniso$isotropic ) ans$fit.param <- c(
+ ans$fit.param, object$initial.objects$fit.aniso
+ )
+
+ ans$terms <- NA
+ ans$scale <- sqrt(object$param["nugget"])
+ ans$control$method <- "TODO: PRINT GLSROB CONTROL PARAMETERS HERE"
+
+ se <- sqrt(diag(expand(object$cov$cov.betahat)))
+ est <- object$coefficients
+ tval <- est/se
+
+ ans$coefficients <- cbind(
+ est, se, tval, 2 * pt( abs(tval), object$df.residual, lower.tail = FALSE )
+ )
+ dimnames( ans$coefficients ) <- list(
+ names(est), c("Estimate", "Std. Error", "t value", "Pr(>|t|)")
+ )
+
+ if( correlation ){
+ ans$correlation <- expand( object$cov$cov.betahat ) / outer(se, se)
+ }
+
+ ans$param <- as.matrix( object$param, ncol = 1 )
+
+ if( !object$aniso$isotropic ) ans$param <- rbind(
+ ans$param,
+ as.matrix( object$aniso$aniso, ncol = 1 ) * c( rep( 1, 2 ), rep( 180/pi, 3 ) )
+ )
+
+ colnames( ans$param ) <- "Estimate"
+
+ ## compute confidence intervals of variogram parameters from observed
+ ## Fisher information matrix (Gaussian REML only)
+
+ if( !is.null( object$hessian ) ){
+
+ ## initialization
+
+ cor.tf.param <- cov.tf.param <- matrix(
+ NA, nrow = nrow( object$hessian ), ncol = nrow( object$hessian ),
+ dimnames = dimnames( object$hessian )
+ )
+
+ se <- rep( NA, nrow( object$hessian ) )
+ names( se ) <- rownames( object$hessian)
+
+ ci <- matrix( NA, nrow = nrow( ans$param ), ncol = 2 )
+ colnames( ci ) <- c( "Lower", "Upper" )
+ rownames( ci ) <- rownames( ans$param )
+
+ ## select parameters that are not on boundary of parameter space
+
+ sr <- !apply( object$hessian, 1, function( x ) all( is.na( x ) ) )
+
+ if( sum( sr ) > 0 ){
+
+ t.chol <- try( chol( object$hessian[sr, sr] ), silent = TRUE )
+
+ if( !identical( class( t.chol ), "try-error" ) ){
+
+ ## compute covariance matrix of fitted transformed parameters
+
+ cov.tf.param[sr, sr] <- chol2inv( t.chol )
+
+ ## correlation matrix and standard errors of fitted transformed
+ ## parameters
+
+ if( sum( sr ) > 1 ){
+ cor.tf.param[sr, sr] <- cov2cor( cov.tf.param[sr, sr] )
+ } else {
+ cor.tf.param[sr, sr] <- 1.
+ }
+
+ se[sr] <- sqrt( diag( cov.tf.param )[sr] )
+
+ ## compute confidence interval on original scale of parameters
+
+ sel.names <- names( object$param[object$initial.objects$fit.param] )
+ if( !object$aniso$isotropic ) sel.names <- c(
+ sel.names,
+ names( object$aniso$aniso[object$initial.objects$fit.aniso] )
+ )
+ sel.names <- sel.names[sr]
+
+ ff <- c( rep( 1, length( object$param ) + 2 ), rep( 180/pi, 3 ) )
+ names( ff ) <- names( c( object$param, object$aniso$aniso ) )
+
+ ci[sel.names, ] <- t(
+ sapply(
+ sel.names,
+ function( x, param, f, se, param.tf, trafo.fct, inv.trafo.fct ){
+ inv.trafo.fct[[param.tf[x]]](
+ trafo.fct[[param.tf[x]]]( param[x] ) +
+ c(-1, 1) * se[x] * qnorm( (1-signif)/2, lower.tail = FALSE )
+ )
+ },
+ param = c( object$param, object$aniso$aniso ),
+ f = ff,
+ se = se,
+ param.tf = object$param.tf,
+ trafo.fct = object$fwd.tf,
+ inv.trafo.fct = object$bwd.tf
+ )
+ )
+ is.angle <- rownames( ci ) %in% c( "omega", "phi", "zeta" )
+ if( sum(is.angle) > 0 ) ci[is.angle, ] <- ci[is.angle, ] * 180/pi
+
+ ans$param <- cbind( ans$param, ci )
+ if( correlation ) ans$cor.tf.param <- cor.tf.param
+
+ } else {
+ warning(
+ "Hessian not positive definite:",
+ "\nconfidence intervals of variogram parameters cannot be computed"
+ )
+ }
+ }
+ }
+
+ ans$se.residuals <- if( !is.null( object$cov$cov.ehat ) ){
+
+ if( is.matrix( object$cov$cov.ehat ) ){
+ sqrt( diag( object$cov$cov.ehat ) )
+ } else {
+ sqrt( object$cov$cov.ehat )
+ }
+
+ } else NULL
+
+ class( ans ) <- c( "summary.georob" )
+
+ ans
+}
+
+## ##############################################################################
+
+print.summary.georob <-
+ function (
+ x, digits = max(3, getOption("digits") - 3),
+ signif.stars = getOption("show.signif.stars"),
+ ...
+ )
+{
+
+ ## ToDos:
+
+ ## - Ausgabe df
+ ## - Wald-Test des Modells y ~ 1
+ ## - ausgewaehlte Angaben zur Fitmethode ausgeben
+
+ ## 2012-01-05 A. Papritz
+ ## 2012-01-31 AP small changes
+ ## 2012-02-07 AP change for anisotropic variograms
+ ## 2012-12-18 AP invisible(x)
+ ## 2013-04-23 AP new names for robustness weights
+
+
+ cat("\nCall:")
+ cat( paste( deparse(x$call), sep = "\n", collapse = "\n"), "\n", sep = "" )
+
+
+ cat("\nTuning constant: ", x[["tuning.psi"]], "\n" )
+
+ if( is.na( x$converged ) ){
+ cat( "\nEstimation with fixed variogram parameters\n" )
+
+ } else {
+
+ if(!(x$converged)) {
+ cat(
+ "\nAlgorithm did not converge, diagnostic code: ",
+ x$convergence.code, "\n"
+ )
+ } else {
+ cat(
+ "\nConvergence in", x$iter[1], "function and",
+ x$iter[2], "Jacobian/gradient evaluations\n"
+ )
+ }
+
+ attr( x$gradient, "eeq.emp" ) <- NULL
+ attr( x$gradient, "eeq.exp" ) <- NULL
+
+ cat( "\nEstimating equations (gradient)\n")
+ print( x$gradient, digits = digits, ... )
+
+ if( x$tuning.psi >=
+ georob.control()$tuning.psi.nr ) cat(
+ "\nMaximized restricted log-likelihood:",
+ x$loglik, "\n"
+ )
+
+ }
+
+ df <- x$df.residual
+
+ bhat <- x$bhat
+ cat( "\nPredicted latent variable (z):\n")
+ if(df > 5){
+ nam <- c("Min", "1Q", "Median", "3Q", "Max")
+ rq <- structure( quantile(bhat), names = nam )
+ print( rq, digits = digits, ...)
+ }
+ else print( bhat, digits = digits, ...)
+
+ resid <- x$residuals
+ cat( "\nResiduals (epsilon):\n")
+ if(df > 5){
+ nam <- c("Min", "1Q", "Median", "3Q", "Max")
+ rq <- structure( quantile(resid), names = nam )
+ print( rq, digits = digits, ...)
+ }
+ else print( resid, digits = digits, ...)
+
+ if( !is.null( x$se.residuals ) ){
+ resid <- x$residuals / x$se.residuals
+ cat( "\nStandardized residuals:\n")
+ if(df > 5){
+ nam <- c("Min", "1Q", "Median", "3Q", "Max")
+ rq <- structure( quantile(resid), names = nam )
+ print( rq, digits = digits, ...)
+ }
+ else print( resid, digits = digits, ...)
+ }
+
+ cat( "\nVariogram: ", x$variogram.model, "\n" )
+ rownames( x$param ) <- ifelse(
+ x$fit.param,
+ rownames( x$param ),
+ paste( rownames( x$param ), "(fixed)", sep = "" )
+ )
+ ## print( format( x$param, digits = digits ), print.gap = 2, quote = FALSE )
+ printCoefmat(
+ x$param, digits = digits, signif.stars = FALSE, ...
+ )
+
+
+ if( !is.null( x$cor.tf.param ) ){
+
+ correl <- x$cor.tf.param
+ p <- NCOL(correl)
+ if( p > 1 ){
+ cat("\nCorrelation of (transformed) variogram parameters:\n")
+ correl <- format(round(correl, 2), nsmall = 2,
+ digits = digits)
+ correl[!lower.tri(correl)] <- ""
+ print(correl[-1, -p, drop = FALSE], quote = FALSE)
+ }
+
+ }
+
+ cat( "\nFixed effects coefficients:\n" )
+ printCoefmat(
+ x$coefficients, digits = digits, signif.stars = signif.stars, ...
+ )
+
+ cat(
+ "\nResidual standard error (sqrt(nugget)):",
+ format(signif(x$scale, digits)), "\n"
+ )
+
+ correl <- x$correlation
+ if( !is.null(correl) ){
+ p <- NCOL(correl)
+ if( p > 1 ){
+ cat("\nCorrelation of fixed effects coefficients:\n")
+ correl <- format(round(correl, 2), nsmall = 2,
+ digits = digits)
+ correl[!lower.tri(correl)] <- ""
+ print(correl[-1, -p, drop = FALSE], quote = FALSE)
+ }
+ }
+
+ cat("\n")
+ summarizeRobWeights(x$rweights, digits = digits, ... )
+
+ invisible( x )
+}
+
+## ##############################################################################
+
+vcov.georob <-
+ function( object, ... )
+{
+
+ ## 2012-11-04 AP handling compressed cov.betahat
+
+ result <- expand( object$cov$cov.betahat )
+ attr( result, "struc" ) <- NULL
+ result
+
+}
+
+## ##############################################################################
+
+waldtest.georob <-
+ function(
+ object, ..., vcov = NULL, test = c("Chisq", "F"), name = NULL,
+ fixed = TRUE, verbose = 1
+ )
+{
+
+ ## 2012-02-08 AP change for anisotropic variograms
+
+ ## refit model with fixed variogram parameters
+
+ test <- match.arg( test )
+
+ if( fixed ) {
+
+ if( verbose > 0 ) cat( "\nWald-Test with fixed variogram parameters of model 1\n\n" )
+
+ object <- update(
+ object,
+ 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 = 0
+ )
+
+ }
+
+ ## Wald-Test
+
+ waldtest.default(
+ object = object, ..., vcov = vcov, test = test, name = name
+ )
+
+}
+
+## ##############################################################################
+
+logLik.georob <-
+ function( object, REML = FALSE, ... )
+{
+
+ ## 2012-12-22 method for extracting (restricted) loglikelihood
+
+ val <- if( REML ){
+ val <- object$loglik
+ } else if( object[["tuning.psi"]] >= georob.control()[["tuning.psi.nr"]] ){
+ D <- deviance( object )
+ -0.5 * (
+ D + attr( D, "log.det.covmat" ) + length( object$residuals ) * log( 2 * pi )
+ )
+ } else NA_real_
+
+ attr(val, "nall") <- length(object$residuals)
+ attr(val, "nobs") <- object$df.residual
+ attr(val, "df") <- length(object$coefficients) +
+ sum( object$initial.objects$fit.param ) +
+ sum( object$initial.objects$fit.aniso)
+
+ class(val) <- "logLik"
+ val
+
+}
+
+## ##############################################################################
+
+deviance.georob <-
+ function(
+ object, ...
+ )
+{
+ ## deviance method for class georob
+
+ ## 2012-12-22 A. Papritz
+
+ if( object[["tuning.psi"]] < georob.control()[["tuning.psi.nr"]] ){
+ result <- NA_real_
+ } else {
+ object[["Valpha.objects"]] <- expand( object[["Valpha.objects"]] )
+ G <- sum( object[["param"]][c("variance", "snugget")] ) *
+ t(object[["Valpha.objects"]][["Valpha.ucf"]]) %*% object[["Valpha.objects"]][["Valpha.ucf"]]
+
+ diag( G ) <- diag( G ) + object[["param"]]["nugget"]
+ object[["Valpha.objects"]] <- compress( object[["Valpha.objects"]] )
+ G <- G[object[["Tmat"]], object[["Tmat"]]]
+ iucf <- try(
+ backsolve(
+ chol( G ),
+ diag( length( object[["Tmat"]] ) ),
+ k = length( object[["Tmat"]] )
+ ),
+ silent = TRUE
+ )
+ if( identical( class( iucf ), "try-error" ) ) {
+ stop( "(generalized) covariance matrix of observations not positive definite" )
+ } else {
+ result <- sum( colSums( residuals( object, level = 0 ) * iucf )^2 )
+ attr( result, "log.det.covmat" ) <- -2 * sum( log( diag( iucf ) ) )
+ }
+ }
+ result
+}
+
+
+
+
Property changes on: pkg/R/georob.S3methods.R
___________________________________________________________________
Added: svn:executable
+
Added: pkg/R/georob.exported.functions.R
===================================================================
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/georob -r 2
More information about the Georob-commits
mailing list