Извлечение байесовского p-значения из модели Лмера в R

Question

Я пытаюсь извлечь байесовские p-значения (т. Е. доля оценок, которые> 0, если точечная оценка отрицательна, или там, где точечная оценка положительна, доляоценки <0) </em>) по модели lmer, которую я сделал.Я понимаю, что "p-значения" по своей природе часты, но я требую, чтобы байесовские p-значения успокаивали рецензента ( аналогично этому пользователю ).

Для целей воспроизводимости яиспользуя набор данных из 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-значений длякаждая оценка?

Есть пример с оригинальным пакетом / кодом , но моя модель не настроена как их модель.

Мне не нужно использовать этот пакет, и я был бы так же рад рассчитать его из моих 1000 симуляций. В этом случае мне нужно знать , как рассчитать, какмногие оценки ниже / выше нуля.Это число / 1000 (общее число оценок) будет байесовским р-значением.

Answer 1

Для извлечения байесовских p-значений (, т. Е. Доли оценок, которые> 0, если точечная оценка отрицательна, или там, где точечная оценка положительна, доля оценок, которая составляет <0 </em>) Вы можете извлечь точечные оценки для каждого моделирования, а затем разделить на количество моделирования.

Чтобы сделать это, используя набор данных 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 

В этом случае говорится, что все оценки выше нуля, поэтому байесовское значение р составляет <0,001. Это подтверждает то, что мы видим, когда смотрим только на точечную оценку и 95% вероятных интервалов (т. Е. <code>Time оценка составляет 8,80, а 95% вероятных интервалов (8,38, 9,06). В обоих случаях мы видим сильную поддержку Time оказывает влияние на weight.