Nothing
#' @title Anova Tables
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Performs Wald tests of the significance for the linear
#' predictor components by response variables. This function is useful
#' for joint hypothesis tests of regression coefficients associated with
#' categorical covariates with more than two levels. It is not designed
#' for model comparison.
#'
#' @param object an object of class \code{mcglm}, usually, a result of a
#' call to \code{mcglm()} function.
#' @param ... additional arguments affecting the summary produced. Note
#' that there is no extra options for mcglm object class.
#'
#' @return A \code{data.frame} with Chi-square statistic to test the
#' null hypothesis of a parameter, or a set of parameters, be
#' zero. Degree of freedom (Df) and p-values.
#' The Wald test based on the observed covariance matrix of
#' the parameters is used.
#' @examples
#' x1 <- seq(-1, 1, l = 100)
#' x2 <- gl(5, 20)
#' beta <- c(5, 0, -2, -1, 1, 2)
#' X <- model.matrix(~ x1 + x2)
#' set.seed(123)
#' y <- rnorm(100, mean = X%*%beta, sd = 1)
#' data = data.frame("y" = y, "x1" = x1, "x2" = x2)
#' fit.anova <- mcglm(c(y ~ x1 + x2), list(mc_id(data)), data = data)
#' anova(fit.anova)
#'
#' @method anova mcglm
#' @export
anova.mcglm <- function(object, ...) {
n_resp <- length(object$mu_list)
n_beta <- lapply(object$list_X, ncol)
idx.list <- list()
for (i in 1:n_resp) {
idx.list[[i]] <- rep(i, n_beta[i])
}
vv <- vcov(object)
n_par <- dim(vv)[1]
idx.vec <- do.call(c, idx.list)
n_cov <- n_par - length(idx.vec)
idx.vec <- c(idx.vec, rep(0, n_cov))
temp.vcov <- list()
temp.beta <- list()
for (i in 1:n_resp) {
idx.id <- idx.vec == i
temp.vcov[[i]] <- vv[idx.id, idx.id]
temp.beta[[i]] <-
coef(object, type = "beta", response = i)$Estimates
}
saida <- list()
for (i in 1:n_resp) {
idx <- attr(object$list_X[[i]], "assign")
names <- colnames(object$list_X[[i]])
if (names[1] == "(Intercept)") {
idx <- idx[-1]
names <- names[-1]
temp.beta[[i]] <- temp.beta[[i]][-1]
temp.vcov[[i]] <- temp.vcov[[i]][-1, -1]
}
n_terms <- length(unique(idx))
X2.resp <- list()
for (j in 1:n_terms) {
idx.TF <- idx == j
temp <- as.numeric(
t(temp.beta[[i]][idx.TF]) %*%
solve(as.matrix(temp.vcov[[i]])[idx.TF, idx.TF]) %*%
temp.beta[[i]][idx.TF])
nbeta.test <- length(temp.beta[[i]][idx.TF])
X2.resp[[j]] <-
data.frame(Covariate = names[idx.TF][1],
Chi.Square = round(temp, 4), Df = nbeta.test,
p.value = round(pchisq(temp, nbeta.test,
lower.tail = FALSE),
4))
}
saida[[i]] <- do.call(rbind, X2.resp)
}
cat("Wald test for fixed effects\n")
for(i in 1:n_resp) {
cat("Call: ")
print(object$linear_pred[[i]])
cat("\n")
print(saida[[i]])
cat("\n")
}
return(invisible(saida))
}
#' @title Model Coefficients
#' @name coef.mcglm
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Extract model coefficients for objects
#' of \code{mcglm} class.
#'
#' @param object an object of \code{mcglm} class.
#' @param std.error logical. If \code{TRUE} returns the standard errors
#' for the estimates. Default is \code{FALSE}.
#' @param response a numeric vector specifyng for which response
#' variable the coefficients should be returned.
#' @param type a string vector (can be 1 element length) specifying which
#' coefficients should be returned. \cr
#' Options are \code{"beta"}, \code{"tau"}, \code{"power"}, \code{"tau"} and
#' \code{"correlation"}.
#' @param ... additional arguments affecting the summary produced. Note
#' that there is no extra options for \code{mcglm} object class.
#'
#' @return A \code{data.frame} with parameters names, estimates,
#' response variable number and parameters type.
#'
#' @method coef mcglm
#' @export
coef.mcglm <- function(object, std.error = FALSE,
response = c(NA, 1:length(object$beta_names)),
type = c("beta", "tau", "power", "correlation"),
...) {
n_resp <- length(object$beta_names)
cod_beta <- list()
cod_power <- list()
cod_tau <- list()
type_beta <- list()
type_power <- list()
type_tau <- list()
resp_beta <- list()
resp_power <- list()
resp_tau <- list()
response_for <- 1:n_resp
for (i in response_for) {
cod_beta[[i]] <- paste0(
paste0("beta", i), 0:c(object$Information$n_betas[[i]] - 1))
type_beta[[i]] <- rep("beta", length(cod_beta[[i]]))
resp_beta[[i]] <- rep(response_for[i], length(cod_beta[[i]]))
if (object$Information$n_power[[i]] != 0 |
object$power_fixed[[i]] == FALSE) {
cod_power[[i]] <- paste0(
paste0("power", i), 1:object$Information$n_power[[i]])
type_power[[i]] <- rep("power",
length(cod_power[[i]]))
resp_power[[i]] <- rep(response_for[i],
length(cod_power[[i]]))
}
if (object$Information$n_power[[i]] == 0) {
cod_power[[i]] <- rep(1, 0)
type_power[[i]] <- rep(1, 0)
resp_power[[i]] <- rep(1, 0)
}
cod_tau[[i]] <- paste0(
paste0("tau", i), 1:object$Information$n_tau[[i]])
type_tau[[i]] <- rep("tau", length(cod_tau[[i]]))
resp_tau[[i]] <- rep(response_for[i], length(cod_tau[[i]]))
}
rho_names <- c()
if (n_resp != 1) {
combination <- combn(n_resp, 2)
for (i in 1:dim(combination)[2]) {
rho_names[i] <- paste0(
paste0("rho", combination[1, i]), combination[2, i])
}
}
type_rho <- rep("correlation", length(rho_names))
resp_rho <- rep(NA, length(rho_names))
cod <- c(do.call(c, cod_beta), rho_names,
do.call(c, Map(c, cod_tau)))
type_cod <- c(do.call(c, type_beta), type_rho,
do.call(c, Map(c, type_tau)))
response_cod <- c(do.call(c, resp_beta), resp_rho,
do.call(c, Map(c, resp_tau)))
if (length(cod_power) != 0) {
cod <- c(do.call(c, cod_beta), rho_names,
do.call(c, Map(c, cod_power, cod_tau)))
type_cod <- c(do.call(c, type_beta), type_rho,
do.call(c, Map(c, type_power, type_tau)))
response_cod <- c(do.call(c, resp_beta), resp_rho,
do.call(c, Map(c, resp_power, resp_tau)))
}
Estimates <- c(object$Regression, object$Covariance)
coef_temp <- data.frame(
Estimates = Estimates,
Parameters = cod,
Type = type_cod,
Response = response_cod)
if (std.error == TRUE) {
coef_temp <- data.frame(
Estimates = Estimates,
Std.error = sqrt(diag(object$vcov)),
Parameters = cod, Type = type_cod,
Response = response_cod)
}
output <- coef_temp[
which(coef_temp$Response %in% response &
coef_temp$Type %in% type), ]
return(output)
}
#' @title Confidence Intervals for Model Parameters
#' @name confint.mcglm
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Computes confidence intervals for parameters in a fitted
#' \code{mcglm} model.
#'
#' @param object a fitted \code{mcglm} object.
#' @param parm specifies for which parameters are to be given
#' confidence intervals, either a vector of numbers or a vector of
#' strings. If missing, all parameters are considered.
#' @param level the nominal confidence level.
#' @param ... additional arguments affecting the confidence interval
#' produced. Note that there is no extra options for \code{mcglm}
#' object class.
#'
#' @return A \code{data.frame} with confidence intervals, parameters
#' names, response variable number and parameters type.
#' @method confint mcglm
#' @export
confint.mcglm <- function(object, parm, level = 0.95, ...) {
temp <- coef(object, std.error = TRUE)
if (missing(parm)) {
parm <- 1:length(temp$Estimates)
}
a <- (1 - level)/2
a <- c(a, 1 - a)
fac <- stats::qnorm(a)
ci <- temp$Estimates + temp$Std.error %o% fac
colnames(ci) <- paste0(format(a, 2), "%")
rownames(ci) <- temp$Parameters
return(ci[parm, ])
}
#' @title Fitted Values
#' @name fitted.mcglm
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Extract fitted values for objects of \code{mcglm} class.
#'
#' @param object an object of \code{mcglm} class.
#' @param ... additional arguments affecting the summary produced. Note
#' that there is no extra options for \code{mcglm} object class.
#'
#' @return A matrix with fitted values.
#'
#' @method fitted mcglm
#' @export
fitted.mcglm <- function(object, ...) {
n_resp <- length(object$beta_names)
output <- Matrix(object$fitted, ncol = n_resp, nrow = object$n_obs)
return(output)
}
#' @title Residuals and algorithm check plots
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Residual and algorithm check analysis for objects of
#' \code{mcglm} class.
#'
#' @param x a fitted \code{mcglm} object.
#' @param type specify which graphical analysis will be performed.
#' Options are: \code{"residuals"} and \code{"algorithm"}.
#' @param ... additional arguments affecting the plot produced. Note
#' that there is no extra options for mcglm object class.
#'
#' @return The function \code{plot.mcglm} was designed to offer a fast
#' residuals analysis based on the Pearson residuals. Current version
#' offers a simple Pearson residuals versus fitted values and a quantile
#' plot. When using \code{algorithm = TRUE} the function will plot a summary of
#' the fitting algorithm shows the trajectory or iterations of the fitting
#' algorithm. The iterations are shown in terms of values of the model
#' parameters and also the actually value of the quasi-score and Pearson
#' estimating functions. Hence, a quickly check of the algorithm
#' convergence is obtained.
#' @seealso \code{residuals} and \code{fitted}.
#' @method plot mcglm
#' @export
plot.mcglm <- function(x, type = "residuals", ...) {
object <- x
n_resp <- length(object$beta_names)
if (type == "residuals") {
for (i in 1:n_resp) {
res <- residuals(object, type = "pearson")[, i]
fit_values <- fitted(object)[, i]
graphics::plot(res ~ fit_values, ylab = "Pearson residuals",
xlab = "Fitted values")
temp <-
stats::loess.smooth(fitted(object)[, i],
residuals(object, type = "pearson")[, i])
graphics::lines(temp$x, temp$y)
stats::qqnorm(res)
stats::qqline(res)
}
}
if (type == "algorithm") {
n_iter <- length(na.exclude(object$IterationCovariance[,
1]))
graphics::matplot(object$IterationRegression[1:c(n_iter + 5), ],
type = "l", lty = 2, ylab = "Regression",
xlab = "Iterations")
graphics::matplot(object$IterationCovariance[1:c(n_iter + 5), ],
type = "l", lty = 2, ylab = "Covariance",
xlab = "Iterations")
graphics::matplot(object$ScoreRegression[1:c(n_iter + 5), ], type = "l",
lty = 2, ylab = "Quasi-score Regression",
xlab = "Iterations")
graphics::matplot(object$ScoreCovariance[1:c(n_iter + 5), ], type = "l",
lty = 2, ylab = "Quasi-score Covariance",
xlab = "Iterations")
}
if (type == "partial_residuals") {
list_beta <- mc_updateBeta(list_initial = object$list_initial,
betas = object$Regression,
n_resp = n_resp,
information = object$Information)
comp_X <- list()
for (i in 1:n_resp) {
comp_X[[i]] <- as.matrix(object$list_X[[i]]) *
as.numeric(list_beta$regression[[i]])
}
for (i in 1:n_resp) {
res <- residuals(object, type = "pearson")[, i]
dev.new()
n_cov <- dim(comp_X[[i]])[2]
for (j in 2:n_cov) {
p1 <- comp_X[[i]][, j] + res
graphics::plot(p1 ~ object$list_X[[i]][, j],
xlab = object$beta_names[[i]][j],
ylab = "Partial residuals ")
}
}
}
}
#' @title Print
#' @name print.mcglm
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description The default print method for an object of \code{mcglm} class.
#'
#' @param x fitted model objects of class \code{mcglm} as produced by \code{mcglm()}.
#' @param ... further arguments passed to or from other methods.
#' @return No return value.
#'
#' @seealso \code{summary}.
#' @rdname print.mcglm
#' @method print mcglm
#' @export
print.mcglm <- function(x, ...) {
object <- x
n_resp <- length(object$beta_names)
regression <- mc_updateBeta(list_initial = list(),
betas = object$Regression,
information = object$Information,
n_resp = n_resp)
for (i in 1:n_resp) {
cat("Call: ")
print(object$linear_pred[[i]])
cat("\n")
cat("Link function:", object$link[[i]])
cat("\n")
cat("Variance function:", object$variance[[i]])
cat("\n")
cat("Covariance function:", object$covariance[[i]])
cat("\n")
names(regression[[1]][[i]]) <- object$beta_names[[i]]
cat("Regression:\n")
print(regression[[1]][[i]])
cat("\n")
cat("Dispersion:\n")
tau_temp <- coef(object, response = i,
type = "tau")$Estimate
names(tau_temp) <- rep("", length(tau_temp))
print(tau_temp)
cat("\n")
power_temp <- coef(object, response = i,
type = "power")$Estimate
if (length(power_temp) != 0) {
names(power_temp) <- ""
cat("Power:\n")
print(power_temp)
cat("\n")
}
}
}
#' @title Residuals
#' @name residuals.mcglm
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Compute residuals for an object of \code{mcglm} class.
#'
#' @param object an object of \code{mcglm} class.
#' @param type the type of residuals which should be returned.
#' Options are: \code{"raw"} (default), \code{"pearson"} and
#' \code{"standardized"}.
#' @param ... additional arguments affecting the residuals
#' produced. Note that there is no extra options for mcglm object
#' class.
#'
#' @return The function \code{residuals.mcglm} returns a matrix of
#' residuals values.
#' @seealso \code{fitted}.
#' @method residuals mcglm
#' @export
residuals.mcglm <- function(object, type = "raw", ...) {
n_resp <- length(object$beta_names)
output <- Matrix(object$residuals, ncol = n_resp, nrow = object$n_obs)
if (type == "standardized") {
output <- Matrix(
as.numeric(object$residuals %*% chol(object$inv_C)),
ncol = n_resp, nrow = object$n_obs)
}
if (type == "pearson") {
output <- Matrix(
as.numeric(object$residuals/sqrt(diag(object$C))),
ncol = n_resp, nrow = object$n_obs)
}
return(output)
}
#' @title Summarizing
#' @name summary.mcglm
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description The default summary method for an object of \code{mcglm} class.
#'
#' @param object an object of \code{mcglm} class.
#' @param verbose logical. Print or not the model summary.
#' @param print print only part of the model summary, options are
#' \code{Regression}, \code{power}, \code{Dispersion} and \code{Correlation}.
#' @param ... additional arguments affecting the summary produced. Note
#' the there is no extra options for mcglm object class.
#'
#' @return Print a \code{mcglm} object.
#' @seealso \code{print}.
#' @method summary mcglm
#' @export
summary.mcglm <- function(object, verbose = TRUE,
print = c("Regression","power","Dispersion","Correlation"),
...) {
n_resp <- length(object$beta_names)
output <- list()
for(i in 1:n_resp) {
tab_beta <- coef(object, std.error = TRUE,
response = i, type = "beta")[, 1:2]
tab_beta$"Z value" <- tab_beta[, 1]/tab_beta[, 2]
tab_beta$"Pr(>|z|)" <- 2*pnorm(-abs(tab_beta[, 1]/tab_beta[, 2]))
rownames(tab_beta) <- object$beta_names[[i]]
output[i][[1]]$Regression <- tab_beta
tab_power <- coef(object, std.error = TRUE,
response = i, type = "power")[, 1:2]
tab_power$"Z value" <- tab_power[, 1]/tab_power[, 2]
tab_power$"Pr(>|z|)" <- 2*pnorm(-abs(tab_power[, 1]/tab_power[, 2]))
rownames(tab_power) <- NULL
if (dim(tab_power)[1] != 0) {
output[i][[1]]$Power <- tab_power
}
tab_tau <- coef(object, std.error = TRUE,
response = i, type = "tau")[, 1:2]
tab_tau$"Z value" <- tab_tau[, 1]/tab_tau[, 2]
tab_tau$"Pr(>|z|)" <- 2*pnorm(-abs(tab_tau[, 1]/tab_tau[, 2]))
rownames(tab_tau) <- NULL
output[i][[1]]$tau <- tab_tau
}
tab_rho <- coef(object, std.error = TRUE,
response = NA, type = "correlation")[, c(3, 1, 2)]
tab_rho$"Z value" <- tab_rho[, 2]/tab_rho[, 3]
tab_rho$"Pr(>|z|)" <- 2*pnorm(-abs(tab_rho[, 2]/tab_rho[, 3]))
if (dim(tab_rho)[1] != 0) {
rownames(tab_rho) <- NULL
output$Correlation <- tab_rho
}
if(verbose == TRUE) {
for (i in 1:n_resp) {
if("Regression" %in% print) {
cat("Call: ")
print(object$linear_pred[[i]])
cat("\n")
cat("Link function:", object$link[[i]])
cat("\n")
cat("Variance function:", object$variance[[i]])
cat("\n")
cat("Covariance function:", object$covariance[[i]])
cat("\n")
cat("Regression:\n")
print(output[i][[1]]$Regression)
cat("\n")
}
if (dim(tab_power)[1] != 0) {
if("power" %in% print) {
cat("Power:\n")
print(output[i][[1]]$Power)
cat("\n")
}
}
if("Dispersion" %in% print) {
cat("Dispersion:\n")
print(output[i][[1]]$tau)
cat("\n")
}
}
if (dim(tab_rho)[1] != 0) {
if("Correlation" %in% print) {
cat("Correlation matrix:\n")
print(tab_rho)
cat("\n")
}
}
names(object$con$correct) <- ""
iteration_cov <- length(na.exclude(object$IterationCovariance[, 1]))
names(iteration_cov) <- ""
names(object$con$method) <- ""
cat("Algorithm:", object$con$method)
cat("\n")
cat("Correction:", object$con$correct)
cat("\n")
cat("Number iterations:", iteration_cov)
cat("\n")
}
if (dim(tab_rho)[1] != 0) {
names(output) <- c(paste("Resp.Variable", 1:n_resp),"Correlation")
}
if (dim(tab_rho)[1] == 0) {
names(output) <- paste("Resp.Variable", 1:n_resp)
}
return(invisible(output))
}
#' @title Variance-Covariance Matrix
#' @name vcov.mcglm
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Returns the variance-covariance matrix for an object
#' of \code{mcglm} class.
#'
#' @param object an object of \code{mcglm} class.
#' @param ... additional arguments affecting the summary produced. Note
#' that there is no extra options for mcglm object class.
#'
#' @return A variance-covariance matrix.
#'
#' @method vcov mcglm
#' @export
vcov.mcglm <- function(object, ...) {
cod <- coef(object)$Parameters
colnames(object$vcov) <- cod
rownames(object$vcov) <- cod
return(object$vcov)
}
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.