Nothing
#######################################################################
# #
# Package: lcc #
# #
# File: lcc.R #
# Contains: generic methods #
# #
# Written by Thiago de Paula Oliveira #
# copyright (c) 2017-18, Thiago P. Oliveira #
# #
# First version: 30/10/2019 #
# Last update: 01/11/2019 #
# License: GNU General Public License version 2 (June, 1991) or later #
# #
#######################################################################
#=======================================================================
# is.lcc function
#=======================================================================
##' Reports whether x is a lcc object
##' @rdname is.lcc
##' @param x An object to test
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##' @keywords internal
##' @return Returns if the object belongs to lcc class
##' @export
is.lcc <- function(x) inherits(x, "lcc")
#=======================================================================
# Print method
#=======================================================================
##' @rdname print.lcc
##' @method print lcc
##' @title Print an \code{lcc} Object
##' @usage \method{print}{lcc}(x, digits, ...)
##' @aliases print.lcc
##' @description Prints information about the longitudinal concordance
##' correlation represented by an object of class
##' \code{\link[lcc]{lcc}}. The returned object has a
##' \code{\link[base]{print}} method.
##'
##' @return an object inheriting from class \code{print.lcc}.
##'
##' @param x an object inheriting from class
##' \code{\link[lcc]{lcc}}, representing a fitted longitudinal
##' concordance correlation function.
##'
##' @param digits a non-null value for \code{digits} specifies the minimum
##' number of significant digits to be printed in values. The
##' default, \code{NULL}.
##'
##' @param ... further arguments passed to \code{{\link{print}}}.
##' @return No return value, called for side effects
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##' @seealso \code{\link[lcc]{lcc}}, \code{\link[lcc]{summary.lcc}}
##' @examples
##' \dontrun{
##' ## Second degree polynomial model with random intercept, slope and
##' ## quadratic term
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' print(fm1)
##' }
##'
##' @importFrom stats AIC BIC
##' @export
print.lcc <- function(x, digits = NULL, ...){
message("Longitudinal concordance correlation model fit by ")
message( if(x[1]$model$method == "REML") "REML\n" else "maximum likelihood\n")
AIC <- AIC(x[1]$model)
BIC <- BIC(x[1]$model)
logLik <- c(x[1]$model$logLik)
print(data.frame(AIC, BIC, logLik, row.names = " "), digits = digits, ...)
message("\n")
print(x$Summary.lcc$fitted, digits = digits, ...)
dd <- x$model$dims
Ngrps <- dd$ngrps[1:dd$Q]
message("\n")
message("Number of Observations:", dd[["N"]])
message("\nNumber of Groups: ")
if ((lNgrps <- length(Ngrps)) == 1) { # single nesting
message(Ngrps,"\n")
} else { # multiple nesting
sNgrps <- 1:lNgrps
aux <- rep(names(Ngrps), sNgrps)
aux <- split(aux, array(rep(sNgrps, lNgrps),
c(lNgrps, lNgrps))[!lower.tri(diag(lNgrps))])
names(Ngrps) <- unlist(lapply(aux, paste, collapse = " %in% "))
message("\n")
print(rev(Ngrps), digits = digits, ...)
}
invisible(x)
}
#=======================================================================
# fitted method
#=======================================================================
##' @rdname fitted.lcc
##' @method fitted lcc
##' @title Extract \code{lcc} Fitted Values
##' @usage
##' \method{fitted}{lcc}(object, type, digits, ...)
##' @aliases fitted.lcc
##' @description Fitted values from object of class \code{lcc}
##' returned by modeling functions.
##'
##' @return A data frame with columns given by methods, time, and
##' fitted values.
##'
##' @param object an object inheriting from class \code{lcc},
##' representing a fitted longitudinal concordance correlation
##' function.
##'
##' @param type an optional character string specifying the type of
##' output to be returned. If \code{type="lcc"}, prints the fitted
##' longitudinal concordance correlation. If
##' \code{type="lpc"}, prints the fitted longitudinal Pearson
##' correlation. If \code{type="la"}, prints the fitted longitudinal
##' accuracy. Defaults to \code{type="lcc"}.
##'
##' @param digits a non-null value for \code{digits} specifies the minimum
##' number of significant digits to be printed in values. The
##' default, \code{NULL}.
##'
##' @param ... not used.
##'
##' @return No return value, called for side effects
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##' @seealso \code{\link[lcc]{lcc}}, \code{\link[lcc]{summary.lcc}},
##' \code{\link[lcc]{lccPlot}}
##' @examples
##' data(hue)
##' ## Second degree polynomial model with random intercept, slope and
##' ## quadratic term
##' \dontrun{
##' fm1 <- lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2,
##' components = TRUE)
##' fitted(fm1)
##' fitted(fm1, type="lpc")
##' fitted(fm1, type="la")
##' }
##'
##' @importFrom stats AIC BIC
##' @export
fitted.lcc <- function(object, type = "lcc", digits = NULL, ...){
if (!inherits(object, "lcc")) stop("Object must inherit from class \"lcc\"",
call.=FALSE)
if(missing(type)) type="lcc"
if(object$plot_info$components == FALSE && type == "lpc" ||
object$plot_info$components == FALSE && type == "la"){
stop(paste0("It is necessary to include components = TRUE in the
lcc() function to calculate the fitted values for type ", type))
}
if(type == "lcc" || type == "lpc" || type == "la"){
pr <- switch(type,
"lcc" = cat( "Fitted longitudinal concordance correlation function", "\n"),
"lpc" = cat( "Fitted longitudinal Pearson correlation function", "\n"),
"la" = cat( "Fitted longitudinal accuracy function", "\n"))
message("\n")
fittedBuilder(object = object, type = type)
}else{
stop("Available 'type' are 'lcc', 'lpc', or 'la'", call.=FALSE)
}
}
#=======================================================================
# Summary method
#=======================================================================
#-----------------------------------------------------------------------
# Print method for summary.lcc
#-----------------------------------------------------------------------
##' @rdname print.summary.lcc
##' @title Print the Summary of an \code{lcc} Object
##' @usage
##' \method{print}{summary.lcc}(x, verbose, digits, ...)
##' @method print summary.lcc
##' @aliases print.summary.lcc
##' @description Information summarizing the fitted longitudinal
##' concordance correlation is printed. This includes the AIC, BIC,
##' and log-likelihood at convergence. If \code{type = "lcc"}, prints
##' the fitted values while \code{type = "model"} prints the fixed
##' effects estimates and their standard errors, standard deviations,
##' correlations for the random effects, within-group correlation, and
##' variance function parameters.
##'
##' @param x an object inheriting from class
##' \code{\link[lcc]{summary.lcc}}, representing a fitted longitudinal
##' concordance correlation function.
##'
##' @param verbose an optional logical value used to control the amount
##' of printed output when \code{type = "model"}. Defaults to
##' \code{FALSE}
##'
##' @param digits a non-null value for \code{digits} specifies the
##' minimum number of significant digits to be printed in values. The
##' default, \code{NULL}.
##'
##' @param ... further arguments passed to \code{\link{print}}.
##'
##' @importFrom stats AIC BIC asOneSidedFormula
##'
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##' @return No return value, called for side effects
##' @seealso \code{\link{summary.lcc}}, \code{\link{lccPlot}},
##' \code{\link[lcc]{lcc}}
##'
##' @examples
##' \dontrun{
##' ## Second degree polynomial model with random intercept, slope and
##' ## quadratic term
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' print(summary(fm1, type="model"))
##' }
##'
##' @export
print.summary.lcc <- function(x, verbose = FALSE, digits = NULL, ...){
if (inherits(x, "lcc")) {
message("Longitudinal concordance correlation model fit by ")
message( if(x$model$method == "REML") "REML\n" else "maximum likelihood\n")
AIC <- AIC(x$model)
BIC <- BIC(x$model)
logLik <- c(x$model$logLik)
print(data.frame(AIC, BIC, logLik, row.names = " "), digits = digits, ...)
message("\n")
gof <- x$gof
message(paste0(" gof: ", round(gof, 4)), "\n")
message("\n")
if (inherits(x$comp, "character")) {
if (is.null(x$info$ENV.LCC)) {
message(x$comp, "\n")
fitted <- x$fitted
print(fitted, digits = digits, ...)
}else{
message(paste0(" Lower and upper bound of ", (1-x$plot_info$alpha)*100,"%"), "bootstrap confidence interval", "\n")
message(" Number of bootstrap samples: ", x$plot_info$nboot, "\n")
message("\n")
message(x$comp, "\n")
fitted <- x$fitted
print(fitted, digits = digits, ...)
}
}else{
summ <- sum(sapply(x$comp, length))
if(is.null(x$info$ENV.LCC)){
for(i in 1:summ){
message(x$comp[[i]], "\n")
fitted <- x$fitted
print(fitted[[i]], digits = digits, ...)
message("\n")
}
}else{
message(paste0(" Lower and upper bound of ", (1-x$info$alpha)*100,"%"), "bootstrap confidence interval", "\n")
message(" Number of bootstrap samples: ", x$info$nboot, "\n")
message("\n")
for(i in 1:summ){
message(x$comp[[i]], ": LCC", "\n")
fitted <- x$fitted
print(fitted$LCC[[i]], digits = digits, ...)
message("\n")
message(x$comp[[i]], ": LPC", "\n")
print(fitted$LPC[[i]], digits = digits, ...)
message("\n")
message(x$comp[[i]], ": LA", "\n")
print(fitted$LA[[i]], digits = digits, ...)
message("\n", "\n")
}
}
}
}else{
if (inherits(x, "model")) {
dd <- x$dims
verbose <- verbose || attr(x, "verbose")
message( "Linear mixed-effects model fit by " )
message( if(x$method == "REML") "REML\n" else "maximum likelihood\n")
## method <- x$method
message(" Data:", deparse( x$call$data ), "\n")
if (!is.null(x$call$subset)) {
message(" Subset:", deparse(asOneSidedFormula(x$call$subset)[[2L]]),"\n")
}
print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = c(x$logLik),
row.names = " "), ...)
if (verbose) { message("Convergence at iteration:",x$numIter,"\n") }
message("\n")
print(summary(x$modelStruct), sigma = x$sigma,
reEstimates = x$coef$random, verbose = verbose, ...)
message("Fixed effects: ")
fixF <- x$call$fixed
if (inherits(fixF, "formula") || is.call(fixF)) {
message(deparse(x$call$fixed), "\n")
} else {
message(deparse(lapply(fixF, function(el) as.name(deparse(el)))), "\n")
}
## fixed effects t-table and correlations
xtTab <- as.data.frame(x$tTable)
wchPval <- match("p-value", names(xtTab))
for(i in names(xtTab)[-wchPval]) {
xtTab[, i] <- format(zapsmall(xtTab[, i]))
}
xtTab[,wchPval] <- format(round(xtTab[,wchPval], 4))
if (any(wchLv <- (as.double(levels(xtTab[, wchPval])) == 0))) {
levels(xtTab[, wchPval])[wchLv] <- "<.0001"
}
row.names(xtTab) <- dimnames(x$tTable)[[1L]]
print(xtTab, ...)
if (nrow(x$tTable) > 1) {
corr <- x$corFixed
class(corr) <- "correlation"
print(corr, title = " Correlation:", ...)
}
message("\nStandardized Within-Group Residuals:\n")
print(x$residuals, ...)
message("\nNumber of Observations:",x$dims[["N"]])
message("\nNumber of Groups: ")
Ngrps <- dd$ngrps[1:dd$Q]
if ((lNgrps <- length(Ngrps)) == 1) { # single nesting
message(Ngrps,"\n")
} else { # multiple nesting
sNgrps <- 1:lNgrps
aux <- rep(names(Ngrps), sNgrps)
aux <- split(aux, array(rep(sNgrps, lNgrps),
c(lNgrps, lNgrps))[!lower.tri(diag(lNgrps))])
names(Ngrps) <- unlist(lapply(aux, paste, collapse = " %in% "))
message("\n")
print(rev(Ngrps), ...)
}
invisible(x)
}else{
stop("Available only for classes summary.lcc or summary.lme", call.=FALSE)
}
}
}
#-----------------------------------------------------------------------
# summary.lcc
#-----------------------------------------------------------------------
##' @rdname summary.lcc
##' @title Summarize an \code{lcc} Object
##' @usage
##' \method{summary}{lcc}(object, type, adjustSigma, verbose, ...)
##'
##' @method summary lcc
##' @aliases summary.lcc
##' @description Additional information about the fit of longitudinal
##' concordance correlation, longitudinal Pearson correlation, and
##' longitudinal accuracy represented by an object of class
##' \code{\link[lcc]{lcc}}. The returned object has a
##' \code{\link[base]{print}} method.
##'
##' @return an object inheriting from class \code{summary.lcc}
##' including: \item{fitted}{the fitted values extracted from the
##' \code{lcc} object.} \item{gof}{the goodness of fit (gof) measurement
##' is calculated using the concordance correlation coefficient between
##' fitted and observed values. Value of 1 denote perfect concordance.}
##' \item{AIC}{the Akaike Information Criterion corresponding to object.}
##' \item{BIC}{the Bayesian Information Criterion corresponding to object.}
##' \item{logLik}{If \code{REML=FALSE}, returns the log-likelihood value
##' of the linear mixed-effects model; otherwise, the restricted
##' log-likelihood is returned}
##'
##' @param object an object inheriting from class
##' \code{\link[lcc]{lcc}}, representing a fitted longitudinal
##' concordance correlation function.
##'
##' @param type an optional character string specifying the type of
##' output to be returned. If \code{type="model"}, prints the summary
##' of the polynomial mixed-effects regression model. If
##' \code{type="lcc"}, prints the summary of the fitted and sampled
##' values for LCC, LPC, and LA as well as the concordance correlation
##' coefficient between fitted values from the model and observed
##' values as goodness of fit (gof) measurement. Defaults to
##' \code{type="model"}.
##'
##' @param adjustSigma an optional logical value used when \code{type =
##' model}. If TRUE and the estimation method used to obtain object
##' was maximum likelihood, the residual standard error is multiplied
##' by \code{sqrt(nobs/(nobs - npar))}. See
##' \code{\link[nlme]{summary.lme}} for more information. Default is
##' TRUE.
##'
##' @param verbose an optional logical value used to control the amount
##' of output in the \code{print.summary.lme} method when
##' \code{type = model} is used. Defaults to FALSE.
##'
##' @param ... not used.
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @importFrom stats AIC BIC asOneSidedFormula pt resid
##'
##' @seealso \code{\link{AIC}}, \code{\link{BIC}},
##' \code{print.summary.lcc}, \code{\link[lcc]{lcc}}
##'
##' @examples
##'
##' ## Second degree polynomial model with random intercept, slope and
##' ## quadratic term
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' summary(fm1, type="model")
##' summary(fm1, type="lcc")
##' @export
summary.lcc <- function(object, type, adjustSigma = TRUE,
verbose = FALSE, ...)
{
if (!inherits(object, "lcc")) stop("Object must inherit from class \"lcc\"",
call.=FALSE)
if(missing(type)) type <- "model"
if(type=="model" || type=="lcc"){
if(type == "lcc"){
# Object lcc
object$fitted <- object$Summary.lcc$fitted
object$sampled <- object$Summary.lcc$sampled
object$gof <- object$Summary.lcc$gof
object$data <- object$data
object$AIC <- AIC(object[1]$model)
object$BIC <- BIC(object[1]$model)
object$logLik <- c(object[1]$model$logLik)
object$info <- object$plot_info
object$comp <- object$Summary.lcc$comp
object <- object[names(object) != "Summary.lcc"]
object <- object[names(object) != "plot_info"]
object <- object[names(object) != "data"]
#-----------------------------------------------------------------
## generating the final object
#-----------------------------------------------------------------
structure(object, type = "lcc", oClass = class(object),
class = c("summary.lcc", type))
}else {
#-------------------------------------------------------------------
# Model
#-------------------------------------------------------------------
obj <- object[1]$model
## variance-covariance estimates for fixed effects
fixed <- fixef(obj)
stdFixed <- sqrt(diag(as.matrix(obj$varFix)))
obj$corFixed <- array(t(obj$varFix/stdFixed)/stdFixed,
dim(obj$varFix), list(names(fixed),names(fixed)))
if (adjustSigma && obj$method == "ML")
stdFixed <- stdFixed *
sqrt(obj$dims$N/(obj$dims$N - length(stdFixed)))
## fixed effects coefficients, std. deviations and t-ratios
fDF <- obj$fixDF[["X"]]
tVal <- fixed/stdFixed
obj$tTable <- cbind(Value = fixed, Std.Error = stdFixed, DF = fDF,
"t-value" = tVal, "p-value" = 2 * pt(-abs(tVal), fDF))
## residuals
resd <- resid(obj, type = "pearson")
if (length(resd) > 5) {
resd <- quantile(resd, na.rm = TRUE) # might have NAs from na.exclude
names(resd) <- c("Min","Q1","Med","Q3","Max")
}
obj$residuals <- resd
## generating the final obj
aux <- logLik(obj)
obj$BIC <- BIC(aux)
obj$AIC <- AIC(aux)
structure(obj, verbose = verbose, oClass = class(obj),
class = c("summary.lcc", "model", class(obj)))
}
}else {
stop("Available 'type' are lcc or model", call.=FALSE)
}
}
#=======================================================================
# Plot method
#=======================================================================
##' @rdname plot.lcc
##' @method plot lcc
##' @title Diagnostic Plots of an \code{lcc} Object.
##'
##' @usage
##' \method{plot}{lcc}(x, which = c(1L:6L),
##' caption = list("Residuals vs Fitted",
##' "Residuals vs Time",
##' "Residuals by Subject",
##' "Observed values vs Fitted values",
##' "Normal Q-Q Plot (Conditional residuals)",
##' "Normal Q-Q Plot (Random effects)"),
##' sub.caption = NULL, main = NULL,
##' panel = if(add.smooth) panel.smooth else points,
##' add.smooth = TRUE, ask = TRUE,
##' id.n = 3, labels.id = names(residuals(x)),
##' label.pos = c(4, 2), cex.id = 0.75, cex.caption = 1,
##' cex.oma.man = 1.25, ...)
##'
##' @aliases plot.lcc
##'
##' @description Diagnostic plots for conditional error and random
##' effects from the linear mixed-effects fit are obtained. Six plots
##' plots (selectable by 'which') are currently available: a plot of
##' residuals against fitted values, a plot of residuals against
##' time variable, a boxplot of residuals by subject,
##' a plot of observerd values against fitted values, a normal Q-Q plot
##' with simulation envelopes based on conditional error, and a normal
##' Q-Q plot with simulation envelopes based on the random effects. By
##' default, all plots are provided.
##'
##' @param x an object inheriting from class \code{\link[lcc]{lcc}},
##' representing a fitted longitudinal concordance correlation
##' function.
##' @param which if a subset of the plots is required, specify a subset
##' of the numbers from 1 to 6.
##' @param caption captions to appear above the plots. Vector or list of
##' valid graphics annotations is required. All captions can be
##' supressed using '""' or \code{NA}.
##' @param sub.caption common sub-title (at bottom). Default to
##' \code{NULL}.
##' @param main The main title (on top) above the caption.
##' @param panel panel function. If \code{add.smooth = TRUE},
##' \code{panel.smooth} is used rather than \code{points}.
##' @param add.smooth logical indicating if smoother should be added to
##' most plots; see also \code{panel} above. Defaults to \code{TRUE}.
##' @param ask logical; if \code{TRUE}, the default, the user is _ask_ed
##' before each plot, see \code{\link[graphics]{par}}.
##' @param id.n number of points to be labelled is the first three
##' plots, starting with the most extreme.
##' @param labels.id vector of labels, from which the labels for extreme
##' points will be chosen. Default to \code{NULL} (uses observation
##' numbers).
##' @param label.pos positioning of labels, for the left half and right
##' half of the graph respectively, for plots 1-3.
##' @param cex.id magnification of point label.
##' @param cex.caption controls the size of \code{caption}.
##' @param cex.oma.man controls the size of the \code{sub.caption} only
##' if that is _above_ the figures when there is more than one.
##' @param ... further graphical parameters from 'par'.
##'
##' @details The Q-Q plot uses the normalized residuals. The
##' standardized residuals is pre-multiplied by the inverse
##' square-root factor of the estimated error correlation matrix while
##' the random effects is pre-multiplied by the inverse square root of
##' the estimated variances obtained from matrix G. The simulate
##' envelopes are obtained from package hnp (Moral et al., 2018).
##'
##' Code partially adapted from \code{\link[stats]{plot.lm}}.
##' @return Return plots for conditional error and random effects from the linear mixed-effects
##' @importFrom hnp hnp
##'
##' @importFrom nlme getVarCov ranef
##'
##' @importFrom grDevices as.graphicsAnnot dev.flush dev.hold dev.interactive devAskNewPage extendrange n2mfrow
##'
##' @importFrom graphics abline boxplot mtext panel.smooth par plot points strheight text title
##'
##' @importFrom stats fitted residuals
##'
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @seealso \code{\link{lccPlot}}, \code{\link[lcc]{lcc}},
##' \code{mtext}, \code{text}, \code{plotmath}
##'
##' @examples
##'
##' ## Second degree polynomial model with random intercept, slope and
##' ## quadratic term
##' fm1 <- lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' plot(fm1)
##' @export
plot.lcc <- function(x, which = c(1L:6L),
caption = list("Residuals vs Fitted",
"Residuals vs Time",
"Residuals by Subject",
"Observed values vs Fitted values",
"Normal Q-Q Plot (Conditional residuals)",
"Normal Q-Q Plot (Random effects)"),
sub.caption = NULL, main = NULL,
panel = if(add.smooth) panel.smooth else points,
add.smooth = TRUE, ask = TRUE,
id.n = 3, labels.id = names(residuals(x)),
label.pos = c(4, 2), cex.id = 0.75, cex.caption = 1,
cex.oma.man = 1.25, ...)
{
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
if (!is.lcc(x))
stop("use only with \"lcc\" objects", call. = FALSE)
if(!is.numeric(which) || any(which < 1) || any(which > 6))
stop("'which' must be in 1:6")
#-------------------------------------------------------------------
show <- rep(FALSE, 6)
show[which] <- TRUE
#-------------------------------------------------------------------
# information from lme model
#-------------------------------------------------------------------
model <- x$model
r <- residuals(model)
r_norm <- residuals(model, type = "normalized")
yh <- fitted(model)
n <- length(r)
time <- model$data$time
#-------------------------------------------------------------------
if(id.n > 0L) { ## label the largest residuals
if(is.null(labels.id))
labels.id <- paste(1L:n)
iid <- 1L:id.n
show.r <- sort.list(abs(r), decreasing = TRUE)[iid]
text.id <- function(x, y, ind, adj.x = TRUE) {
if (is.factor(x)) {
labpos <-
if(adj.x) label.pos[1+as.numeric(y > mean(range(y)))] else 3
}else {
labpos <-
if(adj.x) label.pos[1+as.numeric(x > mean(range(x)))] else 3
}
text(x, y, labels.id[ind], cex = cex.id, xpd = TRUE,
pos = labpos, offset = 0.25)
}
}
#-------------------------------------------------------------------
getCaption <- function(k) # allow caption = "" , plotmath etc
if(length(caption) < k) NA_character_ else as.graphicsAnnot(caption[[k]])
one.fig <- prod(par("mfcol")) == 1
if (ask) {
oask <- devAskNewPage(TRUE)
on.exit(devAskNewPage(oask))
}
#-------------------------------------------------------------------
# Individual plots
#-------------------------------------------------------------------
if (show[1L]) {
l.fit <- "Fitted values"
ylim <- range(r, na.rm=TRUE)
if(id.n > 0)
ylim <- extendrange(r = ylim, f = 0.08)
dev.hold()
plot(yh, r, xlab = l.fit, ylab = "Residuals", main = main,
ylim = ylim, type = "n", ...)
panel(yh, r, ...)
if (one.fig)
title(sub = sub.caption, ...)
mtext(getCaption(1), 3, 0.25, cex = cex.caption)
if(id.n > 0) {
y.id <- r[show.r]
y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
text.id(yh[show.r], y.id, show.r)
}
abline(h = 0, lty = 3, col = "gray")
dev.flush()
}
#===================================================================
if (show[2L]) {
l.fit <- "Time"
ylim <- range(r, na.rm=TRUE)
if(id.n > 0)
ylim <- extendrange(r = ylim, f = 0.08)
dev.hold()
plot(time, r, xlab = l.fit, ylab = "Residuals", main = main,
ylim = ylim, type = "n", ...)
panel(time, r, ...)
if (one.fig)
title(sub = sub.caption, ...)
mtext(getCaption(2), 3, 0.25, cex = cex.caption)
if(id.n > 0) {
y.id <- r[show.r]
y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
text.id(time[show.r], y.id, show.r)
}
abline(h = 0, lty = 3, col = "gray")
dev.flush()
}
#===================================================================
if (show[3L]) {
Subject <- model$data$subject
boxplot(r ~ Subject, ylab = "Residuals", main = main, ...)
if (one.fig)
title(sub = sub.caption, ...)
mtext(getCaption(3), 3, 0.25, cex = cex.caption)
if(id.n > 0) {
y.id <- r[show.r]
y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
text.id(Subject[show.r], y.id, show.r)
}
abline(h = 0, lty = 3, col = "blue")
}
#===================================================================
if (show[4L]) {
Response <- model$data$resp
plot(Response ~ yh, ylab = "Observed Values",
xlab = "Fitted Values", main = main, ...)
if (one.fig)
title(sub = sub.caption, ...)
mtext(getCaption(4), 3, 0.25, cex = cex.caption)
abline(0, 1)
}
#===================================================================
if (show[5L]) {
hnp(r_norm, scale = TRUE, halfnormal = FALSE, print.on = TRUE,
main = main, ...)
if (one.fig)
title(sub = sub.caption, ...)
mtext(getCaption(5), 3, 0.25, cex = cex.caption)
}
#===================================================================
if (show[6L]) {
vars <- sqrt(diag(getVarCov(model)))
ranefs <- re <- as.matrix(ranef(model))
re <- ranefs %*% diag(1 / vars)
ncol.re <- ncol(re)
if (ncol.re == 1) {
if (is.null(main)) {
main <- "Random effect (Intercept)"
}
hnp(re, scale = TRUE, halfnormal = FALSE, print.on = TRUE,
main = main, ...)
}else {
par(mfrow = rev(n2mfrow(ncol.re)))
for (i in 1:ncol.re) {
hnp(re[, i], scale = TRUE, halfnormal = FALSE, print.on = TRUE,
main = main, ...)
if (one.fig)
title(sub = sub.caption, ...)
mtext(getCaption(6), 3, 0.25, cex = cex.caption)
mtext(paste0("b", i - 1, "i"), 1, -1.5, cex = cex.caption)
}
}
}
par(mfrow = c(1,1))
invisible()
}
#=======================================================================
# coef function
# =======================================================================
##' @rdname coef.lcc
##' @title Extract Model Coefficients
##' @usage \method{coef}{lcc}(object, ...)
##' @method coef lcc
##' @aliases coef.lcc
##'
##' @description The fixed effects estimated and corresponding random
##' effects estimates are obtained at subject levels less or equal to
##' i. The resulting estimates are returned as a data frame, with rows
##' corresponding to subject levels and columns to coefficients.
##'
##' @param object an object inheriting from class \code{lcc},
##' representing a fitted longitudinal concordance correlation
##' function.
##' @param ... optional arguments passed to the \code{coef.lme}
##' function.
##'
##' @details See methods for \code{\link{nlme}} objects to get more
##' details.
##' @return Coefficients extracted from the model object.
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @seealso \code{\link[lcc]{lcc}}, \code{\link{summary.lcc}},
##' \code{\link{lccPlot}}, \code{\link{vcov.lcc}}
##'
##' @importFrom stats coef
##'
##' @examples
##'
##' \dontrun{
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' coef(fm1)
##' }
##'
##' @export
coef.lcc <- function(object, ...) {
if (!is.lcc(object))
stop("use only with \"lcc\" objects", call. = FALSE)
x <- coef(object$model)
colnames(x) <- gsub(pattern = "fixed", x = colnames(x),
replacement = "Fixed")
colnames(x) <- gsub(pattern = "Poly", x = colnames(x),
replacement = "Time")
colnames(x) <- gsub(pattern = "fmla.", x = colnames(x),
replacement = "")
colnames(x) <- gsub(pattern = "\\(time, degree = qr, raw = TRUE\\)",
x = colnames(x), replacement = "")
colnames(x) <- gsub(pattern = "rand", x = colnames(x),
replacement = "Random")
colnames(x) <- gsub(pattern = "poly", x = colnames(x),
replacement = "Time")
colnames(x) <- gsub(pattern = "method", x = colnames(x),
replacement = "")
class(x) <- c("coef.lcc", "ranef.lcc", "data.frame")
x
}
#=======================================================================
# vcov function
#=======================================================================
##' @rdname vcov.lcc
##' @title Extract Variance-Covariance Matrix of the Fixed Effects
##' @usage \method{vcov}{lcc}(object, ...)
##' @method vcov lcc
##' @aliases vcov.lcc
##'
##' @return Returns the variance-covariance matrix of a fitted
##' \code{lcc} model object.
##'
##' @param object an object inheriting from class \code{lcc},
##' representing a fitted longitudinal concordance correlation
##' function.
##' @param ... optional arguments passed to the \code{vcov.lme}
##' function.
##'
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @details See methods for \code{\link{nlme}} objects to get more
##' details.
##'
##' @seealso \code{\link{summary.lcc}}, \code{\link{lccPlot}},
##' \code{\link[lcc]{lcc}}, \code{\link{coef.lcc}}
##'
##' @importFrom stats vcov
##'
##' @examples
##' \dontrun{
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' vcov(fm1)
##' }
##'
##' @export
vcov.lcc <- function(object, ...) {
if (!is.lcc(object))
stop("use only with \"lcc\" objects", call. = FALSE)
x <- vcov(object$model, ...)
rownames(x) <- colnames(x) <-
gsub(pattern = "fixed", x = colnames(x),
replacement = "Fixed")
rownames(x) <- colnames(x) <-
gsub(pattern = "Poly", x = colnames(x),
replacement = "Time")
rownames(x) <- colnames(x) <-
gsub(pattern = "method", x = colnames(x),
replacement = "")
x
}
#=======================================================================
# getVarCov function
#=======================================================================
##' @rdname getVarCov.lcc
##' @title Extract Variance Components from a Fitted Model
##' @usage \method{getVarCov}{lcc}(obj, type, ...)
##' @method getVarCov lcc
##' @aliases getVarCov.lcc
##'
##' @return Returns the variance-covariance matrix of a fitted
##' \code{lcc} model object.
##'
##' @param obj an object inheriting from class \code{lcc}, representing
##' a fitted longitudinal concordance correlation function.
##' @param type specifies the type of variance covariance matrix. If
##' \code{type = "random.effects"}, the default, extract the
##' random-effects variance-covariance; if \code{type = "conditional"}
##' extract the conditional variance-covariance of the responses; and
##' if \code{type = "marginal"} extracts the the marginal
##' variance-covariance of the responses.
##' @param ... optional arguments passed to the \code{getVarCov}
##' function.
##'
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @seealso \code{\link[lcc]{lcc}}, \code{\link{summary.lcc}},
##' \code{\link{coef.lcc}}, \code{\link{vcov.lcc}}
##'
##' @details See methods for \code{\link{nlme}} objects to get more
##' details.
##'
##' @importFrom nlme getVarCov
##'
##' @examples
##'
##' \dontrun{
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' getVarCov(fm1)
##' }
##'
##' @export
getVarCov.lcc <- function(obj, type = "random.effects", ...) {
if (!is.lcc(obj))
stop("use only with \"lcc\" objects", call. = FALSE)
x <- getVarCov(obj$model, type = type, ...)
rownames(x) <- colnames(x) <-
gsub(pattern = "fmla.", x = colnames(x), replacement = "")
rownames(x) <- colnames(x) <-
gsub(pattern = "\\(time, degree = qr, raw = TRUE\\)",
x = colnames(x), replacement = "")
rownames(x) <- colnames(x) <-
gsub(pattern = "rand", x = colnames(x), replacement = "Random")
rownames(x) <- colnames(x) <-
gsub(pattern = "poly", x = colnames(x),
replacement = "Time")
x
}
#=======================================================================
# Residuals
#=======================================================================
##' @rdname residuals.lcc
##' @title Extract Model Residuals
##' @usage \method{residuals}{lcc}(object, type, ...)
##' @method residuals lcc
##' @aliases residuals.lcc
##'
##' @description Extract the residulas from the model used to estimate
##' the longitudinal concordance correlation function.
##'
##' @param object an object inheriting from class \code{lcc},
##' representing a fitted longitudinal concordance correlation
##' function.
##'
##' @param type an optional character string specifying the type of
##' residulas to be used. If \code{type = "response"}, the default,
##' the residuals at level i are obtained by subtracting the fitted
##' values at that level from the response vector. If \code{type =
##' "pearson"}, the "response" residuals is divided by the estimated
##' within-group standard error. If \code{type = "normalized"}, the
##' normalized residuals are used. Partial matching of arguments is
##' used, so only the first character needs to be provided.
##'
##' @param ... optional arguments passed to the \code{residuals.lme}
##' function.
##'
##' @details See methods for \code{\link{nlme}} objects to get more
##' details.
##' @return Return no value, called for side effects
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @seealso \code{\link[lcc]{lcc}}, \code{\link{summary.lcc}},
##' \code{\link{coef.lcc}}, \code{\link{vcov.lcc}}
##'
##' @examples
##'
##' \dontrun{
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' getVarCov(fm1)
##' }
##'
##' @export
residuals.lcc <- function(object, type = "response", ...) {
if (!is.lcc(object))
stop("use only with \"lcc\" objects", call. = FALSE)
residuals(object$model, type = type, ...)
}
#=======================================================================
# AIC and BIC
#=======================================================================
##' @rdname AIC.lcc
##' @title Akaike and Bayesian Information Criteria for an \code{lcc} Object.
##' @method AIC lcc
##' @aliases AIC.lcc
##' @description Calculate the Akaike's 'An Information Criterion' or
##' the BIC or SBC (Schwarz's Bayesian criterion) for an object of
##' class \code{lcc}.
##'
##' @param object an object inheriting from class \code{lcc},
##' representing a fitted longitudinal concordance correlation
##' function.
##'
##' @param k numeric value, use as penalty coefficient for the number of
##' parameters in the fitted model; the default \code{k = 2} is the
##' classical AIC.
##'
##' @param ... optional arguments passed to the \code{AIC}
##' function.
##'
##' @return A numeric value with the corresponding AIC or BIC
##' value. See methods for \code{\link{AIC}} objects to get more
##' details.
##'
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @seealso \code{\link[lcc]{lcc}}, \code{\link{summary.lcc}},
##' \code{\link{coef.lcc}}, \code{\link{vcov.lcc}}
##'
##' @importFrom stats na.omit nobs
##'
##' @export
AIC.lcc <- function(object, ..., k = 2) {
if (!is.lcc(object))
stop("use only with \"lcc\" objects", call. = FALSE)
ll <- logLik
if (!missing(...)) {
lls <- lapply(list(object$model, ...), ll)
vals <- sapply(lls, function(el) {
no <- attr(el, "nobs")
c(as.numeric(el), attr(el, "df"), if (is.null(no)) NA_integer_ else no)
})
val <- data.frame(df = vals[2L, ], ll = vals[1L, ])
nos <- na.omit(vals[3L, ])
if (length(nos) && any(nos != nos[1L]))
warning("models are not all fitted to the same number of observations")
val <- data.frame(df = val$df, AIC = -2 * val$ll + k *
val$df)
Call <- match.call()
Call$k <- NULL
row.names(val) <- as.character(Call[-1L])
val
}
else {
AIC(object$model)
}
}
##' @rdname AIC.lcc
##' @method BIC lcc
##' @aliases BIC.lcc
##'
##' @importFrom stats na.omit nobs
##'
##' @examples
##' \dontrun{
##' attach(simulated_hue)
##' fm6 <- lcc(data = simulated_hue, subject = "Fruit",
##' resp = "Hue", method = "Method", time = "Time",
##' qf = 2, qr = 1, components = TRUE,
##' time_lcc = list(n=50, from=min(Time), to=max(Time)))
##' AIC(fm6)
##' BIC(fm6)
##' }
##'
##' @export
BIC.lcc <- function (object, ...)
{
if (!is.lcc(object))
stop("use only with \"lcc\" objects", call. = FALSE)
ll <- logLik
Nobs <- nobs
if (!missing(...)) {
lls <- lapply(list(object$model, ...), ll)
vals <- sapply(lls, function(el) {
no <- attr(el, "nobs")
c(as.numeric(el), attr(el, "df"), if (is.null(no)) NA_integer_ else no)
})
val <- data.frame(df = vals[2L, ], ll = vals[1L, ],
nobs = vals[3L, ])
nos <- na.omit(val$nobs)
if (length(nos) && any(nos != nos[1L]))
warning("models are not all fitted to the same number of observations")
unknown <- is.na(val$nobs)
if (any(unknown))
val$nobs[unknown] <- sapply(list(object$model, ...)[unknown],
function(x) tryCatch(Nobs(x), error = function(e) NA_real_))
val <- data.frame(df = val$df, BIC = -2 * val$ll + log(val$nobs) *
val$df)
row.names(val) <- as.character(match.call()[-1L])
val
}
else {
BIC(object$model)
}
}
#=======================================================================
# ranef
# =======================================================================
##' @rdname ranef.lcc
##' @title Extract Model Random Effects
##' @usage \method{ranef}{lcc}(object, ...)
##' @method ranef lcc
##' @aliases ranef.lcc
##'
##' @description Extract the estimated random effects at level i.
##'
##' @return A data frame with rows given by the different groups at that level and
##' columns given by the random effects.
##'
##' @param object an object inheriting from class \code{lcc},
##' representing a fitted longitudinal concordance correlation
##' function.
##' @param ... optional arguments passed to the \code{ranef.lme}
##' function.
##'
##' @details See methods for \code{\link{nlme}} objects to get more
##' details.
##'
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @seealso \code{\link[lcc]{lcc}}, \code{\link{coef.lcc}},
##'
##' @importFrom nlme ranef
##'
##' @examples
##' \dontrun{
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' ranef(fm1)
##' }
##' @export
ranef.lcc <- function(object, ...) {
if (!is.lcc(object))
stop("use only with \"lcc\" objects" , call. = FALSE)
x <- ranef(object$model)
colnames(x) <- gsub(pattern = "fmla.", x = colnames(x),
replacement = "")
colnames(x) <- gsub(pattern = "\\(time, degree = qr, raw = TRUE\\)",
x = colnames(x), replacement = "")
colnames(x) <- gsub(pattern = "rand", x = colnames(x),
replacement = "Random")
colnames(x) <- gsub(pattern = "poly", x = colnames(x),
replacement = "Time")
class(x) <- c("ranef.lcc", "data.frame")
x
}
#=======================================================================
# logLik
#=======================================================================
##' @rdname logLik.lcc
##' @title Extract Log-Likelihood of an \code{lcc} Object
##' @usage \method{logLik}{lcc}(object, ..., REML)
##' @method logLik lcc
##' @aliases logLik.lcc
##'
##' @return If \code{REML=TRUE}, the default, returns the
##' restricted log-likelihood value of the linear mixed-effects model;
##' else the log-likelihood value
##'
##' @param object an object inheriting from class \code{lcc},
##' representing a fitted longitudinal concordance correlation
##' function.
##'
##' @param REML an optional logical value. If \code{TRUE} the
##' restricted log-likelihood is returned, else, if \code{FALSE}, the
##' log-likelihood is returned.
##' @param ... further arguments passed to \code{\link{logLik}}.
##'
##' @details See methods for \code{\link{nlme}} objects to get more
##' details.
##'
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @importFrom stats logLik
##'
##' @seealso \code{\link[lcc]{lcc}}, \code{\link{summary.lcc}}
##'
##' @examples
##'
##' \dontrun{
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' logLik(fm1)
##' }
##'
##' @export
logLik.lcc <- function(object, ..., REML) {
if (!is.lcc(object))
stop("use only with \"lcc\" objects" , call. = FALSE)
logLik(object$model, REML = REML, ...)
}
#=======================================================================
# ANOVA
# =======================================================================
##' @rdname anova.lcc
##' @title Compare Likelihoods of Fitted Models from an \code{lcc}
##' Object
##' @usage \method{anova}{lcc}(object, ..., test, type, adjustSigma,
##' verbose)
##' @method anova lcc
##' @aliases anova.lcc
##'
##' @return If just one \code{lcc} model object is declared, a data
##' frame with the numerator degrees of freedom, denominator degrees
##' of freedom, F-values, and P-values for the fixed terms in the
##' model. Otherwise, when multiple \code{lcc} fitted objects are
##' being compared, a data frame with the degrees of freedom, the
##' (restricted) log-likelihood, the Akaike Information Criterion
##' (AIC), and the Bayesian Information Criterion (BIC) of each object
##' is returned.
##'
##' @param object an object inheriting from class \code{lcc} or \code{lme},
##' representing a fitted longitudinal concordance correlation
##' function.
##'
##' @param ... other optional fitted model objects inheriting from
##' classes "lcc", or "lme".
##
##' @param test an optional logical value controlling whether likelihood
##' ratio tests should be used to compare the fitted models
##' represented by object and the objects in \code{...}. Defaults to
##' TRUE.
##'
##' @param type an optional character string specifying the type of sum
##' of squares to be used in F-tests for the terms in the model. If
##' \code{sequential}, the sequential sum of squares obtained by
##' including the terms in the order they appear in the model is used;
##' else, if \code{marginal}, the marginal sum of squares obtained by
##' deleting a term from the model at a time is used. This argument is
##' only used when a single fitted object is passed to the
##' function. Partial matching of arguments is used, so only the first
##' character needs to be provided. Defaults to \code{sequential}.
##'
##' @param adjustSigma an optional logical value. If \code{TRUE} and the
##' estimation method used to obtain object was maximum likelihood,
##' the residual standard error is multiplied by sqrt(nobs/(nobs -
##' npar)), converting it to a REML-like estimate. This argument is
##' only used when a single fitted object is passed to the
##' function. Default is \code{TRUE}.
##'
##' @param verbose an optional logical value. If \code{TRUE}, the
##' calling sequences for each fitted model object are printed with
##' the rest of the output, being omitted if \code{verbose =
##' FALSE}. Defaults to \code{FALSE}.
##'
##' @details This function is an adaptation from the
##' \code{\link{anova.lme}}. For more details see methods for
##' \code{\link{nlme}}.
##'
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @seealso \code{\link[lcc]{lcc}}, \code{\link{summary.lcc}}
##'
##' @examples
##' \dontrun{
##' ## Testing random effects
##' fm1.aov <- lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 1)
##' fm2.aov <- update(fm1.aov, qr = 2)
##' anova(fm1.aov, fm2.aov)
##' }
##'
##' @examples
##' \dontrun{
##' # Testing fixed effects
##' fm3.aov <- update(fm2.aov, REML = FALSE)
##' fm4.aov <- update(fm2.aov, REML = FALSE, qf = 3)
##' anova(fm3.aov, fm4.aov)
##' }
##'
##' @examples
##' \dontrun{
##' # Comparing the 3 lcc models
##' fm5.aov <- update(fm2.aov, var.class = varExp, weights.form = "time")
##' anova(fm1.aov, fm2.aov, fm5.aov)
##' }
##'
##' @importFrom stats formula pchisq pf terms
##' @importFrom utils head tail
##'
##' @export
anova.lcc <- function (object, ..., test = TRUE, type = c("sequential", "marginal"),
adjustSigma = TRUE, verbose = FALSE)
{
fixSig <- attr(object$model$modelStruct, "fixedSigma")
fixSig <- !is.null(fixSig) && fixSig
dots <- list(...)
if ((rt <- length(dots) + 1L) == 1L) {
if (!inherits(object$model, "lme")) {
stop("object must inherit from class \"lme\" ")
}
vFix <- attr(object$model$fixDF, "varFixFact")
if (adjustSigma && object$model$method == "ML")
vFix <- sqrt(object$model$dims$N/(object$model$dims$N - ncol(vFix))) *
vFix
c0 <- solve(t(vFix), fixef(object$model))
assign <- attr(object$model$fixDF, "assign")
nTerms <- length(assign)
type <- match.arg(type)
Fval <- Pval <- double(nTerms)
nDF <- integer(nTerms)
dDF <- object$model$fixDF$terms
for (i in 1:nTerms) {
nDF[i] <- length(assign[[i]])
if (type == "sequential") {
c0i <- c0[assign[[i]]]
}
else {
c0i <- c(qr.qty(qr(vFix[, assign[[i]], drop = FALSE]),
c0))[1:nDF[i]]
}
Fval[i] <- sum(c0i^2)/nDF[i]
Pval[i] <- 1 - pf(Fval[i], nDF[i], dDF[i])
}
aod <- data.frame(numDF = nDF, denDF = dDF, `F-value` = Fval,
`p-value` = Pval, check.names = FALSE)
rownames(aod) <- names(assign)
attr(aod, "rt") <- rt
}
else {
ancall <- sys.call()
ancall$verbose <- ancall$test <- ancall$type <- NULL
object <- list(object$model, ...)
termsClass <- vapply(object, data.class, "")
valid.cl <- c("lme", "lcc")
if (!all(match(termsClass, valid.cl, 0))) {
valid.cl <- paste0("\"", valid.cl, "\"")
stop(gettextf("objects must inherit from classes %s, or %s",
paste(head(valid.cl, -1), collapse = ", "),
tail(valid.cl, 1)), domain = NA)
}
for (i in 1:length(object)) {
if (is.lcc(object[[i]])) {
object[[i]] <- object[[i]]$model
}else {
object[[i]] <- object[[i]]
}
}
getResponseFormula <-
function(object)
{
## Return the response formula as a one sided formula
form <- formula(object)
if (!(inherits(form, "formula") && (length(form) == 3))) {
stop("'form' must be a two-sided formula")
}
eval(parse(text = paste("~", deparse(form[[2]]))))
}
resp <- vapply(object, function(el) deparse(getResponseFormula(el)[[2L]]),
"")
subs <- as.logical(match(resp, resp[1L], FALSE))
if (!all(subs))
warning("some fitted objects deleted because response differs from the first model")
if (sum(subs) == 1)
stop("first model has a different response from the rest")
object <- object[subs]
rt <- length(object)
termsModel <- lapply(object, function(el) formula(el)[-2])
estMeth <- vapply(object, function(el) if (is.null(val <- el[["method"]]))
NA_character_
else val, "")
if (length(uEst <- unique(estMeth[!is.na(estMeth)])) >
1) {
stop("all fitted objects must have the same estimation method")
}
estMeth[is.na(estMeth)] <- uEst
REML <- uEst == "REML"
if (REML) {
aux <- vapply(termsModel, function(el) {
tt <- terms(el)
val <- paste(sort(attr(tt, "term.labels")),
collapse = "&")
if (attr(tt, "intercept") == 1)
paste(val, "(Intercept)", sep = "&")
else val
}, ".")
if (length(unique(aux)) > 1) {
warning("fitted objects with different fixed effects. REML comparisons are not meaningful.")
}
}
termsCall <- lapply(object, function(el) {
if (is.null(val <- el$call) && is.null(val <- attr(el,
"call")))
stop("objects must have a \"call\" component or attribute")
val
})
termsCall <- vapply(termsCall, function(el) paste(deparse(el),
collapse = ""), "")
aux <- lapply(object, logLik, REML)
if (length(unique(vapply(aux, attr, 1, "nall"))) > 1) {
stop("all fitted objects must use the same number of observations")
}
dfModel <- vapply(aux, attr, 1, "df")
logLik <- vapply(aux, c, 1.1)
aod <- data.frame(call = termsCall, Model = 1:rt, df = dfModel,
AIC = vapply(aux, AIC, 1), BIC = vapply(aux, BIC,
1), logLik = logLik, check.names = FALSE)
if (test) {
ddf <- diff(dfModel)
if (sum(abs(ddf)) > 0) {
effects <- rep("", rt)
for (i in 2:rt) {
if (ddf[i - 1] != 0) {
effects[i] <- paste(i - 1, i, sep = " vs ")
}
}
pval <- rep(NA, rt - 1)
ldf <- as.logical(ddf)
lratio <- 2 * abs(diff(logLik))
lratio[!ldf] <- NA
pval[ldf] <- pchisq(lratio[ldf], abs(ddf[ldf]),
lower.tail = FALSE)
aod <- data.frame(aod, Test = effects, L.Ratio = c(NA,
lratio), `p-value` = c(NA, pval), check.names = FALSE,
stringsAsFactors = TRUE)
}
}
## local function for complete deparsing
c_deparse <- function(...) paste(deparse(..., width.cutoff=500),
collapse="")
row.names(aod) <- vapply(as.list(ancall[-1L]), c_deparse,
"")
attr(aod, "rt") <- rt
attr(aod, "verbose") <- verbose
}
class(aod) <- c("anova.lcc", "data.frame")
aod
}
##' @rdname print.anova.lcc
##' @title Print the Anova of an \code{lcc} Object
##' @usage
##' \method{print}{anova.lcc}(x, verbose, ...)
##' @method print anova.lcc
##' @aliases print.anova.lcc
##' @description Method print for the \code{anova.lcc}.
##'
##' @param x an object inheriting from class
##' \code{\link[lcc]{anova.lcc}}, representing a fitted longitudinal
##' concordance correlation function.
##'
##' @param verbose an optional logical value used to control the amount
##' of printed output. If \code{TRUE}, the calling sequences for each fitted
##' model object are printed with the rest of the output, being omitted
##' if \code{verbose = FALSE}. Defaults to \code{FALSE}.
##'
##' @param ... further arguments passed to \code{\link{print}}.
##'
##' @details Modified from \code{\link{anova.lme}}. For more details see
##' methods for \code{\link{nlme}}.
##' @return Return no value, called for side effects
##' @author Thiago de Paula Oliveira,
##' \email{thiago.paula.oliveira@@alumni.usp.br}
##'
##' @seealso \code{\link{summary.lcc}}, \code{\link{lccPlot}},
##' \code{\link[lcc]{lcc}}
##'
##' @examples
##'
##' \dontrun{
##' ## Second degree polynomial model with random intercept, slope and
##' ## quadratic term
##' fm1<-lcc(data = hue, subject = "Fruit", resp = "H_mean",
##' method = "Method", time = "Time", qf = 2, qr = 2)
##' print(anova(fm1))
##' }
##'
##' @export
print.anova.lcc <- function(x, verbose = attr(x, "verbose"), ...)
{
obj <- x
if ((rt <- attr(x,"rt")) == 1) {
if (!is.null(lab <- attr(x, "label"))) {
cat(lab)
}
pval <- format(round(x[, "p-value"],4))
pval[as.double(pval) == 0] <- "<.0001"
x[, "F-value"] <- format(zapsmall(x[, "F-value"]))
x[, "p-value"] <- pval
print(as.data.frame(x), ...)
} else {
if (verbose) {
cat("Call:\n")
objNams <- row.names(x)
for(i in 1:rt) {
cat(" ",objNams[i],":\n", sep ="")
cat(" ",as.character(x[i,"call"]),"\n")
}
cat("\n")
}
x <- as.data.frame(x[,-1])
for(i in names(x)) {
org <- x[[i]]
if (i == "p-value") {
org <- round(org, 4)
xna <- is.na(org)
org[!xna] <- format(org[!xna])
org[as.double(org) == 0] <- "<.0001"
org[xna] <- ""
} else {
if (match(i, c("AIC", "BIC", "logLik", "L.Ratio"), 0)) {
xna <- is.na(org)
org <- zapsmall(org)
org[xna] <- 0
org <- format(org)
org[xna] <- ""
}
}
x[[i]] <- org
}
print(as.data.frame(x), ...)
}
invisible(obj)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.