名古屋市のHPVVデータの解析

[おことわり] 以下はあくまでもRを使ったデータ解析法の例示で,HPVワクチンについて何らかの結論を主張するものではありません。

データの読み込み

名古屋市のデータのCSV化で作ったCSVファイルを読み込む。

名古屋市のデータのCSV化でも書いたように,CSVに列275は含まれていない。また,列229は006001-024000の範囲がすべてNA(欠損値)になっている。これら自由記述欄以外の部分(半角数字の部分)は正しく変換されているはずである。

Nagoya-HPVV-data-01.zip を展開して kaito.csv を取り出して読む:

kaito = read.csv("kaito.csv", header=FALSE, colClasses="character", fileEncoding="UTF-8")
dim(kaito)
[1] 30793   274
head(kaito)
      V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12  V13 V14 V15 V16 V17 V18 V19  V20
1 000001  0  0  1  0  0  0  1  0   0   0   0 0000   0   0   0   0   0   2 0000
2 000002  0  1  0  0  0  0  0  0   1   0   1 0000   0   0   0   0   0   1 0000
3 000003  0  1  0  0  1  0  0  0   0   0   2 0000   2   0   1   0   0   1 0000
4 000004  0  0  1  0  0  0  0  0   0   1   1 0000   0   0   0   0   0   1 0000
5 000005  0  0  1  0  0  0  0  1   0   0   2 2610   1   1   0   0   0   2 2610
6 000006  1  0  0  0  1  0  0  0   0   0   1 0000   0   0   0   0   0   1 0000
……

同じことが data.table の fread() や readr の read_csv() でもできる:

library(data.table)
kaito_fread = fread("kaito.csv", header=FALSE, colClasses="character")
library(readr)
kaito_readr = read_csv("kaito.csv", col_names=FALSE, col_types=paste0(rep("c",274), collapse=""))

readr の read_csv() では変数名が X1,X2,… となることに注意する。

子宮頸がん予防接種(HPVV)の有無を表にしてみる:

table(kaito$V233)

    0     1     2     3 
  624  9245 20912    12 

チェックなしが624人,なし9245人,あり20912人,両方チェックあり12人である。チェックなしや重複チェックは,他の項目の回答と照らし合わせて吟味するべきかもしれないが,ここでは単に無視する。

例:「身体が自分の意思に反して動く」

症状はいろいろある。ここではHPVV副作用としか考えられない症状(実際には思春期によく見られる)として「身体が自分の意思に反して動く」(V145)を例として調べる。なお,これはあくまでもデータ解析法の検討のための例示で,これだけから何らかの結論を導き出そうとしているわけではない。

table(kaito$V145)

    0     1     2 
  259 30272   262 

チェックなし259人,なし30272人,あり262人である。これをHPVV有無とクロス集計すれば

table(kaito$V233, kaito$V145)

        0     1     2
  0    36   582     6
  1    80  9107    58
  2   143 20571   198
  3     0    12     0

つまり次のようになる:

症状なし症状あり
HPVVなし910758
HPVVあり20571198

Fisherの正確検定をしてみる:

fisher.test(matrix(c(9107,20571,58,198),nrow=2))

	Fisher's Exact Test for Count Data

data:  matrix(c(9107, 20571, 58, 198), nrow = 2)
p-value = 0.005159
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 1.121378 2.063976
sample estimates:
odds ratio 
  1.511313 

非常に有意な関連がありそうである。しかし,他の変数と交絡しているかもしれないので,慎重に調べる必要がある。

まずは年齢が交絡しているかどうか調べる。まず年齢が必要である。この調査では列5〜11にそれぞれ平成6〜12年度生まれをチェックさせている。

table(rowSums(kaito[,5:11] == "1"))

    0     1     2 
  456 30334     3 

チェックなしが456人,二つの年度にわたってチェックしているのが3人いる。これらを除いて年齢を調べる(もっとエレガントな方法があるだろうが,面倒なので…):

birth = ifelse(rowSums(kaito[,5:11]=="1") != 1, NA,
  ifelse(kaito[,5]=="1", 6,
  ifelse(kaito[,6]=="1", 7,
  ifelse(kaito[,7]=="1", 8,
  ifelse(kaito[,8]=="1", 9,
  ifelse(kaito[,9]=="1", 10,
  ifelse(kaito[,10]=="1", 11,
  ifelse(kaito[,11]=="1", 12, NA))))))))

よりエレガントな方法を教えていただいた

birth <- as.numeric(apply(kaito[ ,5:11], 1, paste, collapse=""))
birth <- match(birth, 10^(6:0)) + 5

さらに別解をいただいた:

birth = as.matrix(kaito[,5:11] == "1") %*% 6:12
birth = ifelse(birth < 6 | birth > 12, NA, birth)

これは次のようにしてもよさそう。

birth = (kaito[,5:11] == "1") %*% 6:12
birth = ifelse(birth %in% 6:12, birth, NA)

分布は次の通り:

table(birth, useNA="ifany")
birth
   6    7    8    9   10   11   12 <NA>
4102 4227 4235 4484 4501 4241 4544  459 

生まれた年度とHPVV接種率には強い関係がある:

s = table(kaito$V233, birth)
s
   birth
       6    7    8    9   10   11   12
  0   56   64   90   72   96   90  126
  1  496  428  452  663 1260 2038 3761
  2 3549 3735 3690 3745 3143 2112  656
  3    1    0    3    4    2    1    1

グラフにしてみよう:

plot(6:12, s[3,]/(s[2,]+s[3,])*100, ylim=c(0,100), type="o",
     pch=16, xlab="生まれた年度(平成)", ylab="接種率(%)")
生まれた年度とHPVV接種率

