################################################################################
# Calcula la media y la desviación típica usando un conjunto de valores plausibles
# sdata: dataframe con los datos
# pv: vector con los nombres o índices de las columnas con los valores plausibles
# wght: índice o nombre de la columna con el peso del alumno
# brr: vector con los índices o los nombres de las columnas con los pesos replicados
# Devuelve un vector con la media, la desviación típica y sus errores estandar
################################################################################
# Computes the mean and standard deviation using a set of plausible values
# sdata: dataframe with the data
# pv: vector with the indices or column names of the plausible values
# wght: index or name of the column with the student weight
# brr: vector with the indices or column names with the replicate weights
# Returns a vector with the mean, standard deviation and their standard errors
################################################################################

wght_meansd_pv<-function(sdata,pv,wght,brr) {
	mmeans<-c(0, 0, 0, 0);
	mmeanspv<-rep(0,length(pv));
	stdspv<-rep(0,length(pv));
	mmeansbr<-rep(0,length(pv));
	stdsbr<-rep(0,length(pv));
	names(mmeans)<-c("MEAN","SE-MEAN","STDEV","SE-STDEV");
	swght<-sum(sdata[,wght]);
	for (i in 1:length(pv)) {
		mmeanspv[i]<-sum(sdata[,wght]*sdata[,pv[i]])/swght;
		stdspv[i]<-sqrt((sum(sdata[,wght]*(sdata[,pv[i]]^2))/swght)-mmeanspv[i]^2);
		for (j in 1:length(brr)) {
			sbrr<-sum(sdata[,brr[j]]);
			mbrrj<-sum(sdata[,brr[j]]*sdata[,pv[i]])/sbrr;
			mmeansbr[i]<-mmeansbr[i] + (mbrrj - mmeanspv[i])^2;
			stdsbr[i]<-stdsbr[i] + (sqrt((sum(sdata[,brr[j]]*(sdata[,pv[i]]^2))/sbrr)-mbrrj^2) - stdspv[i])^2;
		}		
	}
	mmeans[1]<-sum(mmeanspv) / length(pv);
	mmeans[2]<-sum((mmeansbr * 4) / length(brr)) / length(pv);
	mmeans[3]<-sum(stdspv) / length(pv);
	mmeans[4]<-sum((stdsbr * 4) / length(brr)) / length(pv);
	ivar <- c(0,0);
	for (i in 1:length(pv)) {
		ivar[1] <- ivar[1] + (mmeanspv[i] - mmeans[1])^2;
		ivar[2] <- ivar[2] + (stdspv[i] - mmeans[3])^2;
	}
	ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
	mmeans[2]<-sqrt(mmeans[2] + ivar[1]);
	mmeans[4]<-sqrt(mmeans[4] + ivar[2]);
	return(mmeans);
}

################################################################################
# Calcula la media y la desviación típica usando un conjunto de valores plausibles
# agrupando los resultados en función de los niveles de una serie de factores
# sdata: dataframe con los datos
# pv: vector con los nombres o índices de las columnas con los valores plausibles
# cfact: vector con los índices o los nombres de las columnas con los factores
# wght: índice o nombre de la columna con el peso del alumno
# brr: vector con los índices o los nombres de las columnas con los pesos replicados
# Devuelve una matriz con la media, la desviación típica y sus errores estandar
# cuyas columnas son los valores para cada nivel de factor
################################################################################
# Computes the mean and standard deviation using a set of plausible values
# grouping by the levels of a series of fators
# sdata: dataframe with the data
# pv: vector with the indices or column names of the plausible values
# cfact: vector with the indices or column names with the factors
# wght: index or name of the column with the student weight
# brr: vector with the indices or column names with the replicate weights
# Returns a matrix with the mean, standard deviation and their standard errors
# whose columns are the values for each level of the factors
################################################################################

