/*
#							paired t-test for agreement coefficients.r
#						 		  (July, 2021)
#Description: This script file contains a series of R functions for testing the difference between 2 agreement
#		 coefficients for statistical significance. The approach used here is based on the linearization method 
#    proposed Gwet (2016).A function is proposed for each agreement coefficient (Gwet's AC2, Fleiss' Kappa, 
#    Conger's Kappa, Krippendorff's Alpha, and Brennan-Prediger Coefficient.
#
#    Each function has only 2 required parameters, which are the 2 data frames containing the 2 sets of ratings for
#    the 2 groups of raters under comparison. Each input data file is in the form of nxr matrix or data frame showing 
#    the actual ratings each rater (column) assigned to each subject (in row). A typical table entry (i,g) represents 
#    the rating associated with subject i and rater g. 
#
#Author: Kilem L. Gwet, Ph.D.  (gwet@agreestat.com)
#
#
#     EXAMPLES OF SIMPLE CALLS OF THE 5 FUNCTIONS
#     ===========================================
# > ttest_fleiss(g1.ratings,g2.ratings)
# > ttest_ac2(g1.ratings,g2.ratings)
# > ttest_alpha(g1.ratings,g2.ratings)
# > ttest_conger(g1.ratings,g2.ratings)
# > ttest_bp(g1.ratings,g2.ratings)
*/

/*
#===========================================
#ttest_fleiss: This function performs the paired t-test for testing the difference  
#              between two correlated Fleiss kappa coefficients for statistical significance_ It
#              implements the linearization method of Gwet (2016).
				(Author: K. Gwet, PhD - e-mail: gwet@agreestat.com)
#-------------
# This function takes 2 required parameters, which are the 2 groups of raters being compared:
# "g1_ratings" and "g2_ratings_" Both datasets must have the exact same number of rows, and each
# column represents one rater and contains its ratings (numeric or alphabetic). 
# The user must exclude all subjects that are not rated by any rater_
# Bibliography:
# ------------
# Fleiss, J. L. (1971). Measuring nominal scale agreement among many raters_
#                        Psychological Bulletin, 76, 378-382_ 
# Gwet, K. L. (2016). "Testing the Difference of Correlated Agreement Coefficients for Statistical Significance,
#                      Educational and Psychological Measurement, Vol 76(4) 609-637
#==========================================================================
*/
start ttest_fleiss(g1_ratings,g2_ratings,weights="unweighted",conflev=0.95,Npop=10**7); 
	n2 = nrow(g2_ratings);
	n1 = nrow(g1_ratings);
	if n2=n1 then do;
		coeff2_i=fleiss_linear_i(g2_ratings,weights,conflev,Npop);
		coeff1_i=fleiss_linear_i(g1_ratings,weights,conflev,Npop);
		di = coeff2_i-coeff1_i;
		fleiss_coeff1 = mean(coeff1_i);
		fleiss_coeff2 = mean(coeff2_i);
		std_err = sqrt(var(di)/n1);
		t_stat = (fleiss_coeff2-fleiss_coeff1)/std_err;
		p_value = 2*(1-cdf('T',abs(t_stat),n1-1));
	end;else do;fleiss_coeff1=-999;fleiss_coeff2=-999;std_err=-999;t_stat=-999;p_value=-999;end;
	diff = fleiss_coeff2-fleiss_coeff1;
	mstats = t(fleiss_coeff1||fleiss_coeff2||diff||std_err||t_stat||p_value);
	mlabels={"Fleiss Kappa1","Fleiss Kappa2","Diff","SE(diff)","t-stat","p-value"};
	return([mstats,mlabels]);