生まれた年度と「身体が自分の意思に反して動く」のクロス集計は次の通り:

u = table(kaito$V145, birth)
u
   birth
       6    7    8    9   10   11   12
  0   22   37   26   35   48   35   33
  1 4037 4141 4162 4414 4423 4176 4487
  2   43   49   47   35   30   30   24

グラフにしてみよう(95%信頼区間も付けた):

plot(6:12, u[3,]/(u[2,]+u[3,])*100, ylim=c(0,2), type="o", pch=16,
     xlab="生まれた年度(平成)", ylab="症状の割合(%)")
ci = sapply(1:7, function(x)binom.test(u[3,x],u[2,x]+u[3,x])$conf.int)
arrows(6:12, 100*ci[1,], 6:12, 100*ci[2,], length=0.03, angle=90, code=3)
生まれた年度と症状の割合

年齢が上(グラフで左側)のほうが少し多い傾向があるように見える。この傾向と,年齢が上のほうが接種率が高いこととが,年齢を交絡因子にしているのかもしれない。

横に生まれた年度(6〜12),縦に「症状なし・接種なし」「症状なし・接種あり」「症状あり・接種なし」「症状あり・接種あり」の度数を並べたクロス集計は次の通り:

t = sapply(1:7, function(x)as.numeric(table(kaito$V233, kaito$V145, birth)[2:3,2:3,x]))
t
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]  486  412  447  651 1232 2014 3722
[2,] 3496 3669 3627 3689 3097 2073  647
[3,]    6    9    1    6    8    9   18
[4,]   37   40   44   29   21   21    3

グラフにしてみよう(オレンジが接種あり,青が接種なし):

# 接種なし
plot(6:12, t[3,]/(t[1,]+t[3,])*100, ylim=c(0,2.5), type="o",
     col="#0068b7", xlab="生まれた年度(平成)", ylab="症状の割合(%)")
ci = sapply(1:7, function(x)binom.test(t[3,x],t[1,x]+t[3,x])$conf.int)
arrows(6:12, 100*ci[1,], 6:12, 100*ci[2,], length=0.03, angle=90, code=3, col="#0068b7")
# 接種あり
points(6:12, t[4,]/(t[2,]+t[4,])*100, type="o", pch=16, col="#f39800")
ci = sapply(1:7, function(x)binom.test(t[4,x],t[2,x]+t[4,x])$conf.int)
arrows(6:12, 100*ci[1,], 6:12, 100*ci[2,], length=0.03, angle=90, code=3, col="#f39800")
生まれた年度と症状の割合

生まれた年度ごとに見ると,「身体が自分の意思に反して動く」症状が見られる割合は,接種群(オレンジ)も非接種群(青)もほぼ1%前後で生じている。両方併せれば,さきほど見たように,年齢が上(グラフで左側)のほうが少し多い傾向があるように見える。ただ,特に高年齢(左)側の非接種群(青)は人数が少ないので,変動が大きい。

メタアナリシスで解説したフォレストプロットを,生まれた年度ごとに描いてみよう:

library(metafor)
dat = escalc(measure="OR", ai=t[4,], bi=t[2,], ci=t[3,], di=t[1,])
res = rma(yi, vi, data=dat)
forest(res, slab=paste0("平成", 6:12, "年度生まれ"))
Log Odds Ratio

95%信頼区間を見ると,平成11年度生まれがかろうじて有意であるが,全体をまとめると対数オッズ比はほぼ 0(95%信頼区間 [-0.41, 0.51])である。

なお,ここでのフォレストプロット(および一つ前のグラフ)の目的は,全体をまとめて一つの信頼区間で表し,それで何らかの結論を出すことではなく,個々の年齢層で違いがあるのかを吟味することである。例えば,接種率が8割を超えていた平成6〜9年度生まれと,急激に減少したそれ以降とで,傾向に違いがあるのか。

例:「なかなか眠れない」

上の例は症例数が少ないのでサンプリング誤差も大きい。そこで,もうちょっと人数の多い(平凡な)症状を調べてみよう。V96(No.13)の「なかなか眠れない」である:

t = sapply(1:7, function(x)as.numeric(table(kaito$V233, kaito$V96, birth)[2:3,2:3,x]))
# 接種あり
plot(6:12, t[3,]/(t[1,]+t[3,])*100, ylim=c(0,21), type="o",
     col="#0068b7", xlab="生まれた年度(平成)", ylab="症状の割合(%)")
ci = sapply(1:7, function(x)binom.test(t[3,x],t[1,x]+t[3,x])$conf.int)
arrows(6:12, 100*ci[1,], 6:12, 100*ci[2,], length=0.03, angle=90, code=3, col="#0068b7")
# 接種あり
points(6:12, t[4,]/(t[2,]+t[4,])*100, type="o", pch=16, col="#f39800")
ci = sapply(1:7, function(x)binom.test(t[4,x],t[2,x]+t[4,x])$conf.int)
arrows(6:12, 100*ci[1,], 6:12, 100*ci[2,], length=0.03, angle=90, code=3, col="#f39800")
plot(6:12, t[3,]/(t[1,]+t[3,])*100, ylim=c(0,21), type="o",
     col="#0068b7", xlab="生まれた年度(平成)", ylab="症状の割合(%)")
「なかなか眠れない」
library(metafor)
dat = escalc(measure="OR", ai=t[4,], bi=t[2,], ci=t[3,], di=t[1,])
res = rma(yi, vi, data=dat)
forest(res, slab=paste0("平成", 6:12, "年度生まれ"))
Log Odds Ratio

これはかなりはっきりした傾向がある。高年齢の非接種群はもともと睡眠障害があったのかもしれない(だから接種を見合わせたのかもしれない)。

(続きはそのうち書くかもしれない。)

リンク