[Rcpp-commits] r1233 - in pkg/RcppGSL: . R man src

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Thu May 13 21:35:18 CEST 2010


Author: edd
Date: 2010-05-13 21:35:17 +0200 (Thu, 13 May 2010)
New Revision: 1233

Modified:
   pkg/RcppGSL/NAMESPACE
   pkg/RcppGSL/R/fastLm.R
   pkg/RcppGSL/man/fastLm.Rd
   pkg/RcppGSL/src/fastLm.cpp
Log:
fastLm is now generic and returns a fastLm object


Modified: pkg/RcppGSL/NAMESPACE
===================================================================
--- pkg/RcppGSL/NAMESPACE	2010-05-13 15:26:07 UTC (rev 1232)
+++ pkg/RcppGSL/NAMESPACE	2010-05-13 19:35:17 UTC (rev 1233)
@@ -1,4 +1,10 @@
 useDynLib(RcppGSL)
 
-#exportPattern("^[[:alpha:]]+")
-export(fastLm)
+##exportPattern("^[[:alpha:]]+")
+export(fastLmCall,
+       fastLm,
+       fastLm.default,
+       fastLm.formula,
+       print.fastLm,
+       summary.fastLm,
+       print.summary.fastLm)

Modified: pkg/RcppGSL/R/fastLm.R
===================================================================
--- pkg/RcppGSL/R/fastLm.R	2010-05-13 15:26:07 UTC (rev 1232)
+++ pkg/RcppGSL/R/fastLm.R	2010-05-13 19:35:17 UTC (rev 1233)
@@ -17,7 +17,7 @@
 ## You should have received a copy of the GNU General Public License
 ## along with RcppGSL.  If not, see <http://www.gnu.org/licenses/>.
 
