新型コロナウイルス感染と、トイレ後の石鹸手洗いの相関関係

2014年に実施された国際調査の資料を見つけた（資料は、投稿末尾にまとめたサイトからダウンロードできる）。
トイレの後に石鹸を用いて手洗いするか、についての国別統計になっている。この調査のポイントは、手洗いをするかどうかではない。手洗いするとしても、それは水洗いかもしれない。石鹸を用いて手洗いすること、がポイントだ。
最下位は中国、その次が日本だ。欧州ではラテン系のランクが低く、50位にフランス、57位にイタリアが並んでいる。

これらの国を見て連想するのがCOVID-19の感染だ。次のサイトで公開されている各国の感染者数との相関分析、クラスタ分析をしてみた。
www.arcgis.com

同時に統計局が公開している高齢者人口の割合を含めて、主成分分析も試みた。

結果として、次の傾向が見つかった。

• 分析対象国をグループ化すると、感染者数に応じて層状にグループ化される。
• 各グループは石鹸手洗い、高齢者人口割合の高低を含んでおり、
  • 石鹸手洗いの割合が低いほど、感染者数が多い傾向がある。
    • しかし日本は例外的に、若干低かった。本格的にPCR検査が実施されれば、同じ傾向を示すのではないだろうか。
  • 高齢者人口割合が高いほど、感染者数が低い傾向がある。

統計処理にはRを用いた。コードの全文は、投稿末尾に掲載すると同時に、GitHubにもアップロードしている。用いたデータも同様にアップロードしている。
ちなみに文章の調子はいたって真面目なのだが、本人は学校の課題感覚、遊び半分で分析している。あまり真に受けないように。

データ収集、整理

各サイトから、次のデータをまとめたExcelファイルを作成する。それぞれのサイトは、投稿末尾にまとめたハイパーリンクから参照してほしい。

石鹸手洗い割合	BVA GroupサイトからダウンロードできるPDFファイルの最終ページから引用
感染者数	CSSEサイト左側に一覧されている感染者数から引用 2020年3月16日のデータ
高齢者人口割合	総務省統計局の情報から引用

手洗いデータの国名はフランス語表記なので、英語表記の感染者数データと合わせる。それに総務省統計局のデータを付け足す。結果として出来上がったのが、data.xlsxだ。
同ファイル中のシート「data」に、次の対象を網羅したデータが記録されている。このデータを分析に用いる。

手洗い、感染者数	対象59か国
手洗い、感染者数、高齢者人口割合	対象14か国

感染者数と、トイレ後の石鹸手洗いとの相関

Excelファイルから対象データを読み込む。
まず感染者数と、石鹸を用いた手洗い割合のデータをプロットする。中国とイタリアの感染者数が外れ値的に突出していることが分かる。
相関を検定すると、相関係数はおよそ-0.49と算出された。負の相関関係があることが分かる。95%の信頼区間でP値は非常に低く、統計的に有意な結果だ。

	Pearson's product-moment correlation

data:  x and y
t = -4.3165, df = 57, p-value = 6.389e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6675695 -0.2752390
sample estimates:
       cor 
-0.4963375 

両データを対数変換してプロットすると、おぼろげな傾向がプロットされる。この場合の相関係数は-0.50、こちらも統計的に優位な結果だ。

	Pearson's product-moment correlation

data:  x and y
t = -4.3876, df = 57, p-value = 5.01e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.6720685 -0.2827639
sample estimates:
       cor 
-0.5024625 

ここでは具体的に、トイレ後の石鹸を用いた手洗い、についてのデータなのだが、実際のところ衛生観念や潔癖さに通じる要因で、それが反映された結果なのだと、私は思う。
トイレ後に石鹸を用いて手洗いする習慣のある国では、感染者数は少ない傾向にある。
勘違いしないでほしいのは、これは因果関係ではないということだ。2014年のデータと、感染者数を突き合せたら、そのような関係が見出せた、というだけの話だ。
トレイ後に石鹸を用いて手洗いしているから感染しない、といった因果関係ではない。

感染者数と、トイレ後の石鹸手洗いとのクラスタ分析

対象59か国、20か国の、いずれの場合でも共通して確認できることは、トイレ後の石鹸手洗いの割合が高い国では、その割合が少ない国に比べて感染者数が少ない傾向が現れた。

対象59か国の場合

まず59か国全てを対象にクラスタ分析してみる。次のクラスタ数が最適であると提案された。

階層的分類による多数決	3
非階層的分類による多数決	3, 14

クラスタ数が3の場合には、クラスタ描画では良好に分類されているように見えるのだが、シルエット描画において、数か国に誤分類の可能性が示唆されている。
ちなみに14クラスタでの分類は、かなり入り交じった様相で全貌がつかみにくい。そこで対象国を絞って、再度クラスタ分析を実施することにした。

対象20か国の場合

対象国を感染者数上位20か国に絞り、同様のクラスタ分析を実施する。

階層的分類による多数決	4
非階層的分類による多数決	15

提案されたクラスタ数の内、対象国20に対して15分類は細かすぎるため、4クラスタに分類することにした。
59か国を対象とした分析では分かりにくかった特徴が明らかになっている。感染者数がそれほど変わらない国が層状にグループ化されている。そして、それぞれのグループには石鹸手洗いの割合が高い国、低い国が混ざっている。

総じて、全体的に右肩下がりの傾向がある。各グループ内の感染者数を比較したとき、石鹸手洗いの割合が少ない国の感染者は、割合の高い国に比べて多くなる傾向が現れている。
例外的なのが日本を含むグループだ。日本の感染者数は、グループ内の他国に比べて若干低位にある。もし日本で本格的にPCR検査が実施されたならば、他のグループと同じ傾向を示すことになるのではないだろうか。

