[Rcpp-commits] r2641 - in pkg/RcppModels: . R src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 1 22:25:21 CET 2010
Author: dmbates
Date: 2010-12-01 22:25:21 +0100 (Wed, 01 Dec 2010)
New Revision: 2641
Modified:
pkg/RcppModels/DESCRIPTION
pkg/RcppModels/R/fastGlm.R
pkg/RcppModels/src/glmFamily.h
pkg/RcppModels/tests/pois.R
Log:
see ChangeLog
Modified: pkg/RcppModels/DESCRIPTION
===================================================================
--- pkg/RcppModels/DESCRIPTION 2010-12-01 21:24:48 UTC (rev 2640)
+++ pkg/RcppModels/DESCRIPTION 2010-12-01 21:25:21 UTC (rev 2641)
@@ -9,6 +9,6 @@
linear and nonlinear models that use linear predictor expressions.
License: GPL (>= 2)
LazyLoad: yes
-Depends: R(>= 2.12.0), Rcpp(>= 0.8.9), RcppArmadillo, methods,
+Depends: R(>= 2.12.0), Rcpp(>= 0.8.9), RcppArmadillo, methods
LinkingTo: RcppArmadillo, Rcpp
Suggests: RUnit
Modified: pkg/RcppModels/R/fastGlm.R
===================================================================
--- pkg/RcppModels/R/fastGlm.R 2010-12-01 21:24:48 UTC (rev 2640)
+++ pkg/RcppModels/R/fastGlm.R 2010-12-01 21:25:21 UTC (rev 2641)
@@ -1,3 +1,12 @@
+##' Create the response module from the model frame
+##'
+##' Create a glmResp object from the model frame and the glm family.
+##' The initialization of the family as applied to the response takes
+##' place here.
+##' @title Create a glmResp object
+##' @param fr a model frame
+##' @param family a glm family object
+##' @return a glmResp object
mkRespMod <- function(fr, family)
{
N <- n <- nrow(fr)
@@ -31,21 +40,16 @@
function(formula, family, data, weights, subset,
na.action, start = NULL, etastart, mustart, offset,
drop.unused.levels = FALSE, doFit = TRUE,
- control = list(...))
+ control = list(...), ...)
{
call <- match.call()
- if (missing(family)) {
- family <- NULL
- } else {
- if(is.character(family))
- family <- get(family, mode = "function", envir = parent.frame())
- if(is.function(family)) family <- family()
- if(is.null(family$family)) {
- print(family)
- stop("'family' not recognized")
- }
- }
- ## extract x, y, etc from the model formula and frame
+ maxit <- if(!is.null(control$maxit)) control$maxit else 8
+ # check the family argument
+ if(is.character(family))
+ family <- get(family, mode = "function", envir = parent.frame())
+ if(is.function(family)) family <- family()
+ if(is.null(family$family)) stop("'family' not recognized")
+ # create the model frame
if(missing(data)) data <- environment(formula)
mf <- match.call(expand.dots = FALSE)
m <- match(c("formula", "data", "subset", "weights", "na.action",
@@ -54,14 +58,22 @@
mf$drop.unused.levels <- TRUE
mf[[1L]] <- as.name("model.frame")
mf <- eval(mf, parent.frame())
+ # create the response and predictor modules
+ rr <- RcppModels:::mkRespMod(mf, family)
+ pp <- new(densePred, model.matrix(formula, mf))
- rr <- mkRespMod(mf, family)
- pp <- new(densePred, model.matrix(formula, mf))
## one iteration of the fixed-point algorithm to establish
## a baseline coefficient vector
wts <- rr$sqrtWrkWt
rr$updateMu(pp$gamma(wts, wts * rr$wrkResp))
- rr$updateWts()
pp$installCoef0()
+ print(pp$coef0)
+ for (i in seq_len(maxit)) {
+ oldwrss <- rr$updateWts()
+ if ((newwrss <- rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres)))
+ < oldwrss) pp$installCoef0()
+ print(pp$coef0)
+ if ((oldwrss - newwrss)/newwrss < 1.e-8) break
+ }
list(rmod = rr, pmod = pp)
}
Modified: pkg/RcppModels/src/glmFamily.h
===================================================================
--- pkg/RcppModels/src/glmFamily.h 2010-12-01 21:24:48 UTC (rev 2640)
+++ pkg/RcppModels/src/glmFamily.h 2010-12-01 21:25:21 UTC (rev 2641)
@@ -110,9 +110,12 @@
BinomialDevRes(double y, double mu, double wt) {
return 2 * wt * (y_log_y(y, mu) + y_log_y(1 - y, 1 - mu));
}
+ static inline double logYMu(double y, double mu) {
+ return y ? log(y/mu) : 0;
+ }
static inline double
GammaDevRes (double y, double mu, double wt) {
- return -2 * wt * y_log_y(y, mu) - (y - mu)/mu;
+ return -2 * wt * (logYMu(y, mu) - (y - mu)/mu);
}
static inline double
GaussianDevRes(double y, double mu, double wt) {
Modified: pkg/RcppModels/tests/pois.R
===================================================================
--- pkg/RcppModels/tests/pois.R 2010-12-01 21:24:48 UTC (rev 2640)
+++ pkg/RcppModels/tests/pois.R 2010-12-01 21:25:21 UTC (rev 2641)
@@ -4,19 +4,22 @@
d.AD <- data.frame(treatment = gl(3,3), outcome = gl(3,1,9),
counts = c(18,17,15,20,10,20,25,13,12))
res <- fastGlm(counts ~ outcome + treatment, poisson, d.AD)
-pp <- res$pmod
-rr <- res$rmod
-rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres))
-(oldwrss <- rr$updateWts())
-pp$installCoef0()
-rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres))
-(oldwrss <- rr$updateWts())
-pp$installCoef0()
-(newwrss <- rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres)))
-(oldwrss <- rr$updateWts())
-pp$installCoef0()
-(newwrss <- rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres)))
-(oldwrss <- rr$updateWts())
-pp$installCoef0()
-(newwrss <- rr$updateMu(pp$gamma(rr$sqrtXwt, rr$wtres)))
+res$pmod$coef0 # estimated coefficients
+res$rmod$wrkResids # working residuals
+res$rmod$wtres # Pearson residuals
+res$rmod$residDeviance # residual deviance
+# A Gamma example, from McCullagh & Nelder (1989, pp. 300-2)
+clotting <-
+ data.frame(u = c(5,10,15,20,30,40,60,80,100),
+ lot1 = c(118,58,42,35,27,25,21,19,18),
+ lot2 = c(69,35,26,21,18,16,13,12,12))
+res1 <- fastGlm(lot1 ~ log(u), Gamma, clotting)
+res1$rmod$wrkResids # working residuals
+res1$rmod$wtres # Pearson residuals
+res1$rmod$residDeviance # residual deviance
+
+res2 <- fastGlm(lot2 ~ log(u), Gamma, clotting)
+res2$rmod$wrkResids # working residuals
+res2$rmod$wtres # Pearson residuals
+res2$rmod$residDeviance # residual deviance
More information about the Rcpp-commits
mailing list