Score information criterion new functions

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)

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)
 }
 

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)
 }