diff --git a/DESCRIPTION b/DESCRIPTION index a84f5f71a1ea966e3709fdfda4e0b78176136da0..b8a49da74eb2bfdc423f63e6bb205cbe44d29bde 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,4 +25,4 @@ LazyData: TRUE URL: http://git.leg.ufpr.br/wbonat/mcglm BugReports: http://git.leg.ufpr.br/wbonat/mcglm/issues Encoding: UTF-8 -VignetteBuilder: knitr \ No newline at end of file +VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 698d0cb331ffb672969454b70a4be0a28f0e0617..bc47691aa40d6df8ba5af6d98a37f548024eced3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,6 +12,7 @@ S3method(vcov,mcglm) export(fit_mcglm) export(mc_bias_corrected_std) export(mc_dexp_gold) +export(mc_dfbetaOij) export(mc_fast_forward) export(mc_influence) export(mc_initial_values) @@ -23,6 +24,7 @@ export(mc_robust_std) export(mc_rw1) export(mc_rw2) export(mc_sic) +export(mc_sic_covariance) export(mc_unstructured) export(mc_variance_function) export(mcglm) diff --git a/R/mc_sic.R b/R/mc_sic.R index 4d91101bb9317aedc9ab16d9537cc5a30529ab38..e7b3ae338211e9a3bfef67cff6593bcc2e6a4cd7 100644 --- a/R/mc_sic.R +++ b/R/mc_sic.R @@ -14,6 +14,9 @@ mc_sic <- function (object, scope, data, response, penalty = 2) { SIC <- c() df <- c() + df_total <- c() + TU <- c() + QQ <- c() for(i in 1:length(scope)){ ini_formula <- object$linear_pred[[response]] ext_formula <- as.formula(paste("~", paste(ini_formula[3], @@ -41,9 +44,14 @@ mc_sic <- function (object, scope, data, response, penalty = 2) { 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] + SIC[i] <- -as.numeric(Tu) + penalty*n_beta + df_total[i] <- n_beta + TU[i] <- as.numeric(Tu) + QQ[i] <- qchisq(0.95, df = df[i]) } - output <- data.frame("SIC" = SIC, "Covariance" = scope, "df" = df) + output <- data.frame("SIC" = SIC, "Covariance" = scope, + "df" = df, "df_total" = df_total, + "Tu" = TU, "Chisq" = QQ) return(output) } diff --git a/R/mc_sic_covariance.R b/R/mc_sic_covariance.R index 3c51823deb53e3092eb97998b30277a5632b9e9d..396eef6d12f618b0fdbc7e974dd526b20797451a 100644 --- a/R/mc_sic_covariance.R +++ b/R/mc_sic_covariance.R @@ -5,6 +5,7 @@ #' @param object an object representing a model of \code{mcglm} class. #' @param scope a list of matrices to be tested in the matrix linear #' predictor. +#' @param idx Indicator of matrices belong to the same effect. #' @param data data frame containing all variables envolved in the model. #' @param penalty penalty term (default = 2). #' @param response Indicate for which response variable SIC is computed. @@ -12,22 +13,30 @@ #' argument. #' @export -mc_sic_covariance <- function(object, scope, data, penalty = 2, +mc_sic_covariance <- function(object, scope, idx, data, penalty = 2, response) { - for (j in 1:length(scope)) { + SIC <- c() + df <- c() + df_total <- c() + TU <- c() + QQ <- c() + n_terms <- length(unique(idx)) + for (j in 1:n_terms) { tau <- coef(object, type = "tau", response = response)$Estimates n_tau <- length(tau) - list_tau_new <- list(c(tau, 0)) - n_tau_new <- n_tau + 1 + n_tau_new <- length(idx[idx == j]) + list_tau_new <- list(c(tau, rep(0, n_tau_new))) + n_tau_total <- n_tau + n_tau_new if(object$power_fixed[[response]]){ list_power <- object$list_initial$power } else { list_power <- list(coef(object, type = "power", response = response)$Estimates) - n_tau_new <- n_tau_new + 1 + n_tau_total <- n_tau_total + 1 + n_tau <- n_tau + 1 } - list_Z_new <- list(c(object$matrix_pred[[response]], scope[[j]])) + list_Z_new <- list(c(object$matrix_pred[[response]], scope[idx == j])) if(length(object$mu_list) == 1){rho = 0} else { rho = coef(object,type = "correlation")$Estimates } @@ -48,19 +57,27 @@ mc_sic_covariance <- function(object, scope, data, penalty = 2, J <- temp_score$Sensitivity Sigma <- temp_score$Variability - Sigma22 <- Sigma[n_tau_new,n_tau_new] - J21 <- J[n_tau_new, 1:n_tau] + Sigma22 <- Sigma[c(n_tau+1):n_tau_total,c(n_tau +1):n_tau_total] + J21 <- J[c(n_tau+1):n_tau_total, 1:n_tau] J11 <- solve(J[1:n_tau,1:n_tau]) - Sigma12 <- Sigma[1:n_tau, n_tau_new] - Sigma21 <- Sigma[n_tau_new, 1:n_tau] - J12 <- J[1:n_tau,n_tau_new] + Sigma12 <- Sigma[1:n_tau, c(n_tau+1):n_tau_total] + Sigma21 <- Sigma[c(n_tau+1):n_tau_total, 1:n_tau] + J12 <- J[1:n_tau,c(n_tau+1):n_tau_total] Sigma11 <- Sigma[1:n_tau,1:n_tau] V2 <- Sigma22 - J21%*%J11%*%Sigma12 - Sigma21%*%J11%*%J12 + J21%*%J11%*%Sigma11%*%J11%*%J12 - Tu <- t(temp_score$Score[n_tau_new]%*%solve(V2)%*%temp_score$Score[n_tau_new]) - sic[j] <- sqrt(as.numeric(Tu)) - 2 + TU[j] <- t(temp_score$Score[c(n_tau+1):n_tau_total]%*%solve(V2)%*% + temp_score$Score[c(n_tau+1):n_tau_total]) + df[j] <- n_tau_new + SIC[j] <- -as.numeric(TU[j]) + penalty*n_tau_total + #TU[j] <- as.numeric(TU[j]) + QQ[j] <- qchisq(0.95, df = n_tau_new) + df_total[j] <- n_tau_total + print(j) } - return(data.frame("SIC" = sic, "Df" = rep(1,length(sic)))) + output <- data.frame("SIC" = SIC, "df" = df, "df_total" = df_total, + "Tu" = TU, "Chisq" = QQ) + return(output) }