% 09/08/2022 -- Code for generating results for the CS ranking study,
% Created by Aniruddha Maiti, edited by Sai Shi

%
% LOAD DATA
%
% Load usnews data
[~, ~, us] = xlsread('all_departments.xlsx','Sheet1');
us = us(2:end,:);
us(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),us)) = {''};

% Load faculty citations data
[~, ~, d] = xlsread('all_professors.xlsx','Sheet1');
d = d(2:end,:);
d(cellfun(@(x) ~isempty(x) && isnumeric(x) && isnan(x),d)) = {''};
d(cellfun(@(x) ~isempty(strmatch('NaN',x)),d)) = {NaN};

%
% BASIC STATISTICS 
% how many professors
disp('Total CS Professors:')
disp(size(d))
disp('Full Professors')
disp(length(strmatch('Full',d(:,10))))
disp('Associate Professors')
disp(length(strmatch('Associate',d(:,10))))
disp('Assistant Professors')
disp(length(strmatch('Assistant',d(:,10))))

% Duplicate d
%
% d(:, 9) = regexprep(d(:, 9),'[^a-zA-Z0-9]','');
d1 = d;
%
% CREATE UNIVERSITY ATTRIBUTES
%
u = unique(d1(:,9));
u = u(~cellfun('isempty',u));
univ = cell(length(u),16);

for i=1:length(u)
    univ{i,1} = u{i};       % UNIVERSITY NAME
    %q = find(strcmp(u(i), d1(:,9)));
    q = strmatch(u(i),d1(:,9),'exact');
    d2 = d1(q,:);
    univ{i,2} = length(q);          % TOTAL PROFESSORS
    %univ{i,3} = length(strcmp('Assistant',d2(:,10))); % TOTAL ASSISTANT
    univ{i,3} = length(strmatch('Assistant',d2(:,10)));
end
%
dtmp=d1;
q=strmatch('Full',dtmp(:,10));
d1 = dtmp(q,:);
q=strmatch('Associate',dtmp(:,10));
d1 = [d1; dtmp(q,:)];   % d1 now has only Associate and Full Professors

%
% Introduce new column in d1 (perc_t10):
%  For each professor, determine the prcentile rank based on t10
%   used formula from https://en.wikipedia.org/wiki/Percentile_rank
d1(~cellfun(@isnumeric,d1(:,8)), 8) = {NaN};
d(~cellfun(@isnumeric,d(:,8)), 8) = {NaN};

tmp = cell2mat(d1(:,8));
q = find(~isnan(tmp));
scores = tmp(q);
d1(q,15) = num2cell((tiedrank(scores) - 0.5)/length(scores));
q = find(isnan(tmp));
d1(q,15) = {NaN};
%
% ### TABLE 1 ####: find t10 thresholds for different percentiles 
tmp = cell2mat(d1(:,[8 15]));

theta = [0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.98 0.99]';
for i=1:length(theta)
    q=find(tmp(:,2)>theta(i));
    theta(i,2)=min(tmp(q,1));
end
disp('Content of Table 1: t10 percentiles for senior faculty:')
disp(theta)

% Faculty with Google Scholar profile and with t10 index
d(~cellfun(@isnumeric,d(:,6)), 6) = {NaN};
d1(~cellfun(@isnumeric,d1(:,6)), 6) = {NaN};


tmp = cell2mat(d(:,6));
q = find(~isnan(tmp));
disp('Number of CS faculty with Google Scholar profiles:')
disp(length(q))
tmp = cell2mat(d1(:,6));
q = find(~isnan(tmp));
disp('Number of senior faculty with Google Scholar profiles:')
disp(length(q))
tmp = cell2mat(d(:,8));
q = find(~isnan(tmp));
disp('Number of CS faculty with t10 index:')
disp(length(q))
tmp = cell2mat(d1(:,8));
q = find(~isnan(tmp));
disp('Number of senior faculty with t10 index:')
disp(length(q));

