Remove the following functions

- mc_fast_forward.R
- mc_influence.R
- mc_qic.R
- mc_qll.R
- mc_rw1.R
- mc_rw2.R
- mc_unstructured.R

R/mc_fast_forward.R

-#' Fast forward selection for multivariate covariance generalized linear
-#' models.
-#'
-#' @description Perform fast forward model selection using the score
-#' information criterion. This function works only for univariate months.
-#'
-#' @param object an object representing a model of \code{mcglm} class.
-#' @param scope a vector specyfing the covariate to be tested.
-#' @param interaction Maximum number of covariates interacting.
-#' @param penalty penalty term (default = 2).
-#' @param n_max Maximum number of models to be fitted.
-#' @return The selected model.
-#' @export
-
-mc_fast_forward <- function(object, scope, interaction = 1,
-                            penalty = 2, n_max = 10) {
-  if (interaction > 1) {
-    int_terms <- list()
-    for (i in 2:interaction) {
-      int_terms[[c(i-1)]] <- combn(length(scope), i)
-    }
-    fun_temp <- function(int_terms) {
-      output <- c()
-      for(i in 1:dim(int_terms)[2]) {
-        output[i] <- paste(scope[int_terms[,i]],collapse = "*")
-        }
-      return(output)
-      }
-  scope <- do.call(c, lapply(int_terms, fun_temp))
-  }
-  for (i in 1:n_max) {
-    if(length(scope) == 0) break
-    sic <- mc_sic(object = object, scope = scope,
-                  data = data, response = 1)
-    if (all(sic$SIC < 0)) break
-    cov_new <- as.numeric(rownames(sic[which.max(sic$SIC),]))
-    cov_enter <- scope[as.numeric(rownames(sic[which.max(sic$SIC),]))]
-    next_cov <- paste("~. +", cov_enter)
-    new_formula <- update.formula(object$linear_pred[[1]], next_cov)
-    temp_models <- mcglm(c(new_formula),
-                       matrix_pred = object$matrix_pred,
-                       link = object$link, variance = object$variance,
-                       covariance = object$covariance, data = data)
-    object <- temp_models
-    scope <- scope[-cov_new]
-    cat(paste("Iteration", i), "\n")
-  }
-  return(object)
-}

R/mc_influence.R

