##################################################################################################
# Calcula la diferencia de medias entre los datos de un conjunto de columnas en función de una conjunto 
# de factores, y su error estandar
# sdata: dataframe con los datos
# cmean: vector con los índices o los nombres de las columnas con los datos
# cfact: vector con los índices o los nombres de las columnas con los factores
# wght: índice o nombre de la columna con el peso de las filas
# brr: vector con los índices o los nombres de las columnas con los pesos replicados
# Devuelve una matriz con una columna para cada diferencia de medias entre dos conjuntos de valores
# de cada una de las columnas de cmean, agrupada por una de las combinaciones diferentes de dos de los 
# valores de cada factor de cfact
# Para cada columna en cmean se crean dos filas, la primera para las diferencias y la segunda para 
# los errores típicos
##################################################################################################
# Calculates the difference of means between the data in a set of columns in function of a set of
# factors, and their standard errors
# sdata: dataframe with the data
# cmean: vector with the indices or names of the columns with the data
# cfact: vector with the indices or names of the columns with the factors
# wght: index or name of the column with the row weight
# brr: vector with the indices or names of the columns with the replicate weights
# Returns a matrix with one column for each difference of means between two sets of values of each one
# of the columns in cmean, grouped by each different combination of two values of each factor in cfact
# There are two rows for each column in cmean, the first for the differences and the second for the
# standard errors

wght_meandiff<-function(sdata,cmean,cfact,wght,brr) {
	nc<-0;
	for (i in 1:length(cfact)) {		
		for (j in 1:(length(levels(as.factor(sdata[,cfact[i]])))-1)) {
			for(k in (j+1):length(levels(as.factor(sdata[,cfact[i]])))) {
				nc <- nc + 1;
			}
		}
	}
	mmeans<-matrix(ncol=nc,nrow=2 * length(cmean));
	mmeans[,]<-0;
	cn<-c();
	for (i in 1:length(cfact)) {
		for (j in 1:(length(levels(as.factor(sdata[,cfact[i]])))-1)) {
			for(k in (j+1):length(levels(as.factor(sdata[,cfact[i]])))) {
				cn<-c(cn, paste(names(sdata)[cfact[i]],levels(as.factor(sdata[,cfact[i]]))[j],levels(as.factor(sdata[,cfact[i]]))[k],sep="-"));
			}
		}
	}
	colnames(mmeans)<-cn;
	rn<-c();
	for (i in 1:length(cmean)) {
		rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF",sep="-"));
		rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF","SE",sep="-"));
	}	
	rownames(mmeans)<-rn;	
	ic<-1;
	for(f in 1:length(cfact)) {
		for (l in 1:(length(levels(as.factor(sdata[,cfact[f]])))-1)) {
			for(k in (l+1):length(levels(as.factor(sdata[,cfact[f]])))) {
				ir<-1;
				rfact1<-sdata[,cfact[f]]==levels(as.factor(sdata[,cfact[f]]))[l];
				rfact2<-sdata[,cfact[f]]==levels(as.factor(sdata[,cfact[f]]))[k];
				swght1<-sum(sdata[rfact1,wght]);
				swght2<-sum(sdata[rfact2,wght]);
				for (i in 1:length(cmean)) {				
					mmeans[ir,ic]<-(sum(sdata[rfact1,wght]*sdata[rfact1,cmean[i]])/swght1) - (sum(sdata[rfact2,wght]*sdata[rfact2,cmean[i]])/swght2);
					for (j in 1:length(brr)) {
						sbrr1<-sum(sdata[rfact1,brr[j]]);
						sbrr2<-sum(sdata[rfact2,brr[j]]);
						mmbrj<-(sum(sdata[rfact1,brr[j]]*sdata[rfact1,cmean[i]])/sbrr1) - (sum(sdata[rfact2,brr[j]]*sdata[rfact2,cmean[i]])/sbrr2);
						mmeans[ir+1,ic]<-mmeans[ir+1,ic] + (mmbrj - mmeans[ir,ic])^2;
					}		
					mmeans[ir+1,ic]<-sqrt((mmeans[ir+1,ic] * 4) / length(brr));
					ir<-ir+2;
				}
				ic<-ic + 1;
			}
		}
	}
	return(mmeans);
}