-fastLm <- function(y, X) {
+fastLmCall <- function(y, X) {
 
     stopifnot(is.matrix(X))
     stopifnot(nrow(y)==nrow(X))
@@ -25,18 +25,77 @@
     res <- .Call("fastLm", y, X, package="RcppGSL")
 }
 
-## What would be nice here:
-##
-##   fastLm <- function(x, ...) UseMethod("fastLm")
-##
-##   fastLm.formula <- ...
-##
-##   fastLm.default <- ...
-##
-##   print.fastLm <- function(x, ...)
-##
-##   summary.fastLm <- function(object, ...
-##
-##   print.summary.fastLm <- ...
-##
-## but on the other hand lm.fit() does not do any of this either
+fastLm <- function(x, ...) UseMethod("fastLm")
+
+fastLm.default <- function(x, y, ...) {
+
+    y <- as.numeric(y)
+    X <- as.matrix(x)
+
+    res <- fastLmCall(y, X)
+
+    res$fitted.values <- as.vector(x %*% res$coefficients)
+    res$residuals <- y - res$fitted.values
+    res$call <- match.call()
+
+    class(res) <- "fastLm"
+    res
+}
+
+print.fastLm <- function(x, ...) {
+    cat("Call:\n")
+    print(x$call)
+    cat("\nCoefficients:\n")
+    print(x$coefficients)
+}
+
+summary.fastLm <- function(object, ...) {
+    se <- object$stderr
+    tval <- coef(object)/se
+
+    TAB <- cbind(Estimate = coef(object),
+                 StdErr = se,
+                 t.value = tval,
+                 p.value = 2*pt(-abs(tval), df=object$df))
+
+    res <- list(call=object$call,
+                coefficients=TAB)
+
+    class(res) <- "summary.fastLm"
+    res
+}
+
+print.summary.fastLm <- function(x, ...) {
+    cat("Call:\n")
+    print(x$call)
+    cat("\n")
+
+    printCoefmat(x$coefficients, P.value=TRUE, has.Pvalue=TRUE)
+}
+
+fastLm.formula <- function(formula, data=list(), ...) {
+    mf <- model.frame(formula=formula, data=data)
+    x <- model.matrix(attr(mf, "terms"), data=mf)
+    y <- model.response(mf)
+
+    res <- fastLm.default(x, y, ...)
+    names(res$coefficients) <- colnames(x)
+    colnames(res$vcov) <- rownames(res$vcov) <- colnames(x)
+    res$call <- match.call()
+    res$formula <- formula
+    res
+}
+
+predict.fastLm <- function(object, newdata=NULL, ...) {
+    if (is.null(newdata)) {
+        y <- fitted(object)
+    } else {
+        if (!is.null(object$formula)) {
+            x <- model.matrix(object$formula, newdata)
+        } else {
+            x <- newdata
+        }
+        y <- as.vector(x %*% coef(object))
+    }
+    y
+}

Modified: pkg/RcppGSL/man/fastLm.Rd
===================================================================
--- pkg/RcppGSL/man/fastLm.Rd	2010-05-13 15:26:07 UTC (rev 1232)
+++ pkg/RcppGSL/man/fastLm.Rd	2010-05-13 19:35:17 UTC (rev 1233)
@@ -1,5 +1,12 @@
 \name{fastLm}
 \alias{fastLm}
+\alias{fastLmCall}
+\alias{fastLm.default}
+\alias{fastLm.formula}
+\alias{print.fastLm}
+\alias{summary.fastLm}
+\alias{print.summary.fastLm}
+\alias{predict.fastLm}
 \concept{regression}
 \title{Bare-bones linear model fitting function}
 \description{
@@ -7,11 +14,25 @@
   function of the \code{GNU GSL} library.
 }
 \usage{
-fastLm(y, X)
+fastLmCall(y, X)
+
+fastLm(x, ...)
+\method{fastLm}{default}(x, y, ...)
+\method{fastLm}{formula}(formula, data = list(), ...)
+\method{print}{fastLm}(x, ...)
+\method{summary}{fastLm}(object, ...)
+\method{predict}{fastLm}(object, newdata=NULL, ...)
 }
 \arguments{
   \item{y}{a vector containing the explained variable.}
   \item{X}{a data.frame or matrix containing the explanatory variables.}
+  \item{x}{a numeric design matrix for the model.}
+  \item{formula}{a symbolic description of the model to be fit.}
+  \item{data}{an optional data frame containing the variables in the model.}
+  \item{newdata}{an optional data frame containing the variables for
+    which a model prediction is desired.}
+  \item{object}{an object of class \code{"fastLm"}, i.e. a fitted model}
+  \item{\ldots}{not used}
 }
 \details{
   Linear models should be estimated using the \code{\link{lm}} function. In
@@ -58,8 +79,17 @@
   Gerard Jungman.  RcppGSL is written by Romain Francois and Dirk Eddelbuettel.
 }
 \examples{
+  require(datasets)
   data(trees)
-  flm <- fastLm( log(trees$Volume),  cbind(rep(1,31), log(trees$Girth)) )
+
+  ## bare-bones direct interface
+  flm <- fastLmCall( log(trees$Volume),  cbind(1, log(trees$Girth)) )
+  print(flm)
+
+  ## standard R interface for formula or data returning object of class fastLm
+  flmmod <- fastLm( log(Volume) ~ log(Girth), data=trees)
+  summary(flmmod)
+
 }
 \keyword{regression}
 

Modified: pkg/RcppGSL/src/fastLm.cpp
===================================================================
--- pkg/RcppGSL/src/fastLm.cpp	2010-05-13 15:26:07 UTC (rev 1232)
+++ pkg/RcppGSL/src/fastLm.cpp	2010-05-13 19:35:17 UTC (rev 1233)
@@ -50,13 +50,15 @@
     std::transform( stderr.begin(), stderr.end(), stderr.begin(), sqrt );
     
     Rcpp::List res = Rcpp::List::create( 
-    	_["coef"] = c, 
-    	_["stderr"] = stderr
+    	_["coefficients"] = coef, 
+    	_["stderr"] = stderr,
+		_["vcov"] = cov,
+		_["df"] = n - k
     	) ;
     
     // free all the GSL vectors and matrices -- as these are really C data structure
 	// we cannot take advantage of automatic memory management
-    c.free() ;
+    coef.free() ;
     cov.free();
     y.free();
     X.free();



More information about the Rcpp-commits mailing list