wght_meansdfact_pv<-function(sdata,pv,cfact,wght,brr) {
	nc<-0;
	for (i in 1:length(cfact)) {
		nc <- nc + length(levels(as.factor(sdata[,cfact[i]])));
	}
	mmeans<-matrix(ncol=nc,nrow=4);
	mmeans[,]<-0;
	cn<-c();
	for (i in 1:length(cfact)) {
		for (j in 1:length(levels(as.factor(sdata[,cfact[i]])))) {
			cn<-c(cn, paste(names(sdata)[cfact[i]],levels(as.factor(sdata[,cfact[i]]))[j],sep="-"));
		}
	}
	colnames(mmeans)<-cn;
	rownames(mmeans)<-c("MEAN","SE-MEAN","STDEV","SE-STDEV");
	ic<-1;
	for(f in 1:length(cfact)) {
		for (l in 1:length(levels(as.factor(sdata[,cfact[f]])))) {
			rfact<-sdata[,cfact[f]]==levels(as.factor(sdata[,cfact[f]]))[l];
			swght<-sum(sdata[rfact,wght]);
			mmeanspv<-rep(0,length(pv));
			stdspv<-rep(0,length(pv));
			mmeansbr<-rep(0,length(pv));
			stdsbr<-rep(0,length(pv));
			for (i in 1:length(pv)) {			
				mmeanspv[i]<-sum(sdata[rfact,wght]*sdata[rfact,pv[i]])/swght;
				stdspv[i]<-sqrt((sum(sdata[rfact,wght]*(sdata[rfact,pv[i]]^2))/swght)-mmeanspv[i]^2);
				for (j in 1:length(brr)) {
					sbrr<-sum(sdata[rfact,brr[j]]);
					mbrrj<-sum(sdata[rfact,brr[j]]*sdata[rfact,pv[i]])/sbrr;
					mmeansbr[i]<-mmeansbr[i] + (mbrrj - mmeanspv[i])^2;		
					stdsbr[i]<-stdsbr[i] + (sqrt((sum(sdata[rfact,brr[j]]*(sdata[rfact,pv[i]]^2))/sbrr)-mbrrj^2) - stdspv[i])^2;					
				}		
			}
			mmeans[1, ic]<- sum(mmeanspv) / length(pv);
			mmeans[2, ic]<-sum((mmeansbr * 4) / length(brr)) / length(pv);
			mmeans[3, ic]<- sum(stdspv) / length(pv);
			mmeans[4, ic]<-sum((stdsbr * 4) / length(brr)) / length(pv);
			ivar <- c(sum((mmeanspv - mmeans[1, ic])^2),sum((stdspv - mmeans[3, ic])^2));
			ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
			mmeans[2, ic]<-sqrt(mmeans[2, ic] + ivar[1]);
			mmeans[4, ic]<-sqrt(mmeans[4, ic] + ivar[2]);
			ic<-ic + 1;
		}
	}
	return(mmeans);
}

################################################################################
# Calcula la diferencia de medias entre paises usando valores plausibles
# sdata: dataframe con los datos
# pv: vector con los nombres o índices de las columnas con los valores plausibles
# cnt: índice o nombre de la columna con el país
# wght: índice o nombre de la columna con el peso del alumno
# brr: vector con los índices o los nombres de las columnas con los pesos replicados
# Devuelve una matriz con la media, y su error estandar
# cuyas columnas son las diferencias entre cada par de países en la muestra
################################################################################
# Computes the mean difference between countries using a set of plausible values
# sdata: dataframe with the data
# pv: vector with the indices or column names of the plausible values
# cnt: index or name of the column with the country
# wght: index or name of the column with the student weight
# brr: vector with the indices or column names with the replicate weights
# Returns a matrix with the mean and their standard error
# whose columns are the differences for each pair of countries
################################################################################

wght_meandiffcnt_pv<-function(sdata,pv,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);
	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(levels(as.factor(sdata[,cnt]))[j],levels(as.factor(sdata[,cnt]))[k],sep="-"));
		}
	}
	colnames(mmeans)<-cn;
	rn<-c("MEANDIFF", "SE");
	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])))) {
			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]);
			mmeanspv<-rep(0,length(pv));
			mmcnt1<-rep(0,length(pv));
			mmcnt2<-rep(0,length(pv));
			mmeansbr1<-rep(0,length(pv));
			mmeansbr2<-rep(0,length(pv));
			for (i in 1:length(pv)) {				
				mmcnt1<-sum(sdata[rcnt1,wght]*sdata[rcnt1,pv[i]])/swght1;
				mmcnt2<-sum(sdata[rcnt2,wght]*sdata[rcnt2,pv[i]])/swght2;
				mmeanspv[i]<- mmcnt1 - mmcnt2;
				for (j in 1:length(brr)) {
					sbrr1<-sum(sdata[rcnt1,brr[j]]);
					sbrr2<-sum(sdata[rcnt2,brr[j]]);
					mmbrj1<-sum(sdata[rcnt1,brr[j]]*sdata[rcnt1,pv[i]])/sbrr1;
					mmbrj2<-sum(sdata[rcnt2,brr[j]]*sdata[rcnt2,pv[i]])/sbrr2;
					mmeansbr1[i]<-mmeansbr1[i] + (mmbrj1 - mmcnt1)^2;
					mmeansbr2[i]<-mmeansbr2[i] + (mmbrj2 - mmcnt2)^2;
				}		
			}
			mmeans[1,ic]<-sum(mmeanspv) / length(pv);
			mmeansbr1<-sum((mmeansbr1 * 4) / length(brr)) / length(pv);
			mmeansbr2<-sum((mmeansbr2 * 4) / length(brr)) / length(pv);
			mmeans[2,ic]<-sqrt(mmeansbr1^2 + mmeansbr2^2);
			ivar <- 0;
			for (i in 1:length(pv)) {
				ivar <- ivar + (mmeanspv[i] - mmeans[1,ic])^2;
			}
			ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
			mmeans[2,ic]<-sqrt(mmeans[2,ic] + ivar);
			ic<-ic + 1;
		}
	}
	return(mmeans);
}

