[Georob-commits] r6 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 6 15:28:14 CEST 2013
Author: papritz
Date: 2013-06-06 15:28:13 +0200 (Thu, 06 Jun 2013)
New Revision: 6
Modified:
pkg/ChangeLog
pkg/NAMESPACE
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/man/cv.georob.Rd
pkg/man/georob.control.Rd
pkg/man/georobObject.Rd
pkg/man/internal.functions.Rd
pkg/man/predict.georob.Rd
Log:
handling design matrices with rank < ncol(x)
changes for solving estimating equations for xi
M pkg/R/georob.cv.R
M pkg/R/georob.S3methods.R
M pkg/R/georob.predict.R
M pkg/R/georob.exported.functions.R
M pkg/R/georob.private.functions.R
M pkg/ChangeLog
M pkg/man/internal.functions.Rd
M pkg/man/cv.georob.Rd
M pkg/man/georob.control.Rd
M pkg/man/georobObject.Rd
M pkg/man/predict.georob.Rd
M pkg/NAMESPACE
Modified: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog 2013-05-23 15:41:57 UTC (rev 5)
+++ pkg/ChangeLog 2013-06-06 13:28:13 UTC (rev 6)
@@ -68,8 +68,25 @@
* georob.exported.functions.R (georob): improved way to handle missing observations and to construct model.frame
* georob.predict.R (predict.georob): correct handling of missing observations
* georob.S3methods.R (georob.residuals): new argument "terms"
-* georob.S3methods.R (ranef.georob, rstandard.georob,deviance.georob): correct handling of missing observations
+* georob.S3methods.R (ranef.georob, residuals.georob, rstandard.georob, deviance.georob): correct handling of missing observations
* variogram.R (plot.georob): correct handling of missing observations
+2013-05-24 Andreas Papritz <papritz at env.ethz.ch>
+* georob.cv.R (cv.georob): separate initial variogram parameters for each cross-validation set
+
+
+2013-05-31 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.S3methods.R (ranef.georob, residuals.georob,rstandard.georob,deviance.georob): correct handling of missing observations
+* georob.S3methods.R (deviance.georob, ranef.georob, rstandard.georob, summary.georob): revised expansion of covariance matrices
+
+
+2013-06-06 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.exported.function.R (georob, georob.control): handling fixed effects model matrices with rank < ncol(x)
+* georob.private.function.R (prepare.likelihood.calculations, compute.estimating.equations, negative.restr.loglikelihood, negative.restr.loglikelihood, gradient.negative.restricted.loglikelihood, georob.fit) : solving estimating equations for xi
+* georob.S3methods.R ranef.georob, rstandard.georob) : solving estimating equations for xi
+* georob.private.function.R (compute.covariances, estimate.xihat, update.xihat, negative.restr.loglikelihood): solving estimating equations for xi
+
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-05-23 15:41:57 UTC (rev 5)
+++ pkg/NAMESPACE 2013-06-06 13:28:13 UTC (rev 6)
@@ -76,16 +76,15 @@
## compute.covariances,
## compute.estimating.equations,
## compute.semivariance,
-## compute.U,
## dcorr.dparam,
-## estimate.betahat.bhat,
+## estimate.xihat,
## gcr,
## georob.fit,
## getCall.georob,
## gradient.negative.restricted.loglikelihood,
## negative.restr.loglikelihood,
## prepare.likelihood.calculations,
-## update.betahat.bhat
+## update.xihat
S3method( deviance, georob )
S3method( fixed.effects, georob )
Modified: pkg/R/georob.S3methods.R
===================================================================
--- pkg/R/georob.S3methods.R 2013-05-23 15:41:57 UTC (rev 5)
+++ pkg/R/georob.S3methods.R 2013-06-06 13:28:13 UTC (rev 6)
@@ -169,21 +169,30 @@
## 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
- ## 2013-05-23 AP correct handling of missing observations
+ ## 2013-05-31 AP correct handling of missing observations
+ ## 2013-05-31 AP revised expansion of covariance matrices
+ ## 2013-05-06 AP changes for solving estimating equations for xi
- object$Valpha.objects <- expand( object$Valpha.objects )
- object$cov <- expand( object$cov )
+ ## temporarily redefine na.action component of object
+ object.na <- object$na.action
+ if( identical( class( object$na.action ), "exclude" ) ){
+ class( object$na.action ) <- "omit"
+ }
+
+ Valpha.objects <- expand( object$Valpha.objects )
+ covmat <- expand( object$cov )
+
bhat <- object$bhat
if( standard ){
- if( is.null( object$cov$cov.bhat ) ){
+ if( is.null( covmat$cov.bhat ) ){
## compute standard errors of residuals
- if( is.null( object$Valpha.objects$Valpha.inverse ) ||
- is.null( object$Valpha.objects$Valpha.ilcf )
+ if( is.null( Valpha.objects$Valpha.inverse ) ||
+ is.null( Valpha.objects$Valpha.ilcf )
) stop(
"'Valpha.objects' incomplete or missing in georob object;\n",
"'Valpha.objects' must include components 'Valpha.inverse' and 'Valpha.ilcf'"
@@ -194,12 +203,12 @@
)
- if( is.null( object$Valpha.objects$Valpha.ucf ) ){
+ if( is.null( 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 ) )
+ Valpha.objects$Valpha.ucf <- t( solve( Valpha.objects$Valpha.ilcf ) )
}
@@ -209,7 +218,9 @@
)[!duplicated( object$Tmat ), , drop = FALSE]
r.cov <- compute.covariances(
- Valpha.objects = object$Valpha.objects,
+ Valpha.objects = Valpha.objects,
+ Aalpha = object[["Aalpha"]],
+ Palpha = expand( object[["Palpha"]] ),
rweights = object$rweights,
XX = X, TT = object$Tmat, names.yy = rownames( model.frame( object ) ),
nugget = object$param["nugget"],
@@ -218,8 +229,8 @@
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.deltabhat = FALSE, full.cov.deltabhat = FALSE,
+ cov.deltabhat.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,
@@ -236,10 +247,10 @@
## extract standard errors of residuals from georob object
- if( is.matrix( object$cov$cov.bhat ) ){
- se <- sqrt( diag( object$cov$cov.bhat ) )
+ if( is.matrix( covmat$cov.bhat ) ){
+ se <- sqrt( diag( covmat$cov.bhat ) )
} else {
- se <- sqrt( object$cov$cov.bhat )
+ se <- sqrt( covmat$cov.bhat )
}
}
@@ -248,8 +259,7 @@
}
- bhat <- naresid( object$na.action, bhat )
- return( bhat )
+ naresid( object.na, bhat )
}
@@ -303,12 +313,19 @@
## 2011-10-13 A. Papritz
## 2011-12-14 AP modified for replicated observations
- ## 2013-05-23 AP modified for computing partial residuals for single terms
+ ## 2013-05-31 AP modified for computing partial residuals for single terms
type <- match.arg( type )
if( !level %in% 1:0 ) stop( "wrong level: must be either 1 or 0" )
+ ## temporarily redefine na.action component of object
+
+ object.na <- object$na.action
+ if( identical( class( object$na.action ), "exclude" ) ){
+ class( object$na.action ) <- "omit"
+ }
+
r <- object$residuals
res <- switch(
type,
@@ -319,7 +336,6 @@
partial = r
)
- res <- naresid(object$na.action, res)
if( level == 0 && any( type %in% c( "working", "response", "partial" ) ) ){
res <- res + ranef( object, standard = FALSE )[object$Tmat]
@@ -328,6 +344,8 @@
if( type == "partial" )
res <- res + predict( object, type = "terms", terms = terms )$fit
drop( res )
+
+ naresid( object.na, res )
}
@@ -352,53 +370,63 @@
## 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
- ## 2013-05-23 AP correct handling of missing observations
+ ## 2013-05-31 AP correct handling of missing observations
+ ## 2013-05-31 AP revised expansion of covariance matrices
+ ## 2013-05-06 AP changes for solving estimating equations for xi
- object <- model
- object$Valpha.objects <- expand( object$Valpha.objects )
- object$cov <- expand( object$cov )
+ ## temporarily redefine na.action component of model
+ model.na <- model$na.action
+ if( identical( class( model$na.action ), "exclude" ) ){
+ class( model$na.action ) <- "omit"
+ }
+
+ Valpha.objects <- expand( model$Valpha.objects )
+ covmat <- expand( model$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 )
+ ( is.null( covmat$cov.ehat ) & level == 1 ) ||
+ ( is.null( covmat$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 )
+ if( is.null( Valpha.objects$Valpha.inverse ) ||
+ is.null( 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(
+ if( is.null( model$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]
+ terms( model),
+ model.frame( model )
+ )[!duplicated( model$Tmat ), , drop = FALSE]
- if( is.null( object$Valpha.objects$Valpha.ucf ) ){
- object$Valpha.objects$Valpha.ucf <- t( solve( object$Valpha.objects$Valpha.ilcf ) )
+ if( is.null( Valpha.objects$Valpha.ucf ) ){
+ Valpha.objects$Valpha.ucf <- t( solve( 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,
+ Valpha.objects = Valpha.objects,
+ Aalpha = model[["Aalpha"]],
+ Palpha = expand( model[["Palpha"]] ),
+ rweights = model$rweights,
+ XX = X, TT = model$Tmat, names.yy = rownames( model.frame( model ) ),
+ nugget = model$param["nugget"],
+ eta = sum( model$param[c( "variance", "snugget")] ) / model$param["nugget"],
+ expectations = model$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.deltabhat = FALSE, full.cov.deltabhat = FALSE,
+ cov.deltabhat.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,
@@ -410,9 +438,9 @@
)
if( level == 1 ){
- object$cov$cov.ehat <-r.cov$cov.ehat
+ covmat$cov.ehat <-r.cov$cov.ehat
} else {
- object$cov$cov.ehat.p.bhat <-r.cov$cov.ehat.p.bhat
+ covmat$cov.ehat.p.bhat <-r.cov$cov.ehat.p.bhat
}
}
@@ -420,9 +448,9 @@
## extract standard errors of residuals from georob object
if( level == 1 ){
- se <- object$cov$cov.ehat
+ se <- covmat$cov.ehat
} else {
- se <- object$cov$cov.ehat.p.bhat
+ se <- covmat$cov.ehat.p.bhat
}
if( is.matrix( se ) ){
se <- sqrt( diag( se ) )
@@ -432,10 +460,8 @@
## compute standardized residuals
- se <- naresid( model$na.action, se )
+ naresid( model.na, residuals( model, level = level ) / se )
- residuals( model, level = level ) / se
-
}
## ##############################################################################
@@ -489,8 +515,9 @@
## 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
+ ## 2013-05-31 AP revised expansion of covariance matrices
- object$cov <- expand( object$cov )
+ covmat <- expand( object$cov )
ans <- object[c(
"call", "residuals", "bhat", "rweights", "converged", "convergence.code",
@@ -507,7 +534,7 @@
ans$scale <- sqrt(object$param["nugget"])
ans$control$method <- "TODO: PRINT GLSROB CONTROL PARAMETERS HERE"
- se <- sqrt(diag(expand(object$cov$cov.betahat)))
+ se <- sqrt(diag(covmat$cov.betahat))
est <- object$coefficients
tval <- est/se
@@ -519,7 +546,7 @@
)
if( correlation ){
- ans$correlation <- expand( object$cov$cov.betahat ) / outer(se, se)
+ ans$correlation <- covmat$cov.betahat / outer(se, se)
}
ans$param <- as.matrix( object$param, ncol = 1 )
@@ -619,12 +646,12 @@
}
}
- ans$se.residuals <- if( !is.null( object$cov$cov.ehat ) ){
+ ans$se.residuals <- if( !is.null( covmat$cov.ehat ) ){
- if( is.matrix( object$cov$cov.ehat ) ){
- sqrt( diag( object$cov$cov.ehat ) )
+ if( is.matrix( covmat$cov.ehat ) ){
+ sqrt( diag( covmat$cov.ehat ) )
} else {
- sqrt( object$cov$cov.ehat )
+ sqrt( covmat$cov.ehat )
}
} else NULL
@@ -876,6 +903,7 @@
## 2012-12-22 A. Papritz
## 2013-05-23 AP correct handling of missing observations
+ ## 2013-05-31 AP revised expansion of covariance matrices
## redefine na.action component of object
@@ -886,12 +914,11 @@
if( object[["tuning.psi"]] < georob.control()[["tuning.psi.nr"]] ){
result <- NA_real_
} else {
- object[["Valpha.objects"]] <- expand( object[["Valpha.objects"]] )
+ Valpha.objects <- expand( object[["Valpha.objects"]] )
G <- sum( object[["param"]][c("variance", "snugget")] ) *
- t(object[["Valpha.objects"]][["Valpha.ucf"]]) %*% object[["Valpha.objects"]][["Valpha.ucf"]]
+ t(Valpha.objects[["Valpha.ucf"]]) %*% 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(
Modified: pkg/R/georob.cv.R
===================================================================
--- pkg/R/georob.cv.R 2013-05-23 15:41:57 UTC (rev 5)
+++ pkg/R/georob.cv.R 2013-06-06 13:28:13 UTC (rev 6)
@@ -71,6 +71,7 @@
## 2012-12-04 AP modifiction for changes in predict.georob
## 2013-04-24 AP changes for parallelization on windows os
## 2013-05-23 AP correct handling of missing observations
+ ## 2013-05-24 AP separate initial variogram parameters for each cross-validation set
## auxiliary function that fits the model and computes the predictions of
## a cross-validation set
@@ -96,6 +97,11 @@
environment( formula ) <- environment()
environment( object$terms ) <- environment()
+
+ ## read-off initial values of variogram parameters
+
+ if( ( is.matrix( param ) || is.data.frame( param ) ) ) param <- param[..i..,]
+ if( ( is.matrix( fit.param ) || is.data.frame( fit.param ) ) ) fit.param <- fit.param[..i..,]
t.georob <- update(
object,
@@ -212,6 +218,7 @@
} else {
+ nset <- length( unique( sets ) )
if( length( sets ) != NROW( data ) ) stop(
"sets must be an integer vector with length equal to the number of observations"
)
@@ -230,6 +237,16 @@
function( x ) x
)
+ ## check dimension of param and fit.param
+
+ if( ( is.matrix( param ) || is.data.frame( param ) ) && nrow( param )!= nset ) stop(
+ "param must have 'nset' rows if it is a matrix or data frame"
+ )
+
+ if( ( is.matrix( fit.param ) || is.data.frame( fit.param ) ) && nrow( param )!= nset ) stop(
+ "fit.param must have 'nset' rows if it is a matrix or data frame"
+ )
+
## loop over all cross-validation sets
if( .Platform$OS.type == "windows" ){
@@ -320,9 +337,9 @@
t.fit <- lapply( t.result, function( x ) return( x$fit ) )
- if( re.estimate && !all( sapply( t.fit, function(x) x$converged ) ) )
- warning(
- "lack of covergence when fitting model to cross-validation sets"
+ if( re.estimate && !all( sapply( t.fit, function(x) x$converged ) ) ) warning(
+ "lack of covergence for ",
+ sum( !sapply( t.fit, function(x) x$converged ) ), " cross-validation sets"
)
result <- list(
Modified: pkg/R/georob.exported.functions.R
===================================================================
--- pkg/R/georob.exported.functions.R 2013-05-23 15:41:57 UTC (rev 5)
+++ pkg/R/georob.exported.functions.R 2013-06-06 13:28:13 UTC (rev 6)
@@ -76,6 +76,7 @@
## 2012-05-28 AP handle missing names of coefficients after calling update
## 2013-04-23 AP new names for robustness weights
## 2013-05-23 AP correct handling of missing observations and to construct model.frame
+ ## 2013-06-03 AP handling design matrices with rank < ncol(x)
## check whether input is complete
@@ -177,12 +178,22 @@
min.max.sv <- range( svd( crossprod( x ) )$d )
condnum <- min.max.sv[1] / min.max.sv[2]
+ if( condnum <= control$min.condnum ){
+ if( initial.param || tuning.psi >= control$tuning.psi.nr ) stop(
+ "singular fixed effects design matrices cannot be handled if 'initial.param = TRUE'",
+ "or for Gaussian REML estimation"
+ )
+ cat(
+ "design matrix has not full column rank (condition number of X^T X: ",
+ signif( condnum, 2 ), ")\ninitial values of regression coefficients are computed by 'lm\n\n'"
+ )
+ control$initial.method <- "lm"
+ warning(
+ "design matrix has not full column rank (condition number of X^T X: ",
+ signif( condnum, 2 ), ")\ninitial values of regression coefficients are computed by 'lm'"
+ )
+ }
- if( condnum <= control$min.condnum ) stop(
- "design matrix has not full column rank (condition number of X^T X: ",
- signif( condnum, 2 ), ")"
- )
-
## subtract offset
yy <- y
@@ -224,7 +235,7 @@
fit$formula <- formula
fit$terms <- mt
fit$xlevels <- .getXlevels(mt, mf)
- fit$call <- call
+ fit$call <- cl
fit$tau <- tau
fit$weights <- w
fit$residuals <- drop( fit$residuals )
@@ -250,6 +261,27 @@
if( !is.null( offset ) ) fit$fitted.values + offset
fit
+ },
+ lm = {
+
+ fit <- if( is.null(w) ){
+ lm.fit(x, y, offset = offset, singular.ok = TRUE, ...)
+ } else {
+ lm.wfit(x, y, w, offset = offset, singular.ok = TRUE, ...)
+ }
+ class(fit) <- c(if (is.matrix(y)) "mlm", "lm")
+ fit$na.action <- attr(mf, "na.action")
+ fit$offset <- offset
+ fit$contrasts <- attr(x, "contrasts")
+ fit$xlevels <- .getXlevels(mt, mf)
+ fit$call <- cl
+ fit$terms <- mt
+ if (model) fit$model <- mf
+ if (ret.x) fit$x <- x
+ if (ret.y) fit$y <- y
+ fit$qr <- NULL
+ fit
+
}
)
@@ -336,12 +368,13 @@
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.deltabhat = control$cov.deltabhat,
+ full.cov.deltabhat = control$full.cov.deltabhat,
+ cov.deltabhat.betahat = control$cov.deltabhat.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,
psi.func = control$psi.func,
tuning.psi.nr = control$tuning.psi.nr,
irwls.initial = control$irwls.initial,
@@ -400,12 +433,13 @@
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.deltabhat = control$cov.deltabhat,
+ full.cov.deltabhat = control$full.cov.deltabhat,
+ cov.deltabhat.betahat = control$cov.deltabhat.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,
psi.func = control$psi.func,
tuning.psi.nr = control$tuning.psi.nr,
irwls.initial = control$irwls.initial,
@@ -469,7 +503,7 @@
georob.control <-
function(
- initial.method = c("lmrob", "rq"),
+ initial.method = c("lmrob", "rq", "lm"),
bhat = NULL,
param.tf = param.transf(),
fwd.tf = fwd.transf(),
@@ -486,8 +520,8 @@
cov.bhat = FALSE, full.cov.bhat = FALSE,
cov.betahat = TRUE,
cov.bhat.betahat = FALSE,
- cov.delta.bhat = TRUE, full.cov.delta.bhat = TRUE,
- cov.delta.bhat.betahat = TRUE,
+ cov.deltabhat = TRUE, full.cov.deltabhat = TRUE,
+ cov.deltabhat.betahat = TRUE,
cov.ehat = TRUE, full.cov.ehat = FALSE,
cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
aux.cov.pred.target = FALSE,
@@ -535,10 +569,10 @@
## should be computed
## cov.bhat.betahat logical, flag controlling whether the covariance matrix of
## bhat and betahat should be computed
- ## cov.delta.bhat logical, flag controlling whether the covariances of z-bhat should be computed
- ## full.cov.delta.bhat logical, flag controlling whether the full covariance matrix of z-bhat
+ ## cov.deltabhat logical, flag controlling whether the covariances of z-bhat should be computed
+ ## full.cov.deltabhat logical, flag controlling whether the full covariance matrix of z-bhat
## is computed (TRUE) or only the diagonal elements (FALSE)
- ## cov.delta.bhat.betahat logical, flag controlling whether the covariance matrix of z-bhat
+ ## cov.deltabhat.betahat logical, flag controlling whether the covariance matrix of z-bhat
## and betahat should be computed
## cov.ehat logical, flag controlling whether the covariances of the resdiuals should be computed
## full.cov.ehat logical, flag controlling whether the full covariance matrix of the residuals
@@ -587,8 +621,8 @@
cov.bhat = cov.bhat, full.cov.bhat = full.cov.bhat,
cov.betahat = cov.betahat,
cov.bhat.betahat = cov.bhat.betahat,
- cov.delta.bhat = cov.delta.bhat, full.cov.delta.bhat = full.cov.delta.bhat,
- cov.delta.bhat.betahat = cov.delta.bhat.betahat,
+ cov.deltabhat = cov.deltabhat, full.cov.deltabhat = full.cov.deltabhat,
+ cov.deltabhat.betahat = cov.deltabhat.betahat,
cov.ehat = cov.ehat, full.cov.ehat = full.cov.ehat,
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,
Modified: pkg/R/georob.predict.R
===================================================================
--- pkg/R/georob.predict.R 2013-05-23 15:41:57 UTC (rev 5)
+++ pkg/R/georob.predict.R 2013-06-06 13:28:13 UTC (rev 6)
@@ -93,7 +93,7 @@
locations.coords, betahat, bhat,
pred.X, pred.coords, newdata,
variogram.model, param, aniso,
- cov.dbhat.betahat.l, cov.betahat.l, cov.bhat.betahat, cov.p.t, Valpha.objects,
+ cov.deltabhat.betahat.l, cov.betahat.l, cov.bhat.betahat, cov.p.t, Valpha.objects,
pwidth, pheight, napp,
signif,
extended.output, full.covmat
@@ -385,8 +385,8 @@
## compute uk variance (= (co-)variance of prediction errors)
aux <- cbind(
- gammaValphai %*% cov.dbhat.betahat.l[1:n, 1:n] - pred.X[!ex, , drop = FALSE ] %*% cov.dbhat.betahat.l[-(1:n), 1:n],
- - pred.X[!ex, , drop = FALSE ] %*% cov.dbhat.betahat.l[-(1:n), -(1:n)]
+ gammaValphai %*% cov.deltabhat.betahat.l[1:n, 1:n] - pred.X[!ex, , drop = FALSE ] %*% cov.deltabhat.betahat.l[-(1:n), 1:n],
+ - pred.X[!ex, , drop = FALSE ] %*% cov.deltabhat.betahat.l[-(1:n), -(1:n)]
)
if( full.covmat ){
@@ -630,16 +630,16 @@
## if needed compute missing covariance matrices
cov.betahat <- is.null( object[["cov"]][["cov.betahat"]] )
- cov.dbhat <- is.null( object[["cov"]][["cov.delta.bhat"]] ) ||
- !is.matrix( object[["cov"]][["cov.delta.bhat"]] )
- cov.dbhat.betahat <- is.null( object[["cov"]][["cov.delta.bhat.betahat"]] )
+ cov.deltabhat <- is.null( object[["cov"]][["cov.deltabhat"]] ) ||
+ !is.matrix( object[["cov"]][["cov.deltabhat"]] )
+ cov.deltabhat.betahat <- is.null( object[["cov"]][["cov.deltabhat.betahat"]] )
cov.bhat <- extended.output & (
is.null( object[["cov"]]$cov.bhat ) || !is.matrix( object[["cov"]]$cov.bhat )
)
cov.bhat.betahat <- extended.output & is.null( object[["cov"]]$cov.bhat.betahat )
cov.p.t <- extended.output & is.null( object[["cov"]]$aux.cov.pred.target )
- if( any( c( cov.betahat, cov.dbhat, cov.dbhat.betahat,
+ if( any( c( cov.betahat, cov.deltabhat, cov.deltabhat.betahat,
extended.output & ( cov.bhat || cov.bhat.betahat || cov.p.t )
)
)
@@ -676,8 +676,8 @@
cov.bhat = cov.bhat, full.cov.bhat = cov.bhat,
cov.betahat = cov.betahat,
cov.bhat.betahat = cov.bhat.betahat,
- cov.delta.bhat = cov.dbhat, full.cov.delta.bhat = cov.dbhat,
- cov.delta.bhat.betahat = cov.dbhat.betahat,
+ cov.deltabhat = cov.deltabhat, full.cov.deltabhat = cov.deltabhat,
+ cov.deltabhat.betahat = cov.deltabhat.betahat,
cov.ehat = FALSE, full.cov.ehat = FALSE,
cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
aux.cov.pred.target = cov.p.t,
@@ -691,9 +691,9 @@
if( is.null( object[["cov"]] ) ) object[["cov"]] <- list()
if( cov.betahat ) object[["cov"]][["cov.betahat"]] <- r.cov[["cov.betahat"]]
- if( cov.dbhat ) object[["cov"]][["cov.delta.bhat"]] <- r.cov[["cov.delta.bhat"]]
- if( cov.dbhat.betahat ) object[["cov"]][["cov.delta.bhat.betahat"]] <-
- r.cov[["cov.delta.bhat.betahat"]]
+ if( cov.deltabhat ) object[["cov"]][["cov.delta"]] <- r.cov[["cov.delta"]]
+ if( cov.deltabhat.betahat ) object[["cov"]][["cov.deltabhat.betahat"]] <-
+ r.cov[["cov.deltabhat.betahat"]]
if( extended.output && cov.bhat ) object[["cov"]][["cov.bhat"]] <- r.cov[["cov.bhat"]]
if( extended.output && cov.bhat.betahat ) object[["cov"]][["cov.bhat.betahat"]] <-
r.cov[["cov.bhat.betahat"]]
@@ -702,26 +702,26 @@
} ## end cov
- ## compute lower cholesky factor of covariance matrix of delta.bhat = (b -
+ ## compute lower cholesky factor of covariance matrix of delta = (b -
## bhat) and betahat - beta
- cov.dbhat.betahat.l <- try(
+ cov.deltabhat.betahat.l <- try(
t(
chol(
rbind(
cbind(
- object[["cov"]][["cov.delta.bhat"]],
- object[["cov"]][["cov.delta.bhat.betahat"]]
+ object[["cov"]][["cov.delta"]],
+ object[["cov"]][["cov.deltabhat.betahat"]]
),
cbind(
- t( object[["cov"]][["cov.delta.bhat.betahat"]] ),
+ t( object[["cov"]][["cov.deltabhat.betahat"]] ),
object[["cov"]][["cov.betahat"]]
)
)
)
), silent = TRUE
)
- if( identical( class( cov.dbhat.betahat.l ), "try-error" ) ) stop(
+ if( identical( class( cov.deltabhat.betahat.l ), "try-error" ) ) stop(
"covariance matrix of kriging errors 'b-bhat' and 'betahat' not positive definite"
)
@@ -890,8 +890,8 @@
},
signal = { ## signal
aux <- cbind(
- cov.dbhat.betahat.l[1:n,1:n] - X %*% cov.dbhat.betahat.l[-(1:n),1:n],
- - X %*% cov.dbhat.betahat.l[-(1:n),-(1:n)]
+ cov.deltabhat.betahat.l[1:n,1:n] - X %*% cov.deltabhat.betahat.l[-(1:n),1:n],
+ - X %*% cov.deltabhat.betahat.l[-(1:n),-(1:n)]
)
aux <- aux[object[["Tmat"]], , drop = FALSE]
if( full.covmat ){
@@ -1137,7 +1137,7 @@
locations.coords, betahat, bhat,
pred.X, pred.coords, newdata,
variogram.model, param, aniso,
- cov.dbhat.betahat.l, cov.betahat.l, cov.bhat.betahat, cov.p.t, Valpha.objects,
+ cov.deltabhat.betahat.l, cov.betahat.l, cov.bhat.betahat, cov.p.t, Valpha.objects,
pwidth, pheight, napp,
signif,
extended.output, full.covmat,
@@ -1164,7 +1164,7 @@
bhat = bhat,
pred.X = pred.X, pred.coords = pred.coords, newdata = newdata,
variogram.model = variogram.model, param = param, aniso = aniso,
- cov.dbhat.betahat.l = cov.dbhat.betahat.l,
+ cov.deltabhat.betahat.l = cov.deltabhat.betahat.l,
cov.betahat.l = cov.betahat.l,
cov.bhat.betahat = cov.bhat.betahat,
cov.p.t = cov.p.t,
@@ -1205,7 +1205,7 @@
variogram.model = object[["variogram.model"]],
param = object[["param"]],
aniso = object[["aniso"]],
- cov.dbhat.betahat.l = cov.dbhat.betahat.l,
+ cov.deltabhat.betahat.l = cov.deltabhat.betahat.l,
cov.betahat.l = cov.betahat.l,
cov.bhat.betahat = cov.bhat.betahat,
cov.p.t = cov.p.t,
@@ -1235,7 +1235,7 @@
variogram.model = object[["variogram.model"]],
param = object[["param"]],
aniso = object[["aniso"]],
- cov.dbhat.betahat.l = cov.dbhat.betahat.l,
+ cov.deltabhat.betahat.l = cov.deltabhat.betahat.l,
cov.betahat.l = cov.betahat.l,
cov.bhat.betahat = cov.bhat.betahat,
cov.p.t = cov.p.t,
@@ -1264,7 +1264,7 @@
variogram.model = object[["variogram.model"]],
param = object[["param"]],
aniso = object[["aniso"]],
- cov.dbhat.betahat.l = cov.dbhat.betahat.l,
+ cov.deltabhat.betahat.l = cov.deltabhat.betahat.l,
cov.betahat.l = cov.betahat.l,
cov.bhat.betahat = cov.bhat.betahat,
cov.p.t = cov.p.t,
Modified: pkg/R/georob.private.functions.R
===================================================================
--- pkg/R/georob.private.functions.R 2013-05-23 15:41:57 UTC (rev 5)
+++ pkg/R/georob.private.functions.R 2013-06-06 13:28:13 UTC (rev 6)
@@ -9,15 +9,15 @@
compute.covariances <-
function(
Valpha.objects,
-# Valpha.inverse.Palpha, Palpha,
+ Aalpha, Palpha,
rweights, XX, TT, names.yy,
nugget, eta,
expectations,
cov.bhat, full.cov.bhat,
cov.betahat,
cov.bhat.betahat,
- cov.delta.bhat, full.cov.delta.bhat,
- cov.delta.bhat.betahat,
+ cov.deltabhat, full.cov.deltabhat,
+ cov.deltabhat.betahat,
cov.ehat, full.cov.ehat,
cov.ehat.p.bhat, full.cov.ehat.p.bhat,
aux.cov.pred.target,
@@ -48,486 +48,698 @@
## 2012-11-04 AP unscaled psi-function
## 2013-02-05 AP covariance matrix of xihat
## 2013-04-23 AP new names for robustness weights
+ ## 2013-05-06 AP changes for solving estimating equations for xi
- n <- nrow( XX )
- sel <- 1:n
- result <- list( error = FALSE )
- a <- expectations["psi2"]
- b <- expectations["dpsi"]
+ ## adjust flags for computing covariance matrices
- TtWiT <- b
- TtDT <- a
- TtT <- as.vector( table( TT ) )
+ cov.bhat.b <- FALSE
+ cov.bhat.e <- FALSE
+ cov.betahat.b <- FALSE
+ cov.betahat.e <- FALSE
- ## ... aggregate elements of Wi and D for replicated observations
-
- if( sum( duplicated( TT ) ) > 0 ){
- TtWiT <- TtWiT * TtT
- TtDT <- TtDT * TtT
+ if( any( c( cov.deltabhat, aux.cov.pred.target ))) cov.bhat.b <- TRUE
+ if( any( c( cov.ehat, aux.cov.pred.target ))) cov.bhat.e <- TRUE
+ if( any( c( cov.deltabhat.betahat, cov.ehat.p.bhat, aux.cov.pred.target ))) cov.betahat.b <- TRUE
+ if( any( c( cov.ehat, cov.ehat.p.bhat, aux.cov.pred.target ))) cov.betahat.e <- TRUE
+ if( any( c( cov.ehat, cov.ehat.p.bhat ) ) ) cov.betahat <- TRUE
+ if( any( c( cov.deltabhat, cov.deltabhat.betahat ))){
+ cov.bhat <- TRUE
+ if( full.cov.deltabhat ) full.cov.bhat <- TRUE
}
+ if( cov.deltabhat.betahat ) cov.bhat.betahat <- TRUE
+ if( cov.ehat ){
+ cov.deltabhat.betahat <- TRUE
+ cov.deltabhat <- TRUE
+ if( full.cov.ehat ) full.cov.deltabhat <- TRUE
+ }
- ## construct matrix M
+ ## compute required auxiliary items
+
+ result.new <- list( error = FALSE )
- TtWiTXX <- TtWiT * XX
- aux <- Valpha.objects$Valpha.inverse / eta
- diag( aux ) <- diag( aux ) + TtWiT
+ a <- expectations["psi2"]
+ b <- expectations["dpsi"]
- M <- rbind(
- cbind( aux, TtWiTXX ),
- cbind( t(TtWiTXX), crossprod( XX , TtWiTXX ) )
+ TtT <- as.vector( table( TT ) )
+
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/georob -r 6
More information about the Georob-commits
mailing list