finish;
/*
#===========================================
#ttest_ac2: This function performs the paired t-test for testing the difference  
#              between two correlated Gwet's AC2 coefficients for statistical significance_ It
#              implements the linearization method of Gwet (2016).
				(Author: K. Gwet, PhD - e-mail: gwet@agreestat.com)
#-------------
# This function takes 2 required parameters, which are the 2 groups of raters being compared:
# "g1_ratings" and "g2_ratings_" Both datasets must have the exact same number of rows, and each
# column represents one rater and contains its ratings (numeric or alphabetic). 
# The user must exclude all subjects that are not rated by any rater_
# Bibliography:
# ------------
# Gwet, K. L. (2008). Computing inter-rater reliability and its variance in the presence of high agreement_
#                     British Journal of Mathematical and Statistical Psychology, 61, 29-48_
# Gwet, K. L. (2016). "Testing the Difference of Correlated Agreement Coefficients for Statistical Significance,
#                      Educational and Psychological Measurement, Vol 76(4) 609-637_
#===========================================
*/
start ttest_ac2(g1_ratings,g2_ratings,weights="unweighted",conflev=0.95,Npop=10**7); 
	n2 = nrow(g2_ratings);
	n1 = nrow(g1_ratings);
	if n2=n1 then do;
		coeff2_i=ac2_linear_i(g2_ratings,weights,conflev,Npop);
		coeff1_i=ac2_linear_i(g1_ratings,weights,conflev,Npop);
		di = coeff2_i-coeff1_i;
		ac2_coeff1 = mean(coeff1_i);
		ac2_coeff2 = mean(coeff2_i);
		std_err = sqrt(var(di)/n1);
		t_stat = (ac2_coeff2-ac2_coeff1)/std_err;
		p_value = 2*(1-cdf('T',abs(t_stat),n1-1));
	end;else do;ac2_coeff1=-999;ac2_coeff2=-999;std_err=-999;t_stat=-999;p_value=-999;end;
	diff = ac2_coeff2-ac2_coeff1;
	mstats = t(ac2_coeff1||ac2_coeff2||diff||std_err||t_stat||p_value);
	mlabels={"AC1/2 Coeff1","AC1/2 Coeff2","Diff","SE(diff)","t-stat","p-value"};
	return([mstats,mlabels]);
finish;
/*
#===========================================
#ttest_alpha: This function performs the paired t-test for testing the difference  
#              between two correlated Krippendorff's alpha coefficients for statistical significance_ It
#              implements the linearization method of Gwet (2016).
				(Author: K. Gwet, PhD - e-mail: gwet@agreestat.com)
#-------------
# This function takes 2 required parameters, which are the 2 groups of raters being compared:
# "g1_ratings" and "g2_ratings_" Both datasets must have the exact same number of rows, and each
# column represents one rater and contains its ratings (numeric or alphabetic). 
# The user must exclude all subjects that are not rated by any rater_
# Bibliography:
# ------------
# Krippendorff, K. (1970). Estimating the reliability, systematic error, and random error of interval data_
#                          Educational and Psychological Measurement, 30, 61-70_
# Gwet, K. L. (2016). "Testing the Difference of Correlated Agreement Coefficients for Statistical Significance,
#                      Educational and Psychological Measurement, Vol 76(4) 609-637_
#===========================================
*/
start ttest_alpha(g1_ratings,g2_ratings,weights="unweighted",conflev=0.95,Npop=10**7); 
	n2 = nrow(g2_ratings);
	n1 = nrow(g1_ratings);
	if n2=n1 then do;
		coeff2_i=alpha_linear_i(g2_ratings,weights,conflev,Npop);
		coeff1_i=alpha_linear_i(g1_ratings,weights,conflev,Npop);
		di = coeff2_i-coeff1_i;
		alpha_coeff1 = mean(coeff1_i);
		alpha_coeff2 = mean(coeff2_i);
		std_err = sqrt(var(di)/n1);
		t_stat = (alpha_coeff2-alpha_coeff1)/std_err;
		p_value = 2*(1-cdf('T',abs(t_stat),n1-1));
	end;else do;alpha_coeff1=-999;alpha_coeff2=-999;std_err=-999;t_stat=-999;p_value=-999;end;
	diff = alpha_coeff2-alpha_coeff1;
	mstats = t(alpha_coeff1||alpha_coeff2||diff||std_err||t_stat||p_value);
	mlabels={"Krippendorff Alpha1","Krippendorff Alpha2","Diff","SE(diff)","t-stat","p-value"};
	return([mstats,mlabels]);
