version 0.0-7

DESCRIPTION

@@ -0,0 +1,29 @@
+Package: robustsur
+Version: 0.0-7
+Date: 2021-10-03
+Title: Robust Estimation for Seemingly Unrelated Regression Models
+Authors@R: 
+  c(person(given = "Claudio",
+           family = "Agostinelli",
+           role = c("aut", "cre"),
+           email = "claudio.agostinelli@unitn.it",
+           comment = c(ORCID = "0000-0001-6702-4312")),
+    person(given = "Giovanni",
+           family = "Saraceno",
+           role = "aut",
+           email = "giovanni.saraceno@unitn.it",
+           comment = c(ORCID = "0000-0002-1753-2367"))
+  )
+Maintainer: Claudio Agostinelli <claudio.agostinelli@unitn.it>
+Description: Data sets are often corrupted by outliers. When data are multivariate outliers can be classified as case-wise or cell-wise. The latters are particularly challenge to handle. We implement a robust estimation procedure for Seemingly Unrelated Regression Models which is able to cope well with both type of outliers. Giovanni Saraceno, Fatemah Alqallaf, Claudio Agostinelli (2021) <arXiv:2107.00975>. 
+Depends: R (>= 3.0.0), robustbase, robreg3S
+Imports: Matrix, GSE
+Suggests: systemfit
+License: GPL (>= 2)
+NeedsCompilation: no
+Packaged: 2021-10-03 14:59:00 UTC; claudio
+Author: Claudio Agostinelli [aut, cre]
+    (<https://orcid.org/0000-0001-6702-4312>),
+  Giovanni Saraceno [aut] (<https://orcid.org/0000-0002-1753-2367>)
+Repository: CRAN
+Date/Publication: 2021-10-04 08:40:02 UTC

MD5

@@ -0,0 +1,11 @@
+c889a622c3e00d1d1e5fb459dfab7fde *DESCRIPTION
+242c16761b156bdea2f86939161455b3 *NAMESPACE
+12e16d40a36dbe7b430ed783c0e0949a *R/eigenkronecker.R
+7f00752e499c6747bf8369dbd2ae7a3a *R/print.surerob.R
+524665f80b0380baaff0e2af9edaeb6c *R/summary.R
+ddeb3e482470f1b473f3df6d408d83a0 *R/surerob.R
+f97ae5f1d44d19bdcccd242d6facc10c *R/surerob2.R
+db01133b16a7576ed72b7442b2df8114 *README
+41301004534268b992a6832f1609a823 *man/eigenkronecker.Rd
+5d705370a273d93e78c39853efd0e601 *man/summary.surerob.Rd
+c90a16427bb7adf5ef7a85febec14fc4 *man/surerob.Rd

NAMESPACE

@@ -0,0 +1,24 @@
+##useDynLib(robustsur, .registration=TRUE) # -> R >= 2.3.0
+
+importFrom("robustbase", "lmrob")
+importFrom("robustbase", "lmrob.control")
+importFrom("robustbase", "Mwgt")
+importFrom("robustbase", ".Mpsi.tuning.default")  
+importFrom("robreg3S", "robreg3S")      
+importFrom("Matrix", "bdiag")
+importFrom("GSE", "TSGS")
+importFrom("stats", "coef", "cor", "formula", "lm", "lm.fit",
+             "model.matrix", "model.response", "pnorm", "printCoefmat",
+             "residuals")
+##importFrom("systemfit", "Kmenta")  
+
+export(
+  eigenkronecker,
+  surerob
+)
+
+
+S3method(summary, surerob)  
+S3method(print, surerob)
+S3method(print, summary.surerob)
+

R/eigenkronecker.R

@@ -0,0 +1,8 @@
+eigenkronecker <- function(x, n) {
+  e <- eigen(x, symmetric=TRUE, only.values=FALSE)
+  d <- e$values
+  v <- e$vectors
+  D <- rep(d, each=n)
+  V <- kronecker(v, diag(n))
+  list(values=D, vectors=V)  
+}

R/print.surerob.R

@@ -0,0 +1,42 @@
+## print a few results of the whole system
+print.surerob <- function( x,
+      digits = max( 3, getOption("digits") - 1 ),... ) {
+
+   cat("\n")
+   cat("surerob results \n")
+   cat("method: ")
+   if(!is.null(x$iter)) if(x$iter>1) cat("iterated ")
+   cat( paste( x$method, "\n\n"))
+   if(!is.null(x$iter)) {
+      if(x$iter>1) {
+         if(x$iter<x$control$maxiter) {
+            cat( paste( "convergence achieved after",x$iter,"iterations\n\n" ) )
+         } else {
+            cat( paste( "warning: convergence not achieved after", x$iter,
+                        "iterations\n\n" ) )
+         }
+      }
+   }
+   cat( "Coefficients:\n" )
+   print( x$coefficients, digits = digits )
+   invisible( x )
+}
+
+## ## print a few results for a single equation
+## print.surerob.equation <- function( x,
+##       digits = max( 3, getOption("digits") - 1 ), ... ) {
+##    cat("\n")
+##    cat( x$method, " estimates for '", x$eqnLabel,
+##             "' (equation ", x$eqnNo, ")\n", sep = "" )
+
+##    cat("Model Formula: ")
+##    print( formula( x$terms ) )
+##    if(!is.null(x$inst)) {
+##       cat("Instruments: ")
+##       print(x$inst)
+##    }
+
+##    cat("\nCoefficients:")
+##    print( x$coefficients, digits = digits )
+##    invisible( x )
+## }

R/summary.R

@@ -0,0 +1,173 @@
+## prepare summary results that belong to the whole system
+summary.surerob <- function(object, residCov=TRUE, equations=TRUE, ...) {
+
+   # if( is.null( useDfSys ) ) {
+   #    useDfSys <- length( coef( object ) ) != object$rank
+   #       # TRUE if there are restrictions imposed
+   # }
+
+   # number of equations
+   nEq <- length( object$eq )
+   # number of observations per equation
+   nObsEq <- rep( NA, nEq )
+   for( i in 1:nEq ) {
+      nObsEq[ i ] <- length( residuals( object$eq[[ i ]], na.rm = TRUE ) )
+   }
+   # total number of observations
+   nObs <- sum( nObsEq )
+
+   # preparing objects that will be returned
+   result <- list()
+   result$method <- object$method
+   result$residuals <- residuals( object )
+   result$residCovEst <- object$residCovEst
+   result$residCov <- object$residCov
+   if( !is.null( result$residCovEst ) ) {
+      dimnames( result$residCovEst ) <- dimnames( result$residCov )
+   }
+   result$residCor <- cor( residuals( object ), use = "na.or.complete" )
+   dimnames( result$residCor ) <- dimnames( result$residCov )
+   result$detResidCov <- det( object$residCov, tol = object$control$solvetol )
+   result$rweights <- object$rweights
+
+   # now prepare summury results for the individual equations
+   result$eq <- list()
+   for( i in 1:length( object$eq ) ) {
+       result$eq[[i]] <- summary(object$eq[[i]])
+       result$eq[[i]]$ssr <- sum( residuals( object$eq[[i]], na.rm = TRUE )^2 )
+       result$eq[[i]]$eqnNo <- object$eq[[i]]$eqnNo
+       result$eq[[i]]$eqnLabel <- object$eq[[i]]$eqnLabel
+   }
+
+   # coefficients, standard errors, ... 
+   p <- object$rank
+   df <- object$df.residual
+   cf.names <- c( "Estimate", "Std. Error", "t value", "Pr(>|t|)" )
+   if(p > 0){
+      n <- p + df
+      se <- sqrt(diag(object$coefCov))
+      est <- object$coefficients
+      tval <- est/se
+      if(!is.null(result$rweights)){
+         result$residual_weighted <- result$residuals * sqrt(object$rweights)
+      }
+      result$df <- c(p, df)
+      result$coefficients <- cbind(est, se, tval, 2*pnorm(abs(tval), df, lower.tail = FALSE))
+      dimnames(result$coefficients) <- list(names(est), cf.names)
+
+      resid <- result$residuals
+      pred <- object$fitted.values
+      resp <- object$y
+      wgt <- object$rweights
+      psi <- object$control$psi
+      correc <- switch(psi,
+                       bisquare = 1.207617,
+                       welsh    = 1.224617, # 1.2246131
+                       optimal  = 1.068939,
+                       hampel   = 1.166891,
+                       lqq      = 1.159232,
+                       stop('unsupported psi function -- should not happen'))
+      resp.mean <-  sum(wgt * resp)/sum(wgt) 
+      yMy <- sum(wgt * (resp - resp.mean)^2)
+      rMr <- sum(wgt * resid^2)
+      result$ssr_weighted <- rMr
+
+      result$r.squared <- r2correc <- (yMy - rMr) / (yMy + rMr * (correc - 1))
+      result$adj.r.squared <- 1 - (1 - r2correc) * ((n - 1) / df)
+      result$coefCov <- object$coefCov
+   } 
+   # else { ## p = 0: "null model"
+   #    result <- object
+   #    result$df <- c(0L, df, length(is.na(coef(object))))
+   #    result$coefficients <- matrix(result$coefficients[0L], 0L, 4L, dimnames = list(NULL, cf.names))
+   #    result$r.squared <- result$adj.r.squared <- 0
+   #    result$coefCov <- object$coefCocov
+   # }
+
+   result$printResidCov  <- residCov
+   result$printEquations <- equations
+   result$control <- object$control
+   result$call <- object$call 
+   class( result ) <- "summary.surerob"
+   return( result )
+}
+
+## print summary results of the whole system
+print.summary.surerob <- function(x,
+      digits = max( 3, getOption("digits") - 1 ),
+      residCov = x$printResidCov, equations = x$printEquations, ...) {
+
+  table <- NULL
+  labels <- NULL
+
+  cat("\n")
+  cat("robustsur results \n")
+  cat("method: ")
+  cat( paste( x$method, "\n\n"))
+
+
+   table.sys <- cbind( round( sum( x$df ),         digits ),
+                       round( x$df[2],             digits ),
+                       round( x$ssr_weighted, digits ),
+                       round( x$detResidCov,       digits ),
+                       round( x$r.squared,     digits ),
+                       round( x$adj.r.squared,     digits ))
+   rownames( table.sys ) <- c( "system" )
+   colnames( table.sys ) <- c( "N", "DF", "SSR", "detRCov", "R2", "Adjusted R2")
+   print( table.sys, quote = FALSE, right = TRUE, digits = digits )
+
+   cat("\n")
+
+  for(i in 1:length( x$eq ) ) {
+    row <- NULL
+    row <- cbind( round( sum( x$eq[[i]]$df ),  digits ),
+                  round( x$eq[[i]]$df[2], digits ),
+                  round( x$eq[[i]]$ssr,   digits ),
+                  round( x$eq[[i]]$r.squared,     digits ),
+                  round( x$eq[[i]]$adj.r.squared, digits ))
+    table  <- rbind( table, row )
+    labels <- rbind( labels, x$eq[[i]]$eqnLabel )
+  }
+  rownames(table) <- c( labels )
+  colnames(table) <- c("N","DF", "SSR", "R2", "Adj R2" )
+  print(table, quote = FALSE, right = TRUE, digits = digits )
+
+  cat("\n")
+
+   if( residCov ){
+      if(!is.null(x$residCovEst)) {
+         cat("The covariance matrix of the residuals used for estimation\n")
+         print( x$residCovEst, digits = digits )
+         cat("\n")
+         if( min(eigen( x$residCov, only.values=TRUE)$values) < 0 ) {
+            cat("warning: this covariance matrix is NOT positive semidefinit!\n")
+            cat("\n")
+         }
+      }
+
+      cat("The covariance matrix of the residuals\n")
+      print( x$residCov, digits = digits )
+      cat("\n")
+
+      cat("The correlations of the residuals\n")
+      print( x$residCor, digits = digits )
+      cat("\n")
+   }
+
+   if( equations ){
+      ## now print the individual equations
+      for(i in 1:length( x$eq ) ) {
+         cat("\n")
+         cat( x$method, " estimates for '", x$eq[[i]]$eqnLabel,
+              "' (equation ", x$eq[[i]]$eqnNo, ")\n", sep = "" )
+
+         cat("Model Formula: ")
+         print( formula( x$eq[[i]]$terms ) )
+         print( x$eq[[i]], digits = digits, showAlgo = FALSE )
+      }
+   } else {
+      cat( "\nCoefficients:\n" )
+      printCoefmat( coef( x ), digits = digits )
+   }
+  invisible( x )
+}