Pythonによる体重・体脂肪データの解析 その3 番外編・欠損データの影響チェック
前回、体重データの周波数領域解析について記述しましたが、そこで1週間および1/3週間周期の微弱な成分が検出されたと言いました。
しかし、元のデータにあった欠損を線形補完していることが原因では無いかという懸念がありました。たとえば毎週日曜日に欠損があったとします。この値はランダムになるはずでした。しかしそれが線形になめらかに補完されるとするなら、毎週日曜日だけ乱雑さが小さくなるわけで、1週間周期の成分が出てきても不思議はなくなります。しかし失われたデータはもう戻ってこないので、実データでは検証のしようがありません。そこでシミュレーションでデータを生成し、欠損の影響を検証する事にしました。
1/fノイズの生成
まず体重データと同じような1/fの周波数特性を持つダミーデータを作成します。このためにcolorednoiseというライブラリを使用しました。 作者のFelix Patzelt氏のgithubはこちら。
使用に先立ち、コマンドプロンプトからpipを使ってインストールします。
pip install colorednoisePython内では以下のように
import colorednoise as cn
beta = 1
samples = 6201
y = cn.powerlaw_psd_gaussian(beta, samples)betaでPSDのスロープを指定します。1/fノイズならbeta=1。sample数は実データと合わせてます。結果はyに代入しておきます。この場合だとyは標準偏差1で1/fの周波数特性を持つガウス雑音になります。
実データ（欠損補間済）とシミュレーションデータのASDを比較します。
シミュレーションデータ（赤）と実データ(青)で概ねASDが一致しました。シミュレーションデータを定数倍するなどして、振幅を揃えないとな、と思っていましたが、偶然にもその必要はありませんでした。もちろんシミュレーションデータには1週間と1/3週間に対応するピークはあらわれません。
データ欠損の再現
つぎに、シミュレーションデータを持ってきて、実データと全く同じ場所に欠損を生じさせます。
a=(df['weight_nan']-df['weight_nan']+1)*(y+62)ここで、df[..]は欠損データをNaNに置き換えたバージョンの（つまり補間直前の）体重データです。前半のdf[..]-df[..]では同じ配列の引き算をして、数値があるばあいは0、NanであればNanを返すようにします。これに1を足しyを掛ければ、データに欠損を生じたようになります。yに62を足したのは体重データっぽく見せるためですが、detrendで平均値は除去してしまうので、本質的には不要な操作ですね。
単なる興味でいくつくらい体重データが欠損しているのかチェックしてみます。
np.sum(np.isnan(a))答えは1211でした。えっ！？思ったより多い！？6201サンプルのうち20%が欠損？うっかりミスに加え、長期出張があるとごっそり抜けるからですかね。仕方ない。
データ欠損のPSD/ASDへの影響をチェック
次に欠損を補間処理したものをbとします。
b=a.interpolate(method='linear', axis=0)下図は欠損の補間後(青)と補間前(赤)の比較です。
このb(青い時系列データ)のパワースペクトルにピークが出現したら、欠損データがピークの原因だったという事になりますが結果はどうか！？
なんとピークは出現しませんでした！1週間・1/3週間変動はリアルだった公算が高いですね！おおむねそのような結論でいいと思いますが、一応最後に欠損データに週間成分が含まれているかどうかだけチェックします。
a=(df['weight_nan']-df['weight_nan'])
b=a.replace(np.nan,1)
b.plot(color='blue',linewidth=.2)
欠損データを1、それ以外を0としたものです。いやほんとに欠損多い。このデータのASDをプロットすると・・・
これをみると、確かに欠損データに週間成分は含まれているようです。しかし、補間データには出ていませんし、いくつか試してみたところ、この欠損データがピークを説明するには、補間によるこの「雑音」が相当大きくないとピークが説明できないという感じでした。ですので、欠損がピークの原因では無かった、といって差し支えないと今は思っています。
以上、検証終わり。
#python3
#データ解析
#データ可視化
#時系列
#matplotlib
#時系列データ
#周波数領域
#スペクトル分析