カテゴリー
統計

Pythonによる分散分析

分散分析(Analysis Of VAriance; ANOVA)は3標本以上の差の検定です。

Pythonの統計パッケージには、scipy.statsやstatsmodelsがあります。最近はpingouinというものもあります。いずれもANOVAをサポートしています。それぞれのパッケージでの使用方法を説明します。

scipy.statsによるANOVA

from scipy import stats

x = [43, 55, 57, 72, 51, 52, 48, 46, 58, 54]
y = [53, 44, 54, 51, 68, 64, 45, 67, 49, 56]
z = [77, 55, 67, 54, 46, 75, 65, 57, 49, 61]
print(stats.f_oneway(x, y, z))

結果

F_onewayResult(statistic=1.6463168290164742, pvalue=0.21152459574247603)

statsmodelsによるANOVA

import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols

x = [43, 55, 57, 72, 51, 52, 48, 46, 58, 54]
y = [53, 44, 54, 51, 68, 64, 45, 67, 49, 56]
z = [77, 55, 67, 54, 46, 75, 65, 57, 49, 61]

values = x + y + z
groups = ['x'] * len(x) + ['y'] * len(y) + ['z'] * len(z)
data = pd.DataFrame({'values': values, 'groups': groups})

lm = ols('values ~ groups', data=data).fit()
sm.stats.anova_lm(lm, typ=2) # Type 2 Anova DataFrame

結果

               sum_sq    df         F    PR(>F)
groups     271.666667   2.0  1.646317  0.211525
Residual  2227.700000  27.0       NaN       NaN

pingouinによるANOVA

import pandas as pd
import pingouin as pg

x = [43, 55, 57, 72, 51, 52, 48, 46, 58, 54]
y = [53, 44, 54, 51, 68, 64, 45, 67, 49, 56]
z = [77, 55, 67, 54, 46, 75, 65, 57, 49, 61]

values = x + y + z
groups = ['x'] * len(x) + ['y'] * len(y) + ['z'] * len(z)
data = pd.DataFrame({'values': values, 'groups': groups})

print(pg.anova(data, dv='values', between='groups'))

結果

   Source  ddof1  ddof2         F     p-unc       np2
0  groups      2     27  1.646317  0.211525  0.108694
カテゴリー
統計

Pythonによる級内相関係数の計算

検者内・検者間の信頼性を表す指標として級内相関係数(ICC; Intraclass correlation)というのがあります[1, 2]。Pythonの統計パッケージにはscipy.statsstatsmodelsがありますが、ICCは実装されていない様です。

最近の統計パッケージpingouinにはありますので紹介します。

pingouinはpipでインストールできます。

pip install pingouin

サンプルコードは以下のとおりです。

import pandas as pd
import pingouin as pg

A = [1,1,3,1,1,2,1,2,1,1]
B = [2,1,3,1,3,2,1,3,3,3]
C = [2,3,3,1,1,1,1,2,3,3]
D = [2,3,3,1,1,2,1,2,3,1]
E = [2,3,3,3,3,2,1,2,3,1]
ratings = A + B + C + D + E
raters = ['A'] * len(A) + ['B'] * len(B) + ['C'] * len(C) + \
         ['D'] * len(D) + ['E'] * len(E)
targets = list('abcdefghij') * 5
data = pd.DataFrame({'targets':targets, 'raters':raters, 'ratings':ratings})
icc = pg.intraclass_corr(data=data, targets='targets', 
                         raters='raters', ratings='ratings')
print(icc.set_index('Type'))

出力結果は以下のとおりです。

                   Description       ICC         F  df1  df2      pval  \
Type                                                                     
ICC1    Single raters absolute  0.266854  2.819923    9   40  0.011519   
ICC2      Single random raters  0.280000  3.221007    9   36  0.005798   
ICC3       Single fixed raters  0.307576  3.221007    9   36  0.005798   
ICC1k  Average raters absolute  0.645380  2.819923    9   40  0.011519   
ICC2k    Average random raters  0.660377  3.221007    9   36  0.005798   
ICC3k     Average fixed raters  0.689538  3.221007    9   36  0.005798   

              CI95%  
Type                 
ICC1   [0.03, 0.64]  
ICC2   [0.05, 0.64]  
ICC3   [0.06, 0.67]  
ICC1k   [0.13, 0.9]  
ICC2k   [0.21, 0.9]  
ICC3k  [0.23, 0.91]  

[1] Shrout, P. E., & Fleiss, J. L. (1979). Intraclass correlations: uses in assessing rater reliability. Psychological bulletin, 86(2), 420. https://pdfs.semanticscholar.org/b8d4/7b0c0b12dd77543e82e6bf6636ddd335cfea.pdf

