从 R 中的 lmer 模型中提取贝叶斯 p 值
Extracting the Bayesian p-value from lmer model in R
我正在尝试提取贝叶斯 p 值（即，如果点估计为负，则估计 >0 的比例，或者点估计是正的，估计的比例 < 0)) 来自我所做的 lmer 模型。我知道 "p-values" 本质上是常客，但我需要贝叶斯 p 值来安抚审稿人 (similar to this user)。
出于可重现性的目的，我使用 R 中的数据集来说明我的问题。数据集：
library(datasets)
data(ChickWeight) #importing data from base R
summary(ChickWeight)
weight Time Chick Diet
Min. : 35.0 Min. : 0.00 13 : 12 1:220
1st Qu.: 63.0 1st Qu.: 4.00 9 : 12 2:120
Median :103.0 Median :10.00 20 : 12 3:120
Mean :121.8 Mean :10.72 10 : 12 4:118
3rd Qu.:163.8 3rd Qu.:16.00 17 : 12
Max. :373.0 Max. :21.00 19 : 12
(Other):506
我的真实数据同时具有连续和离散预测变量以及个体身份的随机效应。
正在创建 lmer 模型：
install.packages("lme4", dependencies=TRUE)
library(lme4)
m1<-lmer(weight ~ Time + Diet+ (1|Chick), data=ChickWeight)
summary(m1)
Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ Time + Diet + (1 | Chick)
Data: ChickWeight
REML criterion at convergence: 5584
Scaled residuals:
Min 1Q Median 3Q Max
-3.0591 -0.5779 -0.1182 0.4962 3.4515
Random effects:
Groups Name Variance Std.Dev.
Chick (Intercept) 525.4 22.92
Residual 799.4 28.27
Number of obs: 578, groups: Chick, 50
Fixed effects:
Estimate Std. Error t value
(Intercept) 11.2438 5.7887 1.942
Time 8.7172 0.1755 49.684
Diet2 16.2100 9.4643 1.713
Diet3 36.5433 9.4643 3.861
Diet4 30.0129 9.4708 3.169
Correlation of Fixed Effects:
(Intr) Time Diet2 Diet3
Time -0.307
Diet2 -0.550 -0.015
Diet3 -0.550 -0.015 0.339
Diet4 -0.550 -0.011 0.339 0.339
与 ChickWeight 数据集不同，我的真实数据集的估计值既有正值也有负值。
然后我想从我的模型中提取 95% 的可信区间 m1:
install.packages(c("MCMCglmm", "arm"), dependencies=TRUE)
library(MCMCglmm)
library(arm)
sm1<-sim(m1,1000)
smfixef=sm1@fixef #fixed effects
smranef=sm1@ranef #random effects
smfixef=as.mcmc(smfixef)
posterior.mode(smfixef) #extract estimates for fixed effects
(Intercept) Time Diet2 Diet3 Diet4
10.489143 8.800899 16.761983 31.684341 28.037318
HPDinterval(smfixef) ##extract 95% credible intervals for fixed effects
lower upper
(Intercept) -0.05392775 21.960966
Time 8.38244319 9.064171
Diet2 -0.46587564 34.061686
Diet3 17.90445947 53.817409
Diet4 11.17259787 48.467258
attr(,"Probability")
[1] 0.95
现在我想得到贝叶斯 p 值：
install.packages("conting", dependencies=TRUE)
library(conting)
bayespval(object=sm1, n.burnin = 0, thin = 1, statistic = "X2")
#this last line is the line I am having trouble with
Error: $ operator not defined for this S4 class
根据我设置模型的方式 m1，为每个估计提取贝叶斯 p 值的正确格式是什么？
original package/code 发布了一个示例，但我的模型没有像他们的模型那样设置。
我不需要使用这个包，并且很乐意从我的 1000 次模拟中计算它。 在那种情况下，我需要知道 如何计算有多少估计值 below/above 为零。该数字 / 1000（估计总数）将是贝叶斯 p 值。
提取贝叶斯 p 值（即，如果点估计为负，则估计 >0 的比例，或者如果点估计为正，则估计 <0 的比例) 您可以提取每个模拟的点估计，然后除以模拟次数。
要使用 ChickWeight 数据集和上述模型执行此操作，您需要：
library(datasets)
data(ChickWeight)
m1<-lmer(weight ~ Time + Diet+ (1|Chick), data=ChickWeight)
sm1<-sim(m1,1000)
smfixef=sm1@fixef
smfixef=as.mcmc(smfixef) #this has the 1000 simulations in it for the fixed effects
as.mcmc(smfixef)
Markov Chain Monte Carlo (MCMC) output:
Start = 1
End = 1000
Thinning interval = 1
(Intercept) Time Diet2 Diet3 Diet4
[1,] 17.52609243 8.381517 7.47169881 46.442343 19.7164997 #simulation 1
[2,] 16.52854430 8.859378 8.83279931 29.017547 25.4610474 #simulation 2
[3,] 4.00702870 8.830302 29.68309621 47.459395 35.1939344 #simulation 3
[4,] 16.44162722 8.599929 15.87393285 31.946265 33.7513144 #simulation 4
[5,] 21.07173579 8.596701 1.81909415 28.934133 19.0499998 #simulation 5
etc.
然后对于每一列，您可以编码哪些模拟高于或低于零：
p_Time=if_else(smfixef[,2]>0, 1,0) #Time variable (i.e., 2nd column)
因为 Time 变量的点估计值为正，您想要计算该变量的估计值低于零的次数：
sum_p_Time=sum(p_Time<1)
> sum_p_Time
0
在这种情况下，它表示所有估计值都在零以上，因此贝叶斯 p 值 < 0.001。这支持我们仅查看点估计和 95% 可信区间（即 Time 估计为 8.80，95% 可信区间为 (8.38, 9.06)）时所看到的结果。在这两种情况下，我们都看到对Time 对 weight 产生影响。
我正在尝试提取贝叶斯 p 值（即，如果点估计为负，则估计 >0 的比例，或者点估计是正的，估计的比例 < 0)) 来自我所做的 lmer 模型。我知道 "p-values" 本质上是常客，但我需要贝叶斯 p 值来安抚审稿人 (similar to this user)。
出于可重现性的目的，我使用 R 中的数据集来说明我的问题。数据集：
library(datasets)
data(ChickWeight) #importing data from base R
summary(ChickWeight)
weight Time Chick Diet
Min. : 35.0 Min. : 0.00 13 : 12 1:220
1st Qu.: 63.0 1st Qu.: 4.00 9 : 12 2:120
Median :103.0 Median :10.00 20 : 12 3:120
Mean :121.8 Mean :10.72 10 : 12 4:118
3rd Qu.:163.8 3rd Qu.:16.00 17 : 12
Max. :373.0 Max. :21.00 19 : 12
(Other):506
我的真实数据同时具有连续和离散预测变量以及个体身份的随机效应。
正在创建 lmer 模型：
install.packages("lme4", dependencies=TRUE)
library(lme4)
m1<-lmer(weight ~ Time + Diet+ (1|Chick), data=ChickWeight)
summary(m1)
Linear mixed model fit by REML ['lmerMod']
Formula: weight ~ Time + Diet + (1 | Chick)
Data: ChickWeight
REML criterion at convergence: 5584
Scaled residuals:
Min 1Q Median 3Q Max
-3.0591 -0.5779 -0.1182 0.4962 3.4515
Random effects:
Groups Name Variance Std.Dev.
Chick (Intercept) 525.4 22.92
Residual 799.4 28.27
Number of obs: 578, groups: Chick, 50
Fixed effects:
Estimate Std. Error t value
(Intercept) 11.2438 5.7887 1.942
Time 8.7172 0.1755 49.684
Diet2 16.2100 9.4643 1.713
Diet3 36.5433 9.4643 3.861
Diet4 30.0129 9.4708 3.169
Correlation of Fixed Effects:
(Intr) Time Diet2 Diet3
Time -0.307
Diet2 -0.550 -0.015
Diet3 -0.550 -0.015 0.339
Diet4 -0.550 -0.011 0.339 0.339
与 ChickWeight 数据集不同，我的真实数据集的估计值既有正值也有负值。
然后我想从我的模型中提取 95% 的可信区间 m1:
install.packages(c("MCMCglmm", "arm"), dependencies=TRUE)
library(MCMCglmm)
library(arm)
sm1<-sim(m1,1000)
smfixef=sm1@fixef #fixed effects
smranef=sm1@ranef #random effects
smfixef=as.mcmc(smfixef)
posterior.mode(smfixef) #extract estimates for fixed effects
(Intercept) Time Diet2 Diet3 Diet4
10.489143 8.800899 16.761983 31.684341 28.037318
HPDinterval(smfixef) ##extract 95% credible intervals for fixed effects
lower upper
(Intercept) -0.05392775 21.960966
Time 8.38244319 9.064171
Diet2 -0.46587564 34.061686
Diet3 17.90445947 53.817409
Diet4 11.17259787 48.467258
attr(,"Probability")
[1] 0.95
现在我想得到贝叶斯 p 值：
install.packages("conting", dependencies=TRUE)
library(conting)
bayespval(object=sm1, n.burnin = 0, thin = 1, statistic = "X2")
#this last line is the line I am having trouble with
Error: $ operator not defined for this S4 class
根据我设置模型的方式 m1，为每个估计提取贝叶斯 p 值的正确格式是什么？
original package/code 发布了一个示例，但我的模型没有像他们的模型那样设置。
我不需要使用这个包，并且很乐意从我的 1000 次模拟中计算它。 在那种情况下，我需要知道 如何计算有多少估计值 below/above 为零。该数字 / 1000（估计总数）将是贝叶斯 p 值。
提取贝叶斯 p 值（即，如果点估计为负，则估计 >0 的比例，或者如果点估计为正，则估计 <0 的比例) 您可以提取每个模拟的点估计，然后除以模拟次数。
要使用 ChickWeight 数据集和上述模型执行此操作，您需要：
library(datasets)
data(ChickWeight)
m1<-lmer(weight ~ Time + Diet+ (1|Chick), data=ChickWeight)
sm1<-sim(m1,1000)
smfixef=sm1@fixef
smfixef=as.mcmc(smfixef) #this has the 1000 simulations in it for the fixed effects
as.mcmc(smfixef)
Markov Chain Monte Carlo (MCMC) output:
Start = 1
End = 1000
Thinning interval = 1
(Intercept) Time Diet2 Diet3 Diet4
[1,] 17.52609243 8.381517 7.47169881 46.442343 19.7164997 #simulation 1
[2,] 16.52854430 8.859378 8.83279931 29.017547 25.4610474 #simulation 2
[3,] 4.00702870 8.830302 29.68309621 47.459395 35.1939344 #simulation 3
[4,] 16.44162722 8.599929 15.87393285 31.946265 33.7513144 #simulation 4
[5,] 21.07173579 8.596701 1.81909415 28.934133 19.0499998 #simulation 5
etc.
然后对于每一列，您可以编码哪些模拟高于或低于零：
p_Time=if_else(smfixef[,2]>0, 1,0) #Time variable (i.e., 2nd column)
因为 Time 变量的点估计值为正，您想要计算该变量的估计值低于零的次数：
sum_p_Time=sum(p_Time<1)
> sum_p_Time
0
在这种情况下，它表示所有估计值都在零以上，因此贝叶斯 p 值 < 0.001。这支持我们仅查看点估计和 95% 可信区间（即 Time 估计为 8.80，95% 可信区间为 (8.38, 9.06)）时所看到的结果。在这两种情况下，我们都看到对Time 对 weight 产生影响。