記事一覧
古典統計学とベイズ統計での2群間の平均の差の検定
今回はRの組み込みのデータセットToothGrowthを用いて
2群間の平均の差を古典統計学とベイズ統計の実装であるRstanをもちいて
検定（推定）していきたいとおもいます。
supp サプリの種類 ”VC” ”OJ”
dose サプリの用量 0.5㎎ 1.0㎎ 2.0㎎
len 歯の長さ ㎜
head(ToothGrowth) len supp dose1 4.2 VC 0.52
平均の差を線形回帰で表現
引き続いて線形回帰で確認します。
t検定と線形回帰は数学的に同値です。
> fit.lm=lm(len~supp,data=ToothGrowth)
結果を確認します。
> summary(fit.lm)
Call:
lm(formula = len ~ supp, data = ToothGrowth)
Residuals:
Min 1Q Median 3Q Max
-12.7633
便利なパッケージ"brms"を使って線形回帰で平均の差を表現
同じことをRのパッケージである"brms"で行います。
このパッケージはベイズ統計をもとに、いろいろな線形回帰が行える
非常に便利なパッケージです。
fit_brms=brm( formula=len~supp, # 線形回帰を指定しています。 family=gaussian(), # 正規分布を指定しています。 data=ToothGrowth)
結果は以下の通りです。
分散分析をしてみよう。
前回までのToothGrowthでは、実は歯の長さlenの因子としては栄養の種類supp(OJとVC)、用量(0.5mg、1.0㎎、2.0㎎)の2要因あります。
交互作用と主効果を評価するためには2元配置分散分析が必要です。
まず、栄養の種類と用量に注目して箱ひげ図を作成します。
boxplot(len ~ dose + supp , data=ToothGrowth)
なんとなく、用量には
共分散分析を重回帰解析で
また同じToothGrowthのデータで考えていきます。
今回は２元配置分散分析を線形回帰(重回帰解析)で検討していきます。
まず、doseは数値型のままです。
lm.fit=lm(len ~ supp * dose , data=ToothGrowth)summary(lm.fit)
Call:
lm(formula = len ~ supp * dose)
Residuals:
Min 1
brmsパッケージで共分散分析をベイズ推定してみた。
library(brms)library(rstan)
用量doseを因子型に変換します。
data=ToothGrowthdata=transform(data,dose=as.factor(dose))
ベイズの線形回帰パッケージbrmsを使用します。
交互作用を考えるモデルを組みます。
brms.fit=brm(formula= len ~ supp*dose,family=gaus
分散分析の悩みどころ
引き続いて分散分析です。
前回はbrmsパッケージで交互作用ありで解析しましたが、交互作用なしで検討する考えもあります。多重共線性を排除するためです。他に今回は当てはまりませんが、中心化したり、標準化することもあります。ToothGrowthは因子×因子なので中心化も標準化も意味がなさそうです。
交互作用なしで解析すると、OJ1mgと2mgで有意差が出てしまいます。おそらく不適切な解析と考えま
２群間の時系列データの解析
友人に２群間の時系列データの解析、検定を尋ねられたので
改めて勉強してみたので、健忘録として。
まず、練習用のデータを作成していきます。
症例は30例ずつの2群で、計測は3時点、処置の有無で層別化します。
id = c(1:60)result0 = rnorm(30,40,10)result1 = rnorm(30,25,7)result2 = rnorm(30,20,5)result3 =