[Georob-commits] r8 - in pkg: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 11 21:44:56 CEST 2013
Author: papritz
Date: 2013-06-11 21:44:56 +0200 (Tue, 11 Jun 2013)
New Revision: 8
Modified:
pkg/ChangeLog
pkg/R/georob.exported.functions.R
pkg/R/georob.private.functions.R
Log:
non-robust estimation with rank deficient fixed effects design matrix
M pkg/R/georob.exported.functions.R
M pkg/R/georob.private.functions.R
M pkg/ChangeLog
Modified: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog 2013-06-10 10:49:41 UTC (rev 7)
+++ pkg/ChangeLog 2013-06-11 19:44:56 UTC (rev 8)
@@ -85,8 +85,13 @@
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.exported.function.R (georob, georob.control): handling fixed effects model matrices with rank < ncol(x) for robust estimation
+* georob.private.function.R (prepare.likelihood.calculations, compute.estimating.equations, 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
+* georob.private.function.R (estimate.xihat, prepare.likelihood.calculations, compute.estimating.equations, negative.restr.loglikelihood, gradient.negative.restricted.loglikelihood, georob.fit): handling fixed effects model matrices with rank < ncol(x) for robust estimation
+
+2013-06-11 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.exported.function.R (georob): handling fixed effects model matrices with rank < ncol(x) for non-robust estimation
+* georob.private.function.R (estimate.xihat, prepare.likelihood.calculations, compute.estimating.equations, negative.restr.loglikelihood, gradient.negative.restricted.loglikelihood, georob.fit): handling fixed effects model matrices with rank < ncol(x) for non-robust estimation
Modified: pkg/R/georob.exported.functions.R
===================================================================
--- pkg/R/georob.exported.functions.R 2013-06-10 10:49:41 UTC (rev 7)
+++ pkg/R/georob.exported.functions.R 2013-06-11 19:44:56 UTC (rev 8)
@@ -175,22 +175,22 @@
## check whether design matrix has full column rank
+ rankdef.x <- FALSE
+
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"
- )
+ rankdef.x <- TRUE
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'"
+ signif( condnum, 2 ), ")\ninitial values of fixed effects coefficients are computed by 'lm'\n\n"
)
control$initial.method <- "lm"
+ initial.param <- FALSE
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'"
+ signif( condnum, 2 ), ")\ninitial values of fixed effects coefficients are computed by 'lm'"
)
}
@@ -375,6 +375,7 @@
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,
@@ -440,6 +441,7 @@
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,
Modified: pkg/R/georob.private.functions.R
===================================================================
--- pkg/R/georob.private.functions.R 2013-06-10 10:49:41 UTC (rev 7)
+++ pkg/R/georob.private.functions.R 2013-06-11 19:44:56 UTC (rev 8)
@@ -798,7 +798,7 @@
estimate.xihat <-
function(
- XX, min.condnum, yy, betahat, TT, xihat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, xihat,
psi.function, tuning.psi, tuning.psi.nr,
maxit, reltol,
nugget, eta, Valpha.inverse,
@@ -834,15 +834,23 @@
## compute projection matrix Palpha and related items
-# browser()
-#
result <- list( error = FALSE )
aux <- t( XX ) %*% Valpha.inverse
- s <- svd( aux %*% XX )
- s$d <- ifelse( s$d / max( s$d ) <= min.condnum, 0., 1. / s$d )
- Palpha <- s$v %*% ( s$d * t( s$u ) )
+ if( rankdef.x ){
+ s <- svd( aux %*% XX )
+ s$d <- ifelse( s$d / max( s$d ) <= min.condnum, 0., 1. / s$d )
+ Palpha <- s$v %*% ( s$d * t( s$u ) ) # Moore-Penrose inverse
+ } else {
+ t.chol <- try( chol( aux %*% XX ), silent = TRUE )
+ if( !identical( class( t.chol ), "try-error" ) ){
+ Palpha <- chol2inv( t.chol )
+ } else {
+ result$error <- TRUE
+ return( result )
+ }
+ }
result$Aalpha <- Palpha %*% aux
dimnames( result$Aalpha ) <- dimnames( t(XX) )
@@ -1055,7 +1063,7 @@
adjustable.param, variogram.model, fixed.param, param.name, aniso.name,
param.tf, bwd.tf, safe.param,
lag.vectors,
- XX, min.condnum, yy, betahat, TT, bhat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, bhat,
psi.function, dpsi.function, tuning.psi, tuning.psi.nr,
irwls.initial, irwls.maxiter, irwls.reltol,
compute.Q,
@@ -1282,7 +1290,7 @@
}
lik.item$effects <- estimate.xihat(
- XX, min.condnum, yy, betahat, TT, bhat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, bhat,
psi.function, tuning.psi, tuning.psi.nr,
irwls.maxiter, irwls.reltol,
lik.item$param["nugget"], lik.item$eta, lik.item$Valpha$Valpha.inverse,
@@ -1318,27 +1326,38 @@
Q <- rbind(
cbind( aux, TtDTX ),
cbind( t(TtDTX), crossprod( XX, TtDTX) )
- )
+ ) / lik.item$param["nugget"]
lik.item$Q <- list( error = TRUE )
- ## compute log(det(Q)) and inverse of Q by cholesky decomposition
-
- t.chol <- try(
- chol( Q / lik.item$param["nugget"] ),
- silent = TRUE
- )
-
- if( !identical( class( t.chol ), "try-error" ) ) {
+ if( rankdef.x ){
+ ## compute log(pseudo.det(Q)) and (Moore-Penrose) pseudo inverse of Q by svd
+
lik.item$Q$error <- FALSE
- lik.item$Q$log.det.Q <- 2 * sum( log( diag( t.chol) ) )
- lik.item$Q$Q.inverse <- chol2inv( t.chol )
+ s <- svd( Q )
+ lik.item$Q$log.det.Q <- sum( log( s$d[s$d / max( s$d ) > min.condnum] ) )
+ s$d <- ifelse( s$d / max( s$d ) <= min.condnum, 0., 1. / s$d )
+ lik.item$Q$Q.inverse <- s$v %*% ( s$d * t( s$u ) )
} else {
- return( lik.item )
+ ## compute log(det(Q)) and inverse of Q by cholesky decomposition
+ t.chol <- try( chol( Q ), silent = TRUE )
+
+ if( !identical( class( t.chol ), "try-error" ) ) {
+
+ lik.item$Q$error <- FALSE
+ lik.item$Q$log.det.Q <- 2 * sum( log( diag( t.chol) ) )
+ lik.item$Q$Q.inverse <- chol2inv( t.chol )
+
+ } else {
+
+ return( lik.item )
+
+ }
+
}
}
@@ -2370,7 +2389,7 @@
variogram.model, fixed.param, param.name, aniso.name,
param.tf, bwd.tf, safe.param,
lag.vectors,
- XX, min.condnum, yy, betahat, TT, bhat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, bhat,
psi.function, dpsi.function,
tuning.psi, tuning.psi.nr,
irwls.initial, irwls.maxiter, irwls.reltol,
@@ -2397,7 +2416,7 @@
adjustable.param, variogram.model, fixed.param, param.name, aniso.name,
param.tf, bwd.tf, safe.param,
lag.vectors,
- XX, min.condnum, yy, betahat, TT, bhat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, bhat,
psi.function, dpsi.function, tuning.psi, tuning.psi.nr,
irwls.initial, irwls.maxiter, irwls.reltol,
compute.Q = FALSE,
@@ -2628,7 +2647,7 @@
variogram.model, fixed.param, param.name, aniso.name,
param.tf, bwd.tf, safe.param,
lag.vectors,
- XX, min.condnum, yy, betahat, TT, bhat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, bhat,
psi.function, dpsi.function,
tuning.psi, tuning.psi.nr,
irwls.initial, irwls.maxiter, irwls.reltol,
@@ -2657,7 +2676,7 @@
adjustable.param, variogram.model, fixed.param, param.name, aniso.name,
param.tf, bwd.tf, safe.param,
lag.vectors,
- XX, min.condnum, yy, betahat, TT, bhat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, bhat,
psi.function, dpsi.function, tuning.psi, tuning.psi.nr,
irwls.initial, irwls.maxiter, irwls.reltol,
compute.Q = TRUE,
@@ -2743,7 +2762,7 @@
variogram.model, fixed.param, param.name, aniso.name,
param.tf, deriv.fwd.tf, bwd.tf, safe.param,
lag.vectors,
- XX, min.condnum, yy, betahat, TT, bhat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, bhat,
psi.function, dpsi.function, d2psi.function,
tuning.psi, tuning.psi.nr,
irwls.initial, irwls.maxiter, irwls.reltol,
@@ -2773,7 +2792,7 @@
adjustable.param, variogram.model, fixed.param, param.name, aniso.name,
param.tf, bwd.tf, safe.param,
lag.vectors,
- XX, min.condnum, yy, betahat, TT, bhat,
+ XX, min.condnum, rankdef.x, yy, betahat, TT, bhat,
psi.function, dpsi.function, tuning.psi, tuning.psi.nr,
irwls.initial, irwls.maxiter, irwls.reltol,
compute.Q = TRUE,
@@ -3188,7 +3207,7 @@
cov.ehat, full.cov.ehat,
cov.ehat.p.bhat, full.cov.ehat.p.bhat,
aux.cov.pred.target,
- min.condnum,
+ min.condnum, rankdef.x,
psi.func,
tuning.psi.nr,
irwls.initial,
@@ -3780,7 +3799,7 @@
bwd.tf = bwd.tf,
safe.param = safe.param,
lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum,
+ 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,
@@ -3831,7 +3850,7 @@
bwd.tf = bwd.tf,
safe.param = safe.param,
lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum,
+ 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,
@@ -3885,7 +3904,7 @@
bwd.tf = bwd.tf,
safe.param = safe.param,
lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum,
+ 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,
@@ -3922,7 +3941,7 @@
bwd.tf = bwd.tf,
safe.param = safe.param,
lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum,
+ 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,
@@ -3960,7 +3979,7 @@
bwd.tf = bwd.tf,
safe.param = safe.param,
lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum,
+ 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,
@@ -3990,7 +4009,7 @@
bwd.tf = bwd.tf,
safe.param = safe.param,
lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum,
+ 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,
More information about the Georob-commits
mailing list