-#' @title Influence measures for McGLMs
-#'
-#' @description Compute influence measures for multivariate covariance
-#' generalized linear models. Leverage, DFBETA and Cook's distance
-#' for unit sample and observations.
-#' @param object An object of \code{mcglm} class.
-#' @param id a vector which identifies the clusters.
-#' The length and order of id should match with the number of
-#' observations.
-#' @return A list with influence measures for cluster and observations.
-#'
-#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
-#' @export
-
-mc_influence <- function(object, id) {
-  temp_data <- data.frame(object$residuals, id)
-  names(temp_data) <- c("res", "id")
-  data_id <- split(temp_data, temp_data$id)
-  n_id <- length(data_id)
-  D <- bdiag(lapply(object$mu_list, function(x) x$D))
-  n_beta <- dim(D)[2]
-  R <- bdiag(lapply(data_id, function(x) {
-    tcrossprod(x[, 1])
-  }))
-  # Leverage for observations
-  inv_M <- -object$inv_S_beta
-  M <- solve(inv_M)
-  uni_id <- as.character(unique(id))
-  Hi <- list()
-  Di <- list()
-  inv_Ci <- list()
-  for (i in 1:n_id) {
-    idTF <- id == uni_id[i]
-    if(sum(idTF) == 1) {
-      Di[[i]] <- Matrix(D[idTF,], nrow = sum(idTF), ncol = dim(D)[2])
-    } else {
-      Di[[i]] <- D[idTF,]
-    }
-    inv_Ci[[i]] <- object$inv_C[idTF,idTF]
-    Hi[[i]] <- Di[[i]]%*%inv_M%*%t(Di[[i]])%*%inv_Ci[[i]]
-  }
-  H <- bdiag(Hi)
-  # Leverage for observations
-  Hobs <- diag(H)
-  # Leverage for clusters
-  Hid <- tapply(Hobs, id, mean)
-  # DFBETAs for clusters
-  I <- Diagonal(dim(temp_data)[1], 1)
-  inv_IH <- solve(I - H)
-  DFBETACi <- list()
-  DFBETAOij <- list()
-  DCLSi <- c()
-  for(i in 1:n_id) {
-    idTF <- id == uni_id[i]
-    if(sum(idTF) == 1) {
-      D_i <- Matrix(D[idTF,], nrow = sum(idTF), ncol = n_beta)
-      inv_C_i <- Matrix(object$inv_C[idTF,idTF],ncol = 1, nrow = 1)
-    } else {
-      D_i <- D[idTF,]
-      inv_C_i <- object$inv_C[idTF,idTF]
-    }
-    inv_IH_i = inv_IH[idTF,idTF]
-    r_i <- object$residuals[idTF]
-    DFBETACi[[i]] <- t(inv_M%*%t(D_i)%*%inv_C_i%*%inv_IH_i%*%r_i)
-    DCLSi[i] <- as.numeric(DFBETACi[[i]]%*%M%*%t(DFBETACi[[i]]))/n_beta
-    if(dim(inv_C_i)[1] == 1) {DFBETAOij[[i]] <- DFBETACi[[i]]
-    } else {
-    DFBETAOij[[i]] <- mc_dfbetaOij(Di = D_i, Ci = object$C[idTF,idTF],
-                                   inv_Ci = inv_C_i, ri = r_i,
-                                   inv_M = inv_M)
-    }
-  }
-  DFBETACi <- do.call(rbind, DFBETACi)
-  DFBETAOij <- do.call(rbind, DFBETAOij)
-  DCOij <- apply(DFBETAOij, MARGIN = 1,
-                function(x,M,n_beta)as.numeric((t(x)%*%M%*%x)/n_beta),
-                M = M, n_beta = n_beta)
-  std.error <- coef(object, type = "beta", std.error = TRUE)$Std.error
-  DFBETACi = data.frame(as.matrix(DFBETACi))/std.error
-  DFBETAOij = data.frame(as.matrix(DFBETAOij))/std.error
-  names(DFBETACi) <- do.call(c, object$beta_names)
-  names(DFBETAOij) <- do.call(c, object$beta_names)
-  DFBETACi$Leverage <- as.numeric(Hid)
-  DFBETACi$Cook <- as.numeric(DCLSi)
-  DFBETAOij$Leverage <- as.numeric(Hobs)
-  DFBETAOij$Cook <- as.numeric(DCOij)
-  DFBETAOij$Id <- rownames(object$data)
-  DFBETACi$Id <- uni_id
-  output <- list("Id" = DFBETACi, "Obs" = DFBETAOij)
-  return(output)
-}
-
-#' @title Influence measures for McGLMs
-#'
-#' @description Compute influence measures for multivariate covariance
-#' generalized linear models. Leverage, DFBETA and Cook's distance
-#' for observations. Auxiliar function for \code{mc_influence}.
-#' @param Di D matrix for the cluster i.
-#' @param Ci C matrix for the cluster i.
-#' @param inv_Ci Inverse of C matrix for the cluster i.
-#' @param ri Residual vector for the cluster i.
-#' @param inv_M Inverse of variance/covariance of regression parameters.
-#' @return Matrix with the DFBETA for observation in the cluster i.
-#'
-#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
-#' @export
-
-mc_dfbetaOij <- function(Di, Ci, inv_Ci, ri, inv_M) {
-  nobs <- dim(inv_Ci)[1]
-  nbeta <- dim(Di)[2]
-  DFBETAOij <- matrix(NA, ncol = nbeta, nrow = nobs)
-  for (i in 1:nobs) {
-    Ci.jj <- Matrix(Ci[-i,i])
-    Cij.j <- t(Ci.jj)
-    inv_Ci.j <- inv_Ci[-i,-i]
-    Di.j <- Di[-i,]
-    Dij_temp <- Di[i,] - Cij.j%*%inv_Ci.j%*%Di.j
-    rij_temp <- as.numeric(ri[i] - Cij.j%*%inv_Ci.j%*%ri[-i])
-    Cij_temp <- as.numeric(Ci[i,i] - Cij.j%*%inv_Ci.j%*%t(Cij.j))
-    Hij_temp <- as.numeric(Dij_temp%*%inv_M%*%t(Dij_temp)%*%solve(Cij_temp))
-    part2 <- Matrix(rij_temp/(Cij_temp*(1-Hij_temp)))
-    DFBETAOij[i,] <- as.numeric(t(inv_M%*%t(Dij_temp)%*%part2))
-  }
-  return(DFBETAOij)
-}
-
-
-

