MATLAB ユーザーコミュニティー

MATLAB & Simulink ユーザーコミュニティー向け日本語ブログ

韓国を事例とする新型コロナウィルス感染のデータ解析

こんにちは、道家です。緊急事態宣言が解除されましたが、まだ安心して外出できていません。 皆さん、どうぞ安全にお過ごしください。

今日は MathWorks の本社で働いている竹内さんからの COVID-19 データ解析に関する ゲスト投稿です。

その前に、 File Exchangeに どの位 COVID-19 に関するファイル がアップされているかみてみましょう。

可視化やアニメーションからデータ解析やシミュレーションまで、結構ありますね。 Authorを見る限り海外のユーザーの方が殆どですが、日本のユーザーの方も見かけます。 (すみません、File Exchange上の名前からしか判断できませんが)

  • covidx by Hisa. こちらは、様々な日本の都市での感染者数の累計と倍加時間を可視化するエントリーです。

興味ある方は是非試してみてください。

それでは、竹内さんの記事です。最初の自己紹介にも書かれていますが、竹内さんは まったくの MATLAB 初心者でしたが、独学でみるみるうちに MATLAB を使いこなせるように なり、マーケティングデータなど様々なデータの解析をされるようになりました。 今では社内でちょっとした有名人になっています。今回はこのような形でゲスト投稿して 頂き、とても嬉しいです。


初めまして、MathWorksの竹内俊明(たけうちとしあき)です。米国採用で、マーケティングに 所属しています。これまで、 Loren on the Art of MATLAB には何度かゲストブロガーをさせて頂いているのですが、 今回は、初めてMATLABユーザーコミュニティーのブログにお邪魔することになりました。 よろしくお願い致します。もともと文系人間で、MATLABは入社する前は全く知らなかった のですが、入社した後、会社のコンピューターに入れて興味本位でいじっているうちにマハって しまいました。

Contents

新型コロナウイルス対策の展望

さて、新型コロナウイルス対策として外出制限が各国で始まって二か月以上経ち、その緩和や 解除がそれぞれの国で始まっています。他国に比べてかなり出遅れている米国から見ると、 東アジア諸国は非常に優等生に見えます。その中でも、韓国は、MERS(中東呼吸器症候群)の 反省から、世界的にも珍しく、かなり早い段階から大規模な検査を徹底したことが知られて います。このため、韓国のデータは精度が高いと期待できます。なので、そのデータがどこかに 公開されていたら、そこから何か学ぶことがあるんじゃないかなと思いつきました。

今回は GitHub 上で KCDC(韓国疾病管理本部)から提供された同国の新型コロナウィルスに 関するデータが公開されているのを見つけました。

上述のページでは、様々なデータが提供されていますが、手始めに、以下の三つのファイルから 始めてみます。

patient = readtable("southKorea\PatientInfo.csv");
route = readtable("southKorea\PatientRoute.csv");
region = readtable("southKorea\Region.csv");

新型コロナウイルス感染者数分布

まずは、 geodensityplot を使って、おおまかな感染者の分布状況を見てみましょう。 R2020a では、ライブスクリプト上のプロット更新速度が速くなっているので、動画として サクサク動きます。動画で見たい方は、スクリプトをR2020a上でライブスクリプトとして 開いて、コメントアウトされたコードを実行してみてください。それが出来ない場合は、 こちらに 動画 があります。

patientR = innerjoin(patient,region,"Keys",{'province','city'});
t = patientR.confirmed_date;
t.Format = "yyyy年MM月dd日";
lat = patientR.latitude;
lon = patientR.longitude;
% dateRange = min(t):max(t);
% cdate = patientR.confirmed_date;
% lookback = days(7);
figure
colormap hot
alphamap(normalize((1:64).^0.2,'range'))
% for ii = 1:length(dateRange)
%     lat2date = lat(cdate <= dateRange(ii) & cdate > ...
%         (dateRange(ii) - lookback));
%     lon2date = lon(cdate <= dateRange(ii) & cdate > ...
%         (dateRange(ii) - lookback));
%     geodensityplot(lat2date,lon2date,"Radius",5*10^4,"FaceColor","interp");
    geodensityplot(lat,lon,"Radius",5*10^4,"FaceColor","interp");
    geolimits([33.3710, 38.3250],[123.9779, 131.7367]);
    gx = gca;
    gx.LatitudeLabel.String = "緯度";
    gx.LongitudeLabel.String = "経度";
    title(compose("韓国の新型コロナウイルス累計感染者数分布 %s 時点",max(t)))