################################################################################
# Calcula la diferencia de medias de valores plausibles agrupados por
# niveles de una serie de factores dentro de cada país y entre paises
# sdata: dataframe con los datos
# pv: vector con los nombres o índices de las columnas con los valores plausibles
# 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 del alumno
# brr: vector con los índices o los nombres de las columnas con los pesos replicados
# Devuelve una lista de matrices con la diferencia de media y su error estandar
# dentro de cada país, mas un elemento con un array de tres dimensiones con
# las diferencias entre cada par de paises
# las columnas son las diferencias entre cada combinación de dos niveles de cada factor
################################################################################
# Computes the mean difference of plausible values grouped by
# levels of a series of factors between each country and between countries
# sdata: dataframe with the data
# pv: vector with the indices or column names of the plausible values
# cnt: index or name of the column with the country
# cfact: vector with the indices or column names with the factors
# wght: index or name of the column with the student weight
# brr: vector with the indices or column names with the replicate weights
# Returns a list of matrixes with the mean difference and their standard error
# into each country, plus an element with a three dimensional array with 
# the differences between each pair of countries
# the columns are the differences between each combination of two levels of each factor
################################################################################

wght_meandifffactcnt_pv<-function(sdata,pv,cnt,cfact,wght,brr) {
	lcntrs<-vector('list',1 + 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];
	}
	names(lcntrs)[1 + length(levels(as.factor(sdata[,cnt])))]<-"BTWNCNT";
	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("MEANDIFF", "SE");
	for (p in 1:length(levels(as.factor(sdata[,cnt])))) {
		mmeans<-matrix(ncol=nc,nrow=2);
		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]])))) {
					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]);
					mmeanspv<-rep(0,length(pv));
					mmeansbr<-rep(0,length(pv));
					for (i in 1:length(pv)) {				
						mmeanspv[i]<-(sum(sdata[rfact1,wght]*sdata[rfact1,pv[i]])/swght1) - (sum(sdata[rfact2,wght]*sdata[rfact2,pv[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,pv[i]])/sbrr1) - (sum(sdata[rfact2,brr[j]]*sdata[rfact2,pv[i]])/sbrr2);
							mmeansbr[i]<-mmeansbr[i] + (mmbrj - mmeanspv[i])^2;
						}		
					}
					mmeans[1,ic]<-sum(mmeanspv) / length(pv);
					mmeans[2,ic]<-sum((mmeansbr * 4) / length(brr)) / length(pv);
					ivar <- 0;
					for (i in 1:length(pv)) {
						ivar <- ivar + (mmeanspv[i] - mmeans[1,ic])^2;
					}
					ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
					mmeans[2,ic]<-sqrt(mmeans[2,ic] + ivar);
					ic<-ic + 1;
				}
			}
		}
		lcntrs[[p]]<-mmeans;
	}
	pn<-c();	
	for (p in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
		for (p2 in (p + 1):length(levels(as.factor(sdata[,cnt])))) {
			pn<-c(pn, paste(levels(as.factor(sdata[,cnt]))[p],levels(as.factor(sdata[,cnt]))[p2],sep="-"));	
		}
	}
	mbtwmeans<-array(0, c(length(rn), length(cn), length(pn)));
	nm <- vector('list',3);
	nm[[1]]<-rn;
	nm[[2]]<-cn;
	nm[[3]]<-pn;
	dimnames(mbtwmeans)<-nm;
	pc<-1;
	for (p in 1:(length(levels(as.factor(sdata[,cnt])))-1)) {
		for (p2 in (p + 1):length(levels(as.factor(sdata[,cnt])))) {
			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]])))) {
						mbtwmeans[1,ic,pc]<-lcntrs[[p]][1,ic]-lcntrs[[p2]][1,ic];
						mbtwmeans[2,ic,pc]<-sqrt((lcntrs[[p]][2,ic]^2)+(lcntrs[[p2]][2,ic]^2));
						ic<-ic + 1;
					}
				}
			}
			pc<-pc+1;			
		}
	}
	lcntrs[[1 + length(levels(as.factor(sdata[,cnt])))]]<-mbtwmeans;
	return(lcntrs);
}