##################################################################################################
# Calcula la diferencia de medias entre países de un conjunto de columnas, y su error estandar
# sdata: dataframe con los datos
# cmean: vector con los índices o los nombres de las columnas con los datos
# cnt: índice o nombre de la columna con el país
# wght: índice o nombre de la columna con el peso de las filas
# brr: vector con los índices o los nombres de las columnas con los pesos replicados
# Devuelve una matriz con una columna para cada diferencia de medias entre dos conjuntos de valores
# de cada una de las columnas de cmean, para cada combinación de dos paises
# Para cada columna en cmean se crean dos filas, la primera para las diferencias y la segunda para 
# los errores típicos
##################################################################################################
# Calculates the difference of means between countries of a set of columns, and their standard errors
# sdata: dataframe with the data
# cmean: vector with the indices or names of the columns with the data
# cnt: index or name of the column with the country
# wght: index or name of the column with the row weight
# brr: vector with the indices or names of the columns with the replicate weights
# Returns a matrix with one column for each difference of means between two sets of values of each one
# of the columns in cmean, for each different combination of two countries
# There are two rows for each column in cmean, the first for the differences and the second for the
# standard errors

wght_meandiffbycnt1<-function(sdata,cmean,cnt,wght,brr) {
	nc<-0;
	for (j in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
		for(k in (j+1):length(levels(as.factor(sdata[,cnt])))) {
			nc <- nc + 1;
		}
	}
	mmeans<-matrix(ncol=nc,nrow=2 * length(cmean));
	mmeans[,]<-0;
	cn<-c();
	for (j in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
		for(k in (j+1):length(levels(as.factor(sdata[,cnt])))) {
			cn<-c(cn, paste(names(sdata)[cnt],levels(as.factor(sdata[,cnt]))[j],levels(as.factor(sdata[,cnt]))[k],sep="-"));
		}
	}
	colnames(mmeans)<-cn;
	rn<-c();
	for (i in 1:length(cmean)) {
		rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF",sep="-"));
		rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF","SE",sep="-"));
	}	
	rownames(mmeans)<-rn;	
	ic<-1;
	for (l in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
		for(k in (l+1):length(levels(as.factor(sdata[,cnt])))) {
			ir<-1;
			rcnt1<-sdata[,cnt]==levels(as.factor(sdata[,cnt]))[l];
			rcnt2<-sdata[,cnt]==levels(as.factor(sdata[,cnt]))[k];
			swght1<-sum(sdata[rcnt1,wght]);
			swght2<-sum(sdata[rcnt2,wght]);
			for (i in 1:length(cmean)) {				
				mcnt1<-sum(sdata[rcnt1,wght]*sdata[rcnt1,cmean[i]])/swght1;
				mcnt2<-sum(sdata[rcnt2,wght]*sdata[rcnt2,cmean[i]])/swght2;
				mmeans[ir,ic]<-mcnt1 - mcnt2;
				secnt1<-0;
				secnt2<-0;
				for (j in 1:length(brr)) {
					sbrr1<-sum(sdata[rcnt1,brr[j]]);
					sbrr2<-sum(sdata[rcnt2,brr[j]]);
					mbcnt1<-sum(sdata[rcnt1,brr[j]]*sdata[rcnt1,cmean[i]])/sbrr1;
					mbcnt2<-sum(sdata[rcnt2,brr[j]]*sdata[rcnt2,cmean[i]])/sbrr2;
					secnt1<-secnt1 + (mcnt1 - mbcnt1)^2;
					secnt2<-secnt2 + (mcnt2 - mbcnt2)^2;
				}		
				secnt1<-sqrt((secnt1 * 4) / length(brr));
				secnt2<-sqrt((secnt2 * 4) / length(brr));
				mmeans[ir+1,ic]<-sqrt(secnt1^2 + secnt2^2);
				ir<-ir+2;
			}
			ic<-ic + 1;
		}
	}
	return(mmeans);
}

