From b4c1fd9335ca419272b5ac77ec1c962364a61cbc Mon Sep 17 00:00:00 2001
From: wbonat <wbonat@gmail.com>
Date: Tue, 10 Nov 2015 16:13:06 +0100
Subject: [PATCH] Score information criterion new functions
---
DESCRIPTION | 2 +-
NAMESPACE | 2 ++
R/mc_sic.R | 12 ++++++++++--
R/mc_sic_covariance.R | 45 +++++++++++++++++++++++++++++--------------
4 files changed, 44 insertions(+), 17 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index a84f5f7..b8a49da 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 698d0cb..bc47691 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 4d91101..e7b3ae3 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 3c51823..396eef6 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)
}
--
GitLab