Nothing
#' @useDynLib jmcm
#' @importFrom Rcpp evalCpp
NULL
#' @import stats
NULL
#' @import graphics
NULL
#' Cattle Data
#'
#' Kenward (1987) reported an experiment in which cattle were assigned randomly
#' to two treatment groups A and B, and their body weights were recorded in
#' kilogram. Thirty animals received treatment A and another 30 received
#' treatment B. The animals were weighted 11 times over a 133-day period; the
#' first 10 measurements for each animal were made at two-week intervals and the
#' last measurement was made one week later. Since no observation was missing,
#' it is considered to be a balanced longitudinal dataset.
#'
#' \itemize{
#' \item id: subject id
#' \item day: measurement time
#' \item group: Treatment A or Treatment B
#' \item weight: cattle weight
#' }
#'
#' @docType data
#' @keywords datasets
#' @name cattle
#' @usage data(cattle)
#' @format A data frame with 660 rows and 4 variables
NULL
#' Aids Data
#'
#' The aids dataset comprises a total of 2376 CD4+ cell counts for 369 HIV
#' infected men with a follow up period of approximately eight and half year.
#' The number of measurements for each individual varies from 1 to 12 and the
#' times are not equally spaced. The CD4+ cell data are highly unbalanced.
#'
#' \itemize{
#' \item id: subject id
#' \item time: measurement time
#' \item cd4: CD4+ cell count
#' }
#'
#' @docType data
#' @keywords datasets
#' @name aids
#' @usage data(aids)
#' @format A data frame with 2376 rows and 8 variables
NULL
#' @title Fit Joint Mean-Covariance Models
#'
#' @description Fit a joint mean-covariance model to longitudinal data, via
#' maximum likelihood.
#'
#' @param formula a two-sided linear formula object describing the covariates
#' for both the mean and covariance matrix part of the model, with the response,
#' the corresponding subject id and measurement time on the left of a operator~,
#' divided by vertical bars ("|").
#' @param data a data frame containing the variables named in formula.
#' @param triple an integer vector of length three containing the degrees of the
#' three polynomial functions for the mean structure, the log innovation
#' -variances and the autoregressive or moving average coefficients when 'mcd'
#' or 'acd' is specified for cov.method. It refers to the degree for the mean
#' structure, variances and angles when 'hpc' is specified for cov.method.
#' @param cov.method covariance structure modelling method,
#' choose 'mcd' (Pourahmadi 1999), 'acd' (Chen and Dunson 2013) or 'hpc'
#' (Zhang et al. 2015).
#' @param optim.method optimization method, choose 'default' or 'BFGS' (vmmin in R)
#' @param control a list (of correct class, resulting from jmcmControl())
#' containing control parameters, see the *jmcmControl documentation for
#' details.
#' @param start starting values for the parameters in the model.
#'
#' @references Pan J, Pan Y (2017). "jmcm: An R Package for Joint Mean-Covariance Modeling of Longitudinal Data." \emph{Journal of Statistical Software}, 82(9), 1--29.
#' @examples
#' cattleA <- cattle[cattle$group=='A', ]
#' fit.mcd <- jmcm(weight | id | I(ceiling(day/14 + 1)) ~ 1 | 1,
#' data=cattleA, triple = c(8, 4, 3), cov.method = 'mcd',
#' control = jmcmControl(trace = TRUE, ignore.const.term = FALSE,
#' original.poly.order = TRUE))
#' @export
jmcm <- function(formula, data = NULL, triple = c(3, 3, 3),
cov.method = c('mcd', 'acd', 'hpc'),
optim.method = c('default','BFGS'),
control = jmcmControl(), start = NULL)
{
mc <- mcout <- match.call()
if (missing(cov.method))
stop("cov.method must be specified")
if (missing(optim.method))
optim.method = "default"
missCtrl <- missing(control)
if (!missCtrl && !inherits(control, "jmcmControl"))
{
if(!is.list(control))
stop("'control' is not a list; use jmcmControl()")
warning("please use jmcmControl() instead", immediate. = TRUE)
control <- do.call(jmcmControl, control)
}
mc[[1]] <- quote(jmcm::ldFormula)
args <- eval(mc, parent.frame(1L))
opt <- do.call(optimizeJmcm,
c(args, cov.method, optim.method, list(control=control, start=start)))
mkJmcmMod(opt=opt, args=args, triple=triple, cov.method=cov.method, optim.method=optim.method, mc=mcout)
}
#' @title Modular Functions for Joint Mean Covariance Model Fits
#'
#' @description Modular Functions for joint mean covariance model fits
#'
#' @param formula a two-sided linear formula object describing the covariates
#' for both the mean and covariance matrix part of the model, with the response,
#' the corresponding subject id and measurement time on the left of a operator~,
#' divided by vertical bars ("|").
#' @param data a data frame containing the variables named in formula.
#' @param triple an integer vector of length three containing the degrees of the
#' three polynomial functions for the mean structure, the log innovation
#' -variances and the autoregressive or moving average coefficients when 'mcd'
#' or 'acd' is specified for cov.method. It refers to the degree for the mean
#' structure, variances and angles when 'hpc' is specified for cov.method.
#' @param cov.method covariance structure modelling method,
#' choose 'mcd' (Pourahmadi 1999), 'acd' (Chen and Dunson 2013) or 'hpc'
#' (Zhang et al. 2015).
#' @param optim.method optimization method, choose 'default' or 'BFGS' (vmmin in R)
#' @param control a list (of correct class, resulting from jmcmControl())
#' containing control parameters, see the *jmcmControl documentation for
#' details.
#' @param start starting values for the parameters in the model.
#' @param m an integer vector of number of measurements for each subject.
#' @param Y a vector of responses for all subjects.
#' @param X model matrix for mean structure model.
#' @param Z model matrix for the diagonal matrix.
#' @param W model matrix for the lower triangular matrix.
#' @param time a vector of time from the data.
#' @param opt optimized results returned by optimizeJmcm.
#' @param args arguments returned by ldFormula.
#' @param mc matched call from the calling function.
#'
#' @name modular
NULL
#> NULL
#' @rdname modular
#' @export
ldFormula <- function(formula, data = NULL, triple = c(3,3,3),
cov.method = c('mcd', 'acd', 'hpc'),
optim.method = c("default", "BFGS"),
control = jmcmControl(), start = NULL)
{
mf <- mc <- match.call()
m <- match(c("formula", "data"), names(mf), 0L)
mf <- mf[c(1, m)]
f <- Formula::Formula(formula)
mf[[1]] <- as.name("model.frame")
mf$formula <- f
mf <- eval(mf, parent.frame())
Y <- Formula::model.part(f, data = mf, lhs = 1)
id <- Formula::model.part(f, data = mf, lhs = 2)
time <- Formula::model.part(f, data = mf, lhs = 3)
X <- model.matrix(f, data = mf,rhs = 1)
Z <- model.matrix(f, data = mf,rhs = 2)
index <- order(id[[1L]], time[[1L]])
Y <- Y[index, ]
id <- id[index, ]
time <- time[index, ]
m <- table(id)
attr(m, "dimnames") <- NULL
if(control$original.poly.order) {
tmp = triple[2]
triple[2] = triple[3]
triple[3] = tmp
}
# covariates from rhs of the formula
Xtmp <- X[index, -1]
Ztmp <- Z[index, -1]
# covariates based on polynomials of time
X <- rep(1, length(time))
Z <- rep(1, length(time))
if (triple[1] != 0)
for (i in 1:triple[1]) X = cbind(X, time^i)
else
X <- as.matrix(X, ncol = 1)
if (triple[2] != 0)
for (i in 1:triple[2]) Z = cbind(Z, time^i)
else
Z <- as.matrix(Z, ncol = 1)
# combine two parts of the covariates
X <- cbind(X, Xtmp)
Z <- cbind(Z, Ztmp)
W <- NULL
for (i in 1:length(m))
{
if (i == 1) {
ti <- time[1:m[1]]
} else {
first_index <- 1+sum(m[1:(i-1)])
last_index <- sum(m[1:i])
ti <- time[first_index:last_index]
}
if(m[i] != 1) {
for (j in 2:m[i])
{
for (k in 1:(j - 1))
{
wijk = (ti[j]-ti[k])^(0:triple[3])
W = rbind(W,wijk)
}
}
}
}
attr(X, "dimnames") <- NULL
attr(Z, "dimnames") <- NULL
attr(W, "dimnames") <- NULL
p <- triple[1]
d <- triple[2]
q <- triple[3]
list(m = m, Y = Y, X = X, Z = Z, W = W, time = time)
}
#' @rdname modular
#' @export
optimizeJmcm <- function(m, Y, X, Z, W, time, cov.method, optim.method, control, start)
{
missStart <- is.null(start)
lbta <- ncol(X)
llmd <- ncol(Z)
lgma <- ncol(W)
if(!missStart && (lbta+llmd+lgma) != length(start))
stop("Incorrect start input")
if (cov.method == 'mcd') {
if (missStart) {
lm.obj <- lm(Y ~ X - 1)
bta0 <- coef(lm.obj)
resid(lm.obj) -> res
lmd0 <- coef(lm(log(res^2) ~ Z - 1))
gma0 <- rep(0, lgma)
start <- c(bta0, lmd0, gma0)
if(anyNA(start)) stop("failed to find an initial value with lm(). NA detected.")
}
est <- mcd_estimation(m, Y, X, Z, W, start, Y, control$trace, control$profile, control$errormsg, FALSE, optim.method)
}
if (cov.method == 'acd') {
if (missStart) {
lm.obj <- lm(Y ~ X - 1)
bta0 <- coef(lm.obj)
resid(lm.obj) -> res
lmd0 <- coef(lm(log(res^2) ~ Z - 1))
gma0 <- rep(0, lgma)
start <- c(bta0, lmd0, gma0)
if(anyNA(start)) stop("failed to find an initial value with lm(). NA detected.")
}
est <- acd_estimation(m, Y, X, Z, W, start, Y, control$trace, control$profile, control$errormsg, FALSE, optim.method)
}
if (cov.method == 'hpc') {
if (missStart) {
lm.obj <- lm(Y ~ X - 1)
bta0 <- coef(lm.obj)
resid(lm.obj) -> res
lmd0 <- coef(lm(log(res^2) ~ Z - 1))
gma0 <- rep(0, lgma); gma0[1] <- pi / 2
start <- c(bta0, lmd0, gma0)
if(anyNA(start)) stop("failed to find an initial value with lm(). NA detected.")
}
est <- hpc_estimation(m, Y, X, Z, W, start, Y, control$trace, control$profile, control$errormsg, FALSE, optim.method)
}
if (!(control$ignore.const.term)) {
const.term = - sum(m) * 0.5 * log(2 * pi)
est$loglik = est$loglik + const.term
est$BIC = est$BIC - 2 / length(m) * const.term
}
est
}
#' @rdname modular
#' @export
mkJmcmMod <- function(opt, args, triple, cov.method, optim.method, mc)
{
if(missing(mc))
mc <- match.call()
isMCD <- (cov.method == "mcd")
isACD <- (cov.method == "acd")
isHPC <- (cov.method == "hpc")
dims <- c(nsub = length(args$m),
max.nobs = max(args$m),
p = triple[1],
d = triple[2],
q = triple[3],
MCD = isMCD,
ACD = isACD,
HPC = isHPC)
new("jmcmMod", call=mc, opt=opt, args= args,
triple=triple, devcomp=list(dims=dims) )
}
###----- Printing etc ----------------------------
methodTitle <- function(object, dims = object@devcomp$dims)
{
MCD <- dims[["MCD"]]
ACD <- dims[["ACD"]]
HPC <- dims[["HPC"]]
kind <- switch(MCD * 1L + ACD * 2L + HPC * 3L, "MCD", "ACD", "HPC")
paste("Joint mean-covariance model based on", kind)
}
cat.f <- function(...) cat(..., fill = TRUE)
.prt.methTit <- function(mtit, class) {
cat(sprintf("%s ['%s']\n", mtit, class))
}
.prt.call <- function(call, long = TRUE) {
if (!is.null(cc <- call$formula))
cat.f("Formula:", deparse(cc))
if (!is.null(cc <- call$triple))
cat.f(" triple:", deparse(cc))
if (!is.null(cc <- call$data))
cat.f(" Data:", deparse(cc))
}
.prt.loglik <- function(n2ll, digits=4)
{
t.4 <- round(n2ll, digits)
cat.f("logLik:", t.4)
}
.prt.bic <- function(bic, digits=4)
{
t.4 <- round(bic, digits)
cat.f(" BIC:", t.4)
}
print.jmcmMod <- function(x, digits=4, ...)
{
dims <- x@devcomp$dims
.prt.methTit(methodTitle(x, dims = dims), class(x))
.prt.call(x@call); cat("\n")
.prt.loglik(x@opt$loglik)
.prt.bic(x@opt$BIC); cat("\n")
cat("Mean Parameters:\n")
print.default(format(drop(x@opt$beta), digits = digits),
print.gap = 2L, quote = FALSE)
if(dims["MCD"] == 1 || dims["ACD"] == 1)
cat("Innovation Variance Parameters:\n")
else if(dims["HPC"] == 1)
cat("Variance Parameters:\n")
print.default(format(drop(x@opt$lambda), digits = digits),
print.gap = 2L, quote = FALSE)
if(dims["MCD"] == 1)
cat("Autoregressive Parameters:\n")
else if(dims["ACD"] == 1)
cat("Moving Average Parameters:\n")
else if(dims["HPC"] == 1)
cat("Angle Parameters:\n")
print.default(format(drop(x@opt$gamma), digits = digits),
print.gap = 2L, quote = FALSE)
invisible(x)
}
#' Print information for jmcmMod-class
#'
#' @param object a fitted joint mean covariance model of class "jmcmMod", i.e.,
#' typically the result of jmcm().
#'
#' @exportMethod show
setMethod("show", "jmcmMod", function(object) print.jmcmMod(object))
summary.jmcmMod <- function(x, digits=4, ...)
{
dims <- x@devcomp$dims
.prt.methTit(methodTitle(x, dims = dims), class(x))
.prt.call(x@call); cat("\n")
.prt.loglik(x@opt$loglik)
.prt.bic(x@opt$BIC); cat("\n")
cat("Mean Parameters:\n")
print.default(format(drop(x@opt$beta), digits = digits),
print.gap = 2L, quote = FALSE)
if(dims["MCD"] == 1 || dims["ACD"] == 1)
cat("Innovation Variance Parameters:\n")
else if(dims["HPC"] == 1)
cat("Variance Parameters:\n")
print.default(format(drop(x@opt$lambda), digits = digits),
print.gap = 2L, quote = FALSE)
if(dims["MCD"] == 1)
cat("Autoregressive Parameters:\n")
else if(dims["ACD"] == 1)
cat("Moving Average Parameters:\n")
else if(dims["HPC"] == 1)
cat("Angle Parameters:\n")
print.default(format(drop(x@opt$gamma), digits = digits),
print.gap = 2L, quote = FALSE)
invisible(x)
}
# bdbind <- function(A, B)
# {
# stopifnot(isSymmetric(A))
# stopifnot(isSymmetric(B))
# dimA <- dim(A)[1]
# dimB <- dim(B)[1]
#
# res <- matrix(0, dimA+dimB, dimA+dimB)
# res[1:dimA, 1:dimA] = A
# res[(dimA+1):(dimA+dimB), (dimA+1):(dimA+dimB)] = B
# res
# }
#
# vcov.merMod <- function(object)
# {
# args <- object@args
# dims <- object@devcomp$dims
# par <- object@opt$par
#
# m <- args[["m"]]
# Y <- args[["Y"]]
# X <- args[["X"]]
# Z <- args[["Z"]]
# W <- args[["W"]]
#
# if(dims[["MCD"]] == 1) {
# mcdobj <- new("MCD", m, Y, X, Z, W)
# res <- solve(mcdobj$meancovi(par, 1)$Sigmai)
# for(i in 2:length(m)) {
# res <- bdbind(res, solve(mcdobj$meancovi(par, i)$Sigmai))
# }
# }
# res <- solve(t(X) %*% res %*% X)
# }
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.