# Copyright 2006-7 Université de Liège, Nicolas Mazziotta, 
# sous les termes de la licence GPL2

# Analyse des spécificités 

# FIXME: la fonction est vraiment lente
spec <- function (df, col, values=NULL,vars=NULL, threshold=.001) {
	# soit un data.frame "df". On explore les modalités de la colonne "col" (des facteurs) en cherchant les spécificités des autres colonnes du "df"
	# pour chaque modalité de "col". On peut restreindre aux modalités "values" et on peut restreindre les colonnes de "df" à "vars". Le "threshold" est 
	# le seuil de probabilité retenu.

	# nombre d'individus
	n<-dim(df)[1]
	# valeurs de la col examinée
	l<-levels(df[,col]);
	if (is.null(vars)) { vars<-colnames(df)	}

	out<-NULL

	# pour chaque valeur de la col
	for (i in 1:length(l)) {
		# pour chaque caractère comparé
		table<-NULL
		# sous-tableau
		subdf<-df[df[,col]==l[i],]
		if (!is.null(values)) {
			subdf<-subdf[subdf[,col] %in% values,]
		}
		
		for (j in 1:length(vars)) {
			if (col == vars[j]) { next }
			vals<-levels(df[,vars[j]])
			# pour chaque modalité
			for (ii in 1:length(vals)) {
		
				# fréquence de la modalité dans le df
				nj<-dim(df[df[,vars[j]]==vals[ii],])[1]
				nk<-dim(subdf)[1]

				nkj<-dim(subdf[subdf[,vars[j]]==vals[ii],])[1]
				
				h<-dhyper(1:nkj, nj, n-nj, nk)
				hh<-cbind(cbind(1:nkj),h)

				mode<-max(hh[hh[,2]==max(hh[,2]),1])
				low<-TRUE;
				
				if (nkj > mode) {low<-FALSE}
				p<-phyper(nkj, nj, n-nj, nk, low=low)
				if (p > threshold) { next }

				int<-round(nkj/nk*100, dig=2)
				glo<-round(nj/n*100, dig=2)
				spec_dir<-sign(int-glo);
				if ( spec_dir == 0 ) { 
					if (low) {spec_dir<-c(1)} 
					else {spec_dir<-c(-1)}
				} 
				line<-rbind(c(int, glo, nkj, nj, round(nkj/nj*100, 2), p, spec_dir))
				rownames(line)<-paste(vars[j],vals[ii],sep=".")
				table<-rbind(table, line)
			}
		}
		table<-data.frame(table)
		
		spec<-NULL
		if (! is.null(table) & dim(table)[2]){
			colnames(table)<-c("%IN","%GL", "FIN", "FGL" ,"RAT", "PRB", "SPC")
			# %age interne, global, fréquence interne, globale, probabilité
			# spécificité positive ou négative
			s<-table;
			attract<-s[s[,7]==1,];
			repuls <-s[s[,7]!=1,];

			if (dim(s)[1]>1) {s<-s[order(-s[,7],-s[,5]*s[,7]),]}
			# s[,"PRB"]<-round(s[,"PRB"], dig=3)
			spec$tab<-s
			class(spec)<-c('spec')
			out[[l[i]]]<-spec

		}
	}
	return(out)
}


#vim: encoding=utf-8 fileencoding=utf-8
