#' @title Print
#' @name print.mglm4twin
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description The default print method for an object of \code{mglm4twin} class.
#'
#' @param x fitted model objects of class \code{mglm4twin} as produced
#' by \code{mglm4twin()}.
#' @param ... further arguments passed to or from other methods.
#'
#' @seealso \code{summary}.
#' @rdname print.mglm4twin
#' @method print mglm4twin
#' @export
print.mglm4twin <- function(x, ...) {
object <- x
n_resp <- length(object$beta_names)
regression <- mt_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")
names(regression[[1]][[i]]) <- object$beta_names[[i]]
cat("Regression:\n")
print(regression[[1]][[i]])
cat("\n")
}
cat("Dispersion:\n")
tau_temp <- object$Covariance
names(tau_temp) <- rep("", length(tau_temp))
print(tau_temp)
cat("\n")
}
#' @title Model Coefficients
#' @name coef.mglm4twin
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Extract model coefficients for objects
#' of \code{mglm4twin} class.
#'
#' @param object an object of \code{mglm4twin} class.
#' @param std.error logical. If \code{TRUE} returns the standard errors
#' for the estimates. Default is \code{FALSE}.
#' @param model Type of the fitted model. Options are \code{E}, \code{AE},
#' \code{CE}, \code{ACE} and \code{ADE}.
#' @param response For which response regression coefficients are extracted.
#' Default = NULL. It means all estimates are returned including
#' dispersion estimates for all responses.
#' @param ... additional arguments affecting the summary produced. Note
#' that there is no extra options for \code{mglm4twin} object class.
#'
#' @return A \code{data.frame} with parameters names, estimates,
#' response variable number and parameters type.
#' @importFrom stats pnorm deriv
#' @method coef mglm4twin
#' @export
coef.mglm4twin <- function(object, std.error = FALSE, model,
response = NULL, ...) {
n_resp <- length(object$beta_names)
Estimates <- c(object$Regression, object$Covariance)
cod_beta <- list()
beta_num <- list()
power_code <- list()
for(i in 1:n_resp) {
cod_beta[[i]] <- rep(i, object$Information$n_betas[[i]])
beta_num[[i]] <- rep(0:c(object$Information$n_betas[[i]]-1))
if(object$Information$n_power[[i]] != 0) {
power_code[[i]] <- rep(i, object$Information$n_power[[i]])
}
}
cod_beta <- unlist(cod_beta)
beta_nul <- unlist(beta_num)
beta_nome <- paste0("beta", paste0(cod_beta, beta_nul))
if(length(power_code) != 0) {
power_nome <- paste0("power", unlist(power_code))
}
# Genetic components
E_main <- paste0("E", 1:n_resp)
A_main <- paste0("A", 1:n_resp)
C_main <- paste0("C", 1:n_resp)
D_main <- paste0("D", 1:n_resp)
if(n_resp > 1) {
idx <- combn(n_resp, 2)
E_cross <- list()
A_cross <- list()
C_cross <- list()
D_cross <- list()
for(i in 1:dim(idx)[2]) {
E_cross[[i]] <- paste0("E", paste0(idx[1,i],idx[2,i]))
A_cross[[i]] <- paste0("A", paste0(idx[1,i],idx[2,i]))
C_cross[[i]] <- paste0("C", paste0(idx[1,i],idx[2,i]))
D_cross[[i]] <- paste0("D", paste0(idx[1,i],idx[2,i]))
}
E_cross <- unlist(E_cross)
A_cross <- unlist(A_cross)
C_cross <- unlist(C_cross)
D_cross <- unlist(D_cross)
}
if(model == "E") {
if(n_resp > 1) {
cod_cov <- c(E_main, E_cross)
}
if(n_resp == 1) {
cod_cov <- E_main
}
}
if(model == "AE") {
if(n_resp > 1) {
cod_cov <- c(E_main, E_cross, A_main, A_cross)
}
if(n_resp == 1) {
cod_cov <- c(E_main, A_main)
}
}
if(model == "CE") {
if(n_resp > 1) {
cod_cov <- c(E_main, E_cross, C_main, C_cross)
}
if(n_resp == 1) {
cod_cov <- c(E_main, C_main)
}
}
if(model == "ACE") {
if(n_resp > 1) {
cod_cov <- c(E_main, E_cross, A_main, A_cross, C_main, C_cross)
}
if(n_resp == 1) {
cod_cov <- c(E_main, A_main, C_main)
}
}
if(model == "ADE") {
if(n_resp > 1) {
cod_cov <- c(E_main, E_cross, A_main, A_cross, D_main, D_cross)
}
if(n_resp == 1) {
cod_cov <- c(E_main, A_main, D_main)
}
}
if(length(power_code) != 0) {
Parameters <- c(beta_nome, power_nome, cod_cov)
}
if(length(power_code) == 0) {
Parameters <- c(beta_nome, cod_cov)
}
coef_temp <- data.frame(Parameters = Parameters,
Estimates = Estimates)
if (std.error == TRUE) {
coef_temp$std.error = sqrt(diag(object$vcov))
}
if(!is.null(response)) {
coef_temp$Response <- c(cod_beta, rep(NA, object$Information$n_cov))
coef_temp <- coef_temp[which(coef_temp$Response %in% response), ]
}
return(coef_temp)
}
#' @title Summarizing
#' @name summary.mglm4twin
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description The default summary method for an object of \code{mglm4twin} class.
#'
#' @param object an object of \code{mglm4twin} class.
#' @param model String indicating twin model.
#' @param biometric Logical. If TRUE biometric measures such as heritability,
#' common environment, genetic correlation etc are printed. (Default = FALSE).
#' @param ... additional arguments affecting the summary produced. Note
#' the there is no extra options for mglm4twin object class.
#'
#' @return Print a \code{mglm4twin} object.
#' @seealso \code{print}.
#' @method summary mglm4twin
#' @export
summary.mglm4twin <- function(object, model, biometric = FALSE, ...) {
n_resp <- length(object$beta_names)
out <- coef(object, std.error = TRUE, model = model)
VCOV <- vcov(object, model = model)
if(n_resp > 1) {
cross_terms <- combn(n_resp, 2)
}
output <- list()
for(i in 1:n_resp) {
tab_beta <- coef(object, std.error = TRUE, model = model,
response = i, type = "beta")[, 2:3]
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
}
out$"z value" <- out[, 2]/out[, 3]
out$"Pr(>|z|)" <- 2*pnorm(-abs(out[, 2]/out[, 3]))
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("Regression:\n")
print(output[i][[1]]$Regression)
cat("\n")
if(object$Information$n_power[[i]] != 0) {
cat("Power:\n")
tt = paste0("power",i)
rownames(out)[as.numeric(rownames(out[which(out$Parameters == tt),-1]))] <- tt
print(out[which(out$Parameters == tt),-1])
cat("\n")
}
}
n_beta <- sum(unlist(object$Information$n_betas))
n_power <- sum(unlist(object$Information$n_power))
n_cov <- object$Information$n_taus
n_comp <- n_resp + n_resp*(n_resp-1)/2
n_total <- dim(out)[1]
out_cov <- out[c(n_beta+n_power+1):n_total,]
rownames(out_cov) <- out_cov$Parameters
out_cov$"z value" <- out_cov[, 2]/out_cov[, 3]
out_cov$"Pr(>|z|)" <- 2*pnorm(-abs(out_cov[, 2]/out_cov[, 3]))
output$Dispersion <- out_cov
cat("Dispersion:\n")
if(model == "E") {
cat("Environment component:\n")
print(out_cov[1:c(n_comp),-1])
cat("\n")
}
if(model == "AE") {
cat("Genetic component:\n")
print(out_cov[c(n_comp+1):c(2*n_comp),-1])
cat("\n")
cat("Environment component:\n")
print(out_cov[1:c(n_comp),-1])
cat("\n")
}
if(model == "CE") {
cat("Common environment component:\n")
print(out_cov[c(n_comp+1):c(2*n_comp),-1])
cat("\n")
cat("Environment component:\n")
print(out_cov[1:c(n_comp),-1])
cat("\n")
}
if(model == "ACE") {
cat("Genetic component:\n")
print(out_cov[c(n_comp+1):c(2*n_comp),-1])
cat("\n")
cat("Common environment component:\n")
print(out_cov[c(2*n_comp+1):c(3*n_comp),-1])
cat("\n")
cat("Environment component:\n")
print(out_cov[1:c(n_comp),-1])
cat("\n")
}
if(model == "ADE") {
cat("Genetic component:\n")
print(out_cov[c(n_comp+1):c(2*n_comp),-1])
cat("\n")
cat("D component:\n")
print(out_cov[c(2*n_comp+1):c(3*n_comp),-1])
cat("\n")
cat("Environment component:\n")
print(out_cov[1:c(n_comp),-1])
cat("\n")
}
if(biometric == TRUE) {
if(model == "E") {
tab_rho_e <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "E", n_resp = n_resp)
if(n_resp > 1) {
output$E_cross <- tab_rho_e
}
cat("Environmental correlation:\n")
print(tab_rho_e)
cat("\n")
}
if(model == "AE") {
tab_rho_a <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "A", n_resp = n_resp)
tab_rho_e <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "E", n_resp = n_resp)
gen <- mt_compute_gen(Estimates = out, vcov = VCOV,
model = model, n_resp = n_resp)
# Environment main effects
Env_main <- data.frame("Estimates" = diag(gen[[1]]),
"std.error" = diag(gen[[3]]))
Env_main$"z value" <- Env_main[, 1]/Env_main[, 2]
Env_main$"Pr(>|z|)" <- 2*pnorm(-abs(Env_main[, 1]/Env_main[, 2]))
rownames(Env_main) <- paste0("e", 1:n_resp)
output$Env_main <- Env_main
# Environment cross effects
Env_cross <- "Bivariate environmentality no available."
if(n_resp > 1) {
Env_cross <- data.frame("Estimates" = gen[[1]][upper.tri(gen[[1]])],
"std.error" = gen[[3]][upper.tri(gen[[3]])])
Env_cross$"z value" <- Env_cross[, 1]/Env_cross[, 2]
Env_cross$"Pr(>|z|)" <- 2*pnorm(-abs(Env_cross[, 1]/Env_cross[, 2]))
rownames(Env_cross) <- paste0("e", cross_terms[1,], cross_terms[2,])
output$Env_cross <- Env_cross
}
# Genetic main effects
A_main <- data.frame("Estimates" = diag(gen[[2]]),
"std.error" = diag(gen[[4]]))
A_main$"z value" <- A_main[, 1]/A_main[, 2]
A_main$"Pr(>|z|)" <- 2*pnorm(-abs(A_main[, 1]/A_main[, 2]))
rownames(A_main) <- paste0("h", 1:n_resp)
output$A_main <- A_main
# Genetic cross effects
A_cross <- "Bivariate heritability no available."
if(n_resp > 1) {
A_cross <- data.frame("Estimates" = gen[[2]][upper.tri(gen[[2]])],
"std.error" = gen[[4]][upper.tri(gen[[4]])])
A_cross$"z value" <- A_cross[, 1]/A_cross[, 2]
A_cross$"Pr(>|z|)" <- 2*pnorm(-abs(A_cross[, 1]/A_cross[, 2]))
rownames(A_cross) <- paste0("h", cross_terms[1,], cross_terms[2,])
output$A_cross <- A_cross
}
cat("Measures associated with genetic structure:\n")
cat("Genetic correlation:\n")
print(tab_rho_a)
cat("\n")
cat("Heritability:\n")
print(A_main)
cat("\n")
cat("Bivariate heritability:\n")
print(A_cross)
cat("\n")
cat("Measures associated with environment structure:\n")
cat("Environment correlation:\n")
print(tab_rho_e)
cat("\n")
cat("Environmentality:\n")
print(Env_main)
cat("\n")
cat("Bivariate environmentality:\n")
print(Env_cross)
cat("\n")
}
if(model == "CE") {
tab_rho_c <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "C", n_resp = n_resp)
tab_rho_e <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "E", n_resp = n_resp)
gen <- mt_compute_gen(Estimates = out, vcov = VCOV,
model = model, n_resp = n_resp)
# Environment main effects
E_main <- data.frame("Estimates" = diag(gen[[1]]),
"std.error" = diag(gen[[3]]))
E_main$"z value" <- E_main[, 1]/E_main[, 2]
E_main$"Pr(>|z|)" <- 2*pnorm(-abs(E_main[, 1]/E_main[, 2]))
rownames(E_main) <- paste0("e", 1:n_resp)
output$E_main <- E_main
# Environment cross effects
E_cross <- "Bivariate environmentality no available."
if(n_resp > 1) {
E_cross <- data.frame("Estimates" = gen[[1]][upper.tri(gen[[1]])],
"std.error" = gen[[3]][upper.tri(gen[[3]])])
E_cross$"z value" <- E_cross[, 1]/E_cross[, 2]
E_cross$"Pr(>|z|)" <- 2*pnorm(-abs(E_cross[, 1]/E_cross[, 2]))
rownames(E_cross) <- paste0("e", cross_terms[1,], cross_terms[2,])
output$E_cross <- E_cross
}
# Common environment main effects
C_main <- data.frame("Estimates" = diag(gen[[2]]),
"std.error" = diag(gen[[4]]))
C_main$"z value" <- C_main[, 1]/C_main[, 2]
C_main$"Pr(>|z|)" <- 2*pnorm(-abs(C_main[, 1]/C_main[, 2]))
rownames(C_main) <- paste0("c", 1:n_resp)
output$C_main <- C_main
# Common environment cross effects
C_cross <- "Bivariate common environmentality no available."
if(n_resp > 1) {
C_cross <- data.frame("Estimates" = gen[[2]][upper.tri(gen[[2]])],
"std.error" = gen[[4]][upper.tri(gen[[4]])])
C_cross$"z value" <- C_cross[, 1]/C_cross[, 2]
C_cross$"Pr(>|z|)" <- 2*pnorm(-abs(C_cross[, 1]/C_cross[, 2]))
rownames(C_cross) <- paste0("c", cross_terms[1,], cross_terms[2,])
output$C_cross <- C_cross
}
cat("Measures associated with common environment structure:\n")
cat("Common environment correlation:\n")
print(tab_rho_c)
cat("\n")
cat("Common environment:\n")
print(C_main)
cat("\n")
cat("Bivariate common environmentality:\n")
print(C_cross)
cat("\n")
cat("Measures associated with environment structure:\n")
cat("Environment correlation:\n")
print(tab_rho_e)
cat("\n")
cat("Environmentality:\n")
print(E_main)
cat("\n")
cat("Bivariate environmentality:\n")
print(E_cross)
cat("\n")
}
if(model == "ACE") {
tab_rho_a <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "A", n_resp = n_resp)
tab_rho_c <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "C", n_resp = n_resp)
tab_rho_e <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "E", n_resp = n_resp)
gen <- mt_compute_gen(Estimates = out, vcov = VCOV,
model = model, n_resp = n_resp)
# Environment main effects
E_main <- data.frame("Estimates" = diag(gen[[1]]),
"std.error" = diag(gen[[4]]))
E_main$"z value" <- E_main[, 1]/E_main[, 2]
E_main$"Pr(>|z|)" <- 2*pnorm(-abs(E_main[, 1]/E_main[, 2]))
rownames(E_main) <- paste0("e", 1:n_resp)
output$E_main <- E_main
# Environment cross effects
E_cross <- "Bivariate environmetality no available."
if(n_resp > 1) {
E_cross <- data.frame("Estimates" = gen[[1]][upper.tri(gen[[1]])],
"std.error" = gen[[4]][upper.tri(gen[[4]])])
E_cross$"z value" <- E_cross[, 1]/E_cross[, 2]
E_cross$"Pr(>|z|)" <- 2*pnorm(-abs(E_cross[, 1]/E_cross[, 2]))
rownames(E_cross) <- paste0("e", cross_terms[1,], cross_terms[2,])
output$E_cross <- E_cross
}
# Genetic main effects
A_main <- data.frame("Estimates" = diag(gen[[2]]),
"std.error" = diag(gen[[5]]))
A_main$"z value" <- A_main[, 1]/A_main[, 2]
A_main$"Pr(>|z|)" <- 2*pnorm(-abs(A_main[, 1]/A_main[, 2]))
rownames(A_main) <- paste0("h", 1:n_resp)
output$A_main <- A_main
# Genetic cross effects
A_cross <- "Bivariate heritability no available."
if(n_resp > 1) {
A_cross <- data.frame("Estimates" = gen[[2]][upper.tri(gen[[2]])],
"std.error" = gen[[5]][upper.tri(gen[[5]])])
A_cross$"z value" <- A_cross[, 1]/A_cross[, 2]
A_cross$"Pr(>|z|)" <- 2*pnorm(-abs(A_cross[, 1]/A_cross[, 2]))
rownames(A_cross) <- paste0("h", cross_terms[1,], cross_terms[2,])
output$A_cross <- A_cross
}
# Commmon environment main effects
C_main <- data.frame("Estimates" = diag(gen[[3]]),
"std.error" = diag(gen[[6]]))
C_main$"z value" <- C_main[, 1]/C_main[, 2]
C_main$"Pr(>|z|)" <- 2*pnorm(-abs(C_main[, 1]/C_main[, 2]))
rownames(C_main) <- paste0("c", 1:n_resp)
output$C_main <- C_main
# Common environment cross effects
C_cross <- "Bivariate common environmentality no available."
if(n_resp > 1) {
C_cross <- data.frame("Estimates" = gen[[3]][upper.tri(gen[[3]])],
"std.error" = gen[[6]][upper.tri(gen[[6]])])
C_cross$"z value" <- C_cross[, 1]/C_cross[, 2]
C_cross$"Pr(>|z|)" <- 2*pnorm(-abs(C_cross[, 1]/C_cross[, 2]))
rownames(C_cross) <- paste0("c", cross_terms[1,], cross_terms[2,])
output$C_cross <- C_cross
}
cat("Measures associated with genetic structure:\n")
cat("Genetic correlation:\n")
print(tab_rho_a)
cat("\n")
cat("Heritability:\n")
print(A_main)
cat("\n")
cat("Bivariate heritability:\n")
print(A_cross)
cat("\n")
cat("Measures associated with common environment structure:\n")
cat("Common environment correlation:\n")
print(tab_rho_c)
cat("\n")
cat("Common environment:\n")
print(C_main)
cat("\n")
cat("Bivariate common environment:\n")
print(C_cross)
cat("\n")
cat("Measures associated with environment structure:\n")
cat("Environment correlation:\n")
print(tab_rho_e)
cat("\n")
cat("Environmentality:\n")
print(E_main)
cat("\n")
cat("Bivariate environmentality:\n")
print(E_cross)
cat("\n")
}
if(model == "ADE") {
tab_rho_a <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "A", n_resp = n_resp)
tab_rho_d <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "D", n_resp = n_resp)
tab_rho_e <- mt_compute_rho(Estimates = out, vcov = VCOV,
component = "E", n_resp = n_resp)
gen <- mt_compute_gen(Estimates = out, vcov = VCOV,
model = model, n_resp = n_resp)
# Environment main effects
E_main <- data.frame("Estimates" = diag(gen[[1]]),
"std.error" = diag(gen[[4]]))
E_main$"z value" <- E_main[, 1]/E_main[, 2]
E_main$"Pr(>|z|)" <- 2*pnorm(-abs(E_main[, 1]/E_main[, 2]))
rownames(E_main) <- paste0("e", 1:n_resp)
output$E_main <- E_main
# Environment cross effects
E_cross <- "Bivariate environmentality no available."
if(n_resp > 1) {
E_cross <- data.frame("Estimates" = gen[[1]][upper.tri(gen[[1]])],
"std.error" = gen[[4]][upper.tri(gen[[4]])])
E_cross$"z value" <- E_cross[, 1]/E_cross[, 2]
E_cross$"Pr(>|z|)" <- 2*pnorm(-abs(E_cross[, 1]/E_cross[, 2]))
rownames(E_cross) <- paste0("e", cross_terms[1,], cross_terms[2,])
output$E_cross <- E_cross
}
# Genetic main effects
A_main <- data.frame("Estimates" = diag(gen[[2]]),
"std.error" = diag(gen[[5]]))
A_main$"z value" <- A_main[, 1]/A_main[, 2]
A_main$"Pr(>|z|)" <- 2*pnorm(-abs(A_main[, 1]/A_main[, 2]))
rownames(A_main) <- paste0("h", 1:n_resp)
output$A_main <- A_main
# Genetic cross effects
A_cross <- "Bivariate heritability no available."
if(n_resp > 1) {
A_cross <- data.frame("Estimates" = gen[[2]][upper.tri(gen[[2]])],
"std.error" = gen[[5]][upper.tri(gen[[5]])])
A_cross$"z value" <- A_cross[, 1]/A_cross[, 2]
A_cross$"Pr(>|z|)" <- 2*pnorm(-abs(A_cross[, 1]/A_cross[, 2]))
rownames(A_cross) <- paste0("h", cross_terms[1,], cross_terms[2,])
output$A_cross <- A_cross
}
# D main effects
D_main <- data.frame("Estimates" = diag(gen[[3]]),
"std.error" = diag(gen[[6]]))
D_main$"z value" <- D_main[, 1]/D_main[, 2]
D_main$"Pr(>|z|)" <- 2*pnorm(-abs(D_main[, 1]/D_main[, 2]))
rownames(D_main) <- paste0("d", 1:n_resp)
output$D_main <- D_main
# D cross effects
D_cross <- "Bivariate D components no available."
if(n_resp > 1) {
D_cross <- data.frame("Estimates" = gen[[3]][upper.tri(gen[[3]])],
"std.error" = gen[[6]][upper.tri(gen[[6]])])
D_cross$"z value" <- D_cross[, 1]/D_cross[, 2]
D_cross$"Pr(>|z|)" <- 2*pnorm(-abs(D_cross[, 1]/D_cross[, 2]))
rownames(D_cross) <- paste0("d", cross_terms[1,], cross_terms[2,])
output$D_cross <- D_cross
}
cat("Measures associated with genetic structure:\n")
cat("Genetic correlation:\n")
print(tab_rho_a)
cat("\n")
cat("Heritability:\n")
print(A_main)
cat("\n")
cat("Bivariate heritability:\n")
print(A_cross)
cat("\n")
cat("Measures associated with D structure:\n")
cat("D correlation:\n")
print(tab_rho_d)
cat("\n")
cat("D:\n")
print(D_main)
cat("\n")
cat("Bivariate D:\n")
print(D_cross)
cat("\n")
cat("Measures associated with environment structure:\n")
cat("Environment correlation:\n")
print(tab_rho_e)
cat("\n")
cat("Environmentality:\n")
print(E_main)
cat("\n")
cat("Bivariate environmentality:\n")
print(E_cross)
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")
return(invisible(output))
}
#' @title Variance-Covariance Matrix
#' @name vcov.mglm4twin
#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
#'
#' @description Returns the variance-covariance matrix for an object
#' of \code{mglm4twin} class.
#'
#' @param object an object of \code{mglm4twin} class.
#' @param model String indicating twin model.
#' @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 mglm4twin
#' @export
vcov.mglm4twin <- function(object, model, ...) {
cod <- coef(object, model = model)$Parameters
colnames(object$vcov) <- cod
rownames(object$vcov) <- cod
return(object$vcov)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.