Nothing
R/clmm.ssr.R
In ordinal: Regression Models for Ordinal Data
Defines functions
update.uC
nllBase.uC
getNLA.ssr
getNAGQ.ssr
getNGHQ.ssr
clmm.fit.ssr
rho.clm2clmm.ssr
#############################################################################
## Copyright (c) 2010-2022 Rune Haubo Bojesen Christensen
##
## This file is part of the ordinal package for R (*ordinal*)
##
## *ordinal* is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 2 of the License, or
## (at your option) any later version.
##
## *ordinal* is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## A copy of the GNU General Public License is available at
## <https://www.r-project.org/Licenses/> and/or
## <http://www.gnu.org/licenses/>.
#############################################################################
## This file contains:
## Functions for fitting CLMMs with a single simple random-effects
## term (ssr).
rho.clm2clmm.ssr <- function(rho, retrms, ctrl)
### Version of rho.clm2clmm that is set up to use the C
### implementations of Laplace, AGQ and GHQ for a single random
### effect.
{
gfList <- retrms$gfList
rho$grFac <- gfList[[1]]
rho$ctrl <- ctrl
rho$sigma <- rep(1, nrow(rho$B1))
rho$lambda <- 0
rho$nlev <- as.vector(sapply(gfList, nlevels))
rho$random.names <- sapply(gfList, levels)
rho$tau.names <- names(gfList)
rho$nrandom <- sum(rho$nlev) ## no. random effects
rho$Niter <- 0L
rho$neval <- 0L
rho$u <- rho$uStart <- rep(0, rho$nrandom)
rho$linkInt <- switch(rho$link,
logit = 1L,
probit = 2L,
cloglog = 3L,
loglog = 4L,
cauchit = 5L)
rho$ST <- lapply(retrms$retrms, `[[`, "ST")
}
## set.AGQ <- function(rho, nAGQ) {
## rho$nAGQ <- nAGQ
## if(nAGQ %in% c(0L, 1L)) return(invisible())
## ghq <- gauss.hermite(abs(nAGQ))
## rho$ghqns <- ghq$nodes
## rho$ghqws <-
## if(nAGQ > 0) ghq$weights ## AGQ
## else log(ghq$weights) + (ghq$nodes^2)/2 ## GHQ
## }
clmm.fit.ssr <-
function(rho, control = list(), method=c("nlminb", "ucminf"),
Hess = FALSE)
### Fit a clmm with a single simple random effects term using AGQ, GHQ
### or Laplace.
{
optim.error <- function(fit, method)
if(inherits(fit, "try-error"))
stop("optimizer ", method, " terminated with an error", call.=FALSE)
### OPTION: Could have an argument c(warn, fail, ignore) to optionally
### return the fitted model despite the optimizer failing.
method <- match.arg(method)
## Set appropriate objective function:
obj.fun <-
if(rho$nAGQ < 0) getNGHQ.ssr
else if(rho$nAGQ > 1) getNAGQ.ssr
else getNLA.ssr ## nAGQ %in% c(0, 1)
init.val <- obj.fun(rho, rho$par)
if(!is.finite(init.val))
stop(gettextf("non-finite likelihood at starting value (%g)",
init.val), call.=FALSE)
## Fit the model:
if(method == "ucminf") {
fit <- try(ucminf(rho$par, function(par) obj.fun(rho, par),
control = control), silent=TRUE)
## Check if optimizer converged without error:
optim.error(fit, method)
## Save return value:
value <- fit$value
} else if(method == "nlminb") {
## hack to remove ucminf control settings:
keep <- !names(control) %in% c("grad", "grtol")
control <- if(length(keep)) control[keep] else list()
fit <- try(nlminb(rho$par, function(par) obj.fun(rho, par),
control = control), silent=TRUE)
## Check if optimizer converged without error:
optim.error(fit, method)
## Save return value:
value <- fit$objective
}
else stop("unkown optimization method: ", method)
## Extract parameters from optimizer results:
rho$par <- fit$par
## Ensure random mode estimation at optimum:
nllBase.uC(rho)
update.uC(rho)
rho$ST <- par2ST(rho$tau, rho$ST)
names(rho$ST) <- names(rho$dims$nlev.re)
## Format ranef modes and condVar:
ranef <- rho$u * rho$tau
condVar <- 1/rho$D * rho$tau^2
## names(ranef) <- names(condVar) <- rho$random.names
## ranef <- list(ranef)
## condVar <- list(condVar)
## names(ranef) <- names(condVar) <- rho$tau.names
## Prepare list of results:
res <- list(coefficients = fit$par[1:rho$dims$nfepar],
ST = rho$ST,
optRes = fit,
logLik = -value,
fitted.values = rho$fitted,
ranef = ranef,
condVar = condVar,
dims = rho$dims,
u = rho$u)
## Add gradient vector and optionally Hessian matrix:
## bound <- as.logical(paratBoundary2(rho))
## optpar <- fit$par[!bound]
if(Hess) {
## gH <- deriv12(function(par) obj.fun(rho, par, which=!bound),
gH <- deriv12(function(par) obj.fun(rho, par),
x=fit$par)
res$gradient <- gH$gradient
res$Hessian <- gH$Hessian
} else {
## res$gradient <- grad.ctr(function(par) getNLA(rho, par, which=!bound),
res$gradient <- grad.ctr(function(par) obj.fun(rho, par),
x=fit$par)
}
## Setting Niter and neval after gradient and Hessian evaluations:
res$Niter <- rho$Niter
res$neval <- rho$neval
return(res)
}
getNGHQ.ssr <- function(rho, par) {
### negative log-likelihood by standard Gauss-Hermite quadrature
### implemented in C:
if(!missing(par)) {
rho$par <- par
if(any(!is.finite(par)))
stop(gettextf(paste(c("Non-finite parameters occured:",
formatC(par, format="g")), collapse=" ")))
}
rho$neval <- rho$neval + 1L
nllBase.uC(rho) ## Update tau, eta1Fix and eta2Fix
with(rho, {
.C("getNGHQ_C",
nll = double(1),
as.integer(grFac),
as.double(tau),
as.double(eta1Fix),
as.double(eta2Fix),
as.double(o1),
as.double(o2),
as.double(sigma),
as.double(wts),
length(sigma),
length(uStart),
as.double(ghqns),
as.double(ghqws),
as.integer(abs(nAGQ)),
as.integer(linkInt),
as.double(ghqns * tau),
as.double(lambda))$nll
})
}
getNAGQ.ssr <- function(rho, par) {
### negative log-likelihood by adaptive Gauss-Hermite quadrature
### implemented in C:
if(!missing(par)) {
rho$par <- par
if(any(!is.finite(par)))
stop(gettextf(paste(c("Non-finite parameters occured:",
formatC(par, format="g")), collapse=" ")))
}
rho$neval <- rho$neval + 1L
if(!update.uC(rho)) return(Inf)
if(any(rho$D < 0)) return(Inf)
with(rho, {
.C("getNAGQ",
nll = double(1),
as.integer(grFac),
as.double(tau),
as.double(eta1Fix),
as.double(eta2Fix),
as.double(o1),
as.double(o2),
as.double(sigma),
as.double(wts),
length(sigma),
length(uStart),
as.double(ghqns),
as.double(log(ghqws)),
as.double(ghqns^2),
as.double(u),
as.double(D),
as.integer(abs(nAGQ)),
as.integer(linkInt),
as.double(lambda))$nll
})
}
getNLA.ssr <- function(rho, par) {
### negative log-likelihood by the Laplace approximation
### (with update.u2 in C or R):
if(!missing(par)) {
rho$par <- par
if(any(!is.finite(par)))
stop(gettextf(paste(c("Non-finite parameters occured:",
formatC(par, format="g")), collapse=" ")))
}
rho$neval <- rho$neval + 1L
if(!update.uC(rho)) return(Inf)
if(any(rho$D <= 0)) return(Inf)
logDetD <- sum(log(rho$D))
rho$negLogLik - rho$nrandom*log(2*pi)/2 + logDetD/2
}
nllBase.uC <- function(rho) {
### updates tau, eta1Fix and eta2Fix given new parameter values
with(rho, {
tau <- exp(par[nalpha + nbeta + 1:ntau])
eta1Fix <- drop(B1 %*% par[1:(nalpha + nbeta)])
eta2Fix <- drop(B2 %*% par[1:(nalpha + nbeta)])
})
return(invisible())
}
update.uC <- function(rho) {
### C-implementation of NR-algorithm.
nllBase.uC(rho) ## update: tau, eta1Fix, eta2Fix
fit <- with(rho, {
.C("NRalgv3",
as.integer(ctrl$trace),
as.integer(ctrl$maxIter),
as.double(ctrl$gradTol),
as.integer(ctrl$maxLineIter),
as.integer(grFac), ## OBS
as.double(tau), # stDev
as.double(o1),
as.double(o2),
as.double(eta1Fix),
as.double(eta2Fix),
as.double(sigma), ## rep(1, n)
as.integer(linkInt), ##
as.double(wts), ## pre. weights
u = as.double(uStart),
fitted = as.double(fitted), ## pre. pr
funValue = double(1),
gradValues = as.double(uStart),
hessValues = as.double(rep(1, length(uStart))),
length(fitted),
length(uStart),
maxGrad = double(1),
conv = 0L,
as.double(lambda), ##
Niter = as.integer(Niter) ## OBS
)[c("u", "fitted", "funValue", "gradValues",
"hessValues", "maxGrad", "conv", "Niter")] })
## Get message:
message <- switch(as.character(fit$conv),
"1" = "max|gradient| < tol, so current iterate is probably solution",
"0" = "Non finite negative log-likelihood",
"-1" = "iteration limit reached when updating the random effects",
"-2" = "step factor reduced below minimum when updating the random effects")
## Check for convergence and report warning/error:
if(rho$ctrl$trace > 0 && fit$conv == 1)
cat("\nOptimizer converged! ", "max|grad|:",
fit$maxGrad, message, fill = TRUE)
if(fit$conv != 1 && rho$ctrl$innerCtrl == "warnOnly")
warning(message, "\n at iteration ", rho$Niter)
else if(fit$conv != 1 && rho$ctrl$innerCtrl == "giveError")
stop(message, "\n at iteration ", rho$Niter)
## Store values and return:
rho$Niter <- fit$Niter
rho$fitted <- fit$fitted
rho$u <- fit$u
rho$D <- fit$hessValues
rho$gradient <- fit$gradValues
if(!is.finite(rho$negLogLik <- fit$funValue))
return(FALSE)
return(TRUE)
}
Try the ordinal package in your browser
Any scripts or data that you put into this service are public.
ordinal documentation built on Sept. 11, 2024, 7:44 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions update.uC nllBase.uC getNLA.ssr getNAGQ.ssr getNGHQ.ssr clmm.fit.ssr rho.clm2clmm.ssr
#############################################################################
## Copyright (c) 2010-2022 Rune Haubo Bojesen Christensen
##
## This file is part of the ordinal package for R (*ordinal*)
##
## *ordinal* is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 2 of the License, or
## (at your option) any later version.
##
## *ordinal* is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## A copy of the GNU General Public License is available at
## <https://www.r-project.org/Licenses/> and/or
## <http://www.gnu.org/licenses/>.
#############################################################################
## This file contains:
## Functions for fitting CLMMs with a single simple random-effects
## term (ssr).
rho.clm2clmm.ssr <- function(rho, retrms, ctrl)
### Version of rho.clm2clmm that is set up to use the C
### implementations of Laplace, AGQ and GHQ for a single random
### effect.
{
gfList <- retrms$gfList
rho$grFac <- gfList[[1]]
rho$ctrl <- ctrl
rho$sigma <- rep(1, nrow(rho$B1))
rho$lambda <- 0
rho$nlev <- as.vector(sapply(gfList, nlevels))
rho$random.names <- sapply(gfList, levels)
rho$tau.names <- names(gfList)
rho$nrandom <- sum(rho$nlev) ## no. random effects
rho$Niter <- 0L
rho$neval <- 0L
rho$u <- rho$uStart <- rep(0, rho$nrandom)
rho$linkInt <- switch(rho$link,
logit = 1L,
probit = 2L,
cloglog = 3L,
loglog = 4L,
cauchit = 5L)
rho$ST <- lapply(retrms$retrms, `[[`, "ST")
}
## set.AGQ <- function(rho, nAGQ) {
## rho$nAGQ <- nAGQ
## if(nAGQ %in% c(0L, 1L)) return(invisible())
## ghq <- gauss.hermite(abs(nAGQ))
## rho$ghqns <- ghq$nodes
## rho$ghqws <-
## if(nAGQ > 0) ghq$weights ## AGQ
## else log(ghq$weights) + (ghq$nodes^2)/2 ## GHQ
## }
clmm.fit.ssr <-
function(rho, control = list(), method=c("nlminb", "ucminf"),
Hess = FALSE)
### Fit a clmm with a single simple random effects term using AGQ, GHQ
### or Laplace.
{
optim.error <- function(fit, method)
if(inherits(fit, "try-error"))
stop("optimizer ", method, " terminated with an error", call.=FALSE)
### OPTION: Could have an argument c(warn, fail, ignore) to optionally
### return the fitted model despite the optimizer failing.
method <- match.arg(method)
## Set appropriate objective function:
obj.fun <-
if(rho$nAGQ < 0) getNGHQ.ssr
else if(rho$nAGQ > 1) getNAGQ.ssr
else getNLA.ssr ## nAGQ %in% c(0, 1)
init.val <- obj.fun(rho, rho$par)
if(!is.finite(init.val))
stop(gettextf("non-finite likelihood at starting value (%g)",
init.val), call.=FALSE)
## Fit the model:
if(method == "ucminf") {
fit <- try(ucminf(rho$par, function(par) obj.fun(rho, par),
control = control), silent=TRUE)
## Check if optimizer converged without error:
optim.error(fit, method)
## Save return value:
value <- fit$value
} else if(method == "nlminb") {
## hack to remove ucminf control settings:
keep <- !names(control) %in% c("grad", "grtol")
control <- if(length(keep)) control[keep] else list()
fit <- try(nlminb(rho$par, function(par) obj.fun(rho, par),
control = control), silent=TRUE)
## Check if optimizer converged without error:
optim.error(fit, method)
## Save return value:
value <- fit$objective
}
else stop("unkown optimization method: ", method)
## Extract parameters from optimizer results:
rho$par <- fit$par
## Ensure random mode estimation at optimum:
nllBase.uC(rho)
update.uC(rho)
rho$ST <- par2ST(rho$tau, rho$ST)
names(rho$ST) <- names(rho$dims$nlev.re)
## Format ranef modes and condVar:
ranef <- rho$u * rho$tau
condVar <- 1/rho$D * rho$tau^2
## names(ranef) <- names(condVar) <- rho$random.names
## ranef <- list(ranef)
## condVar <- list(condVar)
## names(ranef) <- names(condVar) <- rho$tau.names
## Prepare list of results:
res <- list(coefficients = fit$par[1:rho$dims$nfepar],
ST = rho$ST,
optRes = fit,
logLik = -value,
fitted.values = rho$fitted,
ranef = ranef,
condVar = condVar,
dims = rho$dims,
u = rho$u)
## Add gradient vector and optionally Hessian matrix:
## bound <- as.logical(paratBoundary2(rho))
## optpar <- fit$par[!bound]
if(Hess) {
## gH <- deriv12(function(par) obj.fun(rho, par, which=!bound),
gH <- deriv12(function(par) obj.fun(rho, par),
x=fit$par)
res$gradient <- gH$gradient
res$Hessian <- gH$Hessian
} else {
## res$gradient <- grad.ctr(function(par) getNLA(rho, par, which=!bound),
res$gradient <- grad.ctr(function(par) obj.fun(rho, par),
x=fit$par)
}
## Setting Niter and neval after gradient and Hessian evaluations:
res$Niter <- rho$Niter
res$neval <- rho$neval
return(res)
}
getNGHQ.ssr <- function(rho, par) {
### negative log-likelihood by standard Gauss-Hermite quadrature
### implemented in C:
if(!missing(par)) {
rho$par <- par
if(any(!is.finite(par)))
stop(gettextf(paste(c("Non-finite parameters occured:",
formatC(par, format="g")), collapse=" ")))
}
rho$neval <- rho$neval + 1L
nllBase.uC(rho) ## Update tau, eta1Fix and eta2Fix
with(rho, {
.C("getNGHQ_C",
nll = double(1),
as.integer(grFac),
as.double(tau),
as.double(eta1Fix),
as.double(eta2Fix),
as.double(o1),
as.double(o2),
as.double(sigma),
as.double(wts),
length(sigma),
length(uStart),
as.double(ghqns),
as.double(ghqws),
as.integer(abs(nAGQ)),
as.integer(linkInt),
as.double(ghqns * tau),
as.double(lambda))$nll
})
}
getNAGQ.ssr <- function(rho, par) {
### negative log-likelihood by adaptive Gauss-Hermite quadrature
### implemented in C:
if(!missing(par)) {
rho$par <- par
if(any(!is.finite(par)))
stop(gettextf(paste(c("Non-finite parameters occured:",
formatC(par, format="g")), collapse=" ")))
}
rho$neval <- rho$neval + 1L
if(!update.uC(rho)) return(Inf)
if(any(rho$D < 0)) return(Inf)
with(rho, {
.C("getNAGQ",
nll = double(1),
as.integer(grFac),
as.double(tau),
as.double(eta1Fix),
as.double(eta2Fix),
as.double(o1),
as.double(o2),
as.double(sigma),
as.double(wts),
length(sigma),
length(uStart),
as.double(ghqns),
as.double(log(ghqws)),
as.double(ghqns^2),
as.double(u),
as.double(D),
as.integer(abs(nAGQ)),
as.integer(linkInt),
as.double(lambda))$nll
})
}
getNLA.ssr <- function(rho, par) {
### negative log-likelihood by the Laplace approximation
### (with update.u2 in C or R):
if(!missing(par)) {
rho$par <- par
if(any(!is.finite(par)))
stop(gettextf(paste(c("Non-finite parameters occured:",
formatC(par, format="g")), collapse=" ")))
}
rho$neval <- rho$neval + 1L
if(!update.uC(rho)) return(Inf)
if(any(rho$D <= 0)) return(Inf)
logDetD <- sum(log(rho$D))
rho$negLogLik - rho$nrandom*log(2*pi)/2 + logDetD/2
}
nllBase.uC <- function(rho) {
### updates tau, eta1Fix and eta2Fix given new parameter values
with(rho, {
tau <- exp(par[nalpha + nbeta + 1:ntau])
eta1Fix <- drop(B1 %*% par[1:(nalpha + nbeta)])
eta2Fix <- drop(B2 %*% par[1:(nalpha + nbeta)])
})
return(invisible())
}
update.uC <- function(rho) {
### C-implementation of NR-algorithm.
nllBase.uC(rho) ## update: tau, eta1Fix, eta2Fix
fit <- with(rho, {
.C("NRalgv3",
as.integer(ctrl$trace),
as.integer(ctrl$maxIter),
as.double(ctrl$gradTol),
as.integer(ctrl$maxLineIter),
as.integer(grFac), ## OBS
as.double(tau), # stDev
as.double(o1),
as.double(o2),
as.double(eta1Fix),
as.double(eta2Fix),
as.double(sigma), ## rep(1, n)
as.integer(linkInt), ##
as.double(wts), ## pre. weights
u = as.double(uStart),
fitted = as.double(fitted), ## pre. pr
funValue = double(1),
gradValues = as.double(uStart),
hessValues = as.double(rep(1, length(uStart))),
length(fitted),
length(uStart),
maxGrad = double(1),
conv = 0L,
as.double(lambda), ##
Niter = as.integer(Niter) ## OBS
)[c("u", "fitted", "funValue", "gradValues",
"hessValues", "maxGrad", "conv", "Niter")] })
## Get message:
message <- switch(as.character(fit$conv),
"1" = "max|gradient| < tol, so current iterate is probably solution",
"0" = "Non finite negative log-likelihood",
"-1" = "iteration limit reached when updating the random effects",
"-2" = "step factor reduced below minimum when updating the random effects")
## Check for convergence and report warning/error:
if(rho$ctrl$trace > 0 && fit$conv == 1)
cat("\nOptimizer converged! ", "max|grad|:",
fit$maxGrad, message, fill = TRUE)
if(fit$conv != 1 && rho$ctrl$innerCtrl == "warnOnly")
warning(message, "\n at iteration ", rho$Niter)
else if(fit$conv != 1 && rho$ctrl$innerCtrl == "giveError")
stop(message, "\n at iteration ", rho$Niter)
## Store values and return:
rho$Niter <- fit$Niter
rho$fitted <- fit$fitted
rho$u <- fit$u
rho$D <- fit$hessValues
rho$gradient <- fit$gradValues
if(!is.finite(rho$negLogLik <- fit$funValue))
return(FALSE)
return(TRUE)
}
Try the ordinal package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.