R/mc_qic.R

-#' @title Compute Quasi Information Criterion (QIC) for McGLMs.
-#' @name qic.mcglm
-#'
-#' @description \code{qic.mcglm} is a function which computes the QIC
-#'     for McGLMs.
-#'
-#' @param object An object of \code{mcglm} class.
-#' @param object.iid An object of \code{mcglm} class contained the model
-#'     fitted using independent covariance structure.
-#'
-#' @return The QIC value.
-#'
-#' @author Wagner Hugo Bonat, \email{wbonat@@ufpr.br}
-#'
-#' @export
-
-qic.mcglm <- function(object, object.iid) {
-    mu <- fitted(object)
-    obs <- object$observed
-    n_resp <- dim(mu)[2]
-    Q <- matrix(NA, ncol = dim(mu)[2], nrow = dim(mu)[1])
-    for (i in 1:n_resp) {
-        if (object$power_fixed[[i]] == FALSE) {
-            power <- coef(object, type = "power", response = i)$Estimate
-        }
-        if (object$power_fixed[[i]] == TRUE) {
-            power = object$list_initial$power[[i]]
-        }
-        Q[, i] <- mc_qll(y = obs[, i],
-                         mu = mu[, i],
-                         variance = object$variance[[i]],
-                         power = power)
-    }
-    Vbeta <- -object$inv_S_beta
-    Vnull <- solve(-object.iid$inv_S_beta)
-    t1 <- -2 * sum(Q)
-    qic <- t1 + 2 * sum(diag(Vnull %*% Vbeta))
-    return(list(Q = t1, qic = qic))
-}

R/mc_qll.R

-#' Compute quasi-likelihood function.
-#'
-#' Given a variance function mc_qll function computes the quasi-likelihood values.
-#' @param y A vector of observed values.
-#' @param mu A vector of fitted values.
-#' @param variance Variance function (constant, tweedie, poisson_tweedie, binomial).
-#' @param power Power parameter value.
-#' @return The quasi-likelihood values.
-#' @export
-
-mc_qll <- function(y, mu, variance, power) {
-    if (variance == "constant") {
-        qll <- -((y - mu)^2)/2  # Gaussian case
-    }
-    if (variance == "tweedie" & power == 1) {
-        qll <- y * log(mu) - mu  # Poisson case
-    }
-    if (variance == "tweedie" & power == 2) {
-        -y/mu - log(mu)  # Gamma case
-    }
-    if (variance == "tweedie" & power != 1 & power != 2) {
-        qll <- (mu^-power) * ((mu * y)/(1 - power) - (mu^2)/(2 - power))  # General Tweedie case
-    }
-    if (variance == "poisson_tweedie" & power == 1) {
-        qll <- (y * log(mu) - mu) + (y * log(mu) - mu)
-    }
-    if (variance == "poisson_tweedie" & power == 2) {
-        qll <- (y * log(mu) - mu) + (-y/mu - log(mu))
-    }
-    if (variance == "poisson_tweedie" & power != 1 & power != 2) {
-        qll <- (y * log(mu) - mu) + (mu^-power) * ((mu * y)/(1 - power) - (mu^2)/(2 - power))
-    }
-    if (variance == "binomial") {
-        qll <- y * log(mu/(1 - mu)) + log(1 - mu)  # Binomial case
-    }
-    return(qll)
-}

