
# Estimation CCS SIxSI


################ FUNCTIONS
################################################################################################
#DESCRIPION: estimates a population total
#USAGE: EstTccsSISI(ECH,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECH	matrix of the interest variable, of size [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a value for the estimated total

EstTccsSISI <- function(ECH,n_m,n_d,N_m,N_d)
	{w=(N_m/n_m)*(N_d/n_d); Tpi=sum(w*c(ECH)) ; return(Tpi) }

################################################################################################
#DESCRIPION: estimates the variance of an estimated total for a CCS SIxSI
#USAGE: EstVARTccsSISI(ECH,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECH	matrix of the interest variable, of size  [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a list of three values (the unbiased estimator Vccs and its approximations V2d and V2d,a)

EstVARTccsSISI <- function(ECH,n_m,n_d,N_m,N_d)
	{
	w_m=N_m/n_m; w_d=N_d/n_d
	VarM=var(apply(ECH*w_d,2,sum));  VM=(N_m^2)*(1/n_m - 1/N_m)*VarM	
	VarD=var(apply(ECH*w_m,1,sum));  VD=(N_d^2)*(1/n_d - 1/N_d)*VarD
	d=apply(ECH,1,mean); e=apply(ECH,2,mean); f=mean(ECH)
	X = t(t(ECH - d) - e) + f; X2=X*X
	Vco=(1/(n_d - 1))*(1/(n_m - 1))*sum(X2)
	VCORR=(N_d^2)*(1/n_d - 1/N_d)*(N_m^2)*(1/n_m - 1/N_m)*Vco	
	VCCS=VM+VD-VCORR ; return(VCCS)
	}

################################################################################################
#DESCRIPION: estimates a population ratio Y/X
#USAGE: EstRccsSISI(ECHY,ECHX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECHY	matrix of the numerator variable, of size [n_d,n_m]
#		ECHX	matrix of the denominator variable, of size [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a value for the estimated ratio

EstRccsSISI <- function(ECHY,ECHX,n_m,n_d,N_m,N_d)
	{Rpi=EstTccsSISI(ECHY,n_m,n_d,N_m,N_d)/EstTccsSISI(ECHX,n_m,n_d,N_m,N_d) ; return(Rpi) }

################################################################################################
#DESCRIPION: estimates the linearized variable of a ratio Y/X
#USAGE: LinearizedR(ECHY,ECHX,n_m,n_d,N_m,N_d) 
#ARGUMENTS: 
#		ECHY	matrix of the numerator variable, of size [n_d,n_m]
#		ECHX	matrix of the denominator variable, of size [n_d,n_m]
#		n_m 	size of the sample drawn in the population U_m
#		n_d 	size of the sample drawn in the population U_d
#		N_m	size of the population U_m
#		N_d	size of the population U_d
#VALUE: the function returns a matrix of length [n_d,n_m]

LinearizedR <- function(ECHY,ECHX,n_m,n_d,N_m,N_d)
	{
	rpi=EstRccsSISI(ECHY,ECHX,n_m,n_d,N_m,N_d)
	txpi=EstTccsSISI(ECHX,n_m,n_d,N_m,N_d)
	upi = (ECHY - ECHX*rpi) / txpi ; return(upi)
	}




########## READING DATA

tableR=read.csv2(".../Data2stCCS.csv")
#or
tableR=read.csv2(file.choose())


echCCS=tableR[tableR$Dummy_CCS==1 ,]  # CCS sample
attach(echCCS)

n_m=320; n_d=25; N_m=544; N_d=365; N=N_m*N_d

echXCCS=matrix(Xik,nrow=n_d)
EstTccsSISI(ECH=echXCCS,n_m,n_d,N_m,N_d)
EstVARTccsSISI(ECH=echXCCS,n_m,n_d,N_m,N_d)


echYCCS=matrix(Yik,nrow=n_d)
EstRccsSISI(ECHY=echYCCS,ECHX=echXCCS,n_m,n_d,N_m,N_d)
LinR=LinearizedR(ECHY=echYCCS,ECHX=echXCCS,n_m,n_d,N_m,N_d)
EstVARTccsSISI(ECH=LinR,n_m,n_d,N_m,N_d)