% Create the remaining university-level attributes (we ignore total, i10)
for i=1:length(u)
    %q = strcmp(u(i),d1(:,9));
    %q = find(strcmp(u(i), d1(:,9)));
    q = strmatch(u(i),d1(:,9),'exact');

    univ{i,4} = length(q);          % total associate and full professors
    tmp = cell2mat(d1(q,[6 8 15])); % h-index, t10, perc_t10 scores   
    %tmp(:,6) = geomean(1+ tmp(:,4:5),2);

    univ{i,5} = length(find(~isnan(tmp(:,1))));   % total with h-index
    univ{i,6} = length(find(~isnan(tmp(:,2))));   % total with t10
    %
    % AGGREGATED measures
    univ{i,7} = median(tmp(find(~isnan(tmp(:,1))),1));  % median h-index
    univ{i,8} = median(tmp(find(~isnan(tmp(:,2))),2));  % median t10
    univ{i,9} = mean(tmp(find(~isnan(tmp(:,3))),3));  % mean perc_t10
    univ{i,10} = geomean(1+ tmp(find(~isnan(tmp(:,2))),2));  % geomean t10
    %
    % CUMMULATIVE measure that counts assistant professors too
    q = find(strcmp(u(i), d(:,9)));
    tmp = cell2mat(d(q,[6 8])); % h-index, t10, perc_t10 scores       
    univ{i,11} = length(find(tmp(:,2)>=theta(2,2))); % how many have t10** over 20%
    univ{i,12} = length(find(tmp(:,2)>=theta(4,2))); % how many have over 40%
    univ{i,13} = length(find(tmp(:,2)>=theta(6,2))); % how many have over 60%
    univ{i,14} = length(find(tmp(:,2)>=theta(8,2))); % how many have over 80%
    univ{i,15} = length(find(tmp(:,2)>=theta(9,2))); % how many have over 90%
    %
    univ{i,16} = d{q(1),11};   % university ID
        

end

%
% join univ and us data into joindata
%
univ(~cellfun(@isnumeric,univ(:,16)), 16) = {NaN};


a1 = cell2mat(univ(:,16));
a2 = cell2mat(us(:,7));
a3 = intersect(a1,a2);

joindata = cell(length(a3),23);
% disp(length(a1))
% disp(length(a2))
% disp(length(joindata))
% disp(setxor(a1,a2))

for i=1:length(a3)
    q=find(a1 == a3(i));
    joindata(i,1:16) = univ(q(1),:);
    q=find(a2 == a3(i));
    joindata(i,17:23) = us(q(1),1:7);
end

% ### FIGURE 1 ###: FACULTY SIZE AND CORRELAITON WITH THE SCORE
figure(1)
tmp = cell2mat(joindata(:,[2 18]));
[a1,x] = hist(tmp(:,1),20);
bar(x,a1)   
q=find(tmp(:,2)>0);
xlabel('Faculty Size')
ylabel('Count')
ylim([0 50])
disp('Correlation between faculty size and US News CS ranking:')
disp(corrcoef(tmp(q,:)))

% Correlation between US News CS ranking and university ranking
tmp = cell2mat(joindata(:,[18 22]));
q=find(tmp(:,1)>0 & tmp(:,2)>0);
disp('Correlation between US News CS ranking and university ranking:')
disp(corrcoef(tmp(q,:)))

% #### FIGURE 2 ###: Distribution of t10
figure(2)
tmp = cell2mat(d1(:,8));    % t10 index
q = find(~isnan(tmp));     % find all faculty with t10 index
tmp = tmp(q,:);
[a1,x] = hist(log10(1+tmp),20);
bar(x,a1)   
xlabel('log10(t10)')
ylabel('Count')

% FIND CORRELATION BETWEEN H-INDEX AND T10
%
% Select only professors with both h-index and t10
a = cell2mat(d1(:,[6 8]));    % h-index, t10 indices
q = find(~isnan(a(:,1)) & ~isnan(a(:,2)));
disp('Correlation between h-index and t10-index:')
disp(corrcoef(a(q,1:2)))
% take the log transform of both, because the distributions are lognormal
tmp = log10(1+a(q,1:2));
disp('Correlation between logs of h-index and t10-index:')
disp(corrcoef(tmp))
% plot(tmp(:,1),tmp(:,2),'.')
% loglog(a(:,1),a(:,2),'.')

