From 7ccdde87e8f28c03f9d185466488baacd3e8b4da Mon Sep 17 00:00:00 2001
From: wbonat <wbonat@gmail.com>
Date: Wed, 6 Jan 2016 18:58:43 +0100
Subject: [PATCH] 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 | 49 -----------------
R/mc_influence.R | 128 --------------------------------------------
R/mc_qic.R | 39 --------------
R/mc_qll.R | 37 -------------
R/mc_rw1.R | 29 ----------
R/mc_rw2.R | 31 -----------
R/mc_unstructured.R | 22 --------
7 files changed, 335 deletions(-)
delete mode 100644 R/mc_fast_forward.R
delete mode 100644 R/mc_influence.R
delete mode 100644 R/mc_qic.R
delete mode 100644 R/mc_qll.R
delete mode 100644 R/mc_rw1.R
delete mode 100644 R/mc_rw2.R
delete mode 100644 R/mc_unstructured.R
diff --git a/R/mc_fast_forward.R b/R/mc_fast_forward.R
deleted file mode 100644
index 07869e6..0000000
--- a/R/mc_fast_forward.R
+++ /dev/null
@@ -1,49 +0,0 @@
-#' 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)
-}
diff --git a/R/mc_influence.R b/R/mc_influence.R
deleted file mode 100644
index fc5e4c4..0000000
--- a/R/mc_influence.R
+++ /dev/null
@@ -1,128 +0,0 @@
-#' @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)
-}
-
-
-
diff --git a/R/mc_qic.R b/R/mc_qic.R
deleted file mode 100644
index e51bd66..0000000
--- a/R/mc_qic.R
+++ /dev/null
@@ -1,39 +0,0 @@
-#' @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))
-}
diff --git a/R/mc_qll.R b/R/mc_qll.R
deleted file mode 100644
index e9e4c14..0000000
--- a/R/mc_qll.R
+++ /dev/null
@@ -1,37 +0,0 @@
-#' 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)
-}
diff --git a/R/mc_rw1.R b/R/mc_rw1.R
deleted file mode 100644
index 519ca71..0000000
--- a/R/mc_rw1.R
+++ /dev/null
@@ -1,29 +0,0 @@
-#' 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)
-}
diff --git a/R/mc_rw2.R b/R/mc_rw2.R
deleted file mode 100644
index 6788c1c..0000000
--- a/R/mc_rw2.R
+++ /dev/null
@@ -1,31 +0,0 @@
-#' 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)
-}
diff --git a/R/mc_unstructured.R b/R/mc_unstructured.R
deleted file mode 100644
index afc67eb..0000000
--- a/R/mc_unstructured.R
+++ /dev/null
@@ -1,22 +0,0 @@
-#' 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)
-}
--
GitLab