[Lme4-commits] r1623 - in pkg/lme4Eigen: . R inst/tests src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 27 21:09:49 CET 2012
Author: dmbates
Date: 2012-02-27 21:09:49 +0100 (Mon, 27 Feb 2012)
New Revision: 1623
Modified:
pkg/lme4Eigen/DESCRIPTION
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/R/lmer.R
pkg/lme4Eigen/inst/tests/test-glmFamily.R
pkg/lme4Eigen/inst/tests/test-lmer.R
pkg/lme4Eigen/inst/tests/test-lmerResp.R
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/glmFamily.cpp
pkg/lme4Eigen/src/glmFamily.h
pkg/lme4Eigen/src/respModule.cpp
pkg/lme4Eigen/src/respModule.h
Log:
Major overhaul of glmFamily code, New version, added to testthat tests
Modified: pkg/lme4Eigen/DESCRIPTION
===================================================================
--- pkg/lme4Eigen/DESCRIPTION 2012-02-27 11:14:49 UTC (rev 1622)
+++ pkg/lme4Eigen/DESCRIPTION 2012-02-27 20:09:49 UTC (rev 1623)
@@ -1,5 +1,5 @@
Package: lme4Eigen
-Version: 0.9996875-10
+Version: 0.9996875-11
Date: $Date$
Revision: $Revision$
Title: Linear mixed-effects models using Eigen and S4
@@ -31,7 +31,6 @@
MEMSS,
testthat,
ggplot2,
- glmmADMB,
mlmRev,
plyr,
reshape,
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2012-02-27 11:14:49 UTC (rev 1622)
+++ pkg/lme4Eigen/R/AllClass.R 2012-02-27 20:09:49 UTC (rev 1623)
@@ -528,7 +528,7 @@
.Call(glm_wrkResids, ptr())
},
wrkResp = function() {
- 'returns the vector of working residuals'
+ 'returns the vector of working responses'
.Call(glm_wrkResp, ptr())
}
)
@@ -614,15 +614,15 @@
length(dev <- as.numeric(dev)) == 1L)
.Call(glmFamily_aic, ptr(), y, n, mu, wt, dev)
},
- devResid = function(mu, weights, y) {
- 'applies the devResid function to mu, weights and y'
+ devResid = function(y, mu, wt) {
+ 'applies the devResid function to y, mu and wt'
mu <- as.numeric(mu)
- weights <- as.numeric(weights)
- y <- as.numeric(y)
- stopifnot(length(mu) == length(weights),
+ wt <- as.numeric(wt)
+ y <- as.numeric(y)
+ stopifnot(length(mu) == length(wt),
length(mu) == length(y),
- all(weights >= 0))
- .Call(glmFamily_devResid, ptr(), mu, weights, y)
+ all(wt >= 0))
+ .Call(glmFamily_devResid, ptr(), y, mu, wt)
},
link = function(mu) {
'applies the (forward) link function to mu'
Modified: pkg/lme4Eigen/R/lmer.R
===================================================================
--- pkg/lme4Eigen/R/lmer.R 2012-02-27 11:14:49 UTC (rev 1622)
+++ pkg/lme4Eigen/R/lmer.R 2012-02-27 20:09:49 UTC (rev 1623)
@@ -1655,22 +1655,27 @@
"flist" = object at flist,
"beta" = structure(object at beta, names = dimnames(PR$X)[[2]]),
"theta"= {
- tt <- object at theta
- nc <- c(unlist(mapply(function(g,e) {
- mm <- outer(e,e,paste,sep=".")
- diag(mm) <- e
- mm <- mm[lower.tri(mm,diag=TRUE)]
- paste(g,mm,sep=".")
- },
- names(object at cnms),object at cnms)))
- names(tt) <- nc
- tt
+ tt <- object at theta
+ nc <- c(unlist(mapply(function(g,e) {
+ mm <- outer(e,e,paste,sep=".")
+ diag(mm) <- e
+ mm <- mm[lower.tri(mm,diag=TRUE)]
+ paste(g,mm,sep=".")
+ },
+ names(object at cnms),object at cnms)))
+ names(tt) <- nc
+ tt
}, ## *OR* PR $ theta --- which one ?? Should be the same.
"REML" = dims["REML"],
"is_REML" = isREML(object),
"n_rtrms" = length(object at flist), ## should this be length(object at cnms) instead?
+
+ ## Yes, assuming that you want the number of random-effects
+ ## terms in the formula. Multiple terms with the same
+ ## grouping factors are allowed.
+
"devcomp" = dc,
"offset" = rsp$offset,
"..foo.." =# placeholder!
@@ -1680,6 +1685,7 @@
stop(sprintf("Mixed-Effects extraction of '%s' is not available for class \"%s\"",
name, class(object))))
}## {getME}
+
##' @importMethodsFrom Matrix t %*% crossprod diag tcrossprod
##' @importClassesFrom Matrix dgCMatrix dpoMatrix corMatrix
NULL
Modified: pkg/lme4Eigen/inst/tests/test-glmFamily.R
===================================================================
--- pkg/lme4Eigen/inst/tests/test-glmFamily.R 2012-02-27 11:14:49 UTC (rev 1622)
+++ pkg/lme4Eigen/inst/tests/test-glmFamily.R 2012-02-27 20:09:49 UTC (rev 1623)
@@ -1,5 +1,4 @@
library("testthat")
-context("glmFamily")
eps <- .Machine$double.eps
oneMeps <- 1 - eps
set.seed(1)
@@ -22,6 +21,7 @@
pmin(pmax(eps, rbeta(100, 0.1, 3)), oneMeps),
pmin(pmax(eps, rbeta(100, 3, 0.1)), oneMeps))
+context("glmFamily linkInv and muEta")
test_that("inverse link and muEta functions", {
tst.lnki <- function(fam, frm) {
ff <- glmFamily$new(family=fam)
@@ -51,6 +51,7 @@
tst.muEta(inverse.gaussian(), etapos)
})
+context("glmFamily linkFun and variance")
test_that("link and variance functions", {
tst.link <- function(fam, frm) {
ff <- glmFamily$new(family=fam)
@@ -62,19 +63,43 @@
sapply(frm, function(x) expect_that(fam$variance(x), equals(ff$variance(x))))
}
- tst.link( binomial(), mubinom)
- tst.variance(binomial(), mubinom)
- tst.link( binomial("probit"), mubinom)
- tst.variance(binomial("probit"), mubinom)
- tst.link( binomial("cauchit"), mubinom)
- tst.variance(binomial("cauchit"), mubinom)
- tst.link( gaussian(), etas)
- tst.variance(gaussian(), etas)
- tst.link( Gamma(), etapos)
- tst.variance(Gamma(), etapos)
- tst.link( inverse.gaussian(), etapos)
- tst.variance(inverse.gaussian(), etapos)
+ tst.link( binomial(), mubinom)
+ tst.variance(binomial(), mubinom)
+ tst.link( binomial("probit"), mubinom)
+ tst.link( binomial("cauchit"), mubinom)
+ tst.link( gaussian(), etas)
+ tst.variance(gaussian(), etas)
+ tst.link( Gamma(), etapos)
+ tst.variance(Gamma(), etapos)
+ tst.link( inverse.gaussian(), etapos)
+ tst.variance(inverse.gaussian(), etapos)
+ tst.variance(MASS::negative.binomial(1), etapos)
+ tst.variance(MASS::negative.binomial(0.5), etapos)
+ tst.link( poisson(), etapos)
+ tst.variance(poisson(), etapos)
})
+context("glmFamily devResid and aic")
+test_that("devResid and aic", {
+ tst.devres <- function(fam, frm) {
+ ff <- glmFamily$new(family=fam)
+ sapply(frm, function(x) {
+ nn <- length(x)
+ wt <- rep.int(1, nn)
+ n <- wt
+ y <- switch(fam$family,
+ binomial = rbinom(nn, 1L, x),
+ gaussian = rnorm(nn, x),
+ poisson = rpois(nn, x),
+ error("Unknown family"))
+ dev <- ff$devResid(y, x, wt)
+ expect_that(fam$dev.resids(y, x, wt), equals(dev))
+ dd <- sum(dev)
+ expect_that(fam$aic(y, n, x, wt, dd), equals(ff$aic(y, n, x, wt, dd)))
+ })
+ }
-
+ tst.devres(binomial(), mubinom)
+ tst.devres(gaussian(), etas)
+ tst.devres(poisson(), etapos)
+})
Modified: pkg/lme4Eigen/inst/tests/test-lmer.R
===================================================================
--- pkg/lme4Eigen/inst/tests/test-lmer.R 2012-02-27 11:14:49 UTC (rev 1622)
+++ pkg/lme4Eigen/inst/tests/test-lmer.R 2012-02-27 20:09:49 UTC (rev 1623)
@@ -1,6 +1,6 @@
library("testthat")
-context("fitting mixed-effects models")
+context("fitting lmer models")
test_that("lmer", {
expect_that(fm1 <- lmer(Yield ~ 1|Batch, Dyestuff), is_a("lmerMod"))
expect_that(fm1 at resp, is_a("lmerResp"))
@@ -24,14 +24,14 @@
## FIXME: recent version gets 313.09721874266512032 instead?
expect_that(fm2 <- refit(fm1, Dyestuff2$Yield), is_a("lmerMod"))
expect_that(fixef(fm2), is_equivalent_to(5.6656))
- expect_that(VarCorr(fm2)[[1]][1,1], equals(0))
- expect_that(getME(fm2, "theta"), equals(0))
+ expect_that(VarCorr(fm2)[[1]][1,1], is_equivalent_to(0))
+ expect_that(getME(fm2, "theta"), is_equivalent_to(0))
expect_that(X <- getME(fm1, "X"), is_equivalent_to(array(1, c(1, 30))))
expect_that(Zt <- getME(fm1, "Zt"), is_a("dgCMatrix"))
expect_that(dim(Zt), equals(c(6L, 30L)))
expect_that(length(Zt at x), equals(30L))
expect_that(Zt at x, equals(rep.int(1, 30L)))
- expect_that(theta <- getME(fm1, "theta"), equals(0.848330078125))
+ expect_that(theta <- getME(fm1, "theta"), is_equivalent_to(0.848330078125))
expect_that(Lambdat <- getME(fm1, "Lambdat"), is_a("dgCMatrix"))
expect_that(as(Lambdat, "matrix"), is_equivalent_to(diag(theta, 6L, 6L)))
})
Modified: pkg/lme4Eigen/inst/tests/test-lmerResp.R
===================================================================
--- pkg/lme4Eigen/inst/tests/test-lmerResp.R 2012-02-27 11:14:49 UTC (rev 1622)
+++ pkg/lme4Eigen/inst/tests/test-lmerResp.R 2012-02-27 20:09:49 UTC (rev 1623)
@@ -1,5 +1,4 @@
library("testthat")
-context("Response modules")
n <- nrow(Dyestuff)
ones <- rep.int(1, n)
@@ -7,6 +6,7 @@
YY <- Dyestuff$Yield
mYY <- mean(YY)
+context("lmerResp objects")
test_that("lmerResp", {
mres <- YY - mYY
rr <- lmerResp$new(y=YY)
@@ -28,6 +28,7 @@
mlYY <- mean(log(YY))
gmeanYY <- exp(mlYY) # geometric mean
+context("glmResp objects")
test_that("glmResp", {
mres <- YY - gmeanYY
gmean <- rep.int(gmeanYY, n)
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2012-02-27 11:14:49 UTC (rev 1622)
+++ pkg/lme4Eigen/src/external.cpp 2012-02-27 20:09:49 UTC (rev 1623)
@@ -28,7 +28,6 @@
using Rcpp::wrap;
using glm::glmFamily;
- using glm::negativeBinomial;
using lme4Eigen::glmResp;
using lme4Eigen::lmResp;
@@ -177,46 +176,47 @@
SEXP glmFamily_link(SEXP ptr, SEXP mu) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->linkFun(as<MAr1>(mu)));
+ return wrap(XPtr<glmFamily>(ptr)->linkFun(as<MVec>(mu)));
END_RCPP;
}
SEXP glmFamily_linkInv(SEXP ptr, SEXP eta) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->linkInv(as<MAr1>(eta)));
+ return wrap(XPtr<glmFamily>(ptr)->linkInv(as<MVec>(eta)));
END_RCPP;
}
- SEXP glmFamily_devResid(SEXP ptr, SEXP mu, SEXP weights, SEXP y) {
+ SEXP glmFamily_devResid(SEXP ptr, SEXP y, SEXP mu, SEXP wt) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->devResid(as<MAr1>(mu),
- as<MAr1>(weights),
- as<MAr1>(y)));
+ return wrap(XPtr<glmFamily>(ptr)->devResid(as<MVec>(y),
+ as<MVec>(mu),
+ as<MVec>(wt)));
END_RCPP;
}
SEXP glmFamily_aic(SEXP ptr, SEXP y, SEXP n, SEXP mu, SEXP wt, SEXP dev) {
BEGIN_RCPP;
- return ::Rf_ScalarReal(XPtr<glmFamily>(ptr)->aic(as<MAr1>(y),
- as<MAr1>(n),
- as<MAr1>(mu),
- as<MAr1>(wt),
+ return ::Rf_ScalarReal(XPtr<glmFamily>(ptr)->aic(as<MVec>(y),
+ as<MVec>(n),
+ as<MVec>(mu),
+ as<MVec>(wt),
::Rf_asReal(dev)));
END_RCPP;
}
SEXP glmFamily_muEta(SEXP ptr, SEXP eta) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->muEta(as<MAr1>(eta)));
+ return wrap(XPtr<glmFamily>(ptr)->muEta(as<MVec>(eta)));
END_RCPP;
}
SEXP glmFamily_variance(SEXP ptr, SEXP mu) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->variance(as<MAr1>(mu)));
+ return wrap(XPtr<glmFamily>(ptr)->variance(as<MVec>(mu)));
END_RCPP;
}
+#if 0
SEXP negativeBinomial_Create(SEXP fam_) {
BEGIN_RCPP;
negativeBinomial *ans = new negativeBinomial(List(fam_));
@@ -235,6 +235,7 @@
XPtr<negativeBinomial>(ptr)->setTheta(::Rf_asReal(newtheta));
END_RCPP;
}
+#endif
static inline double pwrss(lmResp *rp, merPredD *pp, double fac) {
return rp->wrss() + (fac ? pp->sqrL(fac) : pp->u0().squaredNorm());
@@ -937,9 +938,11 @@
CALLDEF(merPredDupdateRes, 2),
CALLDEF(merPredDupdateXwts, 2),
+#if 0
CALLDEF(negativeBinomial_Create, 1), // generate external pointer
CALLDEF(negativeBinomial_theta, 1),
CALLDEF(negativeBinomial_setTheta, 2),
+#endif
CALLDEF(NelderMead_Create, 5),
CALLDEF(NelderMead_newf, 2),
Modified: pkg/lme4Eigen/src/glmFamily.cpp
===================================================================
--- pkg/lme4Eigen/src/glmFamily.cpp 2012-02-27 11:14:49 UTC (rev 1622)
+++ pkg/lme4Eigen/src/glmFamily.cpp 2012-02-27 20:09:49 UTC (rev 1623)
@@ -1,3 +1,10 @@
+//
+// glmFamily.cpp: implementation of glmFamily and related classes using Eigen
+//
+// Copyright (C) 2011-2012 Douglas Bates, Martin Maechler and Ben Bolker
+//
+// This file is part of lme4Eigen.
+
#include "glmFamily.h"
#include <Rmath.h>
#include <limits>
@@ -7,29 +14,6 @@
using namespace std;
namespace glm {
- /**
- * Evaluate y * log(y/mu) with correct limit at y = 0
- *
- * @param y numerator of log argument, if non-zero
- * @param mu denominator of log argument
- *
- * @return y * log(y/mu) with correct limit at at y = 0
- */
- static inline double y_log_y(const double& y, const double& mu) {
- return y ? y * std::log(y/mu) : 0;
- }
-
- /**
- * Evaluate log(y/mu) returning 0 when y = 0
- *
- * @param y numerator of log argument, if non-zero
- * @param mu denominator of log argument
- *
- * @return log(y/mu) returning 0 when y = 0
- */
- static inline double logYMu(const double& y, const double& mu) {
- return y ? std::log(y/mu) : 0;
- }
/** Cumulative probability function of the complement of the Gumbel distribution
*
* (i.e. pgumbel(q,0.,1.,0) == 1 - pgumbel2(-q,0.,1.,0))
@@ -66,16 +50,20 @@
}
//@{ Templated scalar functors used in links, inverse links, etc.
+ template <typename T>
+ struct logN0 : public std::unary_function<T, T> {
+ const T operator()(const T& x) const {return x ? std::log(x) : T();}
+ };
+
+ static inline ArrayXd Y_log_Y(const ArrayXd& y, const ArrayXd& mu) {
+ return y * (y/mu).unaryExpr(logN0<double>());
+ }
+
template<typename T>
struct Round : public std::unary_function<T, T> {
const T operator()(const T& x) const {return nearbyint(x);}
};
- template <typename T>
- struct logN0 {
- const T operator()(const T& x) const {return x ? std::log(x) : T();}
- };
-
template<typename T>
struct x1mx : public std::unary_function<T, T> {
const T operator() (const T& x) const {
@@ -91,6 +79,27 @@
};
template<typename T>
+ struct cauchitinv : public std::unary_function<T, T> {
+ const T operator() (const T& x) const {
+ return T(::Rf_pcauchy(double(x), 0., 1., 1, 0));
+ }
+ };
+
+ template<typename T>
+ struct cauchit : public std::unary_function<T, T> {
+ const T operator() (const T& x) const {
+ return T(::Rf_qcauchy(double(x), 0., 1., 1, 0));
+ }
+ };
+
+ template<typename T>
+ struct cauchitmueta : public std::unary_function<T, T> {
+ const T operator() (const T& x) const {
+ return T(::Rf_dcauchy(double(x), 0., 1., 0));
+ }
+ };
+
+ template<typename T>
struct logitinv : public std::unary_function<T, T> {
const T operator() (const T& x) const {
return T(::Rf_plogis(double(x), 0., 1., 1, 0));
@@ -147,237 +156,229 @@
};
//@}
-#if 0
+ //@{
+ double binomialDist::aic (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
+ const ArrayXd& wt, double dev) const {
+ ArrayXd m((n > 1).any() ? n : wt);
+ ArrayXd yy((m * y).unaryExpr(Round<double>()));
+ m = m.unaryExpr(Round<double>());
+ double ans(0.);
+ for (int i=0; i < mu.size(); ++i)
+ ans += (m[i] <= 0. ? 0. : wt[i]/m[i]) * ::Rf_dbinom(yy[i], m[i], mu[i], true);
+ return (-2. * ans);
+ }
+ const ArrayXd binomialDist::devResid(const ArrayXd& y, const ArrayXd& mu, const ArrayXd& wt) const {
+ return 2. * wt * (Y_log_Y(y, mu) + Y_log_Y(1. - y, 1. - mu));
+ }
+ const ArrayXd binomialDist::variance(const ArrayXd& mu) const {return mu.unaryExpr(x1mx<double>());}
+ //@}
+
+ //@{
+ double gammaDist::aic (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
+ const ArrayXd& wt, double dev) const {
+ double nn(wt.sum());
+ double disp(dev/nn);
+ double ans(0), invdisp(1./disp);
+ for (int i = 0; i < mu.size(); ++i)
+ ans += wt[i] * ::Rf_dgamma(y[i], invdisp, mu[i] * disp, true);
+ return -2. * ans + 2.;
+ }
+ const ArrayXd gammaDist::devResid(const ArrayXd& y, const ArrayXd& mu, const ArrayXd& wt) const {
+ return -2. * wt * ((y/mu).unaryExpr(logN0<double>()) - (y - mu)/mu);
+ }
+ const ArrayXd gammaDist::variance(const ArrayXd& mu) const {return mu.square();}
+ //@}
+
+ //@{
+ double GaussianDist::aic (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
+ const ArrayXd& wt, double dev) const {
+ double nn(mu.size());
+ return nn * (std::log(2. * M_PI * dev/nn) + 1.) + 2. - wt.log().sum();
+ }
+ const ArrayXd GaussianDist::devResid(const ArrayXd& y, const ArrayXd& mu, const ArrayXd& wt) const {
+ return wt * (y - mu).square();
+ }
+ const ArrayXd GaussianDist::variance(const ArrayXd& mu) const {return ArrayXd::Ones(mu.size());}
+ //@}
+
+ //@{
+ double inverseGaussianDist::aic (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
+ const ArrayXd& wt, double dev) const {
+ double wtsum(wt.sum());
+ return wtsum * (std::log(dev/wtsum * 2. * M_PI) + 1.) + 3. * (y.log() * wt).sum() + 2.;
+ }
+ const ArrayXd inverseGaussianDist::devResid(const ArrayXd& y, const ArrayXd& mu, const ArrayXd& wt) const {
+ return wt * ((y - mu).square())/(y * mu.square());
+ }
+ const ArrayXd inverseGaussianDist::variance(const ArrayXd& mu) const {return mu.cube();}
+ //@}
+
+ //@{
+ double negativeBinomialDist::aic (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
+ const ArrayXd& wt, double dev) const {
+ return 2. * (wt * (y + d_theta) * (mu + d_theta).log() -
+ y * mu.log() + (y + 1).unaryExpr(Lgamma<double>()) -
+ d_theta * std::log(d_theta) + lgamma(d_theta) -
+ (d_theta + y).unaryExpr(Lgamma<double>())).sum();
+ }
+ const ArrayXd negativeBinomialDist::devResid(const ArrayXd &y, const ArrayXd &mu, const ArrayXd &wt) const {
+ return 2. * wt * (Y_log_Y(y, mu) - (y + d_theta) * ((y + d_theta)/(mu + d_theta)).log());
+ }
+ const ArrayXd negativeBinomialDist::variance(const ArrayXd &mu) const {
+ Rcpp::Rcout << "in negativeBinomial::variance with theta = " << d_theta << std::endl;
+ return mu + mu.square()/d_theta;
+ }
+ //@}
+
+ //@{
+ double PoissonDist::aic (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
+ const ArrayXd& wt, double dev) const {
+ double ans(0.);
+ for (int i = 0; i < mu.size(); ++i) ans += ::Rf_dpois(y[i], mu[i], true) * wt[i];
+ return (-2. * ans);
+ }
+ const ArrayXd PoissonDist::devResid(const ArrayXd& y, const ArrayXd& mu, const ArrayXd& wt) const {
+ return 2. * wt * (y * (y/mu).unaryExpr(logN0<double>()) - (y - mu));
+ }
+ const ArrayXd PoissonDist::variance(const ArrayXd& mu) const {return mu;}
+ //@}
+
+ //@{
+ const ArrayXd cauchitLink::linkFun(const ArrayXd& mu) const {return mu.unaryExpr(cauchit<double>());}
+ const ArrayXd cauchitLink::linkInv(const ArrayXd& eta) const {return eta.unaryExpr(cauchitinv<double>());}
+ const ArrayXd cauchitLink::muEta( const ArrayXd& eta) const {return eta.unaryExpr(cauchitmueta<double>());}
+ //@}
+
+ //@{
const ArrayXd logLink::linkFun(const ArrayXd& mu) const {return mu.log();}
const ArrayXd logLink::linkInv(const ArrayXd& eta) const {return eta.exp();}
const ArrayXd logLink::muEta( const ArrayXd& eta) const {return eta.exp();}
+ //@}
+ //@{
const ArrayXd logitLink::linkFun(const ArrayXd& mu) const {return mu.unaryExpr(logit<double>());}
const ArrayXd logitLink::linkInv(const ArrayXd& eta) const {return eta.unaryExpr(logitinv<double>());}
const ArrayXd logitLink::muEta( const ArrayXd& eta) const {return eta.unaryExpr(logitmueta<double>());}
+ //@}
+ //@{
const ArrayXd probitLink::linkFun(const ArrayXd& mu) const {return mu.unaryExpr(probit<double>());}
const ArrayXd probitLink::linkInv(const ArrayXd& eta) const {return eta.unaryExpr(probitinv<double>());}
const ArrayXd probitLink::muEta( const ArrayXd& eta) const {return eta.unaryExpr(probitmueta<double>());}
+ //@}
+ //@{
const ArrayXd identityLink::linkFun(const ArrayXd& mu) const {return mu;}
const ArrayXd identityLink::linkInv(const ArrayXd& eta) const {return eta;}
const ArrayXd identityLink::muEta( const ArrayXd& eta) const {return ArrayXd::Ones(eta.size());}
+ //@}
+ //@{
const ArrayXd inverseLink::linkFun(const ArrayXd& mu) const {return mu.inverse();}
const ArrayXd inverseLink::linkInv(const ArrayXd& eta) const {return eta.inverse();}
const ArrayXd inverseLink::muEta( const ArrayXd& eta) const {return -(eta.inverse().square());}
+ //@}
+ //@{
// const ArrayXd cloglogLink::linkFun(const ArrayXd& mu) const {return mu.unaryExpr(cloglog<double>());}
- const ArrayXd cloglogLink::linkFun(const ArrayXd& mu) const {return mu;}
const ArrayXd cloglogLink::linkInv(const ArrayXd& eta) const {return eta.unaryExpr(clogloginv<double>());}
const ArrayXd cloglogLink::muEta( const ArrayXd& eta) const {return eta.unaryExpr(cloglogmueta<double>());}
-
-#endif
- //@{ Vector functions for links, inverse links, derivatives and variances
- static ArrayXd cubef( const ArrayXd& x) {return x.cube();}
- static ArrayXd expf( const ArrayXd& x) {return x.exp();}
- static ArrayXd identf( const ArrayXd& x) {return x;}
- static ArrayXd invderivf( const ArrayXd& x) {return -(x.inverse().square());}
- static ArrayXd inversef( const ArrayXd& x) {return x.inverse();}
- static ArrayXd logf( const ArrayXd& x) {return x.log();}
- static ArrayXd onef( const ArrayXd& x) {return ArrayXd::Ones(x.size());}
- static ArrayXd sqrf( const ArrayXd& x) {return x.square();}
- static ArrayXd sqrtf( const ArrayXd& x) {return x.sqrt();}
- static ArrayXd twoxf( const ArrayXd& x) {return x * 2.;}
- static ArrayXd x1mxf( const ArrayXd& x) {return x.unaryExpr(x1mx<double>());}
- static ArrayXd logitf( const ArrayXd& x) {return x.unaryExpr(logit<double>());}
- static ArrayXd logitinvf( const ArrayXd& x) {return x.unaryExpr(logitinv<double>());}
- static ArrayXd logitmuf( const ArrayXd& x) {return x.unaryExpr(logitmueta<double>());}
- static ArrayXd probitf( const ArrayXd& x) {return x.unaryExpr(probit<double>());}
- static ArrayXd probitinvf(const ArrayXd& x) {return x.unaryExpr(probitinv<double>());}
- static ArrayXd probitmuf( const ArrayXd& x) {return x.unaryExpr(probitmueta<double>());}
- static ArrayXd clogloginf(const ArrayXd& x) {return x.unaryExpr(clogloginv<double>());}
- static ArrayXd cloglogmuf(const ArrayXd& x) {return x.unaryExpr(cloglogmueta<double>());}
//@}
-
- //@{ deviance residuals functions
- static inline Eigen::ArrayXd Y_log_Y(const Eigen::ArrayXd& y, const Eigen::ArrayXd& mu) {
- return y * (y/mu).unaryExpr(logN0<double>());
- }
- static inline ArrayXd
- BinomialDevRes(const Eigen::ArrayXd& y, const Eigen::ArrayXd& mu, const Eigen::ArrayXd& wt) {
- return 2. * wt * (Y_log_Y(y, mu) + Y_log_Y(1. - y, 1. - mu));
+ glmDist::glmDist(Rcpp::List& ll)
+ : d_devRes (as<SEXP>(ll["dev.resids"])),
+ d_variance(as<SEXP>(ll["variance"])),
+ d_aic( as<SEXP>(ll["aic"])),
+ d_rho( d_aic.environment()) {
}
- static inline ArrayXd
- GammaDevRes (const Eigen::ArrayXd& y, const Eigen::ArrayXd& mu, const Eigen::ArrayXd& wt) {
- return -2. * wt * ((y/mu).unaryExpr(logN0<double>()) - (y - mu)/mu);
- }
- static inline ArrayXd
- GaussianDevRes(const Eigen::ArrayXd& y, const Eigen::ArrayXd& mu, const Eigen::ArrayXd& wt) {
- return wt * (y - mu).square();
- }
- static inline ArrayXd
- PoissonDevRes (const Eigen::ArrayXd& y, const Eigen::ArrayXd& mu, const Eigen::ArrayXd& wt) {
- return 2. * wt * (y * (y/mu).unaryExpr(logN0<double>()) - (y - mu));
- }
- //@}
- //@{ AIC functions (which actually return the deviance, go figure)
- static inline double
- BinomialAIC (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
- const ArrayXd& wt, double dev) {
- Eigen::ArrayXd m((n > 1).any() ? n : wt);
- Eigen::ArrayXd yy((m * y).unaryExpr(Round<double>()));
- m = m.unaryExpr(Round<double>());
- double ans(0.);
- for (int i=0; i < mu.size(); ++i)
- ans += (m[i] <= 0. ? 0. : (wt[i]/m[i]) * ::Rf_dbinom(yy[i], m[i], mu[i], true));
- return (-2. * ans);
+ glmLink::glmLink(Rcpp::List& ll)
+ : d_linkFun(as<SEXP>(ll["linkfun"])),
+ d_linkInv(as<SEXP>(ll["linkinv"])),
+ d_muEta( as<SEXP>(ll["mu.eta"])),
+ d_rho( d_linkFun.environment()) {
}
- static inline double
- PoissonAIC (const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
- const ArrayXd& wt, double dev) {
- double ans(0.);
- for (int i=0; i < mu.size(); ++i) ans += ::Rf_dpois(y[i], mu[i], true) * wt[i];
- return (-2. * ans);
- }
- //@}
-
- // initialize the function maps (i.e. associative arrays of
- // functions). Needed because these are static maps.
- drmap glmFamily::devRes = drmap();
- aicmap glmFamily::aics = aicmap();
+ glmFamily::glmFamily(Rcpp::List ll)
+ : d_family( as<std::string>(as<SEXP>(ll["family"]))),
+ d_linknam(as<std::string>(as<SEXP>(ll["link"]))),
+ d_dist( new glmDist(ll)),
+ d_link( new glmLink(ll)) {
+ if (!ll.inherits("family"))
+ throw std::runtime_error("glmFamily requires a list of (S3) class \"family\"");
- fmap glmFamily::linvs = fmap();
- fmap glmFamily::lnks = fmap();
- fmap glmFamily::muEtas = fmap();
- fmap glmFamily::varFuncs = fmap();
-
- /**
- * Initialize the static maps. The identity link is guaranteed to be initialized if
- * any of the maps are initialized. FIXME: Should probably check the size of the map instead?
- *
- */
- void glmFamily::initMaps() {
- lnks ["log"] = &logf;
- linvs ["log"] = &expf;
- muEtas ["log"] = &expf;
-
- lnks ["sqrt"] = &sqrtf;
- linvs ["sqrt"] = &sqrf;
- muEtas ["sqrt"] = &twoxf;
-
- lnks ["identity"] = &identf;
- linvs ["identity"] = &identf;
- muEtas ["identity"] = &onef;
-
- lnks ["inverse"] = &inversef;
- linvs ["inverse"] = &inversef;
- muEtas ["inverse"] = &invderivf;
-
- lnks ["logit"] = &logitf;
- linvs ["logit"] = &logitinvf;
- muEtas ["logit"] = &logitmuf;
-
- lnks ["probit"] = &probitf;
- linvs ["probit"] = &probitinvf;
- muEtas ["probit"] = &probitmuf;
-
- linvs ["cloglog"] = &clogloginf;
- muEtas ["cloglog"] = &cloglogmuf;
-
- devRes ["Gamma"] = &GammaDevRes;
- varFuncs["Gamma"] = &sqrf; // x^2
+ if (d_linknam == "cauchit") {delete d_link; d_link = new cauchitLink(ll);}
+ if (d_linknam == "cloglog") {delete d_link; d_link = new cloglogLink(ll);}
+ if (d_linknam == "identity") {delete d_link; d_link = new identityLink(ll);}
+ if (d_linknam == "inverse") {delete d_link; d_link = new inverseLink(ll);}
+ if (d_linknam == "log") {delete d_link; d_link = new logLink(ll);}
+ if (d_linknam == "logit") {delete d_link; d_link = new logitLink(ll);}
+ if (d_linknam == "probit") {delete d_link; d_link = new probitLink(ll);}
- aics ["binomial"] = &BinomialAIC;
- devRes ["binomial"] = &BinomialDevRes;
- varFuncs["binomial"] = &x1mxf; // x * (1 - x)
+ if (d_family == "binomial") {delete d_dist; d_dist = new binomialDist(ll);}
+ if (d_family == "gamma") {delete d_dist; d_dist = new gammaDist(ll);}
+ if (d_family == "gaussian") {delete d_dist; d_dist = new GaussianDist(ll);}
+ if (d_family == "inverse.gaussian") {delete d_dist; d_dist = new inverseGaussianDist(ll);}
+ if (d_family.substr(0, 17) ==
+ "negative binomial") {delete d_dist; d_dist = new negativeBinomialDist(ll);}
+ if (d_family == "poisson") {delete d_dist; d_dist = new PoissonDist(ll);}
+ }
- devRes ["gaussian"] = &GaussianDevRes;
- varFuncs["gaussian"] = &onef; // 1
+ glmFamily::~glmFamily() {
+ delete d_dist;
+ delete d_link;
+ }
- varFuncs["inverse.gaussian"] = &cubef; // x^3
+ const ArrayXd glmFamily::devResid(const ArrayXd& y, const ArrayXd& mu, const ArrayXd& wt) const {
+ return d_dist->devResid(y, mu, wt);
+ }
- aics ["poisson"] = &PoissonAIC;
- devRes ["poisson"] = &PoissonDevRes;
- varFuncs["poisson"] = &identf; // x
+ double glmFamily::aic(const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
+ const ArrayXd& wt, double dev) const {
+ return d_dist->aic(y, n, mu, wt, dev);
}
-
- glmFamily::glmFamily(List ll)
- : d_family( as<std::string>(as<SEXP>(ll["family"]))),
- d_link( as<std::string>(as<SEXP>(ll["link"]))),
- d_devRes( as<SEXP>(ll["dev.resids"])),
- d_linkfun( as<SEXP>(ll["linkfun"])),
- d_linkinv( as<SEXP>(ll["linkinv"])),
- d_muEta( as<SEXP>(ll["mu.eta"])),
- d_variance(as<SEXP>(ll["variance"])),
- d_aic( as<SEXP>(ll["aic"])),
- d_rho( d_devRes.environment()) {
- if (!ll.inherits("family"))
- throw std::runtime_error("glmFamily requires a list of (S3) class \"family\"");
- if (!lnks.count("identity")) initMaps();
- }
- ArrayXd glmFamily::linkFun(const ArrayXd& mu) const {
- if (lnks.count(d_link)) return lnks[d_link](mu);
- return as<ArrayXd>(d_linkfun(NumericVector(mu.data(), mu.data() + mu.size())));
+ const ArrayXd glmLink::linkFun(const ArrayXd& mu) const {
+ return as<Eigen::VectorXd>(d_linkFun(NumericVector(mu.data(), mu.data() + mu.size())));
+// return as<ArrayXd>(d_linkfun(NumericVector(mu.data(), mu.data() + mu.size())));
}
- ArrayXd glmFamily::linkInv(const ArrayXd& eta) const {
- if (linvs.count(d_link)) return linvs[d_link](eta);
- return as<ArrayXd>(d_linkinv(NumericVector(eta.data(), eta.data() + eta.size())));
+ const ArrayXd glmLink::linkInv(const ArrayXd& eta) const {
+ return as<Eigen::VectorXd>(d_linkInv(NumericVector(eta.data(), eta.data() + eta.size())));
+// return as<ArrayXd>(d_linkinv(NumericVector(eta.data(), eta.data() + eta.size())));
}
- ArrayXd glmFamily::muEta(const ArrayXd &eta) const {
- if (muEtas.count(d_link)) return muEtas[d_link](eta);
- return as<ArrayXd>(as<SEXP>(d_muEta(NumericVector(eta.data(), eta.data() + eta.size()))));
+ const ArrayXd glmLink::muEta(const ArrayXd &eta) const {
+ return as<Eigen::VectorXd>(as<SEXP>(d_muEta(NumericVector(eta.data(), eta.data() + eta.size()))));
+// return as<ArrayXd>(as<SEXP>(d_muEta(NumericVector(eta.data(), eta.data() + eta.size()))));
}
- ArrayXd glmFamily::variance(const ArrayXd &mu) const {
- if (varFuncs.count(d_family)) return varFuncs[d_family](mu);
- return as<ArrayXd>(as<SEXP>(d_variance(NumericVector(mu.data(), mu.data() + mu.size()))));
+ const ArrayXd glmDist::variance(const ArrayXd &mu) const {
+ return as<Eigen::VectorXd>(as<SEXP>(d_variance(NumericVector(mu.data(), mu.data() + mu.size()))));
+// return as<ArrayXd>(as<SEXP>(d_variance(NumericVector(mu.data(), mu.data() + mu.size()))));
}
-/// FIXME: change this so the looping is done in the devResid[d_family] function
- ArrayXd
- glmFamily::devResid(const ArrayXd &mu, const ArrayXd &weights, const ArrayXd &y) const {
- if (devRes.count(d_family)) return devRes[d_family](y, mu, weights);
+ const ArrayXd glmDist::devResid(const ArrayXd &y, const ArrayXd &mu, const ArrayXd &wt) const {
int n = mu.size();
- return as<ArrayXd>(as<SEXP>(d_devRes(NumericVector(y.data(), y.data() + n),
- NumericVector(mu.data(), mu.data() + n),
- NumericVector(weights.data(), weights.data() + n))));
+// return as<ArrayXd>(as<SEXP>(d_devRes(NumericVector(y.data(), y.data() + n),
+ return as<Eigen::VectorXd>(as<SEXP>(d_devRes(NumericVector(y.data(), y.data() + n),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/lme4 -r 1623
More information about the Lme4-commits
mailing list