[Lme4-commits] r1729 - in pkg/lme4: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 15 19:16:14 CEST 2012
Author: dmbates
Date: 2012-05-15 19:16:14 +0200 (Tue, 15 May 2012)
New Revision: 1729
Modified:
pkg/lme4/R/lmer.R
pkg/lme4/src/external.cpp
pkg/lme4/src/predModule.cpp
pkg/lme4/src/predModule.h
Log:
Allow step-halving on glmerPwrssUpdate when compDev=TRUE
Modified: pkg/lme4/R/lmer.R
===================================================================
--- pkg/lme4/R/lmer.R 2012-05-15 17:14:45 UTC (rev 1728)
+++ pkg/lme4/R/lmer.R 2012-05-15 17:16:14 UTC (rev 1729)
@@ -329,8 +329,11 @@
c(reTrms[c("Zt","theta","Lambdat","Lind")],
n=nrow(X), list(X=X)))
rho$resp <- mkRespMod(fr, family=family)
+ if (length(unique(rho$resp$y)) < 2L)
+ stop("Response is constant - cannot fit the model")
+ rho$verbose <- as.integer(verbose)
# initialize (from mustart)
- glmerPwrssUpdate(rho$pp, rho$resp, tolPwrss, GHrule(0L), compDev)
+ .Call(glmerLaplace, pp$ptr(), resp$ptr(), 0L, tolPwrss, verbose)
rho$lp0 <- rho$pp$linPred(1) # each pwrss opt begins at this eta
rho$pwrssUpdate <- glmerPwrssUpdate
rho$compDev <- compDev
@@ -569,23 +572,24 @@
resp$resDev() + pp$sqrL(1)
}
-glmerPwrssUpdate <- function(pp, resp, tol, GQmat, compDev=TRUE, fac) {
+glmerPwrssUpdate <- function(pp, resp, tol, GQmat, compDev=TRUE, grpFac=NULL) {
nAGQ <- nrow(GQmat)
if (compDev) {
if (nAGQ < 2L)
- return(.Call(glmerLaplace, pp$ptr(), resp$ptr(), nAGQ, tol))
- return(.Call(glmerAGQ, pp$ptr(), resp$ptr(), tol, GQmat, fac))
+ return(.Call(glmerLaplace, pp$ptr(), resp$ptr(), nAGQ, tol, verbose))
+ return(.Call(glmerAGQ, pp$ptr(), resp$ptr(), tol, GQmat, grpFac, verbose))
}
oldpdev <- .Machine$double.xmax
+ uOnly <- nAGQ == 0L
repeat {
- ## FIXME:: how should uOnly be set here?
- pdev <- RglmerWrkIter(pp, resp, uOnly=FALSE)
+ oldu <- pp$delu
+ if (uOnly) oldbeta <- pp$delb
+ pdev <- RglmerWrkIter(pp, resp, uOnly=uOnly)
## check convergence first so small increases don't trigger errors
if (abs((oldpdev - pdev) / pdev) < tol)
break
-### FIXME: Replace the stop by step-halving
if (pdev > oldpdev) {
- stop("PIRLS step failed")
+ stop("PIRLS update failed")
}
oldpdev <- pdev
}
Modified: pkg/lme4/src/external.cpp
===================================================================
--- pkg/lme4/src/external.cpp 2012-05-15 17:14:45 UTC (rev 1728)
+++ pkg/lme4/src/external.cpp 2012-05-15 17:16:14 UTC (rev 1729)
@@ -263,22 +263,34 @@
return rp->resDev() + pp->sqrL(1.);
}
- static void pwrssUpdate(glmResp *rp, merPredD *pp, bool uOnly, double tol) {
+ static void pwrssUpdate(glmResp *rp, merPredD *pp, bool uOnly, double tol, int verbose) {
double oldpdev=std::numeric_limits<double>::max();
while (true) {
+ Vec olddelu(pp->delu()), olddelb(pp->delb());
double pdev=internal_glmerWrkIter(pp, rp, uOnly);
if (std::abs((oldpdev - pdev) / pdev) < tol) break;
- // FIXME: Replace the throwing of an error by step-halving
- if (pdev > oldpdev) throw std::invalid_argument("PIRLS step failed");
+ if (pdev > oldpdev) { // try step halving
+ if (verbose > 2)
+ Rcpp::Rcout << "Trying step halving: oldpdev = " << oldpdev
+ << ", pdev = " << pdev << std::endl;
+ for (int k=0; k < 10 && pdev > oldpdev; k++) {
+ pp->setDelu((olddelu + pp->delu())/2.);
+ if (!uOnly) pp->setDelb((olddelb + pp->delb())/2.);
+ pdev = internal_glmerWrkIter(pp, rp, uOnly);
+ if (verbose > 2)
+ Rcpp::Rcout << "k = " << k << ", pdev = " << pdev << std::endl;
+ }
+ if (pdev > oldpdev) throw runtime_error("PIRLS step failed");
+ }
oldpdev = pdev;
}
}
- SEXP glmerLaplace(SEXP pp_, SEXP rp_, SEXP nAGQ_, SEXP tol_) {
+ SEXP glmerLaplace(SEXP pp_, SEXP rp_, SEXP nAGQ_, SEXP tol_, SEXP verbose_) {
BEGIN_RCPP;
XPtr<glmResp> rp(rp_);
XPtr<merPredD> pp(pp_);
- pwrssUpdate(rp, pp, ::Rf_asInteger(nAGQ_), ::Rf_asReal(tol_));
+ pwrssUpdate(rp, pp, ::Rf_asInteger(nAGQ_), ::Rf_asReal(tol_), ::Rf_asInteger(verbose_));
return ::Rf_ScalarReal(rp->Laplace(pp->ldL2(), pp->ldRX2(), pp->sqrL(1.)));
END_RCPP;
}
@@ -291,17 +303,18 @@
static double sqrt2pi = std::sqrt(2. * PI);
- SEXP glmerAGQ(SEXP pp_, SEXP rp_, SEXP tol_, SEXP GQmat_, SEXP fac_) {
+ SEXP glmerAGQ(SEXP pp_, SEXP rp_, SEXP tol_, SEXP GQmat_, SEXP fac_, SEXP verbose_) {
BEGIN_RCPP;
XPtr<glmResp> rp(rp_);
XPtr<merPredD> pp(pp_);
const MiVec fac(as<MiVec>(fac_));
double tol(::Rf_asReal(tol_));
+ double verb(::Rf_asReal(verbose_));
if (fac.size() != rp->mu().size())
throw std::invalid_argument("size of fac must match dimension of response vector");
- pwrssUpdate(rp, pp, true, tol); // should be a no-op
+ pwrssUpdate(rp, pp, true, tol, verb); // should be a no-op
const Ar1 devc0(devcCol(fac, pp->u(1.), rp->devResid()));
const unsigned int q(pp->u0().size());
@@ -330,20 +343,6 @@
END_RCPP;
}
- // SEXP glmerLaplace(SEXP pp_, SEXP rp_, SEXP theta_, SEXP lp0_, SEXP uOnly_, SEXP tol_) {
- // BEGIN_RCPP;
-
- // XPtr<glmResp> rp(rp_);
- // XPtr<merPredD> pp(pp_);
-
- // pp->setTheta(as<MVec>(theta_));
- // rp->updateMu(as<MVec>(lp0_));
- // pwrssUpdate(rp, pp, ::Rf_asLogical(uOnly_), ::Rf_asReal(tol_));
- // return ::Rf_ScalarReal(rp->Laplace(pp->ldL2(), pp->ldRX2(), pp->sqrL(1.)));
-
- // END_RCPP;
- // }
-
void nstepFac(nlsResp *rp, merPredD *pp, int verb) {
double prss0(pwrss(rp, pp, 0.));
@@ -855,7 +854,7 @@
CALLDEF(glmFamily_variance, 2),
CALLDEF(glmerAGQ, 5),
- CALLDEF(glmerLaplace, 4),
+ CALLDEF(glmerLaplace, 6),
CALLDEF(golden_Create, 2),
CALLDEF(golden_newf, 2),
Modified: pkg/lme4/src/predModule.cpp
===================================================================
--- pkg/lme4/src/predModule.cpp 2012-05-15 17:14:45 UTC (rev 1728)
+++ pkg/lme4/src/predModule.cpp 2012-05-15 17:16:14 UTC (rev 1729)
@@ -260,6 +260,18 @@
std::copy(nBeta.data(), nBeta.data() + d_p, d_beta0.data());
}
+ void merPredD::setDelb(const VectorXd& newDelb) {
+ if (newDelb.size() != d_p)
+ throw invalid_argument("setDelb: dimension mismatch");
+ std::copy(newDelb.data(), newDelb.data() + d_p, d_delb.data());
+ }
+
+ void merPredD::setDelu(const VectorXd& newDelu) {
+ if (newDelu.size() != d_q)
+ throw invalid_argument("setDelu: dimension mismatch");
+ std::copy(newDelu.data(), newDelu.data() + d_q, d_delu.data());
+ }
+
void merPredD::setU0(const VectorXd& newU0) {
if (newU0.size() != d_q) throw invalid_argument("setU0: dimension mismatch");
std::copy(newU0.data(), newU0.data() + d_q, d_u0.data());
Modified: pkg/lme4/src/predModule.h
===================================================================
--- pkg/lme4/src/predModule.h 2012-05-15 17:14:45 UTC (rev 1728)
+++ pkg/lme4/src/predModule.h 2012-05-15 17:16:14 UTC (rev 1729)
@@ -87,6 +87,8 @@
void installPars(const Scalar& f);
void MCMC_beta_u(const Scalar& sigma);
void setBeta0(const VectorXd&);
+ void setDelb(const VectorXd&);
+ void setDelu(const VectorXd&);
void setTheta(const VectorXd&);
void setU0(const VectorXd&);
void updateDecomp();
More information about the Lme4-commits
mailing list