%simulate correlations with various risk rates
n=10000;
clear correlation
clear cohend
clear actualprop
clear variance
for run=1:99
    prop=run/100;
    clear data
    data=randn(n,1); %first column is a random normal variable
    data(:,2)=binornd(1,prop,[n,1]); %second column is group membership, random binary sequency of length n with probability of 1 at prop
    data(data(:,2)==1,1)=data(data(:,2)==1,1)+std(data(:,1)); %if group membership is 1, add 1 SD to the value
    actualprop(run)=sum(data(:,2)==1)./n; %let's see how many members the algorithm actually generated
    cohend(run)=computeCohen_d(data(data(:,2)==1,1),data(data(:,2)==0,1));
    correlation(run)=corr2(data(:,1),data(:,2));
    variance(run)=var(data(:,1));
end

plot(0.01:0.01:0.99,correlation)
hold on
ylabel('Biomarker-disease correlation')
yyaxis right
plot(0.01:0.01:0.99,cohend)
ylabel('Biomarker Cohen d, healthy vs. sick')
xlabel('Proportion of sample with disease')
set(gcf,'Color','white')

plot(0.01:0.01:0.99,variance)
hold on
ylabel('Biomarker variance')
xlabel('Proportion of sample with disease')
set(gcf,'Color','white')


%simulating vaccines based on Vokó 2022

pfizer.persondays=651.71; % *100k


pfizer.deaths_vacc=0.4;
pfizer.deaths_unvacc=1.56;
pfizer.efficacy=1-(pfizer.deaths_vacc/pfizer.deaths_unvacc);
sino.deaths_vacc=0.6;
sino.deaths_unvacc=1.79;
sino.efficacy=1-(sino.deaths_vacc/sino.deaths_unvacc);

clear correlation_pfizer
clear cohend_pfizer
clear cohend_sino
clear correlation_sino
for run=1:9999
    mort=run/10000; %mortality without vaccine
    clear vaccdata_pfizer
    vaccdata_pfizer=[zeros(n/2,1)' ones(n/2,1)']'; %half the sample is vaccinated
    vaccdata_pfizer(vaccdata_pfizer(:,1)==0,2)=binornd(1,mort,[n/2,1]); %unvaccinated are dead in proportion to mortality
    vaccdata_pfizer(vaccdata_pfizer(:,1)==1,2)=binornd(1,mort*(1-pfizer.efficacy),[n/2,1]); %unvaccinated are dead in proportion to mortality
    clear vaccdata_sino
    vaccdata_sino=[zeros(n/2,1)' ones(n/2,1)']'; %half the sample is vaccinated
    vaccdata_sino(vaccdata_sino(:,1)==0,2)=binornd(1,mort,[n/2,1]); %unvaccinated are dead in proportion to mortality
    vaccdata_sino(vaccdata_sino(:,1)==1,2)=binornd(1,mort*(1-sino.efficacy),[n/2,1]); %unvaccinated are dead in proportion to mortality
    cohend_pfizer(run)=computeCohen_d(vaccdata_pfizer(vaccdata_pfizer(:,2)==1,1),vaccdata_pfizer(vaccdata_pfizer(:,2)==0,1));
    correlation_pfizer(run)=corr2(vaccdata_pfizer(:,1),vaccdata_pfizer(:,2));
    cohend_sino(run)=computeCohen_d(vaccdata_sino(vaccdata_sino(:,2)==1,1),vaccdata_sino(vaccdata_sino(:,2)==0,1));
    correlation_sino(run)=corr2(vaccdata_sino(:,1),vaccdata_sino(:,2));
    tetracorr_sino(run)=tetrachoric(vaccdata_sino);
    tetracorr_pfizer(run)=tetrachoric(vaccdata_pfizer);
end

plot(0.01:0.01:99.99,correlation_pfizer)
hold on
plot(0.01:0.01:99.99,correlation_sino)
legend('Pfizer-efficacy vaccine', 'Sinopharm-efficacy vaccine')
% yyaxis right
% plot(0.01:0.01:99.99,cohend_pfizer)
% plot(0.01:0.01:99.99,cohend_sino)
% ylabel('Cohen d')
ylabel('Vaccination-mortality correlation')
xlabel('Disease mortality (%)')
set(gcf,'Color','white')


%Simulating a fixed mortality rate

mort=0.0018;
n2=1000000;
clear vaccdata2_pfizer
vaccdata_pfizer2=[zeros(n2/2,1)' ones(n2/2,1)']'; %half the sample is vaccinated
vaccdata_pfizer2(vaccdata_pfizer2(:,1)==0,2)=binornd(1,mort,[n2/2,1]); %unvaccinated are dead in proportion to mortality
vaccdata_pfizer2(vaccdata_pfizer2(:,1)==1,2)=binornd(1,mort*(1-pfizer.efficacy),[n2/2,1]); %unvaccinated are dead in proportion to mortality
clear vaccdata2_sino
vaccdata_sino2=[zeros(n2/2,1)' ones(n2/2,1)']'; %half the sample is vaccinated
vaccdata_sino2(vaccdata_sino2(:,1)==0,2)=binornd(1,mort,[n2/2,1]); %unvaccinated are dead in proportion to mortality
vaccdata_sino2(vaccdata_sino2(:,1)==1,2)=binornd(1,mort*(1-sino.efficacy),[n2/2,1]); %unvaccinated are dead in proportion to mortality
cohend_pfizer2=computeCohen_d(vaccdata_pfizer2(vaccdata_pfizer2(:,2)==1,1),vaccdata_pfizer2(vaccdata_pfizer2(:,2)==0,1));
correlation_pfizer2=corr2(vaccdata_pfizer2(:,1),vaccdata_pfizer2(:,2));
cohend_sino2=computeCohen_d(vaccdata_sino2(vaccdata_sino2(:,2)==1,1),vaccdata_sino2(vaccdata_sino2(:,2)==0,1));
correlation_sino2=corr2(vaccdata_sino2(:,1),vaccdata_sino2(:,2));

oddstab=crosstab(vaccdata_pfizer2(:,1),vaccdata_pfizer2(:,2));
odds(oddstab)