%
% SHOW THAT FACULTY WITH GOOGLE PROFILES ARE UPWARD BIASED
%
a = cell2mat(d1(:,[6 8]));    % h-index, t10 indices
q = find(~isnan(a(:,2)));       % find all faculty with t10 index
disp('number, geomean, median: faculty with t10')
disp([length(q) geomean(1 + a(q,2))  median(a(q,2))])
q = find(~isnan(a(:,1)) & ~isnan(a(:,2))); % find all faculty with google scholar
disp('number, geomean, median: faculty with Google Scholar')
disp([length(q) geomean(1 + a(q,2))  median(a(q,2))])
q = find(isnan(a(:,1)) & ~isnan(a(:,2))); % find all faculty without google scholar
disp('number, geomean, median: faculty without Google Scholar')
disp([length(q) geomean(1 + a(q,2))  median(a(q,2))])
%
% ### FIGURE 3 ###: FRACTION OF FACULTY WITH SCHOLAR PROFILE IN 10 PERCENTILE BINS:
%
figure(3)
a = cell2mat(d1(:,[6 8 15]));    % h-index, t10 indices
x = 0.05:0.1:0.95;
q = find(~isnan(a(:,1))); % find all faculty with google scholar
[a2,x] = hist(a(q,3),x);
q = find(isnan(a(:,1))); % find all faculty without google scholar
[a3,x] = hist(a(q,3),x);
%
bar(x,[a2;a3]','stacked')   % Stacked bar, shows where are people with and witout google scholar
xlabel('Percentile Range')
ylabel('Count')
%
% ### TABLE 2 ###: FIND CORRELATION BETWEEN THE ATTRIBUTES AND USNEWS SCORES 
%
% create matrices
joindata(~cellfun(@isnumeric,joindata(:,18)), 18) = {NaN};
joindata(cellfun(@isempty, joindata)) = {0};
%disp(joindata);

a1 = cell2mat(joindata(:,18));      % target: us news score
a2 = cell2mat(joindata(:,8:15));    % potential features
% for different thresholds for inclusion of departments based on usnews score
T = [1 2.01 2.21 2.41 2.71 2.91 3.11];
for i=1:length(T)
    q = find(a1>T(i));
    c = corrcoef([a2(q,:) a1(q,1)]);
    res1(i,1:8) = c(end,1:8);  % correlations with the target
    res1(i,8) = length(q);      % how many departments are included
end
disp('Correlations with original measures')
disp(res1)
% 
%  Taking the logarithm:
a3 = log10(1+a2);
T = [1 2.01 2.21 2.41 2.71 2.91 3.11];
for i=1:length(T)
    q = find(a1>T(i));
    c = corrcoef([a3(q,:) a1(q,1)]);
    res2(i,1:8) = c(end,1:8);  % correlations with the target
    res2(i,8) = length(q);      % how many departments are included
end
disp('Correlations with logarithm of the measures')
disp(res2)
%  Taking the square root
a3 = sqrt(a2);
T = [1 2.01 2.21 2.41 2.71 2.91 3.11];
for i=1:length(T)
    q = find(a1>T(i));
    c = corrcoef([a3(q,:) a1(q,1)]);
    res3(i,1:8) = c(end,1:8);  % correlations with the target
    res3(i,8) = length(q);      % how many departments are included
end
disp('Correlations with squared root of the measures')
disp(res3)
% TABLE 2:
disp('Content of Table 2: Correlations with US News CS scores')
disp([res1(1,:); res2(1,:); res3(1,:)])

%
% ### FIGURE 4 ###:
%
figure(4)
a = cell2mat(univ(:,8));    % t10 index
[a1,x] = hist(a,20);
bar(x,a1) 
xlabel('Faculty Size')
ylabel('Count')

%
% REGRESSION MODEL
%
y = cell2mat(joindata(:,18));      % target: us news score
x = cell2mat(joindata(:,[8:10 12:14]));    % potential features
x(:,[1 3 4 5 6]) = sqrt(x(:,[1 3 4 5 6]));

x = [ones(length(x),1) x];
q = find(y > 1);
II = [2 4 3];
JJ = [5 6 7];
yh = [];
for I = 1:length(II)
    for J = 1:length(JJ)
        b{I,J} = regress(y(q)-1,x(q,[II(I) JJ(J)]));
        yhat = 1 + x(:,[II(I) JJ(J)])*b{I,J};    % apply on all data
        yh = [yh yhat]; 
        rsq(I,J) = (1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2));
        tmpp = corrcoef(y(q),yhat(q));
        cc(I,J) = tmpp(1,2);
    end
end
% averaged model
yhat = mean(yh,2);
disp(' R-squared of the averaged model, using 9 individual predictors')
disp(1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2))
disp(' corrcoef of the averged model, using 9 individual predictors')
disp(corrcoef(y(q),yhat(q)))

% ### TABLE 3 ###: R-SQUARED OF INDIVIDUAL MODELS
%
disp('R-squared values for individual models')
disp(rsq)
% ### TABLE 4 ###: correlation coefficients OF INDIVIDUAL MODELS
%
disp('Correlation coefficients values for individual models')
disp(cc)

