Add new functions

- Function for computing the score information criterion (mc_sic).
- Function for automatic model selection (mc_fast_forward).

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).
+#' @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 <- c(scope, do.call(c, lapply(int_terms, fun_temp)))
+  for (i in 1:n_max) {
+  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 = matrix_pred,
+                       link = object$link, variance = object$variance,
+                       covariance = object$covariance, data = data)
+  object <- temp_models
+  scope <- scope[-cov_new]
+  }
+  return(object)
+}

R/mc_sic.R

+#' Compute the score information criterion (SIC) for multivariate
+#' covariance generalized linear models.
+#'
+#' @description Compute the SIC for McGLMS.
+#' @param object an object representing a model of \code{mcglm} class.
+#' @param scope a vector containing all covariate names to be tested.
+#' @param data data frame containing the all variables envolved
+#' @param penalty penalty term (default = 2).
+#' @param response Indicate for which response variable SIC is computed.
+#' @return A data frame with SIC values for each covariate in the scope
+#' argument.
+#' @export
+
+mc_sic <- function (object, scope, data, response, penalty = 2) {
+  SIC <- c()
+  df <- c()
+  for(i in 1:length(scope)){
+  ini_formula <- object$linear_pred[[response]]
+  ext_formula <- as.formula(paste("~", paste(ini_formula[3],
+                                             scope[i], sep = "+")))
+  md <- model.frame(object$linear_pred[[response]], data = data)
+  Y = model.response(md)
+  ini_beta <- coef(object, type = "beta", response = response)$Estimates
+  ext_X <- model.matrix(ext_formula, data = data)
+  n_beta <- dim(ext_X)[2]
+  n_ini_beta <- length(ini_beta)
+  ext_beta <- c(ini_beta, rep(0, n_beta - n_ini_beta))
+  n_total_beta <- length(ext_beta)
+  mu_temp <- mc_link_function(beta = ext_beta, X = ext_X, offset = NULL,
+                              link = object$link[[response]])
+  score_temp <- mc_quasi_score(D = mu_temp$D, inv_C = object$inv_C,
+                               y_vec = Y, mu_vec = mu_temp$mu)
+  S11 <- score_temp$Variability[1:n_ini_beta,1:n_ini_beta]
+  S22 <- score_temp$Variability[c(n_ini_beta+1):n_total_beta,
+                                c(n_ini_beta+1):n_total_beta]
+  S12 <- score_temp$Variability[1:n_ini_beta,
+                                c(n_ini_beta+1):n_total_beta]
+  S21 <- score_temp$Variability[c(n_ini_beta+1):n_total_beta,
+                                1:n_ini_beta]
+  VB <- S22 - S21 %*% solve(S11) %*% S12
+  Tu <- t(score_temp$Score[c(n_ini_beta+1):n_total_beta])%*%
+    solve(VB)%*%score_temp$Score[c(n_ini_beta+1):n_total_beta]
+  df[i] <- n_beta - n_ini_beta
+  SIC[i] <- as.numeric(sqrt(Tu)) - penalty*df[i]
+  }
+  output <- data.frame("SIC" = SIC, "Covariance" = scope, "df" = df)
+  return(output)
+}
+

man/mc_fast_forward.Rd

+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/mc_fast_forward.R
+\name{mc_fast_forward}
+\alias{mc_fast_forward}
+\title{Fast forward selection for multivariate covariance generalized linear
+models.}
+\usage{
+mc_fast_forward(object, scope, interaction = 1, penalty = 2, n_max = 10)
+}
+\arguments{
+\item{object}{an object representing a model of \code{mcglm} class.}
+
+\item{scope}{a vector specyfing the covariate to be tested.}
+
+\item{interaction}{Maximum number of covariates interacting.}
+
+\item{penalty}{penalty term (default = 2).}
+}
+\value{
+The selected model.
+}
+\description{
+Perform fast forward model selection using the score
+information criterion. This function works only for univariate months.
+}
+

man/mc_sic.Rd

+% Generated by roxygen2 (4.1.1): do not edit by hand
+% Please edit documentation in R/mc_sic.R
+\name{mc_sic}
+\alias{mc_sic}
+\title{Compute the score information criterion (SIC) for multivariate
+covariance generalized linear models.}
+\usage{
+mc_sic(object, scope, data, response, penalty = 2)
+}
+\arguments{
+\item{object}{an object representing a model of \code{mcglm} class.}
+
+\item{scope}{a vector containing all covariate names to be tested.}
+
+\item{data}{data frame containing the all variables envolved}
+
+\item{response}{Indicate for which response variable SIC is computed.}
+
+\item{penalty}{penalty term (default = 2).}
+}
+\value{
+A data frame with SIC values for each covariate in the scope
+argument.
+}
+\description{
+Compute the SIC for McGLMS.
+}
+