finish;
/*
#===========================================
#ttest_conger: This function performs the paired t-test for testing the difference  
#              between two correlated Conger's Kappa coefficients for statistical significance_ It
#              implements the linearization method of Gwet (2016).
				(Author: K. Gwet, PhD - e-mail: gwet@agreestat.com)
#-------------
# This function takes 2 required parameters, which are the 2 groups of raters being compared:
# "g1_ratings" and "g2_ratings_" Both datasets must have the exact same number of rows, and each
# column represents one rater and contains its ratings (numeric or alphabetic). 
# The user must exclude all subjects that are not rated by any rater_
# Bibliography:
# ------------
# Conger, A. J. (1980). Integration and generalization of kappas for multiple raters_ 
#                       Psychological Bulletin, 88, 322-328_ 
# Gwet, K. L. (2016). "Testing the Difference of Correlated Agreement Coefficients for Statistical Significance,
#                      Educational and Psychological Measurement, Vol 76(4) 609-637_
#===========================================
*/
start ttest_conger(g1_ratings,g2_ratings,weights="unweighted",conflev=0.95,Npop=10**7); 
	n2 = nrow(g2_ratings);
	n1 = nrow(g1_ratings);
	if n2=n1 then do;
		coeff2_i=conger_linear_i(g2_ratings,weights,conflev,Npop);
		coeff1_i=conger_linear_i(g1_ratings,weights,conflev,Npop);
		di = coeff2_i-coeff1_i;
		conger_coeff1 = mean(coeff1_i);
		conger_coeff2 = mean(coeff2_i);
		std_err = sqrt(var(di)/n1);
		t_stat = (conger_coeff2-conger_coeff1)/std_err;
		p_value = 2*(1-cdf('T',abs(t_stat),n1-1));
	end;else do;conger_coeff1=-999;conger_coeff2=-999;std_err=-999;t_stat=-999;p_value=-999;end;
	diff = conger_coeff2-conger_coeff1;
	mstats = t(conger_coeff1||conger_coeff2||diff||std_err||t_stat||p_value);
	mlabels={"Conger Kappa1","Conger Kappa2","Diff","SE(diff)","t-stat","p-value"};
	return([mstats,mlabels]);
finish;
/*
#=========================================================================================================
#ttest_bp: This function performs the paired t-test for testing the difference  
#              between two correlated Brennan-Prediger coefficients for statistical significance_ It
#              implements the linearization method of Gwet (2016).
				(Author: K. Gwet, PhD - e-mail: gwet@agreestat.com)
#-------------
# This function takes 2 required parameters, which are the 2 groups of raters being compared:
# "g1_ratings" and "g2_ratings_" Both datasets must have the exact same number of rows, and each
# column represents one rater and contains its ratings (numeric or alphabetic). 
# The user must exclude all subjects that are not rated by any rater_
# Bibliography:
# ------------
# Brennan, R_L. & Prediger, D. J. (1981). Coefficient kappa: some uses, misuses, and alternatives_
#                       Educational and Psychological Measurement, 41, 687-699_ 
# Gwet, K. L. (2016). "Testing the Difference of Correlated Agreement Coefficients for Statistical Significance,
#                      Educational and Psychological Measurement, Vol 76(4) 609-637
#==============================================================================================
*/
start ttest_bp(g1_ratings,g2_ratings,weights="unweighted",conflev=0.95,Npop=10**7); 
	n2 = nrow(g2_ratings);
	n1 = nrow(g1_ratings);
	if n2=n1 then do;
		coeff2_i=bp_linear_i(g2_ratings,weights,conflev,Npop);
		coeff1_i=bp_linear_i(g1_ratings,weights,conflev,Npop);
		di = coeff2_i-coeff1_i;
		bp_coeff1 = mean(coeff1_i);
		bp_coeff2 = mean(coeff2_i);
		std_err = sqrt(var(di)/n1);
		t_stat = (bp_coeff2-bp_coeff1)/std_err;
		p_value = 2*(1-cdf('T',abs(t_stat),n1-1));
	end;else do;bp_coeff1=-999;bp_coeff2=-999;std_err=-999;t_stat=-999;p_value=-999;end;
	diff = bp_coeff2-bp_coeff1;
	mstats = t(bp_coeff1||bp_coeff2||diff||std_err||t_stat||p_value);
	mlabels={"BP Coeff1","BP Coeff2","Diff","SE(diff)","t-stat","p-value"};
	return([mstats,mlabels]);
finish;