%     title(compose("韓国の新型コロナウイルス感染者数分布 %s 時点",dateRange(ii)))
%     pause(0.1)
%     drawnow;
% end

時系列で見ていくと、こんな感じです。

  • 一月末にソウル近辺で最初の感染例が出たが、二月中旬までは、拡大は比較的緩やかだった
  • 二月中旬に大邱(テグ)で大規模な集団感染が発生した後、急速に全国に拡大した
  • しかし、三月末までにはかなり沈静化していき、四月中旬以降は大きな変化はなかった

動画を見て感じたのは、やはり大規模な集団感染が発生すると、一挙に感染が広がるという ことです。外出制限の緩和も、集団感染の防止を考えて実施しないといけませんね。三月末に MathWorks上海オフィスの同僚とビデオミーティングした時に、「外出制限は解除されたけど、 やはり自主的に外出は控えている」と言っていましたが、正しい選択ですね。

感染経路の分類

次に、感染例のうち、どの感染経路が多いのかヒストグラムにして見てみましょう。

figure
histogram(categorical(patient.infection_case), ...
    "Orientation","horizontal","DisplayOrder","ascend")
xlabel("感染者数")
title("韓国の新型コロナウイルス感染経路分類")

感染経路として海外からの入国者、国内感染者の二つが最も大きいことが分かります。 それ以外は、etc(その他)を除くと、集団感染ですね。集団感染がどこで発生したかと見ると、 教会、コールセンター、ナイトクラブ、カラオケ、介護施設、ジムなどが目立ちます。共通して いるのは、以下の点です。

  • 人が集まる屋内施設であること
  • 話をしたり歌ったり、運動したりと、活動内容が肺活量を増やす傾向にあること
  • マスク着用に不向きな活動内容であること
  • 長時間、同じ場所にいること

屋内では、時間が経過するにつれてウィルスの濃度が上がるので、感染する可能性が高まりますが、 コールセンターのように、マスクの着用ができない場合、他の職場に比べてさらにリスクが 高まるのではないでしょうか。

では、数の多い感染経路を選んでの推移を時系列で見てみます。

gs1 = groupsummary(patient,{'confirmed_date'});
gs2 = groupsummary(patient,{'confirmed_date','infection_case'});
gs2.infection_case = categorical(gs2.infection_case);
gs2 = renamevars(gs2,{'confirmed_date','infection_case','GroupCount'}, ...
    {'Date','Case','Count'});
figure
cases = ["Total","overseas inflow","Guro-gu Call Center", ...
    "Shincheonji Church","Itaewon Clubs"];
plot(gs1.confirmed_date,gs1.GroupCount)
hold on
plot(gs2.Date(gs2.Case == cases(2)),gs2.Count(gs2.Case == cases(2)), ...
    "LineWidth",1.5)
plot(gs2.Date(gs2.Case == cases(3)),gs2.Count(gs2.Case == cases(3)), ...
    "LineWidth",1.5)
plot(gs2.Date(gs2.Case == cases(4)),gs2.Count(gs2.Case == cases(4)), ...
    "LineWidth",1.5)
plot(gs2.Date(gs2.Case == cases(5)),gs2.Count(gs2.Case == cases(5)), ...
    "LineWidth",1.5)
legend("総感染者数","海外からの入国者","九老区のコールセンター", ...
    "大邱の新天地イエス教会","梨泰院のクラブ","Location","northeast")
xtickformat("MM月")
ylabel("感染者数")
title("韓国の新型コロナウイルス感染経路別推移")

動画で見た通り、ここでも大邱(テグ)で発生した新天地イエス 教会での集団感染以降、 国内感染が急に悪化したことも見て取れます。

興味深いのは、三月中旬以降、国内の感染例が減っているのに、海外からの入国者の感染例が 急に増えていることです。感染が欧州や米国など全世界に広まった 結果、かえって海外経路の リスクが急に高まったのでしょう。しかし、現時点では、入国制限の実施により沈静化している ようです。水際作戦、やはり大事ですね。

