R/MFA.R

MFA <- function (base, group, type = rep("s",length(group)), excl = NULL, ind.sup = NULL, ncp = 5, name.group = NULL, num.group.sup = NULL, graph = TRUE, weight.col.mfa = NULL, row.w = NULL, axes=c(1,2),tab.comp=NULL){

    moy.p <- function(V, poids) {
        res <- sum(V * poids,na.rm=TRUE)/sum(poids[!is.na(V)])
    }
    ec <- function(V, poids) {
        res <- sqrt(sum(V^2 * poids,na.rm=TRUE)/sum(poids[!is.na(V)]))
    }
	funcLg <- function (x, y, ponderation.x, ponderation.y, wt = rep(1/nrow(x), nrow(x)), cor = FALSE) {
      if (is.data.frame(x)) x <- as.matrix(x)
      else if (!is.matrix(x)) stop("'x' must be a matrix or a data frame")
      if (is.data.frame(y)) y <- as.matrix(y)
      else if (!is.matrix(y)) stop("'y' must be a matrix or a data frame")
      if (!all(is.finite(x))) stop("'x' must contain finite values only")
      if (!all(is.finite(y))) stop("'y' must contain finite values only")
      s <- sum(wt)
	  wt <- wt/s
      center <- colSums(wt * x)
      x <- sqrt(wt) * t(t(sweep(x, 2, center, check.margin = FALSE))*sqrt(ponderation.x))
      center <- colSums(wt * y)
      y <- sqrt(wt) * t(t(sweep(y, 2, center, check.margin = FALSE))*sqrt(ponderation.y))
	  Lg <- 0
	  for (i in 1:ncol(x)) Lg <- Lg+sum(crossprod(x[,i],y)^2)
	  Lg
    }
	if (!is.null(excl) & "m"%in%type) stop("Excluding categories is not allowed when some groups are mixed")
if (!is.null(tab.comp)){
  if (!is.null(weight.col.mfa)) stop("Weightings on the variables are not allowed with the tab.comp argument")
}

    if (is.null(rownames(base))) rownames(base) <- 1:nrow(base)
    if (is.null(colnames(base))) colnames(base) <- paste("V",1:ncol(base),sep="")
    base <- as.data.frame(base)
    is.quali <- which(!unlist(lapply(base,is.numeric)))
    base[,is.quali] <- lapply(base[,is.quali,drop=FALSE],as.factor)
	base <- droplevels(base)
    nature.group <- NULL
    nbre.group <- length(group)
for (i in 1:nbre.group){
  if ((type[i] == "n")&&(!(i%in%num.group.sup))) nature.group <- c(nature.group,"quali")
  if ((type[i] == "n")&&(i%in%num.group.sup)) nature.group <- c(nature.group,"quali.sup")
  if (((type[i] == "s")||(type[i] == "c"))&&(!(i%in%num.group.sup))) nature.group <- c(nature.group,"quanti")
  if (((type[i] == "s")||(type[i] == "c"))&&(i%in%num.group.sup)) nature.group <- c(nature.group,"quanti.sup")
  if (((type[i] == "f")||(type[i] == "f2"))&&(!(i%in%num.group.sup))) nature.group <- c(nature.group,"contingency")
  if (((type[i] == "f")||(type[i] == "f2"))&&(i%in%num.group.sup)) nature.group <- c(nature.group,"contingency.sup")
  if ((type[i] == "m")&&(!(i%in%num.group.sup))) nature.group <- c(nature.group,"mixed")
  if ((type[i] == "m")&&(i%in%num.group.sup)) nature.group <- c(nature.group,"mixed.sup")
}
# nature.var <- rep(nature.group,times=group)

list.type.var <- vector(mode = "list", length = nbre.group)
for (i in 1:nbre.group){
  if ((type[i] == "n")&&(!(i%in%num.group.sup))) list.type.var[[i]] <- rep("quali",group[i])
  if ((type[i] == "n")&&(i%in%num.group.sup)) list.type.var[[i]] <- rep("quali.sup",group[i])
  if (((type[i] == "s")||(type[i] == "c"))&&(!(i%in%num.group.sup))) list.type.var[[i]] <- rep("quanti",group[i])
  if (((type[i] == "s")||(type[i] == "c"))&&(i%in%num.group.sup)) list.type.var[[i]] <- rep("quanti.sup",group[i])
  if (((type[i] == "f")||(type[i] == "f2")||(type[i] == "f3"))&&(!(i%in%num.group.sup))) list.type.var[[i]] <- rep(type[i],group[i])
  if (((type[i] == "f")||(type[i] == "f2")||(type[i] == "f3"))&&(i%in%num.group.sup)) list.type.var[[i]] <- rep(paste(type[i],"sup",sep="_"),group[i])
  if ((type[i] == "m")&&(!(i%in%num.group.sup))){
    if (i==1) list.type.var[[i]] <- ifelse(sapply(base[,1:group[1]],is.numeric),"quanti","quali")
	else list.type.var[[i]] <- ifelse(sapply(base[,(sum(group[1:(i-1)])+1):sum(group[1:i])],is.numeric),"quanti","quali")
  }
  if ((type[i] == "m")&&(i%in%num.group.sup)){
    if (i==1) list.type.var[[i]] <- ifelse(sapply(base[,1:group[1]],is.numeric),"quanti.sup","quali.sup")
    else list.type.var[[i]] <- ifelse(sapply(base[,(sum(group[1:(i-1)])+1):sum(group[1:i])],is.numeric),"quanti.sup","quali.sup")
  }
  if (i==1) names(list.type.var[[1]]) <- colnames(base)[1:group[1]]
  else names(list.type.var[[i]]) <- colnames(base)[(sum(group[1:(i-1)])+1):sum(group[1:i])]
}
 nature.var <- type.var <- unlist(list.type.var)
 names(nature.var)=NULL
 
 if (!is.null(ind.sup)) {
      base <- rbind.data.frame(base[-ind.sup,],base[ind.sup,,drop=FALSE])
      ind.sup <- (nrow(base)-length(ind.sup)+1) : nrow(base)
    }
    nbre.var <- ncol(base)
    group.actif <- NULL
    if ("n"%in%type || "m"%in%type){
      niveau <- NULL
      for (j in 1:ncol(base)){
        if (!is.numeric(base[,j])) niveau <- c(niveau,levels(base[,j]))
      }
      for (j in 1:ncol(base)) {
        if (!is.numeric(base[,j])){
          if (sum(niveau%in%levels(base[,j]))!=nlevels(base[,j])) levels(base[,j]) <- paste(colnames(base)[j],levels(base[,j]),sep="_")
        }
      }
    }
    for (i in 1:nbre.group) if (!(i%in%num.group.sup)) group.actif <- c(group.actif,i)
    group.mod <- group
    nbre.ind <- nrow(base)
    nb.actif <- nbre.ind-length(ind.sup)
    if (nbre.var != sum(group)) stop("not convenient group definition")
    if (nbre.group != length(type)) stop("not convenient type definition")
    base <- as.data.frame(base)
    if (!inherits(base, "data.frame")) stop("base should be a data.frame")
    res.separe <- vector(mode = "list", length = nbre.group)
    if (is.null(name.group)) name.group <- paste0("Gr", 1:nbre.group)
    names(res.separe) <- name.group
    ind.grpe <- 0
## modif 2023/11/29
    if (any(is.na(base)) & is.null(tab.comp)){
       for (j in 1:ncol(base)){
	     if (is.numeric(base[,j])) base[,j] <- replace(base[,j],is.na(base[,j]),mean(base[,j],na.rm=TRUE))
		 else base[,j] <- as.factor(replace(as.character(base[,j]),is.na(base[,j]),paste(colnames(base)[j],".NA",sep="")))
	   }
	 }
    # if (any(is.na(base))){
      # if (!("n"%in%type)) for (j in 1:ncol(base)) base[,j] <- replace(base[,j],is.na(base[,j]),mean(base[,j],na.rm=TRUE))
      # else{
        # if (type[1]!="n") for (j in 1:group[1]) base[,j] <- replace(base[,j],is.na(base[,j]),mean(base[,j],na.rm=TRUE))
        # for (g in 2:nbre.group){
         # if (type[g]!="n") for (j in (sum(group[1:(g-1)])+1):sum(group[1:g])) base[,j] <- replace(base[,j],is.na(base[,j]),mean(base[,j],na.rm=TRUE))
        # }
		# if (is.null(tab.comp)){
		  # if (type[1]=="n") for (j in 1:group[1]) base[,j] <- as.factor(replace(as.character(base[,j]),is.na(base[,j]),paste(colnames(base)[j],".NA",sep="")))
          # for (g in 2:nbre.group){
            # if (type[g]=="n") for (j in (sum(group[1:(g-1)])+1):sum(group[1:g])) base[,j] <- as.factor(replace(as.character(base[,j]),is.na(base[,j]),paste(colnames(base)[j],".NA",sep="")))
          # }
		# }
	   #}
    # }

    if (is.null(row.w)) row.w <- rep(1,nb.actif)

    if (any("f" %in% type)+any("f2" %in% type)+any("f3" %in% type)>1) stop("For the contingency tables, the type must the the same")
    if (("f" %in% type)||("f2" %in% type)||("f3" %in% type)) {
		grfrec<-c(which(type=="f"),which(type=="f2"),which(type=="f3"))

## pour avoir individus actifs, que ind.sup soit NULL ou non
##		ind.actif <- !((1:nrow(base))%in%intersect(ind.sup,(1:nrow(base))))
		for (i in grfrec){
			if ((type[i]=="f2")||(type[i]=="f3")||(i%in%num.group.sup)){
				if (i==1) base[,1:group[1]]<- base[,1:group[1]]/sum(base[1:nb.actif,1:group[1]])
				else base[,(sum(group[1:(i-1)])+1):sum(group[1:i])]<-base[,(sum(group[1:(i-1)])+1):sum(group[1:i])]/sum(base[1:nb.actif,(sum(group[1:(i-1)])+1):sum(group[1:i])])
			}
		}
		type.var=="f"
        if(!any(type.var=="f")) sumT <-1
		else sumT <- sum(base[1:nb.actif,as.logical((type.var=="f")+(type.var=="f2")+(type.var=="f3"))])
		if (sumT==0) sumT <- 1
		base[,as.logical((type.var=="f")+(type.var=="f_sup")+(type.var=="f2")+(type.var=="f2_sup")+(type.var=="f3")+(type.var=="f3_sup"))]<-base[,as.logical((type.var=="f")+(type.var=="f_sup")+(type.var=="f2")+(type.var=="f2_sup")+(type.var=="f3")+(type.var=="f3_sup"))]/sumT
		F.jt<-list()
		Fi.t<-list()
		for (j in grfrec){
			if (j==1){
				F.jt[[j]]<-apply(base[1:nb.actif,1:group[1]],2,sum)
				Fi.t[[j]]<-apply(base[,1:group[1]],1,sum)
			}else{
				F.jt[[j]]<-apply(base[1:nb.actif,(sum(group[1:(j-1)])+1):(sum(group[1:(j-1)])+group[j])],2,sum)
				Fi.t[[j]]<-apply(base[,(sum(group[1:(j-1)])+1):(sum(group[1:(j-1)])+group[j])],1,sum)
			}
		}
		if ("f"%in%type.var){
    		row.w[1:nrow(base)]<-0
			for (j in grfrec){
				if (!(j%in%num.group.sup)) row.w<-row.w+Fi.t[[j]]
			}
		}

		F..t<-numeric()
		for (j in grfrec)	F..t[j]<-sum(Fi.t[[j]][1:nb.actif])
	
		for (t in grfrec){
			if (t==1) {
			    base[,1:group[t]]<-sweep(base[,1:group[t]],2,F.jt[[t]],FUN="/")
			    base[,1:group[t]] <- sweep(base[,1:group[t]],1,Fi.t[[t]]/F..t[t],FUN="-")
			    base[,1:group[t]] <- sweep(base[,1:group[t]],1,row.w,FUN="/")
			}
			else {
			    base[,(sum(group[1:(t-1)])+1):sum(group[1:t])]<-sweep(base[,(sum(group[1:(t-1)])+1):sum(group[1:t])],2,F.jt[[t]],FUN="/")
			    base[,(sum(group[1:(t-1)])+1):sum(group[1:t])]<-sweep(base[,(sum(group[1:(t-1)])+1):sum(group[1:t])],1,Fi.t[[t]]/F..t[t],FUN="-")
			    base[,(sum(group[1:(t-1)])+1):sum(group[1:t])]<-sweep(base[,(sum(group[1:(t-1)])+1):sum(group[1:t])],1,row.w,FUN="/")
			}
		}
		row.w <- row.w[1:nb.actif]
	}		

    if (!is.null(ind.sup))  row.w.moy.ec <- c(row.w,rep(0,length(ind.sup)))
    else row.w.moy.ec <- row.w

if (is.null(weight.col.mfa)) weight.col.mfa <- rep(1,sum(group.mod))
### Begin handle missing values
if (!is.null(tab.comp)){
  group.mod <- tab.comp$call$group.mod
  ind.var.group <- tab.comp$call$ind.var
  tab.comp <- tab.comp$tab.disj
}
### End  handle missing values
   for (g in 1:nbre.group) {
        aux.base <- as.data.frame(base[, (ind.grpe + 1):(ind.grpe + group[g])])
        dimnames(aux.base) <- list(rownames(base),colnames(base)[(ind.grpe + 1):(ind.grpe + group[g])])
        if (type[g] == "s") res.separe[[g]] <- PCA(aux.base, ind.sup =ind.sup, scale.unit = TRUE, ncp = ncp, row.w=row.w, graph = FALSE, col.w = weight.col.mfa[(ind.grpe + 1):(ind.grpe + group[g])])
        if (type[g] == "c") res.separe[[g]] <- PCA(aux.base, ind.sup =ind.sup, scale.unit = FALSE, ncp = ncp, row.w=row.w,graph = FALSE, col.w = weight.col.mfa[(ind.grpe + 1):(ind.grpe + group[g])])
        if (type[g]=="f"||type[g]=="f2"||type[g]=="f3")  res.separe[[g]] <- PCA(aux.base, ind.sup =ind.sup, scale.unit = FALSE, ncp = ncp, row.w=row.w,graph = FALSE, col.w = F.jt[[g]]*weight.col.mfa[(ind.grpe + 1):(ind.grpe + group[g])])
        if (type[g] == "n") {
          for (v in (ind.grpe + 1):(ind.grpe + group[g])) {
            if (!is.factor(base[, v])) stop("factors are not defined in the qualitative groups")
          }
          res.separe[[g]] <- MCA(aux.base, excl = excl[[g]], ind.sup = ind.sup, ncp=ncp, graph = FALSE, row.w=row.w)
        }
        if (type[g] == "m") res.separe[[g]] <- FAMD(aux.base, ind.sup = ind.sup, ncp=ncp, graph = FALSE, row.w=row.w)
###  Begin handle missing values
if (!is.null(tab.comp)){
 if (type[g] == "s") res.separe[[g]] <- PCA(tab.comp[,ind.var.group[[g]]],scale.unit=TRUE,row.w=row.w,ind.sup=ind.sup,col.w=weight.col.mfa[(ind.grpe + 1):(ind.grpe + group[g])],graph=FALSE)
 if (type[g] == "c") res.separe[[g]] <- PCA(tab.comp[,ind.var.group[[g]]],scale.unit=FALSE,row.w=row.w,ind.sup=ind.sup,col.w=weight.col.mfa[(ind.grpe + 1):(ind.grpe + group[g])],graph=FALSE)
 if (type[g] == "n") res.separe[[g]] <- MCA(aux.base, ind.sup = ind.sup, ncp=ncp, graph = FALSE, row.w=row.w,tab.disj=tab.comp[,ind.var.group[[g]]])
 if (type[g] == "m") res.separe[[g]] <- FAMD(aux.base, ind.sup = ind.sup, ncp=ncp, graph = FALSE, row.w=row.w,tab.disj=tab.comp[,ind.var.group[[g]]])
}
###  End handle missing values

        ind.grpe <- ind.grpe + group[g]
    }
    data <- matrix(0,nbre.ind,0)
    ind.grpe <- ind.grpe.mod <- 0
    ponderation <- vector(length = sum(group.mod))
    poids.bary <- NULL
    for (g in 1:nbre.group) {
        aux.base <- base[, (ind.grpe + 1):(ind.grpe + group[g]),drop=FALSE]
###  Begin handle missing values
        if (!is.null(tab.comp)){
          if (g==1) aux.base <- tab.comp[,1:group.mod[1]]
          else aux.base <- tab.comp[,(cumsum(group.mod)[g-1]+1):cumsum(group.mod)[g],drop=FALSE]
        }
###  End handle missing values
        aux.base <- as.data.frame(aux.base)
        colnames(aux.base) <- colnames(base)[(ind.grpe + 1):(ind.grpe + group[g])]
        if (type[g] == "s") {
          centre.aux.base <- apply(as.data.frame(aux.base), 2, moy.p, row.w.moy.ec)
          aux.base <- t(t(as.matrix(aux.base))-centre.aux.base)
          ecart.type.aux.base <- apply(as.data.frame(aux.base), 2, ec, row.w.moy.ec)
          ecart.type.aux.base[ecart.type.aux.base <= 1e-08] <- 1
          aux.base <- t(t(aux.base)/ecart.type.aux.base)
          type[g] <- "c"
        }
        if (type[g] == "c") {
          data <- cbind.data.frame(data, aux.base)
          ponderation[(ind.grpe.mod + 1):(ind.grpe.mod + group.mod[g])] <- 1/res.separe[[g]]$eig[1,1]
        }
        if (type[g] == "f"||type[g] == "f2") {
			data <- cbind.data.frame(data, aux.base)
			ponderation[(ind.grpe.mod+1):(ind.grpe.mod+group[g])]<-F.jt[[g]]/res.separe[[g]]$eig[1,1]
        }
        if (type[g] == "n") {
###  Begin handle missing values
            if (!is.null(tab.comp)){
              if (g==1) tmp <- tab.comp[,1:group.mod[1]]
              else tmp <- tab.comp[,(cumsum(group.mod)[g-1]+1):cumsum(group.mod)[g],drop=FALSE]
            } else {
			  tmp <- tab.disjonctif(aux.base)
              group.mod[g] <- ncol(tmp)
			}
###  End handle missing values
            centre.tmp <- apply(tmp, 2, moy.p, row.w.moy.ec)
            centre.tmp <- centre.tmp/sum(row.w.moy.ec)
            tmp2 <- tmp*(row.w.moy.ec/sum(row.w.moy.ec))
            poids.bary <- c(poids.bary,colSums(tmp2))
            poids.tmp <- 1-apply(tmp2, 2, sum)
            if(!is.null(excl[[g]])) poids.tmp[excl[[g]]] <- 0
            ponderation[(ind.grpe.mod + 1):(ind.grpe.mod + group.mod[g])] <- poids.tmp/(res.separe[[g]]$eig[1,1] * group[g])
            tmp <- tmp/sum(row.w.moy.ec)
            tmp <- t(t(as.matrix(tmp))-centre.tmp)
            ecart.type.tmp <- apply(tmp, 2, ec, row.w.moy.ec)
### Pb if the disjunctive table doesn't have only 0 and 1
			if (!is.null(tab.comp)) ecart.type.tmp <- sqrt(centre.tmp*sum(row.w.moy.ec) * (1-centre.tmp*sum(row.w.moy.ec) ))/sum(row.w.moy.ec)
### End pb if the disjunctive table doesn't have only 0 and 1
            ecart.type.tmp[ecart.type.tmp <= 1e-08] <- 1
			tmp <- t(t(as.matrix(tmp))/ecart.type.tmp)
            data <- cbind.data.frame(data, as.data.frame(tmp))
        }
        if (type[g] == "m") {
###  Begin handle missing values
            if (is.null(tab.comp)){
	          numAct <- which(sapply(aux.base,is.numeric))
	          facAct <- which(!sapply(aux.base,is.numeric))
              QuantiAct <- as.matrix(aux.base[,numAct,drop=FALSE])
              QualiAct <- tab.disjonctif(aux.base[,facAct,drop=FALSE])
              group.mod[g] <- ncol(QuantiAct)+ncol(QualiAct)
			} else {
              if (g==1){
 			    aux <- tab.comp[,1:group.mod[1],drop=FALSE]
			  } else {
			    aux <- tab.comp[,(cumsum(group.mod)[g-1]+1):cumsum(group.mod)[g],drop=FALSE]
			  }
			  QuantiAct <- aux[,which(colnames(aux)==colnames(QuantiAct))]
			  QualiAct <- aux[,which(colnames(aux)!=colnames(QuantiAct))]
			}
###  End handle missing values

	        centre <- apply(QuantiAct,2,moy.p,row.w.moy.ec)
	        QuantiAct <- t(t(QuantiAct)-centre)
	        ecart.type.tmp <- apply(QuantiAct,2,ec,row.w.moy.ec)
	        QuantiAct <- t(t(QuantiAct)/ecart.type.tmp)
            
            tmp2 <- QualiAct*(row.w.moy.ec/sum(row.w.moy.ec))
            poids.bary <- c(poids.bary,colSums(tmp2))
			
			nk_sur_n <- apply(QualiAct, 2, moy.p, row.w.moy.ec)
            QualiAct <- t(t(as.matrix(QualiAct))/sqrt(nk_sur_n))
            centre.tmp <- apply(QualiAct, 2, moy.p, row.w.moy.ec)
            QualiAct <- t(t(as.matrix(QualiAct))-centre.tmp)

            ponderation[(ind.grpe.mod + 1):(ind.grpe.mod + group.mod[g])] <- 1/(res.separe[[g]]$eig[1,1])
### Pb if the disjunctive table doesn't have only 0 and 1
#			if (!is.null(tab.comp)) ecart.type.tmp <- sqrt(centre*sum(row.w.moy.ec) * (1-centre*sum(row.w.moy.ec)))/sum(row.w.moy.ec)
### End pb if the disjunctive table doesn't have only 0 and 1
            # ecart.type.tmp[ecart.type.tmp <= 1e-08] <- 1
			data <- cbind.data.frame(data,QuantiAct,QualiAct)
        }

## Fin modif
        if(!is.null(excl)) ponderation[ponderation == 0] <- 1e-15
        ind.grpe <- ind.grpe + group[g]
        ind.grpe.mod <- ind.grpe.mod + group.mod[g]
    }
    data.group.sup.indice <- data.group.sup <- NULL
    data.pca <- data
    rownames(data.pca) <- rownames(base)
    if (!is.null(num.group.sup)){
      ponderation.tot <- ponderation
      ponderation.group.sup <- NULL
      nb.of.var <- 0
      supp.quanti <- supp.quali <- NULL
      colnames.data.group.sup <- NULL
      for (i in 1:nbre.group) {
        if (i%in%num.group.sup){
          if ((type[i]=="c")||(type[i]=="f")) supp.quanti <- c(supp.quanti,(nb.of.var+1):(nb.of.var+group.mod[i]))
          if (type[i]=="n") supp.quali <- c(supp.quali,(1+nb.of.var):(nb.of.var+group.mod[i]))
          if (type[i]=="m"){
		    supp.quanti <- c(supp.quanti,nb.of.var+(1:sum(sapply(base[, (sum(group[1:(i-1)])+1):(sum(group[1:i])),drop=FALSE],is.numeric))))
			supp.quali <- c(supp.quali,nb.of.var+(sum(sapply(base[, (sum(group[1:(i-1)])+1):(sum(group[1:i])),drop=FALSE],is.numeric)):(group.mod[i])))
          }
		  if (is.null(data.group.sup)) data.group.sup <- as.data.frame(data[,(1+nb.of.var):(nb.of.var+group.mod[i])])
          else data.group.sup <- cbind.data.frame(data.group.sup,data[,(1+nb.of.var):(nb.of.var+group.mod[i])])
          if (ncol(data.group.sup)>1) colnames.data.group.sup <- c(colnames.data.group.sup,colnames(data)[(1+nb.of.var):(nb.of.var+group.mod[i])])
          else colnames.data.group.sup <- colnames(data)[1+nb.of.var]
          ponderation.group.sup <- c(ponderation.group.sup,ponderation[(1+nb.of.var):(nb.of.var+group.mod[i])])
        }
        nb.of.var <- nb.of.var + group.mod[i]
      }
      colnames(data.group.sup) <- colnames.data.group.sup
      ponderation <- ponderation.tot[-c(supp.quanti,supp.quali)]
      data <- data[,-c(supp.quanti,supp.quali)]
      data.group.sup.indice <- (ncol(data)+1):(ncol(data)+ncol(data.group.sup))
      data.pca <- cbind.data.frame(data,data.group.sup)
    }
# modif 2023-01-16
#    ncp.tmp <- min(nb.actif-1, ncol(data))
    ncp.tmp <- min(nb.actif-1, ncol(data)-sum((group[group.actif])[type[group.actif]=="n"]))
    ind.var <- 0
    ind.quali <- NULL
    for (g in 1:nbre.group) {
        if (type[g] == "n")  ind.quali <- c(ind.quali, c((ind.var + 1):(ind.var + group[g])))
        if (type[g] == "m"){
		  if (g==1) ind.quali <- which(!sapply(base[, 1:group[1],drop=FALSE],is.numeric))
		  else ind.quali <- c(ind.quali, ind.var + which(!sapply(base[, (1+sum(group[1:(g-1)])):sum(group[1:g]),drop=FALSE],is.numeric)))
		}
        ind.var <- ind.var + group[g]
    }
    aux_quali_sup_indice <- aux_quali_sup <- data.sup <- NULL
    if (!is.null(ind.quali)){
      aux_quali_sup <- as.data.frame(base[, ind.quali,drop=FALSE])
      if (is.null(data.group.sup)) aux_quali_sup_indice <- (ncol(data)+1):(ncol(data)+ncol(aux_quali_sup))
      else aux_quali_sup_indice <- (ncol(data)+ncol(data.group.sup)+1):(ncol(data)+ncol(data.group.sup)+ncol(aux_quali_sup))
      data.pca <- cbind.data.frame(data.pca,aux_quali_sup)
    }
row.w <- row.w[1:nb.actif]
###  Begin handle missing values
if ((!is.null(tab.comp))&(any("n"%in%type) || any("m"%in%type))){
  data.pca <- data.pca[,-aux_quali_sup_indice]
  aux_quali_sup_indice <- NULL
}
###  End handle missing values
 res.globale <- PCA(data.pca, scale.unit = FALSE, col.w = ponderation, row.w=row.w,ncp = ncp, ind.sup = ind.sup, quali.sup = aux_quali_sup_indice, quanti.sup = data.group.sup.indice, graph = FALSE)
###  Begin handle missing values
if ((!is.null(tab.comp))&(any("n"%in%type) || any("m"%in%type))){
  if (!is.null(type.var%in%"quali")) {
    if (any(type.var%in%c("freq","quanti"))){
	  res.globale$quali.var$coord <- res.globale$var$coord[-(1:sum(type.var%in%c("freq","quanti"))),]
      res.globale$quali.var$cos2 <- res.globale$var$cos2[-(1:sum(type.var%in%c("freq","quanti"))),]
      res.globale$quali.var$contrib <- res.globale$var$contrib[-(1:sum(type.var%in%c("freq","quanti"))),]
	} else{
	  res.globale$quali.var$coord <- res.globale$var$coord
      res.globale$quali.var$cos2 <- res.globale$var$cos2
      res.globale$quali.var$contrib <- res.globale$var$contrib
	}
  }
  res.globale$call$quali.sup$barycentre <- sweep(crossprod(tab.comp[,unlist(ind.var.group[type%in%"n" || type%in%"m"])],as.matrix(data.pca)),1,apply(tab.comp[,unlist(ind.var.group[type%in%"n" || type%in%"m"])],2,sum),FUN="/")
  res.globale$quali.sup$coord <- sweep(crossprod(tab.comp[,unlist(ind.var.group[type%in%"n" || type%in%"m"])],res.globale$ind$coord),1,apply(tab.comp[,unlist(ind.var.group[type%in%"n" || type%in%"m"])],2,sum),FUN="/")
}
###  End handle missing values
    ncp <- min(ncp, nrow(res.globale$eig))
    call <- res.globale$call
    call$group <- group
    call$type <- type
    call$ncp <- ncp
    call$group.mod <- group.mod
    call$num.group.sup <- num.group.sup
    call$name.group <- name.group
    call$X <- base
    call$XTDC <- data
	call$nature.group <- nature.group
	call$nature.var <- nature.var
	call$list.type.var <- list.type.var
    contrib.group <- matrix(NA, length(group.actif), ncp)
    dimnames(contrib.group) <- list(name.group[group.actif], paste("Dim", c(1:ncp), sep = "."))
    dist2.group <- vector(length = length(group.actif))
    ind.var <- ind.var.sup <- 0
    for (g in 1:length(group.actif)) {
      if (group.mod[group.actif[g]]!=1) contrib.group[g, ] <- apply(res.globale$var$contrib[(ind.var + 1):(ind.var + group.mod[group.actif[g]]), 1:ncp]/100, 2, sum)
      else contrib.group[g, ] <- res.globale$var$contrib[ind.var + 1, 1:ncp]/100
      ind.var <- ind.var + group.mod[group.actif[g]]
      dist2.group[g] <- sum((res.separe[[group.actif[g]]]$eig[,1]/res.separe[[group.actif[g]]]$eig[1,1])^2)
    }
    coord.group <- t(t(contrib.group)*res.globale$eig[1:ncol(contrib.group),1])
    cos2.group <- coord.group^2/dist2.group
	
    if (!is.null(num.group.sup)){
      coord.group.sup <- matrix(NA, length(num.group.sup), ncp)
      dimnames(coord.group.sup) <- list(name.group[num.group.sup], paste("Dim", c(1:ncp), sep = "."))
      ind.gc <- 0
      for (gc in 1:length(num.group.sup)) {
        for (k in 1:ncp){
          if (is.null(ind.sup)) coord.group.sup[gc,k] <- funcLg(res.globale$ind$coord[,k,drop=FALSE],data.group.sup[,(ind.gc+1):(ind.gc+group.mod[num.group.sup[gc]]),drop=FALSE],ponderation.x=1/res.globale$eig[k,1],ponderation.y=ponderation.group.sup[(ind.gc+1):(ind.gc+group.mod[num.group.sup[gc]])],wt=row.w/sum(row.w))
          else coord.group.sup[gc,k] <- funcLg(res.globale$ind$coord[-ind.sup,k,drop=FALSE],data.group.sup[-ind.sup,(ind.gc+1):(ind.gc+group.mod[num.group.sup[gc]]),drop=FALSE],ponderation.x=1/res.globale$eig[k,1],ponderation.y=ponderation.group.sup[(ind.gc+1):(ind.gc+group.mod[num.group.sup[gc]])],wt=row.w/sum(row.w))
        }
        ind.gc <- ind.gc + group.mod[num.group.sup[gc]]
      }
    }
    Lg <- matrix(0, nbre.group+1, nbre.group+1)
    ind.gl <- 0
    for (gl in c(group.actif,num.group.sup)) {
        ind.gc <- 0
        for (gc in c(group.actif,num.group.sup)) {
            if (gc>=gl){
			  if (is.null(num.group.sup)) {
			    if (is.null(ind.sup)) Lg[gl, gc] <- Lg[gc, gl] <- funcLg(x=data[,ind.gl + (1:group.mod[gl]),drop=FALSE],y=data[,ind.gc + (1:group.mod[gc]),drop=FALSE],ponderation.x=ponderation[ind.gl + (1:group.mod[gl])],ponderation.y=ponderation[ind.gc + (1:group.mod[gc])],wt=row.w/sum(row.w))
				else Lg[gl, gc] <- Lg[gc, gl] <- funcLg(x=data[-ind.sup,ind.gl + (1:group.mod[gl]),drop=FALSE],y=data[-ind.sup,ind.gc + (1:group.mod[gc]),drop=FALSE],ponderation.x=ponderation[ind.gl + (1:group.mod[gl])],ponderation.y=ponderation[ind.gc + (1:group.mod[gc])],wt=row.w/sum(row.w))
			  } else {
			    if (is.null(ind.sup)) Lg[gl, gc] <- Lg[gc, gl] <- funcLg(x=cbind.data.frame(data,data.group.sup)[,ind.gl + (1:group.mod[gl]),drop=FALSE],y=cbind.data.frame(data,data.group.sup)[,ind.gc + (1:group.mod[gc]),drop=FALSE],ponderation.x=c(ponderation,ponderation.group.sup)[ind.gl + (1:group.mod[gl])],ponderation.y=c(ponderation,ponderation.group.sup)[ind.gc + (1:group.mod[gc])],wt=row.w/sum(row.w))
			    else Lg[gl, gc] <- Lg[gc, gl] <- funcLg(x=cbind.data.frame(data,data.group.sup)[-ind.sup,ind.gl +(1:group.mod[gl]),drop=FALSE],y=cbind.data.frame(data,data.group.sup)[-ind.sup,ind.gc + (1:group.mod[gc]),drop=FALSE],ponderation.x=c(ponderation,ponderation.group.sup)[ind.gl +(1:group.mod[gl])],ponderation.y=c(ponderation,ponderation.group.sup)[ind.gc + (1:group.mod[gc])],wt=row.w/sum(row.w))
			  }
			}
            ind.gc <- ind.gc + group.mod[gc]
        }
        ind.gl <- ind.gl + group.mod[gl]
    }	
   Lg[nbre.group+1,] <- Lg[,nbre.group+1] <- apply(Lg[group.actif,],2,sum)/res.globale$eig[1,1]
   Lg[nbre.group+1,nbre.group+1] <- sum(Lg[group.actif,nbre.group+1])/res.globale$eig[1,1]
    dist2.group <- diag(Lg)
    if (!is.null(num.group.sup)){
      dist2.group.sup <- dist2.group[num.group.sup]
      dist2.group <- dist2.group[-num.group.sup]
    }
    RV <- sweep(Lg, 2, sqrt(diag(Lg)), "/")
    RV <- sweep(RV, 1, sqrt(diag(Lg)), "/")
    rownames(Lg) <- colnames(Lg) <- rownames(RV) <- colnames(RV) <- c(name.group,"MFA")
    data.partiel <- vector(mode = "list", length = nbre.group)
    names(data.partiel) <- name.group
    ind.col <- 0
    for (g in 1:nbre.group) {
      if (g%in%group.actif){
        data.partiel[[g]] <- as.data.frame(matrix(res.globale$call$centre, nrow(data), ncol(data), byrow = TRUE, dimnames = dimnames(data)))
        data.partiel[[g]][, (ind.col + 1):(ind.col + group.mod[g])] <- data[, (ind.col + 1):(ind.col + group.mod[g])]
        ind.col <- ind.col + group.mod[g]
      }
    }
    res.ind.partiel <- vector(mode = "list", length = nbre.group)
    names(res.ind.partiel) <- name.group
    for (g in group.actif){
      Xis <- t(t(as.matrix(data.partiel[[g]]))-res.globale$call$centre) 
      Xis <- t(t(Xis)/res.globale$call$ecart.type)
      coord.ind.sup <- length(group.actif) * as.matrix(Xis)
      coord.ind.sup <- t(t(coord.ind.sup)*res.globale$call$col.w)
      coord.ind.sup <- crossprod(t(coord.ind.sup),res.globale$svd$V)
      res.ind.partiel[[g]]$coord.sup <- coord.ind.sup
    }
    cor.grpe.fact <- as.matrix(matrix(NA, length(group.actif), ncp))
    colnames(cor.grpe.fact) <- paste("Dim", c(1:ncp), sep = ".")
    rownames(cor.grpe.fact) <- name.group[group.actif]
    for (f in 1:ncp) {
        for (g in 1:length(group.actif))  cor.grpe.fact[g, f] <- cov.wt(cbind.data.frame(res.ind.partiel[[group.actif[g]]]$coord.sup[1:nb.actif, f], res.globale$ind$coord[, f]),wt=row.w/sum(row.w),method="ML",cor=TRUE)$cor[1,2]
    }
    It <- vector(length = ncp)
    for (g in group.actif)  It <- It + apply(res.ind.partiel[[g]]$coord.sup[1:nb.actif,]^2*row.w,2,sum)
    rap.inertie <- apply(res.globale$ind$coord^2*row.w,2,sum) * length(group.actif) / It 

    res.groupes <- list(Lg = Lg, RV = RV, coord = coord.group[, 1:ncp], contrib = contrib.group[, 1:ncp] * 100,  cos2 = cos2.group[, 1:ncp], dist2 = dist2.group[-length(dist2.group)], correlation = cor.grpe.fact[, 1:ncp])
    if (!is.null(num.group.sup)){
      res.groupes$coord.sup <- coord.group.sup[,1:ncp,drop=FALSE]
      res.groupes$cos2.sup <- coord.group.sup[,1:ncp,drop=FALSE]^2/dist2.group.sup
      res.groupes$dist2.sup <- dist2.group.sup
    }
  ####CHANGED THIS!!!! ------------------
  # OLD#
  # nom.ligne <- NULL
  # for (i in 1:nb.actif) nom.ligne <- c(nom.ligne, paste(rownames(base)[i], name.group[group.actif], sep = "."))
  # NEW#
  nom.ligne <- c(sapply( rownames(base)[1:nb.actif],  paste, name.group[group.actif], sep = "."))

    tmp <- array(0,dim=c(nrow(res.globale$ind$coord),ncp,length(group.actif)))
    for (g in 1:length(group.actif)) tmp[,,g] <- (res.ind.partiel[[group.actif[g]]]$coord.sup[1:nb.actif,1:ncp,drop=FALSE]-res.globale$ind$coord[,1:ncp,drop=FALSE])^2/length(group.actif)
## Ajour Avril 2011
    tmp <- tmp*row.w
    variab.auxil <- apply(tmp,2,sum)   ## attention, array
    tmp <- sweep(tmp,2,variab.auxil,FUN="/") * 100  ## attention, array
    inertie.intra.ind <- apply(tmp,c(1,2),sum)
  ####CHANGED THIS!!!! --------------------
  #OLD#
  #  inertie.intra.ind.partiel <- as.data.frame(matrix(NA, (nb.actif * length(group.actif)), ncp.tmp))
#  for (i in 1:nb.actif) inertie.intra.ind.partiel[((i - 1) * length(group.actif)  + 1):(i * length(group.actif)), ] <- t(tmp[i,1:ncp,])
  #NEW#
  inertie.intra.ind.partiel <- data.frame(do.call( rbind, lapply(1:dim(tmp)[1], function(x) t(tmp[x, 1:ncp, ]))))
  inertie.intra.ind.partiel <- as.matrix(inertie.intra.ind.partiel)
  # xxx <- as.data.frame(matrix(NA, (nb.actif * length(group.actif)), ncp))
  # colnames(inertie.intra.ind.partiel) <- colnames(xxx)

    rownames(inertie.intra.ind) <- rownames(res.globale$ind$coord)
    rownames(inertie.intra.ind.partiel) <- nom.ligne
    colnames(inertie.intra.ind) <- colnames(inertie.intra.ind.partiel) <- paste("Dim", c(1:ncp), sep = ".")
    tab.partial.axes <- matrix(NA, nb.actif, ncp * nbre.group)
    rownames(tab.partial.axes) <- rownames(data)[1:nb.actif]
    nom.axes <- paste("Dim", c(1:ncp), sep = "")
    nom.col <- NULL
    debut <- 0
    for (g in 1:nbre.group) {
      nom.col <- c(nom.col, paste(nom.axes, name.group[g],sep="."))
      nbcol <- min(ncp, ncol(res.separe[[g]]$ind$coord))
      tab.partial.axes[, (debut + 1):(debut + nbcol)] <- res.separe[[g]]$ind$coord[,1:nbcol]
      debut <- debut + ncp
    }
   colnames(tab.partial.axes) <- nom.col
    indice.col.NA <- which(!is.na(tab.partial.axes[1, ]))
    tab.partial.axes <- tab.partial.axes[, indice.col.NA]
    centre <- apply(tab.partial.axes, 2, moy.p, res.globale$call$row.w)
    tab.partial.axes <- t(t(tab.partial.axes)-centre)
    ecart.type <- apply(tab.partial.axes, 2, ec, res.globale$call$row.w)
    ecart.type[ecart.type <= 1e-08] <- 1
    tab.partial.axes <- t(t(tab.partial.axes)/ecart.type)
    coord.res.partial.axes <- t(tab.partial.axes*res.globale$call$row.w)
    coord.res.partial.axes <- crossprod(t(coord.res.partial.axes),res.globale$svd$U[,1:ncp])
	contrib.res.partial.axes <- coord.res.partial.axes*0
    debut <- 0
	for (g in 1:nbre.group) {
	  nbcol <- min(ncp, ncol(res.separe[[g]]$ind$coord))
	  if (g %in% group.actif) contrib.res.partial.axes[(debut + 1):(debut + nbcol),] <- coord.res.partial.axes[(debut + 1):(debut + nbcol),]^2*res.separe[[g]]$eig[1:nbcol,1]/res.separe[[g]]$eig[1,1]
	  debut <- debut + nbcol
    }
	contrib.res.partial.axes <- t(t(contrib.res.partial.axes)/apply(contrib.res.partial.axes,2,sum)) *100
	
    sigma <- apply(tab.partial.axes, 2, ec, res.globale$call$row.w)
    cor.res.partial.axes <- coord.res.partial.axes/sigma
    colnames(coord.res.partial.axes) <- paste("Dim", c(1:ncol(coord.res.partial.axes)), sep = ".")
    dimnames(contrib.res.partial.axes) <- dimnames(cor.res.partial.axes) <- dimnames (coord.res.partial.axes)
    summary.n <- as.data.frame(matrix(NA, 0, 4))
    colnames(summary.n) <- c("group", "variable", "modalite", "effectif")
    summary.c <- as.data.frame(matrix(NA, 0, 6))
    colnames(summary.c) <- c("group", "variable", "moyenne", "ecart.type", "minimum", "maximum")
   for (g in 1:nbre.group) {
        if ((type[g] == "c")||(type[g]=="f")) {
            statg <- as.data.frame(matrix(NA, ncol(res.separe[[g]]$call$X), 6))
            colnames(statg) <- c("group", "variable", "moyenne", "ecart.type", "minimum", "maximum")
            statg[, "group"] <- rep(g, nrow(statg))
            statg[, "variable"] <- colnames(res.separe[[g]]$call$X)
            statg[, "moyenne"] <- res.separe[[g]]$call$centre
            if (!is.null(res.separe[[g]]$call$ecart.type)) statg[, "ecart.type"] <- res.separe[[g]]$call$ecart.type
            statg[, "minimum"] <- apply(res.separe[[g]]$call$X, 2, min)
            statg[, "maximum"] <- apply(res.separe[[g]]$call$X, 2, max)
            if (!is.null(res.separe[[g]]$call$ecart.type)) statg[, -c(1, 2)] <- round(statg[, -c(1, 2)], digits = 2)
            else statg[, -c(1, 2,4)] <- round(statg[, -c(1, 2,4)], digits = 2)
            summary.c <- rbind(summary.c, statg)
        }
        if (type[g]=="n") {
            if(is.null(excl[[g]])) statg <- as.data.frame(matrix(NA, length(res.separe[[g]]$call$marge.col), 4))
            else statg <- as.data.frame(matrix(NA, length(res.separe[[g]]$call$marge.col[-excl[[g]]]), 4))
            colnames(statg) <- c("group", "variable", "modalite", "effectif")
            statg[, "group"] <- rep(g, nrow(statg))
            res.separe[[g]]$call$X <- as.data.frame(res.separe[[g]]$call$X)
            nb.var <- ncol(res.separe[[g]]$call$X)
            nb.mod <- NULL
            nom.mod <- NULL
            nb.mod.orig <- NULL
            nb.mod.high <- 0
            for (v in 1:nb.var) {
              if(is.null(excl[[g]])) {
                nb.mod <- c(nb.mod, nlevels(res.separe[[g]]$call$X[, v]))
              } else {
                nb.mod.orig <- c(nb.mod.orig, nlevels(res.separe[[g]]$call$X[, v]))
                nb.mod.low <- nb.mod.high
                nb.mod.high <- nb.mod.high + nlevels(res.separe[[g]]$call$X[, v])
                nb.mod.rm.select <- (nb.mod.low < res.separe[[g]]$call$excl & res.separe[[g]]$call$excl <= nb.mod.high)
                nb.mod.rm <- res.separe[[g]]$call$excl[nb.mod.rm.select]
                nb.mod <- c(nb.mod, (nb.mod.orig[v]-length(nb.mod.rm)))
              }
              nom.mod <- c(nom.mod, levels(res.separe[[g]]$call$X[, v]))
            }
            if(!is.null(excl[[g]])) nom.mod <- nom.mod[-excl[[g]]]
            statg[, "variable"] <- rep(colnames(res.separe[[g]]$call$X), nb.mod)
            statg[, "modalite"] <- nom.mod
            if(is.null(excl[[g]])) statg[, "effectif"] <- res.separe[[g]]$call$marge.col * nbre.ind * nb.var
            else statg[, "effectif"] <- res.separe[[g]]$call$marge.col[-excl[[g]]] * nbre.ind * nb.var
            summary.n <- rbind(summary.n, statg)			
        }
###not managed
        if (type[g]=="m") {
		    nbvar.num <- sum(sapply(res.separe[[g]]$call$X,is.numeric))
            statg <- as.data.frame(matrix(NA, nbvar.num, 6))
            colnames(statg) <- c("group", "variable", "moyenne", "ecart.type", "minimum", "maximum")
            statg[, "group"] <- rep(g, nrow(statg))
            statg[, "variable"] <- colnames(res.separe[[g]]$call$X[,which(sapply(res.separe[[g]]$call$X,is.numeric)),drop=FALSE])
            statg[, "moyenne"] <- res.separe[[g]]$call$centre[1:nbvar.num]
            if (!is.null(res.separe[[g]]$call$ecart.type)) statg[, "ecart.type"] <- res.separe[[g]]$call$ecart.type[1:nbvar.num]
            statg[, "minimum"] <- apply(res.separe[[g]]$call$X[,which(sapply(res.separe[[g]]$call$X,is.numeric)),drop=FALSE], 2, min)
            statg[, "maximum"] <- apply(res.separe[[g]]$call$X[,which(sapply(res.separe[[g]]$call$X,is.numeric)),drop=FALSE], 2, max)
            if (!is.null(res.separe[[g]]$call$ecart.type)) statg[, -c(1, 2)] <- round(statg[, -c(1, 2)], digits = 2)
            else statg[, -c(1, 2,4)] <- round(statg[, -c(1, 2,4)], digits = 2)
            summary.c <- rbind(summary.c, statg)
            if(is.null(excl[[g]])) statg <- as.data.frame(matrix(NA, length(res.separe[[g]]$call$prop), 4))
            else statg <- as.data.frame(matrix(NA, length(res.separe[[g]]$call$prop[-excl[[g]]]), 4))
            colnames(statg) <- c("group", "variable", "modalite", "effectif")
            statg[, "group"] <- rep(g, nrow(statg))
            nb.var <- sum(!sapply(res.separe[[g]]$call$X,is.numeric))
            nb.mod <- unlist(sapply(res.separe[[g]]$call$X[,which(!sapply(res.separe[[g]]$call$X,is.numeric)),drop=FALSE],nlevels))
            nom.mod <- unlist(sapply(res.separe[[g]]$call$X[,which(!sapply(res.separe[[g]]$call$X,is.numeric)),drop=FALSE],levels))
            if(!is.null(excl[[g]])){
			  nb.mod <- nb.mod - length(excl[[g]])
			  nom.mod <- nom.mod[-excl[[g]]]
			}  
            statg[, "variable"] <- rep(colnames(res.separe[[g]]$call$X[,which(!sapply(res.separe[[g]]$call$X,is.numeric)),drop=FALSE]), nb.mod)
            statg[, "modalite"] <- nom.mod
            if(is.null(excl[[g]])) statg[, "effectif"] <- res.separe[[g]]$call$prop * nbre.ind 
            else statg[, "effectif"] <- res.separe[[g]]$call$marge.col[-excl[[g]]] * nbre.ind 
            summary.n <- rbind(summary.n, statg)			
        }
    }
   eig <- res.globale$eig[1:ncp.tmp,]
  ## CHANGED THIS!!!! ---------------
  # nom.ligne <- NULL
  # for (i in 1:nbre.ind) {
  #   ind.tmp <- rownames(base)[i]
  #   nom.ligne <- c(nom.ligne, paste(ind.tmp, name.group[group.actif], sep = "."))
  # }
  
  nom.ligne <- c( sapply( rownames(base), paste, name.group[group.actif], sep = "."))
    coord.ind.partiel <- matrix(NA, (nbre.ind * length(group.actif)), ncp)
    rownames(coord.ind.partiel) <- nom.ligne
    colnames(coord.ind.partiel) <- paste("Dim", c(1:ncp), sep = ".")
    coord.ind <- rbind(res.globale$ind$coord[, 1:ncp,drop=FALSE],res.globale$ind.sup$coord[, 1:ncp,drop=FALSE])
    cos2.ind <- rbind(res.globale$ind$cos2[, 1:ncp,drop=FALSE],res.globale$ind.sup$cos2[, 1:ncp,drop=FALSE])
    contrib.ind <- res.globale$ind$contrib[, 1:ncp,drop=FALSE]
    liste.ligne <- seq(1, nbre.ind * length(group.actif), by = length(group.actif))
    for (g in 1:length(group.actif)) coord.ind.partiel[liste.ligne+g-1, ] <- res.ind.partiel[[group.actif[g]]]$coord.sup[, 1:ncp,drop=FALSE]
    if (!is.null(ind.sup)) {
      res.ind.sup <- list(coord = coord.ind[(nb.actif+1):nrow(coord.ind),,drop=FALSE], cos2 = cos2.ind[(nb.actif+1):nrow(coord.ind),,drop=FALSE], coord.partiel = coord.ind.partiel[(length(group.actif)*nb.actif+1):nrow(coord.ind.partiel),,drop=FALSE])
      res.ind <- list(coord = coord.ind[1:nb.actif,,drop=FALSE], contrib = contrib.ind, cos2 = cos2.ind[1:nb.actif,,drop=FALSE], within.inertia = inertie.intra.ind[1:nb.actif,1:ncp,drop=FALSE], coord.partiel = coord.ind.partiel[1:(length(group.actif)*nb.actif),,drop=FALSE], within.partial.inertia = inertie.intra.ind.partiel[1:(length(group.actif)*nb.actif),1:ncp,drop=FALSE] )
    }
    else res.ind <- list(coord = coord.ind, contrib = contrib.ind, cos2 = cos2.ind, within.inertia = inertie.intra.ind[,1:ncp,drop=FALSE], coord.partiel = coord.ind.partiel, within.partial.inertia = inertie.intra.ind.partiel[,1:ncp,drop=FALSE])
    res.quali.var <- res.quali.var.sup <- NULL
    bool.act <- FALSE
    bool.sup <- FALSE
    if (!is.null(ind.quali)) {
      coord.quali <- res.globale$quali.sup$coord[, 1:ncp,drop=FALSE]
      cos2.quali <- res.globale$quali.sup$cos2[, 1:ncp,drop=FALSE]
      val.test.quali <- res.globale$quali.sup$v.test[, 1:ncp,drop=FALSE]
      contrib.quali <- coord.quali * 0
	  commun <- intersect(rownames(res.globale$var$contrib),rownames(contrib.quali))
      if (!is.null(commun)) contrib.quali[commun,] <- res.globale$var$contrib[commun,1:ncp,drop=FALSE]
      barycentre <- res.globale$call$quali.sup$barycentre
      coord.quali.partiel <- matrix(NA, (nrow(barycentre) * length(group.actif)), ncp)
      nom.ligne.bary <- NULL
      for (q in 1:nrow(barycentre)) {
        ind.tmp <- rownames(barycentre)[q]
        nom.ligne.bary <- c(nom.ligne.bary, paste(ind.tmp, name.group[group.actif], sep = "."))
      }
      rownames(coord.quali.partiel) <- nom.ligne.bary
      liste.ligne <- seq(1, (nrow(barycentre) * length(group.actif) ), by = length(group.actif))
      inertie.intra.cg.partiel <- matrix(NA, (nrow(barycentre) * length(group.actif) ), ncp)
      tmp <- array(0,dim=c(nrow(res.globale$quali.sup$coord),ncp,length(group.actif)))
      ind.col <- 0
for (g in 1:length(group.actif)) {
        cg.partiel <- as.data.frame(matrix(res.globale$call$centre[sum(group.mod[group.actif[1:length(group.actif)]])], nrow(barycentre), sum(group.mod[group.actif[1:length(group.actif)]]), byrow = TRUE, dimnames = dimnames(barycentre[,1:sum(group.mod[group.actif[1:length(group.actif)]]),drop=FALSE])))
        cg.partiel[, (ind.col + 1):(ind.col + group.mod[group.actif[g]])] <- barycentre[, (ind.col + 1):(ind.col + group.mod[group.actif[g]])]
        ind.col <- ind.col + group.mod[group.actif[g]]
        Xis <- t((t(cg.partiel)-res.globale$call$centre)/res.globale$call$ecart.type)
        coord.quali.sup <- length(group.actif) * as.matrix(Xis)
        coord.quali.sup <- t(t(coord.quali.sup)*res.globale$call$col.w)
        coord.quali.sup <- crossprod(t(coord.quali.sup),res.globale$svd$V)
        coord.quali.partiel[liste.ligne + g - 1, ] <- coord.quali.sup[,1:ncp]
        tmp[,,g] <- (coord.quali.sup[,1:ncp,drop=FALSE] - res.globale$quali.sup$coord[,1:ncp,drop=FALSE])^2 / length(group.actif) 
      }
      colnames(coord.quali.partiel) <-  paste("Dim", 1:ncp, sep = ".")
      tmp <- sweep(tmp,2,variab.auxil,FUN="/") * 100   ### attention array
### modif mais attention, si changement dans PCA, remettre ?
##      tmp <- sweep(tmp,1,res.globale$call$quali.sup$nombre,FUN="*")
      tmp <- sweep(tmp,1,poids.bary*sum(row.w),FUN="*")  ### attention array
      inertie.intra.cg <- apply(tmp,c(1,2),sum)
      for (i in 1:nrow(barycentre)) inertie.intra.cg.partiel[((i - 1) * length(group.actif)  + 1):(i * length(group.actif)), ] <- t(tmp[i,1:ncp,])
      rownames(inertie.intra.cg) <- rownames(res.globale$quali.sup$coord)
      rownames(inertie.intra.cg.partiel) <- nom.ligne.bary
      colnames(inertie.intra.cg) <- colnames(inertie.intra.cg.partiel) <- paste("Dim", c(1:ncp), sep = ".")

      ind.col <- 0
      ind.col.act <- NULL
      ind.col.sup <- NULL
      ind.excl <- NULL
      ind.excl.act <- NULL
###not managed
      for (g in 1:nbre.group) {
        if (type[g] =="n" | type[g]=="m"){
          if (g%in%num.group.sup) ind.col.sup <- c(ind.col.sup, ind.col+(1:(group.mod[g]-sum(list.type.var[[g]]=="quanti.sup"))))
          else ind.col.act <- c(ind.col.act, ind.col+(1:(group.mod[g]-sum(list.type.var[[g]]=="quanti"))))
          if(!is.null(excl[[g]])) {
            ind.excl <- ind.col + excl[[g]]
            ind.excl.act <- c(ind.excl.act, ind.excl)
          }
          ind.col <- ind.col + group.mod[g]-sum(list.type.var[[g]]=="quanti")-sum(list.type.var[[g]]=="quanti.sup")
        }
      }
	  
      if(!is.null(ind.excl.act)) ind.col.act <- ind.col.act[-ind.excl.act]

      if (!is.null(ind.col.sup)) {
            coord.quali.sup <- coord.quali[ind.col.sup,,drop=FALSE]
            cos2.quali.sup <- cos2.quali[ind.col.sup,,drop=FALSE]
            val.test.quali.sup <- val.test.quali[ind.col.sup,,drop=FALSE]
            coord.quali.partiel.sup <- coord.quali.partiel[unlist(lapply(ind.col.sup, function(k) seq(length(group.actif)*(k-1)+1,length=length(group.actif)))),]
            inertie.intra.cg.sup <- inertie.intra.cg[ind.col.sup,1:ncp]
            inertie.intra.cg.partiel.sup <- inertie.intra.cg.partiel[unlist(lapply(ind.col.sup, function(k) seq(length(group.actif)*(k-1)+1,length=length(group.actif)))),1:ncp]
            bool.sup <- TRUE
       }
       if (!is.null(ind.col.act)) {
            coord.quali.act <- coord.quali[ind.col.act,,drop=FALSE]
            contrib.quali.act <- contrib.quali[ind.col.act,,drop=FALSE]
			val.test.quali.act <- NULL
            if (is.null(tab.comp)){
			  cos2.quali.act <- cos2.quali[ind.col.act,,drop=FALSE]
              val.test.quali.act <- val.test.quali[ind.col.act,,drop=FALSE]
			} else {
			  cos2.quali.act <- res.globale$quali.var$cos2
			}
            coord.quali.partiel.act <- coord.quali.partiel[unlist(lapply(ind.col.act, function(k) seq(length(group.actif)*(k-1)+1,length=length(group.actif)))),]
            inertie.intra.cg.act <- inertie.intra.cg[ind.col.act,1:ncp]
            inertie.intra.cg.partiel.act <- inertie.intra.cg.partiel[unlist(lapply(ind.col.act, function(k) seq(length(group.actif)*(k-1)+1,length=length(group.actif)))),1:ncp]
            bool.act <- TRUE
        }
      if (bool.act) res.quali.var <- list(coord = coord.quali.act, contrib = contrib.quali.act, cos2 = cos2.quali.act, v.test = val.test.quali.act, coord.partiel = coord.quali.partiel.act, within.inertia = inertie.intra.cg.act, within.partial.inertia = inertie.intra.cg.partiel.act)
      if (bool.sup) res.quali.var.sup <- list(coord = coord.quali.sup, cos2 = cos2.quali.sup, v.test = val.test.quali.sup, coord.partiel = coord.quali.partiel.sup, within.inertia = inertie.intra.cg.sup, within.partial.inertia = inertie.intra.cg.partiel.sup)
    }
    indice.quanti <- NULL
    indice.freq <- NULL
    num.tmp <- 0
    for (g in group.actif) {
        if (type[g] == "c")  indice.quanti <- c(indice.quanti, c((num.tmp + 1):(num.tmp + group.mod[g])))
        if (type[g] == "m"){
		  if (g==1) indice.quanti <- indice.quanti <- c(indice.quanti,c((num.tmp+1):(num.tmp+sum(sapply(base[,1:group[1],drop=FALSE],is.numeric)))))
		  else indice.quanti <- c(indice.quanti, c((num.tmp+1):(num.tmp+sum(sapply(base[,(sum(group[1:(g-1)])+1):sum(group[1:g]),drop=FALSE],is.numeric)))))
		}
        if (type[g] == "f")  indice.freq <- c(indice.freq, c((num.tmp + 1):(num.tmp + group.mod[g])))
        num.tmp <- num.tmp + group.mod[g]
    }
    res.quanti.var <- NULL
    if (!is.null(indice.quanti)){
      coord.quanti.var <- res.globale$var$coord[indice.quanti,1:ncp,drop=FALSE]
      cos2.quanti.var <- res.globale$var$cos2[indice.quanti, 1:ncp,drop=FALSE]
      contrib.quanti.var <- res.globale$var$contrib[indice.quanti, 1:ncp,drop=FALSE]
      cor.quanti.var <- res.globale$var$cor[indice.quanti, 1:ncp,drop=FALSE]
      res.quanti.var <- list(coord = coord.quanti.var, contrib = contrib.quanti.var, cos2 = cos2.quanti.var, cor = cor.quanti.var)
    }
    res.freq <- NULL
    if (!is.null(indice.freq)){
      coord.freq <- res.globale$var$coord[indice.freq,1:ncp,drop=FALSE]
      cos2.freq <- res.globale$var$cos2[indice.freq, 1:ncp,drop=FALSE]
      contrib.freq <- res.globale$var$contrib[indice.freq, 1:ncp,drop=FALSE]
      res.freq <- list(coord = coord.freq, contrib = contrib.freq, cos2 = cos2.freq)
    }

    res.quanti.var.sup <- NULL
    res.freq.sup <- NULL
    if (!is.null(num.group.sup)){
      num.tmp <- 0
      indice.quanti <- NULL
      indice.freq <- NULL
      for (g in num.group.sup) {
        if (type[g] == "c") indice.quanti <- c(indice.quanti, num.tmp + (1:group.mod[g]))
        if (type[g] == "m") indice.quanti <- c(indice.quanti, num.tmp + (1:sum(list.type.var[[g]]=="quanti.sup")))
        if (type[g]=="f") indice.freq <- c(indice.freq, c((num.tmp + 1):(num.tmp + group.mod[g])))
        num.tmp <- num.tmp + group.mod[g]
      }
      if (!is.null(indice.quanti)){
        coord.quanti.var.sup <- res.globale$quanti.sup$coord[indice.quanti,1:ncp,drop=FALSE]
        cos2.quanti.var.sup <- res.globale$quanti.sup$cos2[indice.quanti, 1:ncp,drop=FALSE]
        cor.quanti.var.sup <- res.globale$quanti.sup$cor[indice.quanti, 1:ncp,drop=FALSE]
        res.quanti.var.sup <- list(coord = coord.quanti.var.sup, cos2 = cos2.quanti.var.sup, cor = cor.quanti.var.sup)
      }
      if (!is.null(indice.freq)){
        coord.freq.sup <- res.globale$quanti.sup$coord[indice.freq,1:ncp,drop=FALSE]
        cos2.freq.sup <- res.globale$quanti.sup$cos2[indice.freq, 1:ncp,drop=FALSE]
        res.freq.sup <- list(coord = coord.freq.sup, cos2 = cos2.freq.sup)
      }
    }

    aux <- res.separe[[1]]$ind$coord
    name.aux <- paste(colnames(res.separe[[1]]$ind$coord),name.group[1],sep=".")
    for (g in 2:nbre.group) {
      aux <- cbind(aux,res.separe[[g]]$ind$coord)
      name.aux <- c(name.aux,paste(colnames(res.separe[[g]]$ind$coord),name.group[g],sep="."))
    }
cor.partial.axes <- cov.wt(aux,wt=row.w/sum(row.w),method="ML",cor=TRUE)$cor
    dimnames(cor.partial.axes) <- list(name.aux,name.aux)
    res.partial.axes <- list(coord = coord.res.partial.axes[, 1:ncp], cor = cor.res.partial.axes[, 1:ncp], contrib = contrib.res.partial.axes[, 1:ncp], cor.between = cor.partial.axes)
    resultats <- list(separate.analyses = res.separe, eig = eig, group = res.groupes, 
        inertia.ratio = rap.inertie[1:ncp], ind = res.ind)
    if (!is.null(ind.sup)) resultats$ind.sup <- res.ind.sup
    if (!is.null(c(res.quanti.var,res.quanti.var.sup))) resultats$summary.quanti <- summary.c
    if (!is.null(c(bool.act,bool.sup))) resultats$summary.quali <- summary.n
    if (!is.null(res.quanti.var)) resultats$quanti.var <- res.quanti.var
    if (!is.null(res.quanti.var.sup)) resultats$quanti.var.sup <- res.quanti.var.sup
    if (!is.null(res.freq)) resultats$freq <- res.freq
    if (!is.null(res.freq.sup)) resultats$freq.sup <- res.freq.sup
    if (bool.act) resultats$quali.var <- res.quali.var
    if (bool.sup) resultats$quali.var.sup <- res.quali.var.sup
    resultats$partial.axes <- res.partial.axes
    resultats$call <- call
	resultats$call$call <- match.call()
    resultats$global.pca <- res.globale
    class(resultats) <- c("MFA", "list")

    if (graph & (ncp>1)){
      if (bool.act | bool.sup){
        cg.plot.partial <- NULL
        if (!is.null(resultats["quali.var"]$quali.var)){
          max.inertia <- order(apply(resultats$quali.var$within.inertia[,1:2],1,sum))
          cg.plot.partial <- rownames(resultats$quali.var$coord)[max.inertia[1:length(max.inertia)]]
        }
        if (!is.null(resultats$quali.var.sup)){
          max.inertia <- order(apply(resultats$quali.var.sup$within.inertia[,1:2],1,sum))
          cg.plot.partial <- c(cg.plot.partial,rownames(resultats$quali.var.sup$coord)[max.inertia[1:length(max.inertia)]])
        }
 print(plot.MFA(resultats,choix="ind",invisible="ind",partial=cg.plot.partial,habillage="group",axes=axes,new.plot=TRUE))
      }
      max.inertia <- order(apply(resultats$ind$within.inertia[,1:2],1,sum))
      print(plot.MFA(resultats,choix="axes",habillage="group",axes=axes,new.plot=TRUE,shadowtext=TRUE))
      print(plot.MFA(resultats,choix="ind",invisible="quali",partial=rownames(resultats$ind$coord)[max.inertia[c(1:2,nrow(resultats$ind$coord)-1,nrow(resultats$ind$coord))]],habillage="group",axes=axes,new.plot=TRUE))
      if ("c"%in%type) print(plot.MFA(resultats,choix="var",habillage="group",axes=axes,new.plot=TRUE,shadowtext=TRUE))
      if ("f"%in%type) print(plot.MFA(resultats,choix="freq",habillage="group",axes=axes,new.plot=TRUE))
      print(plot.MFA(resultats,choix="ind",invisible="quali",habillage = "ind",axes=axes,new.plot=TRUE,col.hab=1+3*(1:nbre.ind)%in%ind.sup))
      print(plot.MFA(resultats,choix="group",axes=axes,new.plot=TRUE))
    }
    return(resultats)
}

Try the FactoMineR package in your browser

Any scripts or data that you put into this service are public.

FactoMineR documentation built on May 29, 2024, 3:36 a.m.