/*
#===========================================================================================
#ac2_linear_i: This function produces all n linear components associated with the n subjects, 
	and the linearized version of the ac2 coefficient.
#-------------
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings from n subjects and 
	r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Gwet, K. L. (2008). ``Computing inter-rater reliability and its variance in the presence of high
#		agreement." British Journal of Mathematical and Statistical Psychology, 61, 29-48.
#============================================================================================
*/
start ac2_linear_i(ratings,weights="unweighted",conflev=0.95,Npop=10**7);
	ratings_mat = ratings; 
	if type(ratings_mat)="C" then do;
		ratings_mat = strip(upcase(ratings_mat));
	end;
	n  =  nrow(ratings_mat); * number of subjects;
	r  =  ncol(ratings_mat); * number of raters;
	f  =  n/Npop; * finite population correction;
  
	/* creating a vector containing all categories used by the raters */
	categ  =  unique_vec(shape(ratings_mat,n*r,1));
	call sort(categ);
	q = nrow(categ);

	/* creating the weights matrix */
	weights_mat = weight_mat(weights,categ);

	/* creating the nxq agreement matrix representing the distribution of raters by subjects and category */  
	agree_mat = j(n,q,0);
	do k=1 to q;
		if type(ratings_mat)="N" then do;
			in_categ_k =(ratings_mat=categ[k]);
			agree_mat[,k] = in_categ_k*repeat(1,r); 
		end;else agree_mat[,k] = (ratings_mat=categ[k])*repeat(1,r);
	end;
	agree_mat_w = t(weights_mat*t(agree_mat));
  
	/* calculating gwet's ac1 coefficient */
	ri_vec = agree_mat*repeat(1,q);
	sum_q = (agree_mat#(agree_mat_w-1))*repeat(1,q);
	n2more = sum(ri_vec>=2);
	pa = sum(sum_q[loc(ri_vec>=2)]/((ri_vec#(ri_vec-1))[loc(ri_vec>=2)]))/n2more;
	pi_vec = t(t(repeat(1/n,n))*(agree_mat/(ri_vec*t(repeat(1,q)))));
	pe = sum(weights_mat) # sum(pi_vec#(1-pi_vec)) / (q*(q-1));
	gwet_ac1 = (pa-pe)/(1-pe);

	/* calculating variance, stderr & p-value of gwet's ac1 coefficient */
	den_ivec = ri_vec#(ri_vec-1);
	den_ivec = den_ivec - (den_ivec=0); /* this operation replaces each 0 value with -1 to make the next ratio calculation always possible*/
	pa_ivec = sum_q/den_ivec;

	pe_r2 = pe*(ri_vec>=2);
	ac1_ivec = (n/n2more)*(pa_ivec-pe_r2)/(1-pe);
	pe_ivec = (sum(weights_mat)/(q*(q-1))) * (agree_mat*(1-pi_vec))/ri_vec;
	ac1_ivec_x = ac1_ivec - 2*(1-gwet_ac1) * (pe_ivec-pe)/(1-pe);
	return(ac1_ivec_x);
finish;

/*
#=============================================================================================
#alpha_linear_i: This function produces all n linear components associated with the n subjects, 
				 and the linearized version of Krippendorff's alpha coefficient.
#-------------
#The algorithm used to compute krippendorff's alpha is very different from anything that was published 
	on this topic. Instead, it follows the equations presented by K. Gwet (2021, vol 1 p. 79).
The input data "ratings" is a nxr data frame of raw alphanumeric ratings from n subjects and r raters. 
Exclude all subjects that are not rated by any rater.
#Bibliography:
#Gwet, K. (2021). Handbook of Inter-Rater Reliability: The Definitive Guide to Measuring the Extent of 
	Agreement Among Multiple Raters, 5th Edition. Vol 1: Analysis of Categorical Ratings,
	AgreeStat Analytics, LLC;
#Krippendorff (1970). "Bivariate agreement coefficients for reliability of data." Sociological Methodology,2,139-150
#Krippendorff (1980). Content analysis: An introduction to its methodology (2nd ed.), New-bury Park, CA: Sage.
#======================================================================================
*/
start alpha_linear_i(ratings,weights="unweighted",conflev=0.95,Npop=10**7);
	ratings_mat = ratings; 
	if type(ratings_mat)="C" then do;
		ratings_mat = strip(upcase(ratings_mat));
	end;
	n  =  nrow(ratings_mat); * number of subjects;
	r  =  ncol(ratings_mat); * number of raters;
	f  =  n/Npop; * finite population correction;
  
	/* creating a vector containing all categories used by the raters */
	categ  =  unique_vec(shape(ratings_mat,n*r,1));
	call sort(categ);
	q = nrow(categ);

	/* creating the weights matrix */
	weights_mat = weight_mat(weights,categ);

	/* creating the nxq agreement matrix representing the distribution of raters by subjects and category */  
	agree_mat = j(n,q,0);
	do k=1 to q;
		if type(ratings_mat)="N" then do;
			in_categ_k =(ratings_mat=categ[k]);
			agree_mat[,k] = in_categ_k*repeat(1,r); 
		end;else agree_mat[,k] = (ratings_mat=categ[k])*repeat(1,r);
	end;
	agree_mat_w = t(weights_mat*t(agree_mat));
  
	/*calculating krippendorff's alpha coefficient*/
	ri_vec = agree_mat*repeat(1,q);
	agree_mat = agree_mat[loc(ri_vec>=2),];
	agree_mat_w = agree_mat_w[loc(ri_vec>=2),];
	ri_vec = ri_vec[loc(ri_vec>=2)];
	ri_mean = ri_vec[:];
	n = nrow(agree_mat);
	epsi = 1/ri_vec[+];
	sum_q = (agree_mat#(agree_mat_w-1))*repeat(1,q);
	pa = (1-epsi)* sum(sum_q/(ri_mean#(ri_vec-1)))/n + epsi;
	pi_vec = t(t(repeat(1/n,n))*(agree_mat/ri_mean));
	pe = sum(weights_mat # (pi_vec*t(pi_vec)));
	krippen_alpha = (pa-pe)/(1-pe);
  
	/* calculating variance, stderr & p-value of gwet's ac1 coefficient */
	den_ivec = ri_mean#(ri_vec-1);
	pa_ivec = sum_q/den_ivec;
	pa_v = pa_ivec[:];
	pa_ivec = (1-epsi)#(pa_ivec-pa_v#(ri_vec-ri_mean)/ri_mean) + epsi;

	krippen_ivec = (pa_ivec-pe)/(1-pe);
	pi_vec_wk_ = weights_mat*pi_vec;
	pi_vec_w_k = t(weights_mat)*pi_vec;
	pi_vec_w = (pi_vec_wk_ + pi_vec_w_k)/2;
	pe_ivec = (agree_mat*pi_vec_w)/ri_mean - sum(pi_vec) # (ri_vec-ri_mean)/ri_mean;
	krippen_ivec_x = krippen_ivec - 2*(1-krippen_alpha) * (pe_ivec-pe)/(1-pe);
	return(krippen_ivec_x);
finish;


/*
#===========================================================================================
#conger_linear_i: This function produces all n linear components associated with  
#		           the n subjects, and the linearized version of Conger's kappa coefficient.
#-------------
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings
#from n subjects and r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Conger, A. J. (1980), ``Integration and Generalization of Kappas for Multiple Raters,"
#		Psychological Bulletin, 88, 322-328.
#======================================================================================
*/
start conger_linear_i(ratings,weights="unweighted",conflev=0.95,Npop=10**7);
	ratings_mat = ratings; 
	if type(ratings_mat)="C" then do;
		ratings_mat = strip(upcase(ratings_mat));
	end;
	n  =  nrow(ratings_mat); * number of subjects;
	r  =  ncol(ratings_mat); * number of raters;
	f  =  n/Npop; * finite population correction;
  
	/* creating a vector containing all categories used by the raters */
	categ  =  unique_vec(shape(ratings_mat,n*r,1));
	call sort(categ);
	q = nrow(categ);

	/* creating the weights matrix */
	weights_mat = weight_mat(weights,categ);

	/* creating the nxq agreement matrix representing the distribution of raters by subjects and category */  
	agree_mat = j(n,q,0);
	do k=1 to q;
		if type(ratings_mat)="N" then do;
			in_categ_k =(ratings_mat=categ[k]);
			agree_mat[,k] = in_categ_k*repeat(1,r); 
		end;else agree_mat[,k] = (ratings_mat=categ[k])*repeat(1,r);
	end;
	agree_mat_w = t(weights_mat*t(agree_mat));
  
	/* creating the rxq rater-category matrix representing the distribution of subjects by rater and category */
	classif_mat = j(r,q,0);
	do k=1 to q;
		if type(ratings_mat)="N" then do;
			without_mis =(t(ratings_mat)=categ[k]);
			classif_mat[,k] = without_mis*repeat(1,n);
		end;else classif_mat[,k] = (t(ratings_mat)=categ[k])*repeat(1,n);
	end;
  
	/* calculating conger's kappa coefficient */  
	ri_vec = agree_mat*repeat(1,q);
	sum_q = (agree_mat#(agree_mat_w-1))*repeat(1,q);
	n2more = sum(ri_vec>=2);
	pa = sum(sum_q[loc(ri_vec>=2)]/((ri_vec#(ri_vec-1))[loc(ri_vec>=2)]))/n2more;
	ng_vec = classif_mat*repeat(1,q);
	pgk_mat = classif_mat/(ng_vec*t(repeat(1,q)));
	p_mean_k = (t(pgk_mat)*repeat(1,r))/r; 
	s2kl_mat = (t(pgk_mat)*pgk_mat - r * p_mean_k*t(p_mean_k))/(r-1);
	pe = sum(weights_mat # (p_mean_k*t(p_mean_k) -  s2kl_mat/r));
	conger_kappa = (pa-pe)/(1-pe);
  
	/* calculating variance, stderr & p-value of conger's kappa coefficient */ 
	bkl_mat = (weights_mat+t(weights_mat))/2;
	pe_ivec1 = r*(agree_mat*t(t(p_mean_k)*bkl_mat));
	pe_ivec2 = repeat(0,n);
	if type(ratings_mat)="N" then do;
		do l=1 to q;
		  delta_ig_mat =(ratings_mat=categ[l]);
		  pe_ivec2 = pe_ivec2 + delta_ig_mat*(pgk_mat*bkl_mat[,l]);
		end;
	end;else do;
		do k=1 to q;
		  delta_ig_mat = (ratings_mat=categ[k]);
		  pe_ivec2 = pe_ivec2 + delta_ig_mat*(pgk_mat*bkl_mat[,k]);
		end;
	end;
	pe_ivec = (pe_ivec1-pe_ivec2)/(r*(r-1)); 
	den_ivec = ri_vec#(ri_vec-1);
	den_ivec = den_ivec - (den_ivec=0); /* this operation replaces each 0 value with -1 to make the next ratio calculation always possible*/
	pa_ivec = sum_q/den_ivec;
	pe_r2 = pe*(ri_vec>=2);
	conger_ivec = (n/n2more)*(pa_ivec-pe_r2)/(1-pe); 
	conger_ivec_x = conger_ivec - 2*(1-conger_kappa) * (pe_ivec-pe)/(1-pe);
	return(conger_ivec_x);
finish;

/*
#===========================================================================================
#bp_linear_i: This function produces all n linear components associated with  
#		      the n subjects, and the linearized version of B-P's coefficient.
#-------------
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings
#from n subjects and r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Brennan, R.L., and Prediger, D. J. (1981). ``Coefficient Kappa: some uses, misuses, and alternatives."
#    Educational and Psychological Measurement, 41, 687-699.
#======================================================================================
*/
start bp_linear_i(ratings,weights="unweighted",conflev=0.95,Npop=10**7);
	ratings_mat = ratings; 
	if type(ratings_mat)="C" then do;
		ratings_mat = strip(upcase(ratings_mat));
	end;
	n  =  nrow(ratings_mat); * number of subjects;
	r  =  ncol(ratings_mat); * number of raters;
	f  =  n/Npop; * finite population correction;
  
	/* creating a vector containing all categories used by the raters */
	categ  =  unique_vec(shape(ratings_mat,n*r,1));
	call sort(categ);
	q = nrow(categ);

	/* creating the weights matrix */
	weights_mat = weight_mat(weights,categ);

	/* creating the nxq agreement matrix representing the distribution of raters by subjects and category */  
	agree_mat = j(n,q,0);
	do k=1 to q;
		if type(ratings_mat)="N" then do;
			in_categ_k =(ratings_mat=categ[k]);
			agree_mat[,k] = in_categ_k*repeat(1,r); 
		end;else agree_mat[,k] = (ratings_mat=categ[k])*repeat(1,r);
	end;
	agree_mat_w = t(weights_mat*t(agree_mat));
  
	/* calculating gwet's ac1 coefficient */  
	ri_vec = agree_mat*repeat(1,q);
	sum_q = (agree_mat#(agree_mat_w-1))*repeat(1,q);
	n2more = sum(ri_vec>=2);
	pa = sum(sum_q[loc(ri_vec>=2)]/((ri_vec#(ri_vec-1))[loc(ri_vec>=2)]))/n2more;
	pi_vec = t(t(repeat(1/n,n))*(agree_mat/(ri_vec*t(repeat(1,q)))));
	pe = sum(weights_mat) / (q**2);
	bp_coeff = (pa-pe)/(1-pe);
	/* calculating variance, stderr & p-value of gwet's ac1 coefficient */  
	den_ivec = ri_vec#(ri_vec-1);
	den_ivec = den_ivec - (den_ivec=0); /* this operation replaces each 0 value with -1 to make the next ratio calculation always possible */
	pa_ivec = sum_q/den_ivec;
	pe_r2 = pe*(ri_vec>=2);
	bp_ivec = (n/n2more)*(pa_ivec-pe_r2)/(1-pe);
	return(bp_ivec);
finish;

/*
#=====================================================================================
#fleiss_linear_i: This function produces all n linear components associated with the n 
	subjects, and the linearized version of the Fleiss' kappa coefficient.
			(Author: K. Gwet, PhD. E-mail: gwet@agreestat.com)
#-------------
#The input data "ratings" is a nxr data frame of raw alphanumeric ratings from n subjects 
	and r raters. Exclude all subjects that are not rated by any rater.
#Bibliography:
#Fleiss, J. L. (1981). Statistical Methods for Rates and Proportions. John Wiley & Sons.
#======================================================================================
*/
start fleiss_linear_i(ratings,weights="unweighted",conflev=0.95,Npop=10**7);
	ratings_mat = ratings; 
	if type(ratings_mat)="C" then do;
		ratings_mat = strip(upcase(ratings_mat));
	end;
	n  =  nrow(ratings_mat); * number of subjects;
	r  =  ncol(ratings_mat); * number of raters;
	f  =  n/Npop; * finite population correction;
  
	/* creating a vector containing all categories used by the raters */
	categ  =  unique_vec(shape(ratings_mat,n*r,1));
	call sort(categ);
	q = nrow(categ);

	/* creating the weights matrix */
	weights_mat = weight_mat(weights,categ);

	/* creating the nxq agreement matrix representing the distribution of raters by subjects and category */  
	agree_mat = j(n,q,0);
	do k=1 to q;
		if type(ratings_mat)="N" then do;
			in_categ_k =(ratings_mat=categ[k]);
			agree_mat[,k] = in_categ_k*repeat(1,r); 
		end;else agree_mat[,k] = (ratings_mat=categ[k])*repeat(1,r);
	end;
	agree_mat_w = t(weights_mat*t(agree_mat));

	/* calculating fleiss's generalized kappa coefficient */
	ri_vec = agree_mat*repeat(1,q);
	sum_q = (agree_mat#(agree_mat_w-1))*repeat(1,q);
	n2more = sum(ri_vec>=2);
	pa = sum(sum_q[loc(ri_vec>=2)]/((ri_vec#(ri_vec-1))[loc(ri_vec>=2)]))/n2more;
	pi_vec = t(t(repeat(1/n,n))*(agree_mat/(ri_vec*t(repeat(1,q)))));
	pe = sum(weights_mat # (pi_vec*t(pi_vec)));
	fleiss_kappa = (pa-pe)/(1-pe);

	/* calculating variance, stderr & p-value of gwet's ac1 coefficient */
	den_ivec = ri_vec#(ri_vec-1);
	den_ivec = den_ivec - (den_ivec=0); /* this operation replaces each 0 value with -1 to make the next ratio calculation always possible*/
	pa_ivec = sum_q/den_ivec;
	pe_r2 = pe*(ri_vec>=2);
	kappa_ivec = (n/n2more)*(pa_ivec-pe_r2)/(1-pe);
	pi_vec_wk_ = weights_mat*pi_vec;
	pi_vec_w_k = t(weights_mat)*pi_vec;
	pi_vec_w = (pi_vec_wk_ + pi_vec_w_k)/2;
	pe_ivec = (agree_mat*pi_vec_w)/ri_vec;
	kappa_ivec_x = kappa_ivec - 2*(1-fleiss_kappa) * (pe_ivec-pe)/(1-pe);
	return(kappa_ivec_x);
finish;