感染者ごとの被感染者数

データには、感染源がIDで明記されているケースがあるので、それを集計することにより、 感染者ごとの被感染者数を得ることができます。それが明記されていない場合は、感染経路が 確認できななったか、国内感染ではなかったと思われます。なお、これは感染者単位の集計です ので、集団感染の場合でも、一人あたりの被感染者数が低ければ、プロット上は大きな数字に なりません。

infEdge = patient(~isnan(patient.infected_by),{'infected_by','patient_id'});
infEdge = renamevars(infEdge,{'infected_by','patient_id'},{'Src','Tgt'});
infectors = groupsummary(infEdge,"Src");
infectors = innerjoin(infectors,patient,"LeftKeys","Src","RightKeys","patient_id");
infectors = renamevars(infectors,"GroupCount","Count");
infectors = sortrows(infectors,"Count","descend");
figure
scatter(infectors.Src,infectors.Count)
text(infectors.Src([1:3,5,7]),infectors.Count([1:3,5,7]), ...
    infectors.infection_case([1:3,5,7]))
xlabel("感染者ID")
ylabel("被感染者数")
title("感染者ごとの被感染者数")

接触追跡の対象になった人のうち、感染してなかった人数はこれに含まれていませんので、 接触追跡そのものの規模はもっと大きかったと思われます。幸いに、ほとんどの場合、一人 当たりの被感染者数は数名です。韓国の接触追跡体制が概ねよく機能していることがうかがわれます。 しかし、二十名以上に感染を広めてしまっている事例も十件ほどあります。こうしたケースは、 スーパースプレッダーと呼ばれます。

誰がスーパースプレッダーになるのか

データから、感染源―被感染者の関係が分かるので、その情報を使って、最も一人当たりの 被感染者数が多かった集団感染三件を、 graph を使ってネットワークとして可視化してみましょう。ちなみに、 findInfectedContactsは、ページの最後に定義されているローカル関数です。

[g,~] = findInfectedContacts(infEdge,infectors.Src);
[bins,binsizes] = conncomp(g);
smlBins = find(binsizes < 40);
g = rmnode(g,g.Nodes.Name(ismember(bins,smlBins)));
figure
h = plot(g,"Layout","force");
spreaders = infectors.Src(1:10);
infSrc = infectors.infected_by(1:10);
infSrc(ismissing(infSrc)) = [];
infType = infectors(ismember(infectors.Src,spreaders(~ismissing(infSrc))), ...
    {'Src','infection_case'});
labelnode(h,string(infType.Src),["感染者との接触","天安市内ジム施設", ...
    "感染者との接触"])
highlight(h,g.Nodes.Name(ismember(g.Nodes.Name,string(spreaders))), ...
    "NodeColor","red")
labelnode(h,g.Nodes.Name(ismember(g.Nodes.Name,string(infSrc))),"発端症例")
title(["スーパースプレッダーが関与する";"発端症例と被感染者のネットワーク"])
xlabel("スーパースプレッダーは赤で表示")

これで見ると、スーパースプレッダーが発端症例であるケース(右下の事例)以外に、 二次感染からスーパースプレッダーになるケース(左上の事例)や、発端症例と被感染者の両方が スーパースプレッダーになる事例(左下の事例)があることが分かります。幸いにも、三次感染が 少ないのは、それが発生するまでの時点に情報が伝達されて、接触追跡により濃厚接触者への 対応が追い付いているためでしょう。

さらに、データには、性別と年齢の情報が含まれていますので、スーパースプレッダーの属性を 見てみましょう。

tiledlayout(2,2);
nexttile
pie(categorical(infectors.sex(1:10)))
title("スーパースプレッダーの性別")
nexttile
pie(categorical(infectors.age(1:10)))
title("スーパースプレッダーの年齢")
nexttile
pie(categorical(patient.sex))
title("感染者全体の性別")
nexttile
pie(categorical(patient.age))
title("感染者全体の年齢")

感染者全体では、男女比はやや女性が多めですが、それほど差はありません。また、感染者の 年齢も五十代以下が七割以上です。しかし、スーパースプレッダーは、圧倒的に女性、また 世代も四十代以上の中高年層に集中しています。女性は長寿ですので、介護施設などで 集団感染が発生すると、中高年層の女性がスーパースプレッダーになりやすいのかも知れません。

