R/clmm.methods.R In ordinal: Regression Models for Ordinal Data

Documented in VarCorrVarCorr.clmm

```## This file contains:
## Implementation of various methods for clmm objects.

formatVC <- function(varc, digits = max(3, getOption("digits") - 2))
### "format()" the 'VarCorr' matrix of the random effects -- for
### show()ing
### Borrowed from lme4/R/lmer.R with minor modifications.
{
recorr <- lapply(varc, attr, "correlation")
reStdDev <- lapply(varc, attr, "stddev")
reLens <- unlist(lapply(reStdDev, length))
nr <- sum(reLens)
reMat <- array('', c(nr, 4),
list(rep.int('', nr),
c("Groups", "Name", "Variance", "Std.Dev.")))
reMat[1+cumsum(reLens)-reLens, 1] <- names(reLens)
reMat[,2] <- unlist(lapply(varc, colnames))
reMat[,3] <- format(unlist(reStdDev)^2, digits = digits)
reMat[,4] <- format(unlist(reStdDev), digits = digits)
if(any(reLens > 1)) {
maxlen <- max(reLens)
corr <- do.call("rBind", lapply(recorr, function(x, maxlen)
{
if(is.null(x)) return("")
x <- as(x, "matrix")
cc <- format(round(x, 3), nsmall = 3)
cc[!lower.tri(cc)] <- ""
nr <- dim(cc)[1]
if (nr >= maxlen) return(cc)
cbind(cc, matrix("", nr, maxlen-nr))
},
maxlen))
colnames(corr) <- c("Corr", rep.int("", maxlen - 1))
cbind(reMat, corr)
} else reMat
}

varcov <-
function(object, format=FALSE,
digits=max(3, getOption("digits") - 2), ...)
### VarCorr method for model environments - should be the same for
### fitted model objects.
{
## Compute variance-covariance matrices of the random effects.
res <- lapply(object\$ST, function(st) {
## Variance-covariance matrix for the random effects:
VC <- tcrossprod(st)
## Standard deviations:
stddev <- sqrt(diag(VC))
corr <- t(VC / stddev)/stddev
attr(VC, "stddev") <- stddev
## correlation:
if(NCOL(st) > 1) {
diag(corr) <- 1
attr(VC, "correlation") <- corr
}
VC
})
names(res) <- names(object\$dims\$nlev.re)
if(format) noquote(formatVC(res, digits=digits)) else res
}

VarCorr <- function(x, ...) UseMethod("VarCorr")
VarCorr.clmm <- function(x, ...) varcov(x, ...)

print.clmm <-
function(x, digits = max(3, getOption("digits") - 3), ...)
{
if(x\$nAGQ >= 2)
else if(x\$nAGQ <= -1)
cat(paste("Cumulative Link Mixed Model fitted with the",
else
cat("Cumulative Link Mixed Model fitted with the Laplace approximation\n",
fill=TRUE)
cat("formula:", deparse(x\$formula), fill=TRUE)
if(!is.null(data.name <- x\$call\$data))
cat("data:   ", deparse(data.name), fill=TRUE)
if(!is.null(x\$call\$subset))
cat("subset: ", deparse(x\$call\$subset), fill=TRUE)
cat("\n")

print(x\$info, row.names=FALSE, right=FALSE)

cat("\nRandom effects:\n")
print(formatVC(varcov(x), digits=digits), quote=FALSE, ...)
nlev.char <- paste(names(x\$dims\$nlev.gf), " ", x\$dims\$nlev.gf, sep="", collapse=",  ")
cat("Number of groups: ", nlev.char, "\n")

if(length(x\$beta)) {
cat("\nCoefficients:\n")
print(x\$beta, digits=digits, ...)
} else {
cat("\nNo Coefficients\n")
}
if(length(x\$alpha) > 0) {
cat("\nThresholds:\n")
print(x\$alpha, digits=digits, ...)
}

if(nzchar(mess <- naprint(x\$na.action))) cat("(", mess, ")\n", sep="")
return(invisible(x))
}

vcov.clmm <- function(object, ...) vcov.clm(object, method="Cholesky")

summary.clmm <- function(object, correlation = FALSE, ...)
{
if(is.null(object\$Hessian))
stop("Model needs to be fitted with Hess = TRUE")

nfepar <- object\$dims\$nfepar
coef <- matrix(0, nfepar, 4,
dimnames = list(names(object\$coefficients[1:nfepar]),
c("Estimate", "Std. Error", "z value", "Pr(>|z|)")))
coef[, 1] <- object\$coefficients[1:nfepar]
vc <- try(vcov(object), silent = TRUE)
if(class(vc) == "try-error") {
warning("Variance-covariance matrix of the parameters is not defined")
coef[, 2:4] <- NaN
if(correlation) warning("Correlation matrix is unavailable")
object\$condHess <- NaN
}
else {
coef[, 2] <- sd <- sqrt(diag(vc)[1:nfepar])
## Cond is Inf if Hessian contains NaNs:
object\$condHess <-
if(any(is.na(object\$Hessian))) Inf
else with(eigen(object\$Hessian, only.values = TRUE),
abs(max(values) / min(values)))
coef[, 3] <- coef[, 1]/coef[, 2]
coef[, 4] <- 2 * pnorm(abs(coef[, 3]), lower.tail=FALSE)
if(correlation) ## {
## sd <- sqrt(diag(vc))
object\$correlation <- cov2cor(vc)
## (vc / sd) / rep(sd, rep(object\$edf, object\$edf))
}
object\$info\$cond.H <- formatC(object\$condHess, digits=1, format="e")
object\$coefficients <- coef
class(object) <- "summary.clmm"
return(object)
}

print.summary.clmm <-
function(x, digits = max(3, getOption("digits") - 3),
signif.stars = getOption("show.signif.stars"), ...)
{
if(x\$nAGQ >= 2)
else if(x\$nAGQ <= -1)
cat(paste("Cumulative Link Mixed Model fitted with the",
else
cat("Cumulative Link Mixed Model fitted with the Laplace approximation\n",
fill=TRUE)
cat("formula:", deparse(x\$formula), fill=TRUE)
if(!is.null(data.name <- x\$call\$data))
cat("data:   ", deparse(data.name), fill=TRUE)
if(!is.null(x\$call\$subset))
cat("subset: ", deparse(x\$call\$subset), fill=TRUE)
cat("\n")

print(x\$info, row.names=FALSE, right=FALSE)

cat("\nRandom effects:\n")
print(formatVC(varcov(x), digits=digits), quote=FALSE, ...)
nlev.char <- paste(names(x\$dims\$nlev.gf), " ", x\$dims\$nlev.gf, sep="", collapse=",  ")
cat("Number of groups: ", nlev.char, "\n")

nbeta <- length(x\$beta)
nalpha <- length(x\$alpha)
if(nbeta > 0) {
cat("\nCoefficients:\n")
printCoefmat(x\$coefficients[nalpha + 1:nbeta, , drop=FALSE],
digits=digits, signif.stars=signif.stars,
has.Pvalue=TRUE, ...)
} else {
cat("\nNo Coefficients\n")
}
if(nalpha > 0) { ## always true
cat("\nThreshold coefficients:\n")
printCoefmat(x\$coefficients[seq_len(nalpha), -4, drop=FALSE],
digits=digits, has.Pvalue=FALSE, signif.stars=FALSE,
...)
}

if(nzchar(mess <- naprint(x\$na.action))) cat("(", mess, ")\n", sep="")
if(!is.null(correl <- x\$correlation)) {
cat("\nCorrelation of Coefficients:\n")
ll <- lower.tri(correl)
correl[ll] <- format(round(correl[ll], digits))
correl[!ll] <- ""
print(correl[-1, -ncol(correl)], quote = FALSE, ...)
}
return(invisible(x))
}

## anova.clmm <- function(object, ...)
##   anova.clm(object, ...)

anova.clmm <- function(object, ...) {
### This essentially calls anova.clm(object, ...), but the names of
### the models were not displayed correctly in the printed output
### unless the following dodge is enforced.
mc <- match.call()
arg.list <- as.list(mc)
arg.list[[1]] <- NULL
return(do.call(anova.clm, arg.list))
}

logLik.clmm <- function(object, ...)
structure(object\$logLik, df = object\$edf, class = "logLik")

extractAIC.clmm <- function(fit, scale = 0, k = 2, ...) {
edf <- fit\$edf
c(edf, -2*fit\$logLik + k * edf)
}

nobs.clmm <- function(object, ...) object\$dims\$nobs

### FIXME: define edf method
edf.clmm <- function(object, ...) object\$dims\$edf

## anova.clmm <- function(object, ...)
##   anova.clm(object, ...)

anova.clmm <- function(object, ...) {
### This essentially calls anova.clm(object, ...), but the names of
### the models were not displayed correctly in the printed output
### unless the following dodge is enforced.
mc <- match.call()
arg.list <- as.list(mc)
arg.list[[1]] <- NULL
return(do.call(anova.clm, arg.list))
}

logLik.clmm <- function(object, ...)
structure(object\$logLik, df = object\$edf, class = "logLik")

extractAIC.clmm <- function(fit, scale = 0, k = 2, ...) {
edf <- fit\$edf
c(edf, -2*fit\$logLik + k * edf)
}
```

Try the ordinal package in your browser

Any scripts or data that you put into this service are public.

ordinal documentation built on May 2, 2019, 5:47 p.m.