R/mc_rw1.R

-#' Random walk first order model
-#'
-#' @description Builds a random walk first order model matrix.
-#'
-#' @param n_time Number observations time.
-#' @param intrinsic Logical indicating if the models is intrinsic (rho = 1) or not.
-#' @return A matrix. Note that the function assumes that the data are in the correct order.
-#' @export
-mc_rw1 <- function(n_time, intrinsic = TRUE) {
-    R <- Matrix(0, nrow = n_time, ncol = n_time, sparse = TRUE)
-    ## Border restriction
-    ncol = n_time
-    R[1, c(1, 2)] <- c(1, -1)
-    R[ncol, c(ncol - 1, ncol)] <- c(-1, 1)
-    ## Body of matrix
-    n <- ncol - 1
-    for (i in 2:n) {
-        R[i, c(i - 1, i, i + 1)] <- c(-1, 2, -1)
-    }
-    if (intrinsic == TRUE) {
-        output <- list(R)
-    }
-    if (intrinsic == FALSE) {
-        R1 <- Diagonal(n_time, diag(R))
-        diag(R) <- 0
-        output <- list(Z1 = R1, Z2 = R)
-    }
-    return(output)
-} 

R/mc_rw2.R

-#' Random walk second order model
-#'
-#' @description Builds a random walk second order model matrix.
-#'
-#' @param n_time Number observations time.
-#' @param intrinsic Logical indicating if the models is intrinsic (rho = 1) or not.
-#' @return A matrix. Note that the function assumes that the data are in the correct order.
-#' @export
-mc_rw2 <- function(n_time, intrinsic = TRUE) {
-    R <- Matrix(0, nrow = n_time, ncol = n_time, sparse = TRUE)
-    ## Border restriction
-    ncol = n_time
-    R[1, c(1, 2, 3)] <- c(1, -2, 1)
-    R[ncol, c(c(ncol - 2):ncol)] <- c(1, -2, 1)
-    R[2, c(1:4)] <- c(-2, 5, -4, 1)
-    R[c(ncol - 1), c(c(ncol - 3):c(ncol))] <- c(1, -4, 5, -2)
-    ## Body of matrix
-    n <- ncol - 2
-    for (i in 3:n) {
-        R[i, c(i - 2, i - 1, i, i + 1, i + 2)] <- c(1, -4, 6, -4, 1)
-    }
-    if (intrinsic == TRUE) {
-        output <- list(R)
-    }
-    if (intrinsic == FALSE) {
-        R1 <- Diagonal(n_time, diag(R))
-        diag(R) <- 0
-        output <- list(Z1 = R1, Z2 = R)
-    }
-    return(output)
-} 

R/mc_unstructured.R

-#' Unstructured model
-#'
-#' @description Builds a unstructured model matrix.
-#'
-#' @param n_time Number of observations per unit sample.
-#' @return A matrix. Note that the function assumes that the data are in the correct order.
-#' @export
-
-mc_unstructured <- function(n_time) {
-    mat.temp <- Matrix(0, ncol = n_time, nrow = n_time, sparse = TRUE)
-    non.diagonal.terms <- list()
-    non.diagonal <- t(combn(n_time, 2))
-    n.cor.par <- dim(non.diagonal)[1]
-    ## Covariance elementary matrices
-    for (i in 1:n.cor.par) {
-        non.diagonal.terms[i][[1]] <- mat.temp
-        non.diagonal.terms[i][[1]][non.diagonal[i, 1], non.diagonal[i, 2]] <- non.diagonal.terms[i][[1]][non.diagonal[i,
-            2], non.diagonal[i, 1]] <- 1
-    }
-    ## Output
-    return(non.diagonal.terms)
-}