スーパースプレッダーの動きと被感染者との接触

スーパースプレッダーの移動経路データと、被感染者の移動経路を重ね合わせ、重複する地点を 特定すれば、どこで集団感染が発生したか推定できそうです。ネットワークとして可視化した 集団感染のうち、右下の天安(チョナン)市内ジム施設で発生した集団感染を例として詳細を geoscattergeoplot を使って見ていきましょう。ちなみに、この事例のスーパースプレッダーは、 四十から五十代の女性でした。また、被感染者も全員、二十代から六十代までの女性でした。

spreaders2 = str2double(g.Nodes.Name(ismember(g.Nodes.Name, ...
    string(spreaders))));
spreaders2 = intersect(spreaders2,route.patient_id);
indexCase = route(route.patient_id == spreaders2(2),:);
indexCaseGS = groupsummary(indexCase,{'latitude','longitude','type'});
infected = str2double(neighbors(g,string(spreaders2(2))));
infectedCase = route(ismember(route.patient_id,infected),:);
infectedCaseGS = groupsummary(infectedCase,{'latitude','longitude','type'});
overlap = intersect(indexCaseGS(:,1:end-1),infectedCaseGS(:,1:end-1));
figure
geoscatter(indexCaseGS.latitude,indexCaseGS.longitude, ...
    indexCaseGS.GroupCount*10,"r","filled")
hold on
geoscatter(infectedCaseGS.latitude,infectedCaseGS.longitude, ...
    infectedCaseGS.GroupCount*10,"m")
geoplot(indexCaseGS.latitude,indexCaseGS.longitude,"r-","LineWidth",1.5)
geoplot(infectedCase.latitude,infectedCase.longitude,"m:")
legend("スーパースプレッダー","被感染者","Location","northeast")
text(37.5665,126.9780,"ソウル")
text(36.8151,127.1139,"天安(チョナン)市","FontWeight","bold")
text(35.1796,129.0756,"釜山(プサン)")
title([compose("スーパースプレッダー %s の動きと被感染者", ...
    string(spreaders2(2)));"チョナン市内ジム施設の集団感染"])
gx = gca;
gx.LatitudeLabel.String = "緯度";
gx.LongitudeLabel.String = "経度";

この集団感染の場合、スーパースプレッダーはソウル近郊まで足を延ばしているいるものの、 被感染者の動きはほぼ天安(チョナン)市に集中しています。しかし、一部の被感染者は ソウルや釜山(プサン)ンなどの他の都市に移動しています。感染経路はジム施設と記されて いるけど、全国から人が集まるジムなんでしょうかね。やはり関係者の動きを見ても、感染が 発生した場所は天安(チョナン)市以外にはありえなさそうです。天安(チョナン)市付近を 拡大してみましょう。

figure
geoscatter(indexCaseGS.latitude,indexCaseGS.longitude, ...
    indexCaseGS.GroupCount*10,"r","filled")
hold on
geoscatter(infectedCaseGS.latitude,infectedCaseGS.longitude, ...
    infectedCaseGS.GroupCount*10,"m")
geoscatter(overlap.latitude,overlap.longitude,500,"LineWidth",2)
geoplot(indexCaseGS.latitude,indexCaseGS.longitude,"r-")
geoplot(infectedCase.latitude,infectedCase.longitude,"m:")
text(36.8088093,127.1063,"フィットネスセンター")
legend("スーパー・スプレッダー","被感染者","重複する移動経路", ...
    "Location","northwest")
text(overlap.latitude,overlap.longitude,overlap.type,"Interpreter","none")
text(36.817,127.12,"天安(チョナン)市","FontWeight","bold")
gx = gca;
gx.LatitudeLabel.String = "緯度";
gx.LongitudeLabel.String = "経度";
gx.MapCenter = [36.8157 127.125];
gx.ZoomLevel = 14;
title([compose("スーパー・スプレッダー %s の動きと感染経路", ...
    string(spreaders2(2)));"天安(チョナン)市付近を拡大"])