##################################################################################################
# Calcula la diferencia de medias por país entre los datos de un conjunto de columnas en función de una conjunto 
# de factores, y su error estandar
# sdata: dataframe con los datos
# cmean: vector con los índices o los nombres de las columnas con los datos
# cnt: índice o nombre de la columna con el país
# cfact: vector con los índices o los nombres de las columnas con los factores
# wght: índice o nombre de la columna con el peso de las filas
# brr: vector con los índices o los nombres de las columnas con los pesos replicados
# Devuelve una lista de matrices, con una matriz por país, con una columna para cada diferencia de medias entre dos conjuntos de valores
# de cada una de las columnas de cmean, agrupada por una de las combinaciones diferentes de dos de los 
# valores de cada factor de cfact
# Para cada columna en cmean se crean dos filas, la primera para las diferencias y la segunda para 
# los errores típicos
##################################################################################################
# Calculate the difference of means by country between the data in a set of columns in function of a set of
# factors, and their standard errors
# sdata: dataframe with the data
# cmean: vector with the indices or names of the columns with the data
# cnt: index or name of the column with the country
# cfact: vector with the indices or names of the columns with the factors
# wght: index or name of the column with the row weight
# brr: vector with the indices or names of the columns with the replicate weights
# Returns a list of matrixes, with a matrix by country, with one column for each difference of means between two sets of values of each one
# of the columns in cmean, grouped by each different combination of two values of each factor in cfact
# There are two rows for each column in cmean, the first for the differences and the second for the
# standard errors

wght_meandiffbycnt2<-function(sdata,cmean,cnt,cfact,wght,brr) {
	lcntrs<-vector('list',length(levels(as.factor(sdata[,cnt]))));
	for (p in 1:length(levels(as.factor(sdata[,cnt])))) {
		names(lcntrs)[p]<-levels(as.factor(sdata[,cnt]))[p];
	}
	nc<-0;
	for (i in 1:length(cfact)) {		
		for (j in 1:(length(levels(as.factor(sdata[,cfact[i]])))-1)) {
			for(k in (j+1):length(levels(as.factor(sdata[,cfact[i]])))) {
				nc <- nc + 1;
			}
		}
	}
	cn<-c();
	for (i in 1:length(cfact)) {
		for (j in 1:(length(levels(as.factor(sdata[,cfact[i]])))-1)) {
			for(k in (j+1):length(levels(as.factor(sdata[,cfact[i]])))) {
				cn<-c(cn, paste(names(sdata)[cfact[i]],levels(as.factor(sdata[,cfact[i]]))[j],levels(as.factor(sdata[,cfact[i]]))[k],sep="-"));
			}
		}
	}
	rn<-c();
	for (i in 1:length(cmean)) {
		rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF",sep="-"));
		rn<-c(rn, paste(names(sdata)[cmean[i]],"MEANDIFF","SE",sep="-"));
	}	
	for (p in 1:length(levels(as.factor(sdata[,cnt])))) {
		mmeans<-matrix(ncol=nc,nrow=2 * length(cmean));
		mmeans[,]<-0;
		colnames(mmeans)<-cn;
		rownames(mmeans)<-rn;	
		ic<-1;
		for(f in 1:length(cfact)) {
			for (l in 1:(length(levels(as.factor(sdata[,cfact[f]])))-1)) {
				for(k in (l+1):length(levels(as.factor(sdata[,cfact[f]])))) {
					ir<-1;
					rfact1<-(sdata[,cfact[f]]==levels(as.factor(sdata[,cfact[f]]))[l])&(sdata[,cnt]==levels(as.factor(sdata[,cnt]))[p]);
					rfact2<-(sdata[,cfact[f]]==levels(as.factor(sdata[,cfact[f]]))[k])&(sdata[,cnt]==levels(as.factor(sdata[,cnt]))[p]);
					swght1<-sum(sdata[rfact1,wght]);
					swght2<-sum(sdata[rfact2,wght]);
					for (i in 1:length(cmean)) {				
						mmeans[ir,ic]<-(sum(sdata[rfact1,wght]*sdata[rfact1,cmean[i]])/swght1) - (sum(sdata[rfact2,wght]*sdata[rfact2,cmean[i]])/swght2);
						for (j in 1:length(brr)) {
							sbrr1<-sum(sdata[rfact1,brr[j]]);
							sbrr2<-sum(sdata[rfact2,brr[j]]);
							mmbrj<-(sum(sdata[rfact1,brr[j]]*sdata[rfact1,cmean[i]])/sbrr1) - (sum(sdata[rfact2,brr[j]]*sdata[rfact2,cmean[i]])/sbrr2);
							mmeans[ir+1,ic]<-mmeans[ir+1,ic] + (mmbrj - mmeans[ir,ic])^2;
						}		
						mmeans[ir+1,ic]<-sqrt((mmeans[ir+1,ic] * 4) / length(brr));
						ir<-ir+2;
					}
					ic<-ic + 1;
				}
			}
		}
		lcntrs[[p]]<-mmeans;
	}
	return(lcntrs);
}