[2] 医療系データのとり方・まとめ方、対馬栄輝・石田水里、東京図書 (2013年). https://www.amazon.co.jp/dp/4489021437

カテゴリー
プログラミング

M5StickCでパルスオキシメーターを作る

M5StickCでパルスオキシメーターを作りました。センサにはM5Stack用心拍セEンサユニット(MAX30100)を使用しました。新型コロナウイルス感染症(COVID-19)では自覚症状がないうちに肺炎が悪化するとの報告もありますので、日頃から血中酸素飽和度(SpO2)を測定するのもいいかもしれません。

計測の様子。親指をセンサ部分に軽く押し当てる。

使用したパーツは以下のとおりです。

合計金額3,000円以下です。

M5StickCとM5Stack用心拍センサユニットの組み合わせ。Grove用4ピンケーブルで接続する。ケーブルは心拍センサユニットに付属する。

ArduinoのスケッチはMAX30100libライブラリに付属するサンプルスケッチMAX30100_Minimal.icoを改変しました。

変更した点は、以下のとおりです。

M5StickCのGroveポートをI2Cとして用いる事を宣言するコマンドを追加しました。これをしないと通信ができませんでした。

Wire.begin(32, 33);

LEDの電流量を50mAから7.6mAに変更しました。これをしないと測定値が飽和してしまい正しいデータが得られません。下記のサイトを参考にしました。

pox.setIRLedCurrent(MAX30100_LED_CURR_7_6MA);

M5StickCのLCD表示に関わる部分を追加しました。また、表示の周期を2秒にしました。

M5.Lcd.fillScreen(BLACK);

M5.Lcd.setTextSize(2);
M5.Lcd.setCursor(10, 10);
M5.Lcd.print("SpO2%");
M5.Lcd.setCursor(10, 100);
M5.Lcd.print("PRbpm");
M5.Lcd.setTextSize(5);
M5.Lcd.setCursor(15, 40);
M5.Lcd.printf("%2d", spo2);
M5.Lcd.setTextSize(3);
M5.Lcd.setCursor(20, 130);
M5.Lcd.printf("%3.0f", pulse_rate);
#define REPORTING_PERIOD_MS 2000

センサ部分には透明なセロハンテープを貼りました。これには下記のサイトを参考にしました。

センサに親指の腹を置いて10秒くらい待つとそれらしい数値が得られます。

作成したスケッチは以下のところに置きました。どうぞご利用ください。

https://github.com/attosci/M5StickC_PulseOximeter

カテゴリー
プログラミング

左利きのためのvimキーバインド

vimでインサートモードからコマンドモードに戻る時に便利なキーバインドとして、下記が紹介されている。

inoremap <silent> jj <ESC>

これは、右の人差し指でjを連打するだけでいいので重宝するキーバインドだ。ただ、左利きの人には右の指を使うより左の指を使う方が馴染む。よって、下記のキーバインドを提案する。

inoremap <silent> vv <ESC>

vが二回続く英単語はほとんどないのでおすすめ。

カテゴリー
プログラミング

Paspberry Pi Zero Wにtensorflowとkerasをインストールする

Raspberry Pi Zero W (Wireless)は、無線LAN、BLE(Bluetooth Low Energy)機能を搭載した小型のRaspberry Pi Zeroです。

これにtensorflowとkerasを入れます。
SDカードの容量は8GBです。これにOSとしてRaspbian Stretch Liteを入れてあります。以下のサイトからダウンロードしました。

本題とは関係ありませんが、ロボット制御に使うつもりなのでROS Kineticが入っています。

このボードはメモリ容量512MBと小さいため、メモリ不足が原因でインストールに失敗します。そのため、あらかじめスワップ領域を増やしておきます。

エディタで/etc/dphys-swapfileを開き、CONF_SWAPSIZE という行の値を変更します。デフォルトでは100MBとなっていますが、これを1GBに変更します(もう少し小さくてもいいかもしれません)。

sudo vi /etc/dphys-swapfile

CONF_SWAPSIZEが定義されている行を見つけ、数値を変更します。

CONF_SWAPSIZE=1024

viを終了し、以下のコマンドを実行。

sudo /etc/init.d/dphys-swapfile restart

以下のコマンドで変更されたか確認。

free -h

スワップ領域の変更に関する手続きは以下のサイトを参考にしました。

この後、下記コマンドを実行し、tensorflowをインストールします。とてつもなく時間がかかります。

sudo apt update
sudo apt install libatlas-base-dev
pip install tensorflow

次にkerasです。これもかなり時間がかかります。

sudo apt install libhdf5-serial-dev
pip install h5py
pip install keras

tensorflowとkerasのインストールは下記サイトを参考にしました。

はじめ、h5pyのインストールに失敗したので、下記サイトを参考に libhdf5-serial-devをインストールしています。

以上です。