【Rで時系列解析】cosinorパッケージでコサイナー解析，cosinor2も必要

最近はYoutubeでマリマリマリーの動画を見ています．その影響で，ひ○ゆきのものマネができたらいいなーと思うようになりました．皆さんも，次のセミナーではひ○ゆき風のしゃべりでお願いします．

　ということで，今回は前に紹介したコサイナー解析
chaos-kiyono.hatenablog.com
の続編です．以前の私の記事ではRのパッケージは使いませんでしたが，今回はcosinorパッケージを使います．ただし，cosinorパッケージだけでは，使い勝手が良くないので，cosinor2パッケージも一緒に使います．

周期的リズムのモデル

　ここでは，周期 $T$ を既知として，以下のモデルを仮定します．

　 $\displaystyle y(t) = \mu + A \cos \left(\frac{2 \pi t}{T} + \phi \right) + \varepsilon(t)$

ここで， $\mu$ は，MESOR (midline statistic of rhythm)， $A$ は振幅， $\phi$ はacrophase (頂点位相)と呼ばれるパラメタです． $\varepsilon(t)$ は，ランダムな要因やその他の変動を表します．前の記事では $\phi$ を $-\phi$ としましたが，今回は [Bingham et al. (1982)] に合わせて $+\phi$ )としてあります．

　コサイナー解析では，このモデルのパラメタを推定します．

cosinorとcosinor2パッケージでコサイナー解析

　cosinorパッケージを使ったコサイナー解析の問題点は，acrophase (頂点位相) $\phi$ が自動的に計算できないことです．これは， $\phi$ の推定にアークタンジェントを使っているからだと思います．数学では，アークタンジェントの値域は， $-\pi/2 < \phi < \pi/2$ (-90°から90°)です．しかし，これでは，225°が45°になる問題があるということです．

　cosinorパッケージの関数cosinor.lm(...)の結果をみれば，角度を正しく修正できますが，それでは面倒です．なので，私はcosinorパッケージを使うことはありませんでした．しかし，2022年にcosinor2というパッケージが登場して，この面倒をなくしてくれました．

　使い方の例として，Rスクリプトを載せておきます．

# cosinorパッケージを使った解析
#　【注意】　事前にcosinor, cosinor2をインストール
# install.packages("cosinor")
# install.packages("cosinor2")
#############################
# パッケージを読み込み
require(cosinor)
require(cosinor2)
#############################
# 数値実験
###
# 周期
T <- 24
# テスト用の数値データの生成
time <-  seq(0,24,1)
mesor <- 1
amp <- 3
phase <- 2.4
eps <- rnorm(length(time))
y.obs <- mesor+amp*cos(2*pi*time/T+phase) + eps
#####################
# プロット
par(mfrow=c(1,1),pty="m")
plot(time,y.obs,xlab="time",col=4,cex=1.2,xaxs="i")
curve(mesor+amp*cos(2*pi*x/T+phase),add=TRUE,col=8,xlim=c(0,24),lwd=2)
fit.cos <- cosinor.lm(y~time(t), period = 24, data=data.frame(t=time,y=y.obs))
para <- coef(fit.cos)
# paraの中身は
# mesor: para[[1]]
# amp: para[[2]]
# acr: para[[3]]
########################################
# 【重要】 cosinor2でアクロフェイズを補正
para[3] <- correct.acrophase(fit.cos) 
########################################
curve(para[[2]]*cos(2*pi*x/24+para[[3]])+para[[1]],xlim=c(min(time-12),max(time+12)),add=TRUE,col=2,lwd=3,lty=2)
para