天安(チョナン)市を拡大してみると、病院あるいは「etc」(その他)としか表記されて いない場所で スーパースプレッダーと被感染者が接触している可能性があることがわかります。 地元の患者が発症後に近辺の病院に行くのは当然の経路となるので、「etc」(その他)の うちの一つが集団感染が発生したジムのようです。 いろいろと調べてみると 、天安(チョナン)市で開催されたズンバのワークショップで 集団感染が発生したのです。Google Mapsで見ると、該当付近にフィットネスセンターが ありました。どうやら、講師や参加者が他の都市からも集まっていたため、かなり大規模な 集団感染に発展してしまったようです。ズンバは女性に人気があるので、感染者が全員、 四十代を中心とした女性であるのも納得できます。

梨泰院(イテウォン)のナイトクラブでの集団感染

5月上旬は、韓国も日本と同様に連休があり、新型コロナウイルスの感染が沈静化してたことも あって、同時期に外出自粛要請も解除されました。その後間もななく、ソウルの六本木に 相当する梨泰院(イテウォン)で、また集団感染が発生してしまいました。発生後、まだ日が 浅いため、発端症例がどの感染者なのかも含め、あまり詳細なデータがありません。

itae = patientR(patientR.infection_case == "Itaewon Clubs",:);
g = findgroups(itae.latitude,itae.longitude);
cdate = splitapply(@min,itae.confirmed_date,g);
itae.confirmed_date = cdate(g);
itaeGS = groupsummary(itae,{'latitude','longitude','confirmed_date'});
itaeGS = sortrows(itaeGS,"confirmed_date");
itaeGS.confirmed_date.Format = "yyyy年MM月dd日";
itaeR = route(ismember(route.patient_id,itae.patient_id),:);
figure
geoscatter(itaeGS.latitude,itaeGS.longitude,itaeGS.GroupCount*30,"magenta")
hold on
geoplot(itaeR.latitude,itaeR.longitude,"m:")
text(37.54,126.9921,"梨泰院(イテウォン)","FontWeight","bold")
text(37.5324,126.9904,"病院","HorizontalAlignment","right")
text(37.5336,126.9958,"バー")
text(37.5746,127.0402,"病院")
text(37.573,126.9794,"ネットカフェ")
text(37.5172,127.0473,"売店")
text(37.567,127.03,"教会")
text(37.5611,127.0355,"公共交通機関")
geolimits([37.4988, 37.5812],[126.9267, 127.0575]);
gx = gca;
gx.LatitudeLabel.String = "緯度";
gx.LongitudeLabel.String = "経度";
title("ソウル市内梨泰院(イテウォン)のナイトクラブでの集団感染")

感染者の属性を見てみると、ナイトクラブで発生した集団感染なので、当然ながら、全員が 三十代未満の若年層でした。また圧倒的に男性の割合が大きいのがこの集団感染の特徴に なります。

やはり、人が密集する屋内での活動は避けるべきですね。

tiledlayout(1,2);
t = tiledlayout(1,2);
nexttile
pie(categorical(itae.sex))
title("スーパースプレッダーの性別")
nexttile
pie(categorical(itae.age))
title("スーパースプレッダーの年齢")

おわりに

今回は、韓国のデータを使って、外出制限の緩和や解除になった場合に、どのようなリスクが あるのか、考えてみました。提供されているデータの一部しか使っていないので、まだまだ 面白い知見があるかもしれません。興味があれば、ぜひ、データ解析にチャレンジしてみて ください。参考になれば幸いです。

ローカル関数

function [graphObj,target] = findInfectedContacts(data,source,graphObj)
    srcIds = cell(length(source),1);
    tgtIds = cell(length(source),1);

    for ii = 1:length(source)
        tgtIds{ii} = data.Tgt(ismember(data.Src,source(ii)));
        srcIds{ii} = repmat(source(ii),size(tgtIds{ii},1),1);
    end
    source = vertcat(srcIds{:});
    target = vertcat(tgtIds{:});
    if ~ exist("graphObj","var")
        graphObj = graph(string(source),string(target));
    else
        graphObj = addedge(graphObj,string(source),string(target));
        if ismultigraph(graphObj)
            graphObj = simplify(graphObj);
        end
    end
end




Published with MATLAB® R2020a

|

Comments

To leave a comment, please click here to sign in to your MathWorks Account or create a new one.