R/BICZmiss.R

In CorReg: Linear Regression Based on Linear Structure Between Variables

# ' BIC of Z with missing values in X
# ' 
# ' 

BICZmiss <- function(X = X, Z = Z, M = NULL, Bic_vide_vect = NULL, BicOld = NULL, methode = 1, Zold = NULL, 
                     intercept = TRUE, mixmod = NULL, nbit = 100, nbclustmax = 10, mixmodpenalty = TRUE) {
  
  #   res=.Call( "BicZmissRcpp", X, Z, M, Bic_vide_vect, BicOld, methode, Zold, intercept, mixmod, nbit, nbclustmax, PACKAGE = "CorReg")
  #   return(res)
  #
  if (is.null(M)) {
    M = 0 * X
    M[is.na(M)] = 1
    X[is.na(X)] = 0
  }
  p = ncol(Z)
  if (is.null(Bic_vide_vect)) {
    if (!is.null(mixmod)) {
      Bic_vide_vect = rep(0, times = p)
      for (i in 1:p) {
        Bic_vide_vect[i] = mixmod$BIC_vect[i]
      }
    } else {
      Bic_vide_vect = density_estimation(X = X)$BIC_vect
    }
  }
  if (is.null(Zold) | is.null(BicOld)) {
    Zold = 0 * Z
    BicOld = Bic_vide_vect
  }
  BIC_vect = BicOld
  for (j in 1:p) { # pour chaque variable (colonne)
    if (length(which(Z[, j] - Zold[, j] != 0)) > 0) { # differences donc on recalcule Bicloc
      if (sum(Z[, j]) == 0) { # colonne devenue vide
        BIC_vect[j] = Bic_vide_vect[j]
      } else { # Bic doit etre recalcule
        quiloc = (Z[, j] != 0)
        Xloc = as.matrix(X[, quiloc])
        Yloc = X[, j]
        Mloc = as.matrix(M[, quiloc])
        quibon = (M[, j] == 0)
        Yloc = Yloc[quibon]
        Xloc = as.matrix(Xloc[quibon, ])
        Mloc = as.matrix(Mloc[quibon, ])
        mixmodloc = list()
        mixmodloc$nbclust = mixmod$nbclust[quiloc]
        mixmodloc$details = matrix(ncol = 4, nrow = sum(mixmodloc$nbclust))
        k = 1
        l = 1
        m = quiloc[1]
        for (i in 1:nrow(mixmod$details)) {
          if (Z[mixmod$details[i, 4], j] != 0) { # si la variable concernee est dans Xloc on recopie
            if (mixmod$details[i, 4] == m) {
              # on est sur la meme variable donc le numero ne change pas
            } else { # changement de variable donc il faut reindexer
              l = l + 1
              m = mixmod$details[i, 4]
            }
            mixmodloc$details[k, ] = c(mixmod$details[i, 1:3], l)
            k = k + 1
          }
        }
        #         message(mixmodloc$details)
        res = OLS(X = Xloc, Y = Yloc, M = Mloc, intercept = intercept, methode = methode, miss = (sum(Mloc) > 0), 
                  mixmod = mixmodloc, nbit = nbit, nbclustmax = nbclustmax, sigma = TRUE)
        BICLOC = c()
        k = 0 # complexite du modele
        for (i in 1:nrow(Xloc)) { # pour chaque point
          BICLOC = c(BICLOC, GM_Loglikelihood(Y = Yloc[i], X = Xloc[i, ], M = Mloc[i, ], B = res$beta, 
                                              sigma = res$sigma, mixmod = mixmodloc, log = TRUE, intercept = TRUE))
          k = k + ncol(Xloc) + intercept + 1 # Beta+sigma dans tous les cas
          if (sum(Mloc[i, ]) > 0) { # trous donc on ajoute des parametres mixmod
            classesmank = Z[, j] * mixmodloc$nbclust * Mloc[i, ]
            # k=k+3sum(classmank)-nbfois ou classmank=1
            if (mixmodpenalty) {
              k = k + prod(classesmank[classesmank > 1]) # nombre de proportions estimees avec une par classe quand plusieurs classes
              k = k + 2 * sum(classesmank) # moyennes et ecarts-types
            }
          }
        }
        k = k / nrow(Xloc)
        BIC_vect[j] = -2 * sum(BICLOC) + k * log(nrow(Xloc))
      }
    }
  } # fin de la boucle sur les colonnes
  
  return(BIC_vect)
}