# CA script by Elie Maza and Christian Chervin, AgroToulouse
# Loading packages

library (FactoMineR)

# Loading data

DataForCA=read.table("DataForCA.csv",sep=";",header=TRUE,stringsAsFactors=TRUE)

summary(DataForCA) # to check the raw data (number of replicates, min and max...)

# CA of all data

product=levels(DataForCA$Product)
product # to check the different names of the products

tab=apply(DataForCA[DataForCA$Product==product[1],-c(1,2)],2,sum)
tab # to print the sums for the first product

for (i in 2:length(product)) {
	tab=rbind(tab,apply(DataForCA[DataForCA$Product==product[i],-c(1,2)],2,sum))
}
rownames(tab)=product
tab # to print the full contingency table

res=CA(tab)

pdf("fig_AC.pdf",width=7,height=7)
par(mar=c(4.1,4.1,1.1,1.1)) # edit graph margins
# drawing ellipses at 95% confidence
ellipseCA(res,ellipse="row",autoLab="yes",cex=1.25,cex.axis=1.25,cex.lab=1.25,shadowtext=TRUE,title="",xlim=range(c(res$col$coord[,1],res$row$coord[,1])),ylim=range(c(res$col$coord[,2],res$row$coord[,2])))
dev.off()

# Drawing p-value against sample size

tickInterval=floor(max(DataForCA$NumPanel)/11)
Size=2:floor(max(DataForCA$NumPanel)/tickInterval)*tickInterval # setting the panelist range you want to test, as a function of your panelist number

NbReps=300 # number of random draws for each size

pVal=matrix(nrow=NbReps,ncol=length(Size))

for (i in 1:length(Size)) {
	for (j in 1:NbReps) {
		panelists=sample(1:(nrow(DataForCA)/length(product)),Size[i])
		data.selected.panelists=DataForCA[DataForCA$Product==product[1],-c(1,2)]
		data.selected.panelists=data.selected.panelists[panelists,]
		tab=apply(data.selected.panelists,2,sum)
		for (k in 2:length(product)) {
			data.selected.panelists=DataForCA[DataForCA$Product==product[k],-c(1,2)]
			data.selected.panelists=data.selected.panelists[panelists,]
			tab=rbind(tab,apply(data.selected.panelists,2,sum))
		}
		rownames(tab)=product
		pVal[j,i]=chisq.test(tab)$p.value
	}
}

boxplot(pVal,names=Size,cex.axis=1.25,xlab="Number of panelists",ylab="p-value",cex.lab=1.25,ylim=c(0,1))
abline(h=0.05,col="red",lwd=1.25)

pdf("fig_pVals.pdf",width=7,height=4)
par(mar=c(4.1,4.1,1.1,1.1))
boxplot(pVal,names=Size,cex.axis=1.25,xlab="Number of panelists",ylab="p-value",cex.lab=1.25,ylim=c(0,1))
abline(h=0.05,col="red",lwd=1.25)
dev.off()

#---------
# The End
#---------