#######################################################
## Functions to generate commonly required objects ##
#######################################################
## std.b: Return the spherical random effects or "Standardize" the vector
## vector: \eqn{\Lambda^{-1} vector / \sigma}{Lambda^-1 vector / sigma}
##
## @title Standardized values
## @param object rlmerMod object
## @param sigma to use for standardization
## @param vector vector to standardize
## @rdname std
std.b <- function(object, vector)
object@pp$stdB(vector)
## std.e: Calculate the standardized residuals or "Standardize" the Matrix
## sigma: \eqn{R^{-1} matrix / \sigma}{R^-1 matrix / sigma}
##
## @rdname std
std.e <- function(object, sigma = .sigma(object), matrix, drop=TRUE) {
if (missing(matrix)) return(object@resp$wtres / sigma)
## for the moment: just divide by sigma
if (drop) matrix <- drop(matrix)
matrix/sigma
}
## Calculate scaled squared Mahalanobis distances per group
##
## @title Distances per group
## @param object rlmerMod object
## @param sigma to use for standardization
## @param bs spherical random effects
## @param center whether to return the centered distances.
## @param ... ignored
.dk <- function(object, sigma, center, bs = b.s(object), ...) {
ua <- uArranged(object, bs / sigma)
unlist(lapply(seq_along(object@blocks), function(bt) {
us <- ua[[bt]]
s <- ncol(us)
if (s == 1)
return(us)
## else: square, sum and subtract s
ret <- computeDiagonalOfTCrossproductNumericMatrix(us)
if (center)
ret <- ret - .kappa_b(object)[bt] * s
ret
}))
}
## these helper functions are for the non-centered case, i.e., for estimating
## the random effects themselves. modularize distance function: compute
## sum(bs^2) - s
.d <- function(bs, s = length(bs)) {
if (s == 1)
return(bs)
if (is.matrix(bs))
computeDiagonalOfTCrossproductNumericMatrix(bs)
else
sum(bs * bs)
}
## same function, but assume we've already summed
.d2 <- function(sbs2, s) {
if (s == 1) stop("s must be larger than 1") ## disable this test?
sbs2
}
## inner derivative of .d2:
.Dd2 <- function(sbs2, s) {
##if (s == 1) stop("s must be larger than 1")
2
}
## dist.b: Calculate the distance from 0 standardized by sigma. This is just
## value divided by sigma for uncorrelated observations. For correlated items,
## this is the Mahalanobis distance from 0. If \code{shifted=TRUE} and
## correlated items, the squared distances are centered by -kappa_b*s. This is
## required to compute the weights used for the size of the covariance matrix of
## the random effects.
##
## @title Calculate distance
## @param object object to use
## @param sigma scale for standardization
## @param center whether to use the centered distances
## @param ... passed on to internal functions.
## @rdname dist
dist.b <- function(object, sigma = .sigma(object), center=FALSE, ...) {
db <- .dk(object, sigma, center, ...)
db[object@k]
}
## dist.e: Calculate dist for residuals always assume they are uncorrelated
##
## @rdname dist
dist.e <- function(object, sigma = .sigma(object)) {
std.e(object, sigma) ## just the usual rescaled residuals
}
## wgt.b: Calculate the robustness weights psi(d) / d, standardized by sigma.
## The robustness weights are calculated with d the Mahalanobis distance. Each
## group of correlated items then gets a constant weight. The robustness weights
## for the random effects themselves are different than the ones used for
## estimating the size of the covariance matrix of the random effects. Those are
## additionally centered. That way, inlier can also be downweighted. If
## \code{center=TRUE}, then the centered distances are used to compute the
## robustness weights and the weight function given by rho.sigma.b is used.
##
## @title Calculate robustness weights
## @param object object to use
## @param sigma scale for standardization
## @param center whether return the centered robustness weights, see Details.
## @rdname wgt
## @export
wgt.b <- function(object, sigma = .sigma(object), center = FALSE,
use.rho.sigma = center) {
db <- dist.b(object, sigma, center)
rho <- rho.b(object, if (use.rho.sigma) "sigma" else "default")
ret <- numeric()
for (bt in seq_along(object@blocks)) {
bind <- as.vector(object@idx[[bt]])
ret <- c(ret, rho[[bt]]@wgt(db[bind]))
}
ret
}
## wgt.e: robustness weights of residuals
##
## @param use.rho.sigma return the weights computed using rho.sigma.
## @rdname wgt
## @export
wgt.e <- function(object, sigma = .sigma(object), use.rho.sigma = FALSE)
if (use.rho.sigma) object@rho.sigma.e@wgt(dist.e(object, sigma)) else
object@rho.e@wgt(dist.e(object, sigma))
### Calculate robustness weights * squared effect
### Return sensible result in the infinite case
## Assume: x infinite <=> y infinite
.wgtxy2 <- function(rho, x, y)
.wgtxy(rho, x, y*y)
## wgt(x) * y
.wgtxy <- function(rho, x, y) {
ret <- rho@wgt(x) * y
ret[is.infinite(x) & is.infinite(y)] <- if (rho@psi(Inf) == 0) 0 else Inf
ret
}
## Find blocks of correlated random effects
##
## @title Find blocks in Lambda
## @param obj reTrms-object
## @param Lambdat transpose of matrix Lambda (U_b(theta))
## @param Lind vector of indices, mapping theta to Lambdat
## @return list of blocks
findBlocks <- function(obj, Lambdat=obj$Lambdat(), Lind=obj$Lind) {
LambdaInd <- Lambdat
LambdaInd@x[] <- as.double(Lind)
LambdaInd <- t(LambdaInd)
LambdaInd <- as(LambdaInd, "matrix") ## to avoid attempt to apply non function error
bend <- unique(apply(LambdaInd != 0, 2, function(x) max(which(x))))
nblocks <- length(bend)
bstart <- c(1, bend[-nblocks]+1)
bidx <- lapply(1:nblocks, function(i) seq.int(bstart[i],bend[i]))
blocks <- lapply(bidx, function(idx) LambdaInd[idx,idx])
bind <- match(blocks, ublocks <- unique(blocks))
k <- unlist(lapply(1:nblocks, function(i) rep(i, length(bidx[[i]]))))
bdim <- sapply(ublocks, NCOL)
bidx <- lapply(1:length(ublocks), function (i)
matrix(unlist(bidx[bind == i]),nrow = bdim[i]))
q <- sapply(bidx, length)
list(blocks = ublocks, ind = bind, idx = bidx, dim = bdim, q = q, k = k)
}
lchol <- function(x) {
r <- try(chol.default(x), silent=TRUE)
## if chol fails, return sqrt of diagonal
if (is(r, "try-error")) {
Diagonal(x = sqrt(diag(x)))
} else r
}
#######################################################
## Summary / printing methods ##
#######################################################
.summary.cor.max <- 20
## Print method
## along the lines of printMerenv of lme4
##' @export
print.summary.rlmerMod <- function(x, digits = max(3, getOption("digits") - 3),
correlation = NULL, symbolic.cor = FALSE,
signif.stars = getOption("show.signif.stars"),
ranef.comp = c("Variance", "Std.Dev."),
show.resids = TRUE, ...) {
## check if doFit = FALSE is in call
if (!is.null(x$call$doFit) && !x$call$doFit) {
cat("Unfitted rlmerMod object. Use update(object, doFit=TRUE) to fit it.\n")
return(invisible(x))
}
## title
cat(x$methTitle, "\n")
## these are in fromLme4.R
.prt.call(x$call); cat("\n")
if (show.resids)
## need residuals.merMod() rather than residuals():
## summary.merMod has no residuals method
.prt.resids(x$residuals, digits=digits)
.prt.VC(x$varcor, digits=digits, useScale= x$useScale,
comp = ranef.comp, ...)
.prt.grps(x$ngrps, nobs= x$devcomp$dims[["n"]])
## fixed effecs
## this part is 1:1 from printMerenv
p <- nrow(x$coefficients)
if (p > 0) {
cat("\nFixed effects:\n")
printCoefmat(x$coefficients, zap.ind = 3, #, tst.ind = 4
digits = digits, signif.stars = signif.stars)
if(is.null(correlation)) { # default
correlation <- p <= .summary.cor.max
if(!correlation) {
nam <- deparse(substitute(x))
if(length(nam) > 1 || nchar(nam) >= 32) nam <- "...."
message(sprintf(paste(
"\nCorrelation matrix not shown by default, as p = %d > %d.",
"Use print(%s, correlation=TRUE) or",
" vcov(%s) if you need it\n", sep="\n"),
p, .summary.cor.max, nam, nam))
}
}
else if(!is.logical(correlation)) stop("'correlation' must be NULL or logical")
if(correlation) {
if(is.null(VC <- x$vcov)) VC <- vcov(x, correlation=TRUE)
corF <- VC@factors$correlation
if (is.null(corF)) {
message("\nCorrelation of fixed effects could have been required in summary()")
corF <- cov2cor(VC)
} ## else {
p <- ncol(corF)
if (p > 1) {
rn <- rownames(x$coefficients)
rns <- abbreviate(rn, minlength=11)
cat("\nCorrelation of Fixed Effects:\n")
if (is.logical(symbolic.cor) && symbolic.cor) {
corf <- as(corF, "matrix")
dimnames(corf) <- list(rns,
abbreviate(rn, minlength=1, strict=TRUE))
print(symnum(corf))
} else {
corf <- matrix(format(round(corF@x, 3), nsmall = 3),
ncol = p,
dimnames = list(rns, abbreviate(rn, minlength=6)))
corf[!lower.tri(corf)] <- ""
print(corf[-1, -p, drop=FALSE], quote = FALSE)
} ## !symbolic.cor
} ## if (p > 1)
} ## if (correlation)
} ## if (p>0)
## robustness weights
summarizeRobWeights(x$wgt.e, digits=3,
header="\nRobustness weights for the residuals:")
summarizeRobWeights(x$wgt.b, digits=3,
header="\nRobustness weights for the random effects:")
## rho functions
cat("\nRho functions used for fitting:\n")
cat(" Residuals:\n")
cat(" eff:", .sprintPsiFunc(x$rho.e, short=TRUE), "\n")
cat(" sig:", .sprintPsiFunc(x$rho.sigma.e, short=TRUE), "\n")
for (bt in seq_along(x$rho.b)) {
cat(" Random Effects, variance component ", bt, " (", names(x$rho.b)[bt], "):\n", sep="")
cat(" eff:", .sprintPsiFunc(x$rho.b[[bt]], short=TRUE), "\n")
cat(" vcp:", .sprintPsiFunc(x$rho.sigma.b[[bt]], short=TRUE), "\n")
}
invisible(x)
}
.methTitle <- function(object)
sprintf("Robust linear mixed model fit by %s", object@method)
##' @export
print.rlmerMod <- function(x, digits = max(3, getOption("digits") - 3),
correlation = NULL, symbolic.cor = FALSE,
signif.stars = getOption("show.signif.stars"),
ranef.comp = "Std.Dev.", ...) {
## check if doFit = FALSE is in call
if (!is.null(x@call$doFit) && !x@call$doFit) {
cat("Unfitted rlmerMod object. Use update(object, doFit=TRUE) to fit it.\n")
return(invisible(x))
}
dims <- (devC <- x@devcomp)$dims
## title
cat(.methTitle(x), "\n")
.prt.call(x@call)
useScale <- as.logical(dims[["useSc"]])
varcor <- VarCorr(x)
.prt.VC(varcor, digits=digits, comp = ranef.comp, ...)
ngrps <- sapply(x@flist, function(x) length(levels(x)))
.prt.grps(ngrps, nobs= dims[["n"]])
if(length(cf <- fixef(x)) >= 0) {
cat("Fixed Effects:\n")
print.default(format(cf, digits = digits),
print.gap = 2L, quote = FALSE, ...)
} else cat("No fixed effect coefficients\n")
invisible(x)
}
##' @export
## this follows the lines of summary.merMod of lme4
summary.rlmerMod <- function(object, ...) {
resp <- object@resp
devC <- object@devcomp
dd <- devC$dims
cmp <- devC$cmp
sig <- sigma(object)
coefs <- cbind("Estimate" = fixef(object),
"Std. Error" = sig * sqrt(diag(object@pp$unsc())))
if (nrow(coefs) > 0) {
coefs <- cbind(coefs, coefs[,1]/coefs[,2], deparse.level=0)
colnames(coefs)[3] <- "t value"
}
varcor <- VarCorr(object)
structure(list(methTitle=.methTitle(object), devcomp=devC,
ngrps=sapply(object@flist, function(x) length(levels(x))),
coefficients=coefs, sigma=sig,
vcov=vcov(object, correlation=TRUE, sigm=sig),
varcor=varcor, # and use formatVC(.) for printing.
call=object@call,
wgt.e=wgt.e(object),
wgt.b=wgt.b(object),
rho.e=rho.e(object),
rho.sigma.e=rho.e(object, "sigma"),
rho.b=rho.b(object),
rho.sigma.b=rho.b(object, "sigma"),
residuals=residuals(object, scaled=TRUE)
), class = "summary.rlmerMod")
}
##' Use \code{compare} to quickly compare the estimated parameters of the fits
##' of multiple lmerMod or rlmerMod objects.
##'
##' The functions \code{xtable.comparison.table} and
##' \code{print.xtable.comparison.table} are wrapper functions for the
##' respective \code{\link{xtable}} and \code{\link{print.xtable}} functions.
##'
##' The function \code{getInfo} is internally used to prepare object for
##' producing a comparison chart in \code{compare}.
##'
##' @title Create comparison charts for multiple fits
##' @param ... objects to compare, or, for the \code{\link{xtable}} functions:
##' passed to the respective \code{\link{xtable}} function.
##' @param digits number of digits to show in output
##' @param dnames names of objects given as arguments (optional)
##' @param show.rho.functions whether to show rho functions in output.
##' @keywords models utilities
##' @examples
##' \dontrun{
##' fm1 <- lmer(Yield ~ (1|Batch), Dyestuff)
##' fm2 <- rlmer(Yield ~ (1|Batch), Dyestuff)
##' compare(fm1, fm2)
##' require(xtable)
##' xtable(compare(fm1, fm2))
##' str(getInfo(fm1))
##' }
##' @export
compare <- function(..., digits = 3, dnames = NULL,
show.rho.functions = TRUE) {
linfos <- list(...)
if (!missing(dnames) && !is.null(dnames)) names(linfos) <- dnames
linfos <- lapply(linfos, getInfo)
## check if all methods work at least on the same dataset
if (length(unique(sapply(linfos, function(x) x$data))) > 1)
warning("Comparison for objects not fitted to the same dataset")
## local helper functions
.NULLtoNA <- function(lst)
lapply(lst, function(x) if (is.null(x)) NA else x)
.getComp <- function(lst, comp)
sapply(lst, function(x) .NULLtoNA(x[comp]))
.getComp2 <- function(lst, comp)
sapply(lst, function(x) x[comp])
.dropComp <- function(lst, comp)
lapply(lst, function(x) { x[comp] <- NULL; x })
.getNames <- function(lst)
unique(unlist(lapply(lst, names)))
.combineComp <- function(lst) {
names <- .getNames(lst)
tmp <- .getComp2(lst, names)
if (!is.matrix(tmp))
tmp <- matrix(tmp, ncol = length(tmp))
rownames(tmp) <- names
format(tmp, digits = digits)
}
lnames <- .getNames(linfos)
## header
call <- match.call()
call$digits <- NULL
call$names <- NULL
call$show.rho.functions <- NULL
split <- rep("", length(linfos))
## prepare coef slot:
## add stderr
linfos <- lapply(linfos, function(linfo) {
cf <- paste(format(linfo[["coef"]], digits=digits), " (",
format(linfo[["stderr"]], digits=digits), ")", sep="")
names(cf) <- names(linfo[["coef"]])
linfo[["coef"]] <- cf
linfo[["stderr"]] <- NULL
linfo
})
## combine
## coefficients
ret <- rbind(Coef=split,
.combineComp(.getComp(linfos, "coef")))
## variance components
ret <- rbind(ret, NULL=split, VarComp=split,
.combineComp(.getComp(linfos, "varcomp")))
## correlations if there are any
if (any(sapply(linfos, function(linfo) !is.null(linfo$correlations))))
ret <- rbind(ret, NULL=split, Correlations=split,
.combineComp(.getComp(linfos, "correlations")))
## sigma
ret <- rbind(ret, NULL=split, format(.getComp(linfos, "sigma"), digits = digits))
rownames(ret)[nrow(ret)] <- "sigma"
## drop the items already included
linfos <- .dropComp(linfos, c("data", "coef", "varcomp", "correlations", "sigma"))
## drop rho functions if requested
if (!show.rho.functions)
linfos <- .dropComp(linfos, grep("^rho", .getNames(linfos), value=TRUE))
## show the rest if there is any
if (length(.getNames(linfos)) > 0) {
ret <- rbind(ret, NULL=split)
for (name in .getNames(linfos)) {
ret <- rbind(ret, format(.getComp(linfos, name), digits = digits))
rownames(ret)[nrow(ret)] <- name
}
}
## clean up and finish
ret <- gsub("\\s*(NA|NULL)", "", ret)
colnames(ret) <- if (missing(dnames) && is.null(dnames)) as.character(call)[-1] else dnames
rownames(ret)[rownames(ret) == "NULL"] <- ""
class(ret) <- "comparison.table"
ret
}
##' @export
print.comparison.table <- function(x, ...) {
class(x) <- "matrix"
print(x, ..., quote=FALSE)
}
##' @rdname compare
##' @method getInfo lmerMod
##' @export
getInfo.lmerMod <- function(object, ...) {
lsum <- summary(object)
coefs <- lsum$coefficients
varcor <- lsum$varcor
isREML <- .isREML(object)
.namedVector <- function(mat) {
if (is.vector(mat)) return(mat)
names <- rownames(mat)
ret <- drop(mat)
names(ret) <- names
ret
}
.getVC <- function(varcor) {
vc <- lapply(varcor, function(grp) attr(grp, "stddev"))
ret <- unlist(vc, use.names = FALSE)
names(ret) <-
unlist(lapply(1:length(vc), function(i)
paste(names(vc[[i]]), names(vc)[i], sep=" | ")))
##.namedVector(ret)
ret
}
.getCorr <- function(varcor) {
ret <- lapply(1:length(varcor), function(i) {
grp <- varcor[[i]]
corr <- attr(grp, "correlation")
if (nrow(corr) == 1) return(NULL)
names <- outer(colnames(corr), paste("x", rownames(corr)), paste)
ret <- as.vector(corr[upper.tri(corr)])
names(ret) <- paste(as.vector(names[upper.tri(names)]),
names(varcor)[i], sep = " | ")
ret
})
unlist(ret)
}
ret <- list(data = object@call$data,
coef = .namedVector(coefs[,1,drop=FALSE]),
stderr = .namedVector(coefs[,2,drop=FALSE]),
varcomp = .getVC(varcor),
sigma = sigma(object))
corrs <- .getCorr(varcor)
if (length(corrs) > 0) ret$correlations <- corrs
if (!is(object, "rlmerMod")) {
if (isREML) {
ret$REML <- lme4::REMLcrit(object)
} else {
ret$deviance <- deviance(object, REML=FALSE)
}
}
ret
}
##' @rdname compare
##' @method getInfo rlmerMod
##' @export
getInfo.rlmerMod <- function(object, ...) {
linfo <- getInfo.lmerMod(object)
linfo$REML <- linfo$deviance <- NULL
linfo$rho.e <- .sprintPsiFunc(rho.e(object), TRUE)
linfo$rho.sigma.e <- .sprintPsiFunc(rho.e(object, "sigma"), TRUE)
rho.b <- rho.b(object)
rho.sigma.b <- rho.b(object, "sigma")
for (bt in seq_along(object@blocks)) {
linfo[[paste("rho.b",bt,sep="_")]] <- .sprintPsiFunc(rho.b[[bt]], TRUE)
linfo[[paste("rho.sigma.b",bt,sep="_")]] <- .sprintPsiFunc(rho.sigma.b[[bt]], TRUE)
}
linfo
}
##' The functions \code{xtable.comparison.table} and
##' \code{print.xtable.comparison.table} are wrapper functions for the
##' respective \code{\link{xtable}} and \code{\link{print.xtable}} functions.
##'
##' @rdname compare
##' @param x object of class "comparison.table" or "xtable.comparison.table"
##' @param caption see \code{\link{xtable}}.
##' @param label see \code{\link{xtable}}.
##' @param align see \code{\link{xtable}}.
##' @param display see \code{\link{xtable}}.
##' @seealso \code{\link{xtable}}
##' @importFrom xtable xtable
##' @export
##' @method xtable comparison.table
xtable.comparison.table <- function(x, caption=NULL, label=NULL, align=NULL,
digits=NULL, display=NULL, ...) {
rn <- sapply(rownames(x), function(n) {
switch(n,
Coef="Coefficients (Std. Error)",
VarComp="Variance components",
Correlations="Correlations",
n) })
tbl <- cbind(rn, x)
rownames(tbl) <- NULL
colnames(tbl) <- c(" ", colnames(x))
if (is.null(align)) align <- c("r", "r", rep.int("l", ncol(x)))
xtbl <- xtable(tbl, caption=caption, label=label, align=align,
digits=digits, display=display, ...)
class(xtbl) <- c("xtable.comparison.table", class(xtbl))
xtbl
}
##' @rdname compare
##' @param add.hlines replace empty lines in comparison table by hlines.
##' Supersedes \code{hline.after} argument of \code{print.xtable}.
##' @param latexify.namescol replace \dQuote{sigma} and \dQuote{x} in the first
##' column by latex equivalents.
##' @param include.rownames include row numbers (the object returned by
##' \code{xtable.comparison.table} includes names in the first column)
##' @importFrom xtable print.xtable
##' @seealso \code{\link{print.xtable}}
##' @method print xtable.comparison.table
##' @export
print.xtable.comparison.table <- function(x, add.hlines=TRUE,
latexify.namescol=TRUE,
include.rownames=FALSE, ...) {
args <- list(...)
if (add.hlines) {
rns <- if (is.factor(x[[1]])) levels(x[[1]])[x[[1]]] else x[[1]]
emptyCol <- sapply(rns, nchar) == 0
x <- x[!emptyCol,]
args$hline.after=c(-1, 0, which(emptyCol)-1:sum(emptyCol),nrow(x))
}
if (latexify.namescol) {
if (!is.null(args$sanitize.text.function)) {
sanitize <- args$sanitize.text.function
} else {
## from print.xtable
sanitize <- function(str) {
result <- str
result <- gsub("\\\\","SANITIZE.BACKSLASH",result)
result <- gsub("$","\\$",result,fixed=TRUE)
result <- gsub(">","$>$",result,fixed=TRUE)
result <- gsub("<","$<$",result,fixed=TRUE)
result <- gsub("|","$|$",result,fixed=TRUE)
result <- gsub("{","\\{",result,fixed=TRUE)
result <- gsub("}","\\}",result,fixed=TRUE)
result <- gsub("%","\\%",result,fixed=TRUE)
result <- gsub("&","\\&",result,fixed=TRUE)
result <- gsub("_","\\_",result,fixed=TRUE)
result <- gsub("#","\\#",result,fixed=TRUE)
result <- gsub("^","\\verb|^|",result,fixed=TRUE)
result <- gsub("~","\\~{}",result,fixed=TRUE)
result <- gsub("SANITIZE.BACKSLASH","$\\backslash$",result,fixed=TRUE)
return(result)
}
}
## sanitize text
for (i in 1:ncol(x)) x[[i]] <- sanitize(x[[i]])
args$sanitize.text.function <- identity
x[[1]] <- sapply(x[[1]], function(rn) {
rn <- sub(" x ", " $\\times$ ", rn, fixed=TRUE)
rn <- sub("\\bsigma\\b", "$\\\\sigma$", rn)
rn
})
}
args$include.rownames <- include.rownames
args$x <- x
do.call(print.xtable, args)
}
##' @export
update.rlmerMod <- function(object, formula., ..., evaluate = TRUE) {
## update call
## set old object as init object
## run it
if (is.null(call <- object@call))
stop("object should contain a 'call' component")
extras <- match.call(expand.dots = FALSE)$...
if (!missing(formula.)) {
call$formula <- update.formula(formula(object), formula.)
extras$init <- NULL
} else {
## set init to object, if not explicitly given (and no new data given)
if (is.null(extras[["data"]])) {
if (is.null(extras[["init"]])) {
lcall <- sys.call(sys.parent())
extras$init <- object
}
## copy pp and resp (to really get a new object)
extras$init@pp <- object@pp$copy()
## reset calledInit... fields to FALSE:
fields <- grep("calledInit", names(getRefClass(class(extras$init@pp))$fields()), value=TRUE)
Map(function(field) extras$init@pp$field(field, FALSE), fields)
extras$init@resp <- object@resp$copy()
} else {
extras$init <- NULL
}
}
if (length(extras) > 0) {
existing <- !is.na(match(names(extras), names(call)))
for (a in names(extras)[existing]) call[[a]] <- extras[[a]]
if (any(!existing)) {
call <- c(as.list(call), extras[!existing])
call <- as.call(call)
}
}
if (evaluate)
eval(call, parent.frame())
else call
}
##' @importFrom utils head
createObservationsString <- function(indicator,
max = 6,
obs = which(indicator)) {
n <- length(obs)
if (n == 1) {
return(paste("observation", obs))
}
if (n > max) {
nToShow <- max - 1
} else {
nToShow <- n - 1
}
obsString <- paste("observations",
paste(head(obs, n = nToShow),
collapse = ", "))
if (n > max) {
obsString <- paste(obsString, "and", n - nToShow, "others")
} else {
obsString <- paste(obsString, "and", obs[n])
}
return(obsString)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.