Commit 4a71cbb5 authored by Kristian Soltesz's avatar Kristian Soltesz
Browse files

Adding code and data

parent 98ec054c
#Statistikdatum Totalt_antal_fall Antal_avlidna Antal_intensivvårdade
20-03-01 0 0 0
20-03-02 5 0 0
20-03-03 13 0 0
20-03-04 30 0 0
20-03-05 25 0 0
20-03-06 59 0 1
20-03-07 33 0 1
20-03-08 46 0 1
20-03-09 101 0 0
20-03-10 98 0 2
20-03-11 196 1 1
20-03-12 151 0 0
20-03-13 152 1 2
20-03-14 71 1 6
20-03-15 69 2 5
20-03-16 83 2 5
20-03-17 119 1 3
20-03-18 145 6 15
20-03-19 143 7 8
20-03-20 180 9 17
20-03-21 136 8 13
20-03-22 118 11 27
20-03-23 182 11 37
20-03-24 230 21 32
20-03-25 314 22 31
20-03-26 286 31 40
20-03-27 365 32 29
20-03-28 300 35 28
20-03-29 280 38 41
20-03-30 416 45 29
20-03-31 474 48 34
20-04-01 486 53 48
20-04-02 554 70 47
20-04-03 601 80 39
20-04-04 357 70 36
20-04-05 340 85 44
20-04-06 389 90 44
20-04-07 739 84 45
20-04-08 654 115 47
20-04-09 645 86 36
20-04-10 454 90 36
20-04-11 395 102 45
20-04-12 464 97 36
20-04-13 437 84 43
20-04-14 479 91 41
20-04-15 604 115 32
20-04-16 623 111 33
20-04-17 688 83 39
20-04-18 532 86 29
20-04-19 388 87 33
20-04-20 461 85 28
20-04-21 707 62 34
20-04-22 722 77 49
20-04-23 758 86 27
20-04-24 780 89 46
20-04-25 473 73 28
20-04-26 300 74 26
20-04-27 563 74 28
20-04-28 742 83 32
20-04-29 799 83 25
20-04-30 636 78 33
20-05-01 532 78 15
20-05-02 299 72 28
20-05-03 261 75 27
20-05-04 477 84 25
20-05-05 657 72 20
20-05-06 745 73 26
20-05-07 785 80 28
20-05-08 701 60 26
20-05-09 509 68 14
20-05-10 279 73 17
20-05-11 455 65 15
20-05-12 754 61 16
20-05-13 698 50 19
20-05-14 657 46 15
20-05-15 688 58 21
20-05-16 358 48 18
20-05-17 259 52 20
20-05-18 430 59 23
20-05-19 666 41 13
20-05-20 808 52 14
20-05-21 610 53 13
20-05-22 532 56 15
20-05-23 403 55 16
20-05-24 210 44 17
20-05-25 491 42 28
20-05-26 746 28 14
20-05-27 800 38 15
20-05-28 774 40 19
20-05-29 773 40 13
20-05-30 432 38 19
20-05-31 265 46 14
20-06-01 648 39 20
20-06-02 900 36 15
20-06-03 1046 28 18
20-06-04 1039 43 16
20-06-05 1146 37 20
20-06-06 780 29 20
20-06-07 462 33 14
20-06-08 677 38 16
20-06-09 937 33 16
20-06-10 1437 39 9
20-06-11 1293 35 12
20-06-12 1329 29 12
20-06-13 1032 33 14
20-06-14 418 27 12
20-06-15 685 30 13
20-06-16 1209 28 9
20-06-17 1457 32 13
20-06-18 1494 29 13
20-06-19 1209 30 8
20-06-20 698 29 6
20-06-21 321 22 12
20-06-22 800 20 15
20-06-23 1309 25 7
20-06-24 1698 22 6
20-06-25 1278 23 10
20-06-26 1204 11 5
20-06-27 755 14 5
20-06-28 415 23 8
20-06-29 726 16 4
20-06-30 804 20 5
20-07-01 685 15 6
20-07-02 687 15 10
20-07-03 694 8 3
20-07-04 364 15 1
20-07-05 315 9 5
20-07-06 251 15 2
20-07-07 278 12 2
20-07-08 533 10 1
20-07-09 334 15 2
20-07-10 369 14 2
20-07-11 308 10 2
20-07-12 106 8 3
20-07-13 170 12 3
20-07-14 312 8 2
20-07-15 287 6 0
20-07-16 268 6 4
20-07-17 284 7 3
20-07-18 191 11 0
20-07-19 110 7 2
20-07-20 131 6 3
20-07-21 226 7 5
20-07-22 297 6 1
20-07-23 220 5 1
20-07-24 262 3 1
20-07-25 138 1 0
20-07-26 42 2 1
20-07-27 71 6 1
20-07-28 283 4 0
20-07-29 301 1 2
20-07-30 302 0 0
20-07-31 258 2 1
20-08-01 303 2 3
20-08-02 38 3 4
20-08-03 165 4 1
20-08-04 333 2 1
20-08-05 425 1 3
20-08-06 378 4 2
20-08-07 380 2 2
20-08-08 260 1 5
20-08-09 73 4 2
20-08-10 196 2 1
20-08-11 417 4 3
20-08-12 443 3 2
20-08-13 363 5 4
20-08-14 344 1 0
20-08-15 226 1 2
20-08-16 63 0 1
20-08-17 174 3 1
20-08-18 314 4 2
20-08-19 351 1 2
20-08-20 333 2 1
20-08-21 298 5 0
20-08-22 160 1 1
20-08-23 57 3 1
20-08-24 174 1 3
20-08-25 222 1 1
20-08-26 244 2 1
20-08-27 202 1 1
20-08-28 179 1 0
20-08-29 131 1 0
20-08-30 48 3 1
20-08-31 162 2 0
20-09-01 171 3 1
20-09-02 213 2 3
20-09-03 286 2 3
20-09-04 262 0 0
20-09-05 171 0 1
20-09-06 67 3 0
20-09-07 185 1 0
20-09-08 236 1 0
20-09-09 314 2 1
20-09-10 254 2 0
20-09-11 291 4 2
20-09-12 206 1 4
20-09-13 106 2 1
20-09-14 220 2 1
20-09-15 292 1 2
20-09-16 330 2 1
20-09-17 389 1 0
20-09-18 437 1 1
20-09-19 279 1 2
20-09-20 133 4 0
20-09-21 266 2 0
20-09-22 438 1 2
20-09-23 553 0 3
20-09-24 540 1 2
20-09-25 630 4 1
20-09-26 325 2 0
20-09-27 167 1 0
20-09-28 378 1 3
20-09-29 613 2 0
20-09-30 689 4 4
20-10-01 633 1 3
20-10-02 712 3 0
20-10-03 461 3 3
20-10-04 156 3 1
20-10-05 374 2 5
20-10-06 786 4 2
20-10-07 831 3 3
20-10-08 834 1 3
20-10-09 783 5 3
20-10-10 509 5 5
20-10-11 161 2 1
20-10-12 637 3 3
20-10-13 916 1 3
20-10-14 968 2 4
20-10-15 902 3 1
20-10-16 1179 2 2
20-10-17 697 4 1
20-10-18 321 1 5
20-10-19 771 4 4
20-10-20 1289 4 2
20-10-21 1571 3 4
20-10-22 1666 9 9
20-10-23 1868 7 7
20-10-24 1477 8 4
20-10-25 514 8 7
20-10-26 1069 11 5
20-10-27 2414 9 10
20-10-28 3390 9 12
20-10-29 3262 9 6
20-10-30 4056 9 10
20-10-31 2987 13 7
20-11-01 1297 22 9
20-11-02 1569 20 20
20-11-03 3608 20 14
20-11-04 4482 21 14
20-11-05 4744 22 14
20-11-06 4454 25 15
20-11-07 4451 27 16
20-11-08 2097 22 17
20-11-09 3725 35 11
20-11-10 4497 35 9
20-11-11 5709 27 18
20-11-12 5564 28 13
20-11-13 6730 33 30
20-11-14 3513 38 21
20-11-15 1581 39 18
20-11-16 2544 37 22
20-11-17 4458 44 25
20-11-18 4961 50 20
20-11-19 7609 46 22
20-11-20 5459 45 16
20-11-21 4492 51 19
20-11-22 2422 59 29
20-11-23 3876 55 24
20-11-24 5674 68 30
20-11-25 6066 69 24
20-11-26 6903 60 28
20-11-27 6462 65 15
20-11-28 3825 50 22
20-11-29 2753 50 27
20-11-30 3482 66 23
20-12-01 5820 68 19
20-12-02 6547 71 40
20-12-03 6986 77 23
20-12-04 7338 70 17
20-12-05 4864 43 18
20-12-06 1803 81 26
20-12-07 3782 58 27
20-12-08 7446 54 20
20-12-09 8398 75 22
20-12-10 6693 86 32
20-12-11 7718 88 25
20-12-12 6538 67 17
20-12-13 3057 71 22
20-12-14 3534 94 28
20-12-15 7001 88 32
20-12-16 8825 102 24
20-12-17 9640 117 27
20-12-18 7918 71 29
20-12-19 5510 80 22
20-12-20 3752 92 25
20-12-21 5160 77 31
20-12-22 6597 92 36
20-12-23 11376 87 22
20-12-24 5036 87 19
20-12-25 2793 106 31
20-12-26 2966 94 35
20-12-27 3209 98 33
20-12-28 7095 121 40
20-12-29 8874 72 42
20-12-30 10459 90 30
20-12-31 7002 110 35
21-01-01 2618 86 35
21-01-02 2532 98 18
21-01-03 2756 83 33
21-01-04 6984 86 20
21-01-05 7546 93 36
21-01-06 5031 82 27
21-01-07 7152 102 23
21-01-08 5699 99 30
21-01-09 4737 87 24
21-01-10 2318 90 27
21-01-11 4648 88 24
21-01-12 5342 88 33
21-01-13 6597 89 18
21-01-14 4749 95 26
21-01-15 4209 85 15
21-01-16 2162 81 20
21-01-17 1250 58 17
21-01-18 2118 72 26
21-01-19 4739 86 29
21-01-20 4950 85 15
21-01-21 4204 86 16
21-01-22 3724 57 15
21-01-23 2369 57 11
21-01-24 1148 49 9
21-01-25 1889 42 28
21-01-26 4190 37 18
21-01-27 4109 22 13
21-01-28 2376 22 10
21-01-29 4262 28 14
21-01-30 2695 15 13
21-01-31 1165 11 8
21-02-01 1560 14 12
21-02-02 4325 3 7
21-02-03 3745 2 7
function NyT2021
% Generate all three figures in NyT manuscript
close all
wacker;
sentinel
figure(3),figure(2),figure(1)
end
function sentinel()
data=[...
9 752 13 0 0;...
10 4302 211 0 0;...
11 8990 835 0 0;...
12 10404 911 0 0;...
13 12349 1943 0 0;...
14 17776 3210 0 0;...
15 19880 3711 0 0;...
16 20233 3740 0 0;...
17 24288 4160 0 0;...
18 28997 3728 0 0;...
19 28988 4029 0 0;...
20 32665 3652 0 0;...
21 28796 3466 0 0;...
22 36466 4300 0 0;...
23 49162 6060 0 0;...
24 60678 7229 0 0;...
25 61803 7462 0 0;...
26 75171 7645 0 0;...
27 77642 4935 0 0;...
28 81801 2789 0 0;...
29 69393 2274 0 0;...
30 59143 1376 0 0;...
31 52959 1598 0 0;...
32 53721 2016 0 0;...
33 56487 2058 0 0;...
34 65546 1687 0 0;...
35 85060 1212 0 0;...
36 126219 1335 0 0;...
37 142673 1598 0 0;...
38 139471 2084 0 0;...
39 128852 2926 0 0;...
40 130058 3641 0 0;...
41 137300 4316 36500 3985;...
42 148267 5717 32503 3823;...
43 164742 9223 26726 3885;...
44 189301 18432 25049 3883;...
45 228023 24756 27279 4961;...
46 254295 32705 29550 6077;...
47 260710 30975 30946 7223;...
48 275712 36675 33160 8576;...
49 261229 37544 32280 9317;...
50 270944 46210 37216 11695;...
51 299447 46713 43137 13548;...
52 232114 37276 26885 8912;...
1 198692 41686 22181 8476;...
2 200844 29544 33268 13210;...
];
disp(['Total number of tests: ',num2str(sum(data(:,2)))])
disp(sprintf('Total number of positives : %2.1f (%2.2f)',sum(data(:,3)), sum(data(:,3))/sum(data(:,2))))
[prev,phat,stdphat]=randompcr;
figure(3), plotfix
N=size(data,1);
p1=data(:,3)./data(:,2);
p1s=p1.*(1-p1)/N;
p2=data(:,5)./data(:,4);
p2s=p1.*(1-p1)/N;
h1=plot(8*7+(7:7:7*N),[p1]/20*100,'r-');
hold on
patch(8*7+7*[1:N N:-1:1],[(p1+p1s)' fliplr((p1-p1s)')]/20*100,[1 0.5 0.5],'facealpha',0.7)
for i=1:5
h5=patch((prev.date(i)+prev.days(i)*[0 1 1 0 0]),max(0,phat(i)+stdphat(i)*[1 1 -1 -1 1])*100,[0.5 0.5 1],'facealpha',0.7);
end
legend([h1(1) h5(1)],{'1/20 * Allmänna PCR','Randomiserade PCR'},'location','north')
xdate
ylabel('Skattad prevalens (%)')
end
function [prev,phat, stdphat]=randompcr()
% Raw data
jan1=datenum(2020,1,1);
prev.date=[datenum(2020,04,21) datenum(2020,05,25) datenum(2020,08,24) datenum(2020,09,21) datenum(2020,11,30) ]-jan1;
prev.days=[4 4 5 4 5];
prev.N=[2571 2957 2518 2461 2983];
prev.n=[23 9 0 0 21];
% Prevalence test data, as defined in the preamble
phat=(prev.n+1)./(prev.N+2); % Assuming Beta(1,1)=U(0,1) prior to avoid std=0
stdphat = sqrt(phat.*(1-phat)./prev.N);
end
function [yd,yc,yi]=wacker()
[prev,phat, stdphat]=randompcr();
file='210204.dat';
data=readtable(file);
yc=[zeros(63,1); table2array(data(:, 2))]; % Cases per day
yd=[zeros(63,1); table2array(data(:, 3))]; % ICU per day
yi=[zeros(63,1); table2array(data(:, 4))]; % Deaths per day
b=1/7*ones(1,7);
ycf=filter(b,1,yc);
yif=filter(b,1,yi);
ydf=filter(b,1,yd);
ycf(1:3)=0; % Remove three days to make the MA filter symmetric
N=length(yd);
N=330; % Last day for curve fitting
N=396;
N=365
% yd(t) = b*yi(t-tau) + e
tauvec=1:14;
for tau=1:max(tauvec)
b=yif(1:N-tau)\ydf(tau+1:N);
Vb(tau)=norm(b*yif(1:N-tau)-ydf(tau+1:N))/norm(ydf(tau+1:N));
end
[Vbmin,taubhat]=min(Vb);
bhat=yif(1:N-taubhat)\ydf(taubhat+1:N);
stdbhat=norm(bhat*yif(1:N-taubhat)-ydf(taubhat+1:N))/norm(yif(1:N-taubhat))/sqrt(N-taubhat);
ind=[max(taubhat-3,1):taubhat+3];
figure(11)
[taubhatcont,stdtaubhat,Vbhat,lossb]=DLS(Vb,tauvec,ind);
title('D(t)=b * IVA(t-Delay)')
% yd(t) = c*yc(t-tauc) + e
tauvec=1:20;
Nstart=163; % First day for curve fitting cases
for tau=1:max(tauvec)
c=ycf(Nstart+1:N-tau)\ydf(Nstart+tau+1:N);
Vc(tau)=norm(c*ycf(Nstart+1:N-tau)-ydf(Nstart+tau+1:N))/norm(ycf(Nstart+tau+1:N));
end
[Vcmin,tauchat]=min(Vc);
chat=ycf(Nstart+1:N-tauchat)\ydf(Nstart+tauchat+1:N);
stdchat=norm(chat*ycf(Nstart+1:N-tauchat)-ydf(Nstart+tauchat+1:N))/norm(ycf(Nstart+1:N-tauchat))/sqrt(N-Nstart-tauchat);
ind=[max(tauchat-3,1):tauchat+3];
figure(12)
[tauchatcont,stdtauchat,Vchat,lossc]=DLS(Vc,tauvec,ind);
title('D(t)=b * Cases(t-Delay)')
Vmat{1,1}='Relation';
Vmat{1,2}='Scale';
Vmat{1,3}='Std';
Vmat{1,4}='Delay';
Vmat{1,5}='Std';
Vmat{1,6}='$\bar{V}$';
Vmat{1,7}='LS loss';
Vmat{2,1}='ICU to Death';
Vmat{2,2}=sprintf('%2.3f',bhat);
Vmat{2,3}=sprintf('%2.3f',stdbhat);
Vmat{2,4}=sprintf('%2.3f',taubhatcont);
Vmat{2,5}=sprintf('%2.4f',stdtaubhat);
Vmat{2,6}=sprintf('%2.3f',Vbhat);
Vmat{2,7}=sprintf('%2.3f',lossb);
Vmat{3,1}='Cases to Death';
Vmat{3,2}=sprintf('%2.3f',chat);
Vmat{3,3}=sprintf('%2.3f',stdchat);
Vmat{3,4}=sprintf('%2.3f',tauchatcont);
Vmat{3,5}=sprintf('%2.4f',stdtauchat);
Vmat{3,6}=sprintf('%2.3f',Vchat);
Vmat{3,7}=sprintf('%2.3f',lossc);
Vmat
% Regression model phat = deaths * 1/IFR + stdphat * e
Npop=1e7; % Swedish population
tauvec=1:40; % Possible delays in this interval
for k=tauvec;
th=(ydf(prev.date+round(prev.days)+k)/Npop./stdphat(:).^2)\(phat(:)./stdphat(:).^2);
Vp1(k)=norm(ydf(prev.date+round(prev.days)+k)/Npop-phat(:));
end
[Vp1min,muID1]=min(Vp1);
th1hat=(ydf(prev.date+round(prev.days)+muID1)/Npop./stdphat(:).^2)\(phat(:)./stdphat(:).^2);
varth1=1./sum(ydf(prev.date+round(prev.days)+muID1)/Npop./stdphat(:).^2);
IFR1hat=1/th1hat;
IFR1hatstd=1/th1hat^2*sqrt(varth1);
ind=muID1+[-2:2];
figure(13)
[taup1hat,stdtaup1hat,Vp1hat,lossp1]=DLS(Vp1,tauvec,ind);
title('N * p(t_i)=Deaths(t_i-Delay) / IFR')
% Model 2: Deaths = IFR * phat + e
for k=tauvec;
IFR=(phat(:)*Npop)\ydf(prev.date+round(prev.days)+k);
Vp2(k)=norm(ydf(prev.date+round(prev.days)+k)-(phat(:)*1e7)*IFR);
end
[Vp2min,muID2]=min(Vp2)
IFR2hat=(phat(:)*Npop)\ydf(prev.date+round(prev.days)+muID2);
IFR2hatstd=norm(ydf(prev.date+round(prev.days)+muID2)-(phat(:)*1e7)*IFR)/sqrt((phat*phat'))/1e7
ind=muID2+[-2:2];
figure(14)
[taup2hat,stdtaup2hat,Vp2hat,lossp2]=DLS(Vp2,tauvec,ind);
title('Deaths(t_i)=N * p(t_i+Delay) * IFR')
clear Vmat
Vmat{1,1}='Start day';
Vmat{1,2}='Duration';
Vmat{1,3}='N';
Vmat{1,4}='n';
Vmat{1,5}='$\hat{p}$';
Vmat{1,6}='Std($\hat{p}$)';
for i=1:length(phat);
Vmat{i+1,1}=prev.date(i);
Vmat{i+1,2}=prev.days(i);
Vmat{i+1,3}=prev.N(i);
Vmat{i+1,4}=prev.n(i);
Vmat{i+1,5}=sprintf('%2.4f',phat(i));
Vmat{i+1,6}=sprintf('%2.4f',stdphat(i));
end
Vmat
clear Vmat
Vmat{1,1}='Method';
Vmat{1,2}='Estimated IFR';
Vmat{1,3}='Std of IFR';
Vmat{1,4}='Delay';
Vmat{1,5}='Cont. Delay';
Vmat{1,6}='Std of Delay';
Vmat{2,1}='Measured prevalence';
Vmat{2,2}=sprintf('%2.5f',IFR1hat);
Vmat{2,3}=sprintf('%2.5f',IFR1hatstd);
Vmat{2,4}=muID1;
Vmat{2,5}=sprintf('%2.5f',taup1hat);
Vmat{2,6}=sprintf('%2.5f',stdtaup1hat);
Vmat{3,1}='Measured mortality';
Vmat{3,2}=sprintf('%2.5f',IFR2hat);
Vmat{3,3}=sprintf('%2.5f',IFR2hatstd);
Vmat{3,4}=muID2;
Vmat{3,5}=sprintf('%2.5f',taup2hat);
Vmat{3,6}=sprintf('%2.5f',stdtaup2hat);
Vmat
% Final plot of everything
figure(1), plotfix
Ndays=length(ydf);
b=1/13*ones(1,13);
ycfhat=filter(b,1,ydf);
ycfhat=[ycfhat(7:end);zeros(7,1)];
h1=plot(61:Ndays,ycf(61-tauchat:Ndays-tauchat)*chat,'b-','linewidth',3);
hold on
h2=plot(61:Ndays,yif(61-taubhat:Ndays-taubhat)*bhat,'color',[0 0.5 0],'linewidth',3); % Better green;
h3=plot(61:Ndays,ydf(61:Ndays),'r-','linewidth',3);
legend([h1(1) h2(1) h3(1)],{sprintf('%2.3f * fall %d dagar tidigare',chat,tauchat),sprintf('%2.1f * IVA inlagda %d dagar tidigare',bhat,taubhat),'Döda idag'},'location','north')
xdate
prevbar=ydf(muID2+1:Ndays)/Npop/IFR2hat; % Estimated prevalence
disp(sprintf('P * estimated Itot = %d',sum(prevbar)*Npop))
disp(sprintf('sum(p) = %2.1f',sum(prevbar)))
s=prevbar(prev.date(:))\stdphat(:)
figure(2), plotfix
h1=plot(1:Ndays-muID2,prevbar*100,'r-','linewidth',3);
hold on
h5=patch([1:length(prevbar) length(prevbar):-1:1],[prevbar'*(1+s) fliplr(prevbar'*(1-s))]*100,[1 0.5 0.5],'facealpha',0.7);
for i=1:length(phat)
h5=patch(prev.date(i)+prev.days(i)*[0 1 1 0 0],max(0,phat(i)+stdphat(i)*[1 1 -1 -1 1])*100,[0.5 0.5 1],'facealpha',0.7);
hold on
end
legend([h1(1) h5(1)],{'Skattad prevalens','Uppmätt prevalens'},'location','north')
ylabel('Skattad prevalens (%)')
xdate;
end
function [tauhat,stdtauhat,Vhat,LSloss]=DLS(V,tauvec,ind)
% Discrete Least Squares
N=length(ind);
phi=[ones(size(tauvec)); tauvec; tauvec.^2]';
thhat=phi(ind,:)\V(ind)';
loss=norm(V(ind)'-phi(ind,:)*thhat)/norm(V(ind));
tauhat=-thhat(2)/thhat(3)/2; % Refined estimate
stdtauhat=sqrt(((thhat(1)-thhat(2)^2/4/thhat(3)))/N/thhat(3));
LSloss=1-thhat(2)^2/4/thhat(1)/thhat(3);
taucont=linspace(ind(1),ind(end),100); % Finer resolution in prediction
phicont=[ones(size(taucont)); taucont; taucont.^2]';
Vhat=phicont*thhat;
p