■
昨日のつづき
線形回帰で解析します。
> lm=lm(PRdelta~Medicine,BB)> lm=lm(PRdelta~Medicine,BB)> summary(lm)
Call:lm(formula = PRdelta ~ Medicine, data = BB)
Residuals:
Min 1Q Median 3Q Max
-24.0000 -6.0326 0.8696 5.8696 18.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 11.130 2.078 5.356 2.38e-06 ***
MedicineT 1.870 2.828 0.661 0.512
---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 9.967 on 48 degrees of freedom
Multiple R-squared: 0.009021, Adjusted R-squared: -0.01162 F-statistic: 0.437 on 1 and 48 DF, p-value: 0.5117
古典統計学では95%信頼区間から統計学的有意差なし。p値0.512
次はmcmcでベイズ推定 モデルは Xi1=a+b*Xi1+ei xi2=a+b*Xi2+ei ダミー変数を設定して > model<-" + data{ + int<lower=0> M; + int<lower=0> N; + real x1[M]; + real x2[N]; + } + parameters{ + real<lower=0> s1; + real<lower=0> s2; + real m1; + real m2; + } + model{ + for(i in 1:M) + x1[i] ~ normal(m1, s1); + for(i in 1:N) + x2[i] ~ normal(m2, s2); + }
つぎはbrmsパッケージを用いてマルコフ連鎖モンテカルロ法によるベイズ推定
> brm=brm(
+ formula=PRdelta~Medicine,
+ family=gaussian(),
+ data=BB,
+ )
> print(brm)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: PRdelta ~ Medicine
Data: BB (Number of observations: 50)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept 11.10 2.11 7.06 15.20 1.00
MedicineT 1.88 2.87 -3.76 7.59 1.00
Bulk_ESS Tail_ESS
Intercept 3715 2945
MedicineT 3809 2881
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sigma 10.13 1.05 8.35 12.31 1.00 3516
Tail_ESS
sigma 2552
Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
is a crude measure of effective sample size, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
2群間の平均の差は1.88程度と推定されている。
