circular statistics (角度統計)について

髄膜炎に関する論文

Seasonal dynamics of bacterial meningitis: a time-series analysis. The Lancet Global Health http://dx.doi.org/10.1016/S2214-109X(16)30064-X

以前読んだ論文。前回はWaveletについて述べた。今日は、角度統計(circular statistics)について。

角度統計(Circular statistics)

ある地点における何かの向き（太陽、風向、樹木の樹冠の向き、生物の移動方向など）を記録したデータは、円周上に値をプロットできるデータみなせる。このようなデータを対象とする統計学は角度統計学と呼ばれている。日本語では、方向統計学、円周統計学、英語では、circular statistics, directional statistics, spherical statisticsなどが近い表現である。
角度は円周上の点と1対1対応するため、円周上の分布を用いて解析することができる(0度=360度となる)。
また、周期的な性質を持つデータは角度データとして表現可能であることから、例えば、ある疾病の月毎の発症数データなども角度データとしてみなすことができる。 円周上の分布は多様な方法が提案されてきた。

今回の論文中では、月別の髄膜炎を追っており、真の０や高い低いの評価が適切ではないため、角度統計が採用されている。

検定

Rayleigh test:
角度データには偏りがあるか? 1峰性の分布が想定されるときに有効。
Kuiper test:
角度データはvon Mises分布に従っているのか? 多峰性の分布が想定されてもOK。
Mardia-Watoson-Wheeler test:
2群の角度データは同一分布に従っているのか?

などがある。

信頼区間の推定

信頼区間の推定はブートストラップ信頼区間を用いる事ができる。以下のブログは参考になる。

ブートストラップ法で信頼区間を求めるときの注意点 - ほくそ笑む

参考資料

意外と資料が少ない。日本語だと、以下で外観は掴める。
角度データの統計処理基礎
環境科学における方向統計学の利用
角度分布の概要

R package "circular"

以下で試すことができる。簡単なものを試してみる。
p値の算出はいくつか方法がある

?`circular-package` #circular package help
library(circular) #ciucular
x <- circular(data.frame(runif(50, 0, pi/2))) #trial
par(mfcol=c(1,2))
plot(x, main = "plot") #plot 
rose.diag(x, main="rose diagram") #rose diagram

class(x) 
is.circular(x)

mean.circular(x)
rayleigh.test(x) #Rayleigh test of uniformity, assessing the significance of the mean resultant length.
kuiper.test(x)

sd.circular(x)

mle.vonmises(x) # estimation of mu and kappa 
mle.vonmises(x, mu=circular(0)) # estimation of kappa only


x.bs <- mle.vonmises.bootstrap.ci(x, alpha=.10) 
par(mfcol=c(1,2)) 
rose.diag(x.bs$mu, bins=30, main=expression(mu))
hist(x.bs$kappa, main=expression(kappa))

# Compute the median direction of a random sample of observations. 
median(x) #only the median is returned 
meandeviation(x)  #mean deviation is reported ## Not run:
median.circular(x) #only the median is returned 
medianCircular(x, deviation=TRUE) #both median and deviation are reported

mle.vonmises.bootstrap.ci(x, mu = NULL, bias = FALSE, alpha = 0.05, reps = 1000, control.circular = list())