%
% JOINT MODEL by combining 4 best individual models
%
clear cc rsq b
y = cell2mat(joindata(:,18));      % target: us news score
x = cell2mat(joindata(:,[8:10 12:14]));    % potential features
x(:,[1 3 4 5 6]) = sqrt(x(:,[1 3 4 5 6]));
x = [ones(length(x),1) x];
q = find(y > 1);
II = [2 4];
JJ = [5 6];
yh = [];
for I = 1:length(II)
    for J = 1:length(JJ)
        b{I,J} = regress(y(q)-1,x(q,[II(I) JJ(J)]));
        yhat = 1 + x(:,[II(I) JJ(J)])*b{I,J};    % apply on all data
        yh = [yh yhat]; 
        rsq(I,J) = (1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2));
        tmpp = corrcoef(y(q),yhat(q));
        cc(I,J) = tmpp(1,2);
    end
end
% averaged model
yhat = mean(yh,2);
disp('R-squared of JOINT model')
disp(1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2))
disp(' corrcoef of JOINT model')
disp(corrcoef(y(q),yhat(q)))
%
% JOINT MODEL FORMULA
%
alpha = zeros(4,1);
for I = 1:length(II)
    for J = 1:length(JJ)
        alpha(I) = alpha(I) + b{I,J}(1);
        alpha(length(II)+J) = alpha(length(II)+J) + b{I,J}(2);
    end
end
alpha = alpha / 4; %  => AVERAGE THE COEFFICIENTS
% FINAL FORMULA:
disp(' JOINT MODEL formula')
fprintf('1 + %3.3f*sqrt(m10) + %3.3f*sqrt(g10) +%3.3f*sqrt(c40) + %3.3f*sqrt(c60) \n',alpha)


%text(y,yhat+0.01,joindata(:,1),'HorizontalAlignment','center')
%

%
% INTERESTING INSIGHT
%
% Low correlation for 57 second-tier departments 
disp(length(find(y>0)))
disp(' corrcoef between USN and Joint model ranking, for programs ranked from 0 - 2.7 by USN')
q=find(y>0 & y<2.8); disp(length(q))
corrcoef(y(q),yhat(q))
% High correlation for 73 top-tier departments
disp(' corrcoef between USN and Joint model ranking, for programs ranked above 2.7 by USN')
q=find(y>2.7 ); disp(length(q))
corrcoef(y(q),yhat(q))
disp('R-squared between USN and Joint model ranking, for programs ranked above 2.7 by USN')
disp(1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2))

%
% SCHOLAR MODEL by combining 4 best individual models and training on top
% 73 CS programs
%
clear cc rsq b
y = cell2mat(joindata(:,18));      % target: us news score
x = cell2mat(joindata(:,[8:10 12:14]));    % potential features
x(:,[1 3 4 5 6]) = sqrt(x(:,[1 3 4 5 6]));
x = [ones(length(x),1) x];
q = find(y > 2.7);      % SELECT TOP 73 PROGRAMS BY USN RANKING
II = [2 4];
JJ = [5 6];
yh = [];
for I = 1:length(II)
    for J = 1:length(JJ)
        b{I,J} = regress(y(q)-1,x(q,[II(I) JJ(J)]));
        yhat = 1 + x(:,[II(I) JJ(J)])*b{I,J};    % apply on all data
        yh = [yh yhat]; 
        rsq(I,J) = (1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2));
        tmpp = corrcoef(y(q),yhat(q));
        cc(I,J) = tmpp(1,2);
    end
end
% averaged model
yhat = mean(yh,2);
q = find(y > 0);      % SELECT TOP 73 PROGRAMS BY USN RANKING
disp('R-squared of SCHOLAR model')
disp(1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2))
disp(' corrcoef of SCHOLAR model')
disp(corrcoef(y(q),yhat(q)))
q = find(y > 0 & y < 2.8);
disp('R-squared of SCHOLAR model on programs ranked from 0 to 2.7 by USN')
disp(1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2))
disp(' corrcoef of SCHOLAR model')
disp(corrcoef(y(q),yhat(q)))
q = find(y > 2.7);
disp('R-squared of SCHOLAR model on programs ranked above 2.7 by USN')
disp(1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2))
disp(' corrcoef of SCHOLAR model')
disp(corrcoef(y(q),yhat(q)))
%
% SCHOLAR MODEL FORMULA
%
alpha = zeros(4,1);
for I = 1:length(II)
    for J = 1:length(JJ)
        alpha(I) = alpha(I) + b{I,J}(1);
        alpha(length(II)+J) = alpha(length(II)+J) + b{I,J}(2);
    end
end
alpha = alpha / 4; %  => AVERAGE THE COEFFICIENTS
% FINAL FORMULA:
disp(' SCHOLAR formula');

