[Lme4-commits] r1653 - in pkg/lme4Eigen: . R inst/tests man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Mar 14 23:24:28 CET 2012
Author: dmbates
Date: 2012-03-14 23:24:27 +0100 (Wed, 14 Mar 2012)
New Revision: 1653
Modified:
pkg/lme4Eigen/DESCRIPTION
pkg/lme4Eigen/NAMESPACE
pkg/lme4Eigen/R/AllClass.R
pkg/lme4Eigen/inst/tests/test-glmFamily.R
pkg/lme4Eigen/man/mcmcsamp.Rd
pkg/lme4Eigen/man/profile-methods.Rd
pkg/lme4Eigen/src/external.cpp
pkg/lme4Eigen/src/glmFamily.cpp
pkg/lme4Eigen/src/glmFamily.h
pkg/lme4Eigen/src/respModule.h
Log:
Create theta() and setTheta() methods for negative binomial. Bump the version number.
Modified: pkg/lme4Eigen/DESCRIPTION
===================================================================
--- pkg/lme4Eigen/DESCRIPTION 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/DESCRIPTION 2012-03-14 22:24:27 UTC (rev 1653)
@@ -1,5 +1,5 @@
Package: lme4Eigen
-Version: 0.9996875-14
+Version: 0.9996875-15
Date: $Date$
Revision: $Revision$
Title: Linear mixed-effects models using Eigen and S4
Modified: pkg/lme4Eigen/NAMESPACE
===================================================================
--- pkg/lme4Eigen/NAMESPACE 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/NAMESPACE 2012-03-14 22:24:27 UTC (rev 1653)
@@ -117,7 +117,6 @@
S3method(nobs,merMod)
S3method(plot,coef.mer)
S3method(plot,lmList.confint)
-S3method(plot,merMod)
S3method(plot,ranef.mer)
S3method(predict,merMod)
S3method(print,merMod)
Modified: pkg/lme4Eigen/R/AllClass.R
===================================================================
--- pkg/lme4Eigen/R/AllClass.R 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/R/AllClass.R 2012-03-14 22:24:27 UTC (rev 1653)
@@ -162,7 +162,7 @@
},
unsc = function() {
'the unscaled variance-covariance matrix of the fixed-effects parameters'
- tcrossprod(RXi())
+ .Call(merPredDunsc, ptr())
},
linPred = function(fac) {
'evaluate the linear predictor for step factor fac'
@@ -507,10 +507,18 @@
'returns the sum of the deviance residuals'
.Call(glm_resDev, ptr())
},
+ setTheta = function(theta) {
+ 'sets a new value of theta, for negative binomial distribution only'
+ .Call(glm_setTheta, ptr(), as.numeric(theta))
+ },
sqrtWrkWt = function() {
'returns the square root of the working X weights'
.Call(glm_sqrtWrkWt, ptr())
},
+ theta = function() {
+ 'query the value of theta, for negative binomial distribution only'
+ .Call(glm_theta, ptr())
+ },
updateMu = function(gamma) {
'update mu, residuals, weights, etc. from the linear predictor'
.Call(glm_updateMu, ptr(), as.numeric(gamma))
@@ -642,6 +650,14 @@
Ptr <<- .Call(glmFamily_Create, family)
Ptr
},
+ setTheta = function(theta) {
+ 'sets a new value of theta, for negative binomial distribution only'
+ .Call(glmFamily_setTheta, ptr(), as.numeric(theta))
+ },
+ theta = function() {
+ 'query the value of theta, for negative binomial distribution only'
+ .Call(glmFamily_theta, ptr())
+ },
variance = function(mu) {
'applies the variance function to mu'
.Call(glmFamily_variance, ptr(), as.numeric(mu))
@@ -950,3 +966,45 @@
##' showClass("rePos")
##'
rePos$lock("cnms", "flist", "ncols", "nctot", "nlevs", "terms")
+
+vcRep <-
+ setRefClass("vcRep",
+ fields =
+ list(
+ theta = "numeric",
+ Lambdat = "dgCMatrix",
+ Lind = "integer",
+ Gp = "integer",
+ flist = "list",
+ cnms = "list",
+ ncols = "integer",
+ nctot = "integer",
+ nlevs = "integer",
+ offsets = "integer",
+ terms = "integer"
+ ),
+ methods =
+ list(
+ initialize = function(mer, ...) {
+ stopifnot((ntrms <- length(Cnms <- mer at cnms)) > 0L,
+ (length(Flist <- mer at flist)) > 0L,
+ length(asgn <- as.integer(attr(Flist, "assign"))) == ntrms)
+ theta <<- mer at pp$theta
+ Lambdat <<- mer at pp$dgCMatrix
+ Lind <<- mer at pp$Lind
+ Gp <<- mer at Gp
+ cnms <<- Cnms
+ flist <<- Flist
+ ncols <<- unname(vapply(cnms, length, 0L))
+ nctot <<- unname(as.vector(tapply(ncols, asgn, sum)))
+ nlevs <<- unname(vapply(flist, function(el) length(levels(el)), 0L))
+ offsets <<- c(0L, cumsum(sapply(seq_along(asgn),
+ function(i) ncols[i] * nlevs[asgn[i]])))
+ terms <<- lapply(seq_along(flist), function(i) which(asgn == i))
+ },
+ asCovar = function() {}
+ )
+ )
+
+
+
Modified: pkg/lme4Eigen/inst/tests/test-glmFamily.R
===================================================================
--- pkg/lme4Eigen/inst/tests/test-glmFamily.R 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/inst/tests/test-glmFamily.R 2012-03-14 22:24:27 UTC (rev 1653)
@@ -103,3 +103,25 @@
tst.devres(gaussian(), etas)
tst.devres(poisson(), etapos)
})
+
+context("negative binomial")
+test_that("variance", {
+ tst.variance <- function(fam, frm) {
+ ff <- glmFamily$new(family=fam)
+ sapply(frm, function(x) expect_that(fam$variance(x), equals(ff$variance(x))))
+ }
+ tst.variance(MASS::negative.binomial(1.), etapos)
+ nb1 <- MASS::negative.binomial(1.)
+ cppnb1 <- glmFamily$new(family=nb1)
+ expect_that(cppnb1$theta(), equals(1))
+ nb2 <- MASS::negative.binomial(2.)
+ cppnb1$setTheta(2)
+ sapply(etapos, function(x) expect_that(cppnb1$variance(x), equals(nb2$variance(x))))
+ bfam <- glmFamily$new(family=binomial())
+ expect_that(bfam$theta(), throws_error("theta accessor applies only to negative binomial"))
+ expect_that(bfam$setTheta(2), throws_error("setTheta applies only to negative binomial"))
+})
+
+
+
+
Modified: pkg/lme4Eigen/man/mcmcsamp.Rd
===================================================================
--- pkg/lme4Eigen/man/mcmcsamp.Rd 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/man/mcmcsamp.Rd 2012-03-14 22:24:27 UTC (rev 1653)
@@ -12,7 +12,7 @@
\item{verbose}{should verbose output be given?}
- \item{...}{additional arguments}
+ \item{...}{}
}
\value{
a Markov chain Monte Carlo sample as a matrix
Modified: pkg/lme4Eigen/man/profile-methods.Rd
===================================================================
--- pkg/lme4Eigen/man/profile-methods.Rd 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/man/profile-methods.Rd 2012-03-14 22:24:27 UTC (rev 1653)
@@ -5,21 +5,32 @@
\title{Methods for profile() of [ng]lmer fitted models}
\usage{
\method{profile}{merMod} (fitted, alphamax = 0.01,
- maxpts = 100, delta = cutoff/8, verbose=0, devtol=1e-9,
-startval="prev", ...)
+ maxpts = 100, delta = cutoff/8, verbose = 0,
+ devtol = 1e-09, startval = "prev",
+ optimizer = "bobyqa", ...)
}
-%% FIXME: update from Roxygen!
\arguments{
\item{fitted}{a fitted model, e.g., the result of
\code{\link{lmer}(..)}.}
- \item{alphamax}{used when \code{delta} is unspecified, as
- probability ... to compute \code{delta} ...}
- \item{maxpts}{...}
- \item{delta}{...}
- \item{verbose}{...}
- \item{devtol}{...}
- \item{startval}{...}
- \item{\dots}{potential further arguments for
+
+ \item{alphamax}{maximum alpha value for likelihood ratio
+ confidence regions; used to establish the range of values
+ to be profiled}
+
+ \item{maxpts}{maximum number of points (in each
+ direction, for each parameter) to evaluate in attempting
+ to construct the profile}
+
+ \item{delta}{stepping scale for deciding on next point to
+ profile}
+
+ \item{verbose}{level of output from internal
+ calculations}
+
+ \item{devtol}{tolerance for fitted deviances less than
+ baseline (supposedly minimum) deviance}
+
+ \item{\dots}{potential further arguments for
\code{profile} methods.}
}
\description{
Modified: pkg/lme4Eigen/src/external.cpp
===================================================================
--- pkg/lme4Eigen/src/external.cpp 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/src/external.cpp 2012-03-14 22:24:27 UTC (rev 1653)
@@ -125,12 +125,24 @@
END_RCPP;
}
+ SEXP glm_setTheta(SEXP ptr, SEXP newtheta) {
+ BEGIN_RCPP;
+ XPtr<glmResp>(ptr)->setTheta(::Rf_asReal(newtheta));
+ END_RCPP;
+ }
+
SEXP glm_sqrtWrkWt(SEXP ptr_) {
BEGIN_RCPP;
return wrap(XPtr<glmResp>(ptr_)->sqrtWrkWt());
END_RCPP;
}
+ SEXP glm_theta(SEXP ptr) {
+ BEGIN_RCPP;
+ return ::Rf_ScalarReal(XPtr<glmResp>(ptr)->theta());
+ END_RCPP;
+ }
+
SEXP glm_updateWts(SEXP ptr_) {
BEGIN_RCPP;
return ::Rf_ScalarReal(XPtr<glmResp>(ptr_)->updateWts());
@@ -214,33 +226,24 @@
END_RCPP;
}
- SEXP glmFamily_variance(SEXP ptr, SEXP mu) {
+ SEXP glmFamily_setTheta(SEXP ptr, SEXP ntheta) {
BEGIN_RCPP;
- return wrap(XPtr<glmFamily>(ptr)->variance(as<MVec>(mu)));
+ XPtr<glmFamily>(ptr)->setTheta(::Rf_asReal(ntheta));
END_RCPP;
}
-#if 0
- SEXP negativeBinomial_Create(SEXP fam_) {
+ SEXP glmFamily_theta(SEXP ptr) {
BEGIN_RCPP;
- negativeBinomial *ans = new negativeBinomial(List(fam_));
- return wrap(XPtr<negativeBinomial>(ans, true));
+ return ::Rf_ScalarReal(XPtr<glmFamily>(ptr)->theta());
END_RCPP;
}
- SEXP negativeBinomial_theta(SEXP ptr) {
+ SEXP glmFamily_variance(SEXP ptr, SEXP mu) {
BEGIN_RCPP;
- return ::Rf_ScalarReal(XPtr<negativeBinomial>(ptr)->theta());
+ return wrap(XPtr<glmFamily>(ptr)->variance(as<MVec>(mu)));
END_RCPP;
}
- SEXP negativeBinomial_setTheta(SEXP ptr, SEXP newtheta) {
- BEGIN_RCPP;
- 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());
}
@@ -854,7 +857,9 @@
CALLDEF(glm_link, 1),
CALLDEF(glm_muEta, 1),
CALLDEF(glm_resDev, 1),
+ CALLDEF(glm_setTheta, 2),
CALLDEF(glm_sqrtWrkWt, 1),
+ CALLDEF(glm_theta, 1),
CALLDEF(glm_variance, 1),
CALLDEF(glm_wrkResids, 1),
CALLDEF(glm_wrkResp, 1),
@@ -870,6 +875,8 @@
CALLDEF(glmFamily_linkInv, 2),
CALLDEF(glmFamily_devResid, 4),
CALLDEF(glmFamily_muEta, 2),
+ CALLDEF(glmFamily_setTheta, 2),
+ CALLDEF(glmFamily_theta, 1),
CALLDEF(glmFamily_variance, 2),
CALLDEF(glmerAGQ, 8),
@@ -932,12 +939,6 @@
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),
CALLDEF(NelderMead_setForce_stop, 2),
Modified: pkg/lme4Eigen/src/glmFamily.cpp
===================================================================
--- pkg/lme4Eigen/src/glmFamily.cpp 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/src/glmFamily.cpp 2012-03-14 22:24:27 UTC (rev 1653)
@@ -11,7 +11,6 @@
#include <cmath>
using namespace Rcpp;
-using namespace std;
namespace glm {
/** Cumulative probability function of the complement of the Gumbel distribution
@@ -61,7 +60,7 @@
template<typename T>
struct Round : public std::unary_function<T, T> {
- const T operator()(const T& x) const {return nearbyint(x);}
+ const T operator()(const T& x) const {return std::nearbyint(x);}
};
template<typename T>
@@ -225,7 +224,6 @@
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;
}
//@}
@@ -319,8 +317,8 @@
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.substr(0, 18) ==
+ "Negative Binomial(") {delete d_dist; d_dist = new negativeBinomialDist(ll);}
if (d_family == "poisson") {delete d_dist; d_dist = new PoissonDist(ll);}
}
@@ -339,46 +337,64 @@
}
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())));
+ return as<ArrayXd>(::Rf_eval(::Rf_lang2(as<SEXP>(d_linkFun),
+ as<SEXP>(Rcpp::NumericVector(mu.data(),
+ mu.data() + mu.size()))
+ ), d_rho));
}
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())));
+ return as<ArrayXd>(::Rf_eval(::Rf_lang2(as<SEXP>(d_linkInv),
+ as<SEXP>(Rcpp::NumericVector(eta.data(),
+ eta.data() + eta.size()))
+ ), d_rho));
}
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()))));
+ return as<ArrayXd>(::Rf_eval(::Rf_lang2(as<SEXP>(d_muEta),
+ as<SEXP>(Rcpp::NumericVector(eta.data(),
+ eta.data() + eta.size()))
+ ), d_rho));
}
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()))));
+ return as<ArrayXd>(::Rf_eval(::Rf_lang2(as<SEXP>(d_variance),
+ as<SEXP>(Rcpp::NumericVector(mu.data(),
+ mu.data() + mu.size()))
+ ), d_rho));
}
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),
- return as<Eigen::VectorXd>(as<SEXP>(d_devRes(NumericVector(y.data(), y.data() + n),
- NumericVector(mu.data(), mu.data() + n),
- NumericVector(wt.data(), wt.data() + n))));
+ return as<ArrayXd>(::Rf_eval(::Rf_lang4(as<SEXP>(d_devRes),
+ as<SEXP>(NumericVector(y.data(), y.data() + n)),
+ as<SEXP>(NumericVector(mu.data(), mu.data() + n)),
+ as<SEXP>(NumericVector(wt.data(), wt.data() + n))
+ ), d_rho));
}
double glmDist::aic(const ArrayXd& y, const ArrayXd& n, const ArrayXd& mu,
const ArrayXd& wt, double dev) const {
int nn = mu.size();
- SEXP ans = d_aic(NumericVector(y.data(), y.data() + nn),
- NumericVector(n.data(), n.data() + nn),
- NumericVector(mu.data(), mu.data() + nn),
- NumericVector(wt.data(), wt.data() + nn),
- ::Rf_ScalarReal(dev));
- return ::Rf_asReal(ans);
+ return ::Rf_asReal(::Rf_eval(::Rf_lang6(as<SEXP>(d_aic),
+ as<SEXP>(NumericVector(y.data(), y.data() + nn)),
+ as<SEXP>(NumericVector(n.data(), n.data() + nn)),
+ as<SEXP>(NumericVector(mu.data(), mu.data() + nn)),
+ as<SEXP>(NumericVector(wt.data(), wt.data() + nn)),
+ ::Rf_ScalarReal(dev)
+ ), d_rho));
}
negativeBinomialDist::negativeBinomialDist(Rcpp::List& ll)
: glmDist(ll),
d_theta(::Rf_asReal(as<SEXP>(d_rho[".Theta"]))) {}
+
+ double glmDist::theta() const {
+ throw std::invalid_argument("theta accessor applies only to negative binomial");
+ }
+
+ void glmDist::setTheta(const double& theta) {
+ throw std::invalid_argument("setTheta applies only to negative binomial");
+ }
}
Modified: pkg/lme4Eigen/src/glmFamily.h
===================================================================
--- pkg/lme4Eigen/src/glmFamily.h 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/src/glmFamily.h 2012-03-14 22:24:27 UTC (rev 1653)
@@ -29,6 +29,8 @@
const ArrayXd&, double) const;
/**< in keeping with the botched up nomenclature in the R glm function,
* the value of aic is the deviance */
+ virtual double theta() const;
+ virtual void setTheta(const double&);
};
class binomialDist : public glmDist {
@@ -186,6 +188,8 @@
const ArrayXd variance(const ArrayXd& mu) const {return d_dist->variance(mu);}
double aic(const ArrayXd&, const ArrayXd&, const ArrayXd&,
const ArrayXd&, double) const;
+ double theta() const {return d_dist->theta();}
+ void setTheta(const double& theta) {d_dist->setTheta(theta);}
//@}
};
Modified: pkg/lme4Eigen/src/respModule.h
===================================================================
--- pkg/lme4Eigen/src/respModule.h 2012-03-14 13:24:11 UTC (rev 1652)
+++ pkg/lme4Eigen/src/respModule.h 2012-03-14 22:24:27 UTC (rev 1653)
@@ -104,10 +104,14 @@
double aic() const;
double Laplace(double,double,double) const;
double resDev() const;
+ double theta() const {return d_fam.theta();}
+ //< negative binomial distribution only
double updateMu(const Eigen::VectorXd&);
double updateWts();
void setN(const Eigen::VectorXd&);
+ void setTheta(const double& ntheta) {d_fam.setTheta(ntheta);}
+ // negative binomial distribution only
};
class nlsResp : public lmResp {
More information about the Lme4-commits
mailing list