感染者数、トイレ後の石鹸手洗い、高齢者人口割合の主成分分析

高齢者人口割合のデータを持つ14か国を対象に、主成分分析をする。これまでの分析同様、全てのデータを対数化して分析を実施した。

       eigenvalue percentage of variance cumulative percentage of variance
comp 1  1.9819350               66.06450                          66.06450
comp 2  0.6400748               21.33583                          87.40033
comp 3  0.3779901               12.59967                         100.00000

まず第1主成分、第2主成分で全体の87%を説明できることが分かる。各成分ベクトルをプロットする。
横軸（第1主成分）の左側ほど感染者数が多く、右側ほど石鹸手洗い、高齢者の割合が高くなる傾向が分かる。同時に縦軸（第2主成分）では上側ほど高齢者の割合が高く、下側ほど石鹸手洗いの割合が高くなることが分かる。
また、ベクトルの配置から高齢者人口割合と感染者数との相関は、石鹸手洗いの割合との場合よりも低いことが分かる。ベクトルの角度が直角に近いため、高齢者人口割合と、石鹸手洗いの割合との相関はほぼないに等しい。

ここに対象国をプロットする。左側の中国を先頭として、その他の国々が感染者数の多い順に右へ並ぶのだが、ある国は高齢者人口割合の高い方、またある国は石鹸手洗いの割合が高いほうへ、白葉るように並んでいるのが分かる。どうやら、高齢者人口割合も、石鹸手洗いの割合と同じ傾向を示すかもしれない。

ここで感染者数と高齢者人口の割合のみで、クラスタ分析をしてみる。グループが層状に形成され、それぞれのグループが高齢者人口割合の高い国、低い国を含んでいる。そして、やはり石鹸手洗いの割合の場合と同じような、右肩下がりの傾向が現れる。
意外なことに高齢者人口割合の高い国は、グループ内の他国に比べて感染者数が低くなる傾向が現れた。

参照

www.bva-group.com
www.arcgis.com
www.stat.go.jp
github.com

コード全文

🔎R code

library(readxl)
library(NbClust)
library(factoextra)
library(Rmisc)
library(FactoMineR)
library(corrplot)

#summary and plot
myplot = function(x, y){
  print(summary(x))
  print(summary(y))
  print(cor.test(x, y))
  plot(x, y, xlab = "hand wash", ylab = "infection case")
}

#lead appropriate cluster number
mynb = function(mymx){
  myAHCnum = NbClust(mymx, method = "ward.D", index = "all")
  myNHCnum = NbClust(mymx, method = "kmeans", index = "alllong")
  
  fig1 = fviz_nbclust(myAHCnum, method = "silhouette")
  fig2 = fviz_nbclust(myNHCnum, method = "gap_stat", nboot = 100)

  multiplot(fig1, fig2, cols = 2)
}

#plot cluster
mycl = function(num_cls, dt_cls){
  myAHC = hcut(dt_cls, k = num_cls, stand = TRUE, graph = FALSE)
  myNHC = kmeans(dt_cls, num_cls, iter.max = 100, nstart = nrow(dt_cls))

  fig3 = fviz_silhouette(myAHC, label = TRUE, rotate = TRUE, print.summary = FALSE)
  fig4 = fviz_cluster(myNHC, data = dt_cls, repel = TRUE)
  
  multiplot(fig3)
  multiplot(fig4)
}

#original data
exlData = read_excel("data.xlsx", sheet = "data", 
    col_types = c("text", "numeric", "numeric", "numeric"))
print(exlData)

#correlation
mydata = exlData[,c(2, 3)]
mydata$washlog = log10(mydata$Handwash)
mydata$caselog = log10(mydata$Infection)

myplot(mydata$Handwash, mydata$Infection)
myplot(mydata$washlog, mydata$caselog)

#clustering with 59 countries
clsData60 = as.matrix(mydata[,c(3, 4)])
rownames(clsData60) = exlData$English
mynb(clsData60)

for (i in c(3, 14)){
  mycl(i, clsData60)
}

#clustering with 20 countries
clsData20 = clsData60[order(clsData60[, 2], decreasing = TRUE),]
clsData20 = clsData20[c(1:20),]
mynb(clsData20)

for (i in c(4)){
  mycl(i, clsData20)
}

#data for aged population ratio
df_aged = na.omit(exlData)
df_aged$washlog = log10(df_aged$Handwash)
df_aged$caselog = log10(df_aged$Infection)
df_aged$agedlog = log10(df_aged$Aged)

pcaData = as.matrix(df_aged[,c(5:7)])
row.names(pcaData) = df_aged$English[1:nrow(pcaData)]

#PCA
mypca = PCA(pcaData, scale.unit = TRUE, graph = FALSE)
mypca$eig

fig10 = fviz_pca_var(mypca, col.var = "cos2", gradient.cols = c("blue", "red"), repel = TRUE)
fig11 = fviz_pca_var(mypca, col.var = "contrib", gradient.cols = c("blue", "red"), repel = TRUE)
multiplot(fig10, fig11, cols = 2)

corrplot(mypca$var$cos2, addCoef.col = "gray")
corrplot(mypca$var$contrib, is.corr = FALSE, addCoef.col = "gray")

fviz_pca_biplot(
  mypca,
  col.var = "blue",
  col.ind = "cos2",
  gradient.cols = c("blue", "red"),
  repel = TRUE)

for (i in c(2, 3)){
  mycl(i, agedData[,c(3, 2)])
}