R/SEM.R

In CorReg: Linear Regression Based on Linear Structure Between Variables

SEM <- function(M = M, nbit_gibbs = 1, n = n, nbit_SEM = 50, warm = 10, mixmod = mixmod, X = X, comp_vect = comp_vect, missrow = missrow, quimiss = quimiss, 
                Z = Z, Zc = Zc, vraiX = NULL, alpha = alpha, GibbsIR = TRUE, sigma_IR = sigma_IR, loglikout = FALSE, gibbsfinwarm = 0, nbclust_vect = nbclust_vect,
                Ir = Ir, compout = TRUE, Xout = FALSE, alphaout = TRUE, gibbsfin = 0, verbose = 1, fill = FALSE){
   last = FALSE
   result = list()
   # initialisation
   p = ncol(Z)
   if(is.null(Zc)){
      Zc = colSums(Z)
   }
   if(is.null(Ir)){
      Ir = which(Zc!= 0)
   }
   if(is.null(comp_vect)){
      comp_vect = matrix(rep(1, times = p*n), ncol = p)
      #       for(i in 1:p){
      #          if(nbclust_vect[i]>1){
      #             comp_vect[, i] = which(rmultinom(n  = missrow,  size  =  1, prob = mixmod[[i]][, 1]) == 1, arr.ind  =  TRUE)[, 1]
      #          }
      #       }
      #I mputation des classes####
      for(j in (1:p)[-Ir]){ # pas besoin de classe a gauche car dans Gibbs on connait toujours tout le monde a droite (quitte a le supposer)
         if(nbclust_vect[j] > 1){ # si plusieurs classes (sinon reste a 1)
            for (i in 1:n){
               if(!is.na(X[i, j]) & !fill){ # valeur observee (ou imputee,  peu importe) on estime la classe
                  vect = rmultinom(1,  1,  tik(x = X[i, j], nbclust = nbclust_vect[j], mixmod = mixmod[[j]]))
                  comp_vect[i, j] = match(1, vect)
               }else if(is.na(X[i, j]) | (fill & M[i, j] == 1)){  #manquant donc on impute la classe et la valeur 
                  vect = rmultinom(1,  1,  mixmod[[j]][, 1])
                  comp_vect[i, j] = match(1, vect)
                  X[i, j] = rnorm(n = 1, mean = mixmod[[j]][comp_vect[i, j], 2], sd = sqrt(mixmod[[j]][comp_vect[i, j], 3])) # on tire dans la classe choisie
               }
            }
         }
      }
   }
   if(is.null(alpha)){
      alpha = hatB(Z  =  Z, X = X)
   }
   if(is.null(sigma_IR)){
      sigma_IR = rep(0, times = p)#residus des regressions (fixes avec alpha donc stables pour gibbs)
      for (j in Ir){
         sigma_IR[j] = sd(X[, j]-X%*%alpha[-1, j])
      }
   }
   
   loglik_bool = FALSE
   for (i in 1:(nbit_SEM+warm)){
      if((i == (nbit_SEM+warm)) & (gibbsfin <= 0)){
         last = TRUE
      }
      if((i > warm) & loglikout & (gibbsfin <= 0)){
         loglik_bool = TRUE
      } # on commence a calculer les vraisemblances
      
      # SE step####
      resgibbs2 = Gibbs(M = M, last = last, nbit = nbit_gibbs, mixmod = mixmod, X = X, warm = 0, comp_vect = comp_vect, missrow = missrow, quimiss = quimiss, 
                        Z = Z, Zc = Zc, alpha = alpha, sigma_IR = sigma_IR, nbclust_vect = nbclust_vect, Ir = Ir, GibbsIR = GibbsIR, loglik_bool = loglik_bool)
      comp_vect = resgibbs2$comp_vect
      X = resgibbs2$Xfin
      
      #step M####
      resM = Mstep(Z = Z, X = X, sigma_IR = sigma_IR, Ir = Ir)
      alpha = resM$alpha
      sigma_IR = resM$sigma_IR
      
      if(i>warm){
         if(loglikout){ # on calcule la vraisemblance locale
            loglik = resgibbs2$loglik # scalaire
         }
         if(i == warm+1){
            if(verbose){
               cat("warm-up finished")
            }
            if(alphaout){
               result$alpha = alpha/nbit_SEM
            }
         }else{
            if(alphaout){
               result$alpha = result$alpha+alpha/nbit_SEM
            }# optimiser au format creux
         }   
      }
   }
   
   if(compout){
      result$comp = comp_vect
   }
   
   if(gibbsfin >0 ){
      if(verbose){
         cat("final Gibbs")
      }      
      resgibbs2 = Gibbs(M = M, last = last, warm = gibbsfinwarm, nbit = gibbsfin, mixmod = mixmod, X = X, comp_vect = comp_vect, missrow = missrow, quimiss = quimiss, 
                        Z = Z, Zc = Zc, alpha = result$alpha, GibbsIR = GibbsIR, sigma_IR = sigma_IR, nbclust_vect = nbclust_vect, Ir = Ir, loglik_bool = loglik_bool, Xout = Xout)
      if(Xout){
         result$Xfin = resgibbs2$Xfin
      }
      result$loglik = resgibbs2$loglik
      result$X = resgibbs2$X
      
   }
   return(result)
}