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

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

その前に,文件交换に どの位COVID-19に関するファイルがアップされているかみてみましょう。

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

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

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

---

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

内容

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

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

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

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

患者=可阅读(“韩国\ PatientInfo.csv”)；路由=可读(“韩国\ PatientRoute.csv”）;地区= readtable (“韩国\地区.csv”）;

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

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

patientR = innerjoin(病人、地区“钥匙”，{“省”，“城市”}）;t = patientR.confirmed_date;t.Format =“yyyy年MM月dd日”；纬度=纬度；经度=经度；% daterrange = min(t):max(t);% cdate = patienter .confirmed_date;%回顾=天（7）；数字彩色地图热alphamap(正常化((1:64)。^ 0.2,“范围”））% for ii = 1:length(daterrange)% lat2date = lat(cdate <= dateRange(ii) & cdate >…% (dateRange(ii) - lookback));%lon2date=lon（cdate<=dateRange（ii）&cdate>。。。% (dateRange(ii) - lookback));%geodensityplot（lat2date、lon2date、“半径”、5*10^4、“面色”、“内色”）；地球密度图（纬度、经度、，“半径”5 * 10 ^ 4,“FaceColor”，“插值函数”）;geolimits ([33.3710, 38.3250], [123.9779, 131.7367]);甘氨胆酸gx =;gx.LatitudeLabel.String ="緯度"；gx.LongitudeLabel.String="経度"；标题(组成(“韓国の新型コロナウイルス累計感染者数分布% s時点”，最大（t）））%标题（组成）韓国の新型コロナウイルス感染者数分布 %s時点“，日期范围（ii）））%暂停(0.1)% drawnow;%结束

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

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

感染経路の分類

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

图直方图(分类(patient.infection_case),.．.“方向”，“水平”，“DisplayOrder”，“上升”)包含("感染者数")头衔(“韓国の新型コロナウイルス感染経路分類”）

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

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

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

gs1 = groupsummary(病人,{“confirmed_date”}）;gs2 = groupsummary(病人,{“confirmed_date”，“infection_case”}）;gs2。infection_case =分类(gs2.infection_case);gs2 = renamevars (gs2, {“confirmed_date”，“infection_case”，“GroupCount”},.．.｛“日期”，“案例”，“数”}）;数字案例= [“总”，“海外流入”，“Guro-gu呼叫中心”，.．.“Shincheonji教堂”，梨泰院“俱乐部”];情节(gs1.confirmed_date gs1.GroupCount)在绘图（gs2.Date（gs2.Case==cases（2）），gs2.Count（gs2.Case==cases（2）），.．.“线宽”1.5)情节(gs2.Date (gs2。= =例(3)),gs2。Count(gs2.Case == cases(3)),.．.“线宽”，1.5）绘图（gs2.Date（gs2.Case==cases（4）），gs2.Count（gs2.Case==cases（4）），.．.“线宽”，1.5）绘图（gs2.Date（gs2.Case==cases（5）），gs2.Count（gs2.Case==cases（5）），.．.“线宽”，1.5）图例("総感染者数"，"海外からの入国者"，"九老区のコールセンター"，.．."大邱の新天地イエス教会"，“梨泰院のクラブ”，“位置”，“东北”)xtickformat(“嗯哪月”) ylabel ("感染者数")头衔(“韓国の新型コロナウイルス感染経路別推移”）

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

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

感染者ごとの被感染者数

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

infEdge=患者（~isnan（患者感染者）{“infected_by”，“patient_id”}）;infEdge = renamevars (infEdge, {“infected_by”，“patient_id”},{“Src”，“Tgt”}）;以上= groupsummary (infEdge," Src "）;以上= innerjoin(艾滋病患者、耐心“LeftKeys”，" Src "，“RightKeys”，“patient_id”）;以上= renamevars(艾滋病患者,“GroupCount”，“数”）;以上= sortrows(艾滋病患者,“数”，“下降”）;图散射(infectors.Src infectors.Count)文本(infectors.Src([1:3、5、7]),infectors.Count([1:3、5、7]),.．.infectors.infection_case([1:3、5、7]))包含(“感染者ID”) ylabel ("被感染者数")头衔(“感染者ごとの被感染者数”）

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

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

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

[g，~]=findInfectedContacts（infEdge，infectors.Src）；[bins，binsizes]=conncomp（g）；smlBins=find（binsizes<40）；g=rmnode（g，g.Nodes.Name（ismember（bins，smlBins））；图h=plot（g，“布局”，“力”)；传播者=感染者。Src（1:10）；infSrc=感染者。被感染者（1:10）；infSrc（ismissing（infSrc））=[]；infType=感染者（ismember（infectors.Src，传播者（~ismissing（infSrc））），.．.｛“Src”，“infection_case”})；labelnode（h，字符串（infType.Src）[“感染者との接触”，“天安市内ジム施設”，.．.“感染者との接触”)突出(h, g.Nodes.Name (ismember (g.Nodes.Name字符串(传播者))),.．.“NodeColor”，“红色”)labelnode（h，g.Nodes.Name（ismember（g.Nodes.Name，string（infSrc）），"発端症例")标题(“スーパースプレッダーが関与する”；“発端症例と被感染者のネットワーク”])xlabel(“スーパースプレッダーは赤で表示”）

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

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

平铺布局（2,2）；下一个小馅饼（分类（感染者。性别（1:10）））标题(“スーパースプレッダーの性別”)nexttile pie（分类（感染者、年龄（1:10）））标题(“スーパースプレッダーの年齢”)nexttile pie（分类（患者性别））标题(“感染者全体の性別”) nexttile pie(分类的(患者。年龄))title("感染者全体の年齢"）

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

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

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

spreaders2 = str2double (g.Nodes.Name (ismember (g.Nodes.Name,.．.字符串（吊具））；吊具2=相交（吊具2，route.patient_id）；indexCase=路线（route.patient_id==吊具2（2），：）；indexCaseGS=组摘要（indexCase{“纬度”，“经度”，“类型”}）;感染= str2double(邻居(g,字符串(spreaders2 (2))));infectedCase =路线(ismember (route.patient_id感染):);infectedCaseGS = groupsummary (infectedCase, {“纬度”，“经度”，“类型”}）;重叠=相交(indexCaseGS (: 1: end-1), infectedCaseGS (:, 1: end-1));图geoscatter (indexCaseGS.latitude indexCaseGS.longitude,.．.indexCaseGS。GroupCount * 10,“r”，“填充”)举行在geoscatter (infectedCaseGS.latitude infectedCaseGS.longitude,.．.infectedCaseGS。GroupCount * 10,“m”) geoplot (indexCaseGS.latitude indexCaseGS.longitude,“r -”，“线宽”，1.5）地理地块（感染病例。纬度，感染病例。经度，”男:“)传奇(“スーパースプレッダー”，"被感染者"，“位置”，“东北”)文本(37.5665,126.9780,“ソウル”)文本（36.8151127.1139，“天安(チョナン)市”，“FontWeight”，“大胆”)文本(35.1796,129.0756,“釜山(プサン)”)标题([组成(“スーパースプレッダー% sの動きと被感染者”，.．.字符串(spreaders2 (2)));"チョナン市内ジム施設の集団感染"])gx=gca；gx.LatitudeLabel.String="緯度"；gx.LongitudeLabel.String="経度"；

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

图geoscatter (indexCaseGS.latitude indexCaseGS.longitude,.．.indexCaseGS。GroupCount * 10,“r”，“填充”)举行在geoscatter (infectedCaseGS.latitude infectedCaseGS.longitude,.．.infectedCaseGS。GroupCount * 10,“m”)地球散射（重叠。纬度，重叠。经度，500，“线宽”(2) geoplot indexCaseGS.latitude indexCaseGS.longitude,“r -”)地理地块（传染病例纬度、传染病例经度、，”男:“)文本(36.8088093,127.1063,"フィットネスセンター")传奇(“スーパー・スプレッダー”，"被感染者"，“重複する移動経路”，.．.“位置”，“西北”)文本(overlap.latitude overlap.longitude overlap.type,“口译员”，“没有”)文本(36.817,127.12,“天安(チョナン)市”，“FontWeight”，“大胆”) gx = gca;gx.LatitudeLabel.String ="緯度"；gx.LongitudeLabel.String="経度"；gx。MapCenter = [36.8157 127.125];gx。ZoomLevel = 14;标题([组成(“スーパー・スプレッダー% sの動きと感染経路”，.．.字符串(spreaders2 (2)));"天安（チョナン）市付近を拡大"])

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

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

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

itae = patientR (patientR。infection_case = =梨泰院“俱乐部”，：；g=findgroups（itae.纬度，itae.经度）；cdate=splitapply（@min，itae.已确认的日期，g）；itae.已确认的日期=cdate（g）；itaeGS=groupsummary（itae{“纬度”，“经度”，“confirmed_date”})；itaeGS=索特罗斯（itaeGS，“confirmed_date”）;itaeGS.confirmed_date。格式=“yyyy年MM月dd日”；ita =路线(ismember (route.patient_id itae.patient_id):);图geoscatter (itaeGS.latitude、itaeGS.longitude itaeGS.GroupCount * 30,“红色”)举行在geoplot (itaeR.latitude itaeR.longitude,”男:“)文本（37.54126.9921，"梨泰院（イテウォン）"，“FontWeight”，“大胆”)文本(37.5324,126.9904,"病院"，“水平对齐”，“对”)文本(37.5336,126.9958,“バー”)文本(37.5746,127.0402,"病院")文本(37.573,126.9794,"ネットカフェ")文本(37.5172,127.0473,"売店")文本(37.567,127.03,"教会")文本(37.5611,127.0355,"公共交通機関")地质极限([37.4988,37.5812]，[126.9267,127.0575]);甘氨胆酸gx =;gx.LatitudeLabel.String ="緯度"；gx.LongitudeLabel.String="経度"；标题(“ソウル市内梨泰院(イテウォン)のナイトクラブでの集団感染”）

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

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

tiledlayout(1、2);t = tiledlayout(1、2);nexttile派(分类(itae.sex)标题(“スーパースプレッダーの性別”)nexttile pie（分类（itae.age））标题(“スーパースプレッダーの年齢”）

おわりに

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

ローカル関数

函数[graphhobj,target] = findInfectedContacts(data,source, graphhobj) srcid = cell(length(source)，1); / /获取数据tgtIds =细胞(长度(源),1);为ii = 1:length(source) tgtIds{ii} = data.Tgt(ismember(data.Src,source(ii)));srcIds {2} = repmat(源(ii)、大小(tgtIds {2}, 1), 1);结束source=vertcat（srcIds{：}）；target=vertcat（tgtIds{：}）；如果~ (“graphObj”，“var”) graphhobj = graph(string(source)，string(target));其他的graphObj=addedge（graphObj、字符串（源）、字符串（目标））；如果ismultigraph(graphObj) graphObj = simplify(graphObj);结束结束结束

获取MATLAB代码

发布与MATLAB®R2020a

发布与MATLAB®R2020a