alpha = [0.058; 0.059; 0.121; 0.127];  % Use the same scholar model obtained in 2016
fprintf('1 + %3.3f*sqrt(m10) + %3.3f*sqrt(g10) + %3.3f*sqrt(c40) + %3.3f*sqrt(c60) \n',alpha);
yhat = 1 + x(:,[2 4 5 6])*alpha;


%
% VISUALIZE THE RESULTS
%
% FIGURE 5: Show as a scatterplot
figure(5)
plot([1.5 5.2], [1.5 5.2],'Linewidth',5)
hold
yhat(yhat>5)=5.0;
scatter(y, yhat, 80, 'o', 'black', 'LineWidth', 2)
xticks(0:0.5:5.5)
xtickformat('%.1f')
xlabel('U.S. News Score')
ylabel('Model Score')

disp(' corrcoef of SCHOLAR model using 2017 formula: ')
disp(corrcoef(y(q),yhat(q)))
% Low correlation for 57 second-tier departments 
disp(length(find(y>0)))
disp(' corrcoef between USN and 2017 scholar model ranking, for programs ranked from 0 - 2.7 by USN')
q=find(y>0 & y<2.8); disp(length(q))
corrcoef(y(q),yhat(q))
% High correlation for 73 top-tier departments
disp(' corrcoef between USN and 2017 scholar model ranking, for programs ranked above 2.7 by USN')
q=find(y>2.7 ); disp(length(q))
corrcoef(y(q),yhat(q))
disp('R-squared between USN and 2017 scholar model ranking, for programs ranked above 2.7 by USN')
disp(1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2))

% SCHOLAR Ranking
%  create ranks based on 1 significant digit of the score
yhat(yhat>5) = 5;   % cap all the scores to 5.0
tmp = sort(-round(yhat*10));
[~, rnk] = ismember(-round(yhat*10),tmp);
[~,i] = sort(rnk);
tmp = cell2mat(joindata(i,[2 8 10 9 12 13 14]));
tmp(:,[1:3 5:7]) = round(tmp(:,[1:3 5:7]));
tmp(:,4) = round(tmp(:,4)*100)/100;
% TABLE 4: Scholar Ranking
mycell = [num2cell(rnk(i)) joindata(i,1)  num2cell(tmp) ...
    joindata(i,18) num2cell(round(yhat(i)*10)/10) joindata(i,23)];
[nrows,ncols]= size(mycell);
filename = 'scholar_ranking_data_2022.csv';
fid = fopen(filename, 'w');
% Semi-column separated file
fprintf(fid, 'Rank,University,Size,M10,G10,P10,C40,C60,C80,USN,Scholar,ID\n'); 
for row=1:nrows
    fprintf(fid, '%d,%s,%d,%d,%d,%3.2f,%d,%d,%d,%2.1f,%2.1f,%d\n', mycell{row,:});
end
fclose(fid);

%
% ADDING USN UNIVERSITY RANKING AS AN ADDITIONAL INPUT TO SCHOLAR RANKING MODEL
%
y = cell2mat(joindata(:,18));      % target: us news score
x = cell2mat(joindata(:,[8:10 12:14 22]));    % potential features
x(:,[1 3 4 5 6]) = sqrt(x(:,[1 3 4 5 6]));
x = [ones(length(x),1) x];
q = find(y > 1);
II = [2 4];
JJ = [5 6];
yh = [];
for I = 1:length(II)
    for J = 1:length(JJ)
        b{I,J} = regress(y(q)-1,x(q,[II(I) JJ(J) 8]));
        yhat = 1 + x(:,[II(I) JJ(J) 8])*b{I,J};    % apply on all data
        yh = [yh yhat]; 
        rsq(I,J) = (1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2));
        tmpp = corrcoef(y(q),yhat(q));
        cc(I,J) = tmpp(1,2);
    end
end
% averaged model
yhat = mean(yh,2);
disp('R-squared of enhanced Scholar model that includes USN university ranking')
disp(1 - mean((y(q)-yhat(q)).^2)/mean((y(q)-mean(y(q))).^2))
disp('corrcoef of enhanced Scholar model that includes USN university ranking')
disp(corrcoef(y(q),yhat(q)))

% R-SQUARED OF INDIVIDUAL MODELS for the enhanced mode
%
disp('R-squared values for individual models, enhanced')
disp(rsq)
% ### TABLE 4 ###: correlation coefficients OF INDIVIDUAL MODELS
%
disp('Correlation coefficients values for individual models, enhanced')
disp(cc)