################################################################################
# Realiza una regresión lineal con valores plausibles como variables dependientes
# sdata: dataframe con los datos
# frml: texto con la parte derecha de la fórmula
# pv: vector con los nombres o índices de las columnas con los valores plausibles
# wght: índice o nombre de la columna con el peso del alumno
# brr: vector con los índices o los nombres de las columnas con los pesos replicados
# Devuelve una lista con el modelo lineal para cada valor plausible
# y el modelo lineal final correspondiente al conjunto de valores plausibles
# junto con los errores estandar de los coeficientes
################################################################################
# Performs a linear regression with plausible values as dependent variables
# sdata: dataframe with the data
# frml: text with the right side of the formula
# pv: vector with the indices or column names of the plausible values
# wght: index or name of the column with the student weight
# brr: vector with the indices or column names with the replicate weights
# Returns a list with the linear model for each plausible value
# and the final linear model corresponding of the set of plausible values,
# along with the standard error of the coefficients
################################################################################

wght_lmpv<-function(sdata,frml,pv,wght,brr) {
	listlm <- vector('list', 2 + length(pv));
      listbr <- vector('list', length(pv));	
	for (i in 1:length(pv)) {				
		if (is.numeric(pv[i])) {
			names(listlm)[i] <- colnames(sdata)[pv[i]];
			frmlpv <- as.formula(paste(colnames(sdata)[pv[i]],frml,sep="~"));			
		}
		else {
			names(listlm)[i]<-pv[i];
			frmlpv <- as.formula(paste(pv[i],frml,sep="~"));
		}
		listlm[[i]] <- lm(frmlpv, data=sdata, weights=sdata[,wght]);
		listbr[[i]] <- rep(0,2 + length(listlm[[i]]$coefficients));
		for (j in 1:length(brr)) {
			lmb <- lm(frmlpv, data=sdata, weights=sdata[,brr[j]]);
			listbr[[i]]<-listbr[[i]] + c((listlm[[i]]$coefficients - lmb$coefficients)^2,(summary(listlm[[i]])$r.squared-summary(lmb)$r.squared)^2,(summary(listlm[[i]])$adj.r.squared-summary(lmb)$adj.r.squared)^2);
		}
		listbr[[i]] <- (listbr[[i]] * 4) / length(brr);
	}
	cf <- c(listlm[[1]]$coefficients,0,0);
	names(cf)[length(cf)-1]<-"R2";
	names(cf)[length(cf)]<-"ADJ.R2";
	for (i in 1:length(cf)) {
		cf[i] <- 0;
	}
	for (i in 1:length(pv)) {
		cf<-(cf + c(listlm[[i]]$coefficients,summary(listlm[[i]])$r.squared,summary(listlm[[i]])$adj.r.squared));
	}
	names(listlm)[1 + length(pv)]<-"RESULT";
	listlm[[1 + length(pv)]]<- cf / length(pv); 
	names(listlm)[2 + length(pv)]<-"SE";
	listlm[[2 + length(pv)]] <- rep(0, length(cf));
	names(listlm[[2 + length(pv)]])<-names(cf);
	for (i in 1:length(pv)) {
		listlm[[2 + length(pv)]] <- listlm[[2 + length(pv)]] + listbr[[i]];
	}
	ivar <- rep(0,length(cf));
	for (i in 1:length(pv)) {
		ivar <- ivar + c((listlm[[i]]$coefficients - listlm[[1 + length(pv)]][1:(length(cf)-2)])^2,(summary(listlm[[i]])$r.squared - listlm[[1 + length(pv)]][length(cf)-1])^2,(summary(listlm[[i]])$adj.r.squared - listlm[[1 + length(pv)]][length(cf)])^2);
	}
	ivar = (1 + (1 / length(pv))) * (ivar / (length(pv) - 1));
	listlm[[2 + length(pv)]] <- sqrt((listlm[[2 + length(pv)]] / length(pv)) + ivar);
	return(listlm);
}