Estimate from model not showing at the same scale

Modeling
, ,
1

I am using brms package in R. I have two models, an interaction and a simple model, that roughly have similar distributions of data by a factor of 3 levels (TrialRep) model1 <- brm(Inter_div_abs~TrialRep + (1|ID/AbsSuj), data=totaldf1, family=student) , and an interaction with grouping factor in the second model (TrialRep*Group) model2 <- brm(Inter_div_abs~TrialRep*Group + (1|ID/AbsSuj), data=totaldf2, family=student) . I would expect the population-effect estimates in both models to be around the same scale, but this is not the case. When I plot the conditional_effects I have the following results:

conditional_effects
conditional_effects1003×401 31.1 KB

In the y axis scale it is obvious something is not well. I post also the summary of both models:

This is population effects for model1 (TrialRep)

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     1.46      0.13     1.20     1.73 1.00    11034    21440
TrialRep2    -1.03      0.07    -1.18    -0.90 1.00    57843    51081
TrialRep3    -1.15      0.07    -1.29    -1.01 1.00    56529    50206

And this is population effects for model2 (interaction TrialRep*Group)

Population-Level Effects: 
                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept            3.77      0.20     3.39     4.16 1.00    11466    20753
TrialRep2           -1.83      0.06    -1.94    -1.72 1.00    58790    57464
TrialRep3           -2.21      0.06    -2.33    -2.10 1.00    58552    56927
GroupI              -0.08      0.28    -0.63     0.46 1.00    10058    18532
TrialRep2:GroupI     0.98      0.08     0.82     1.14 1.00    60255    58130
TrialRep3:GroupI     1.04      0.08     0.88     1.20 1.00    60912    57605

If I do the simple model from a frequentist approach, I get the following summary:

> m <- lmer(Inter_div_abs~ TrialRep + (1|ID/AbsSuj), data=totaldf)
> m
Linear mixed model fit by REML ['lmerModLmerTest']
Formula: Inter_div_abs ~ TrialRep + (1 | ID/AbsSuj)
   Data: totaldf
REML criterion at convergence: 63112.62
Random effects:
 Groups    Name        Std.Dev.
 AbsSuj:ID (Intercept) 0.0000  
 ID        (Intercept) 0.9875  
 Residual              4.6132  
Number of obs: 10694, groups:  AbsSuj:ID, 36; ID, 18
Fixed Effects:
(Intercept)    TrialRep2    TrialRep3  
      5.872       -3.458       -4.570  
optimizer (nloptwrap) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

Here, the fixed effects make sense now, they are closer to what I would expect. I am attaching the code I’m using so it can be reproducible and I am uploading the data I am working with. I have been struggling to figure out why this discrepancy is happening between a bayesian/frequentist approach, I also don’t understand why this is not happening with the interaction model.

totaldata2reduced.csv (142.2 KB)

totaldf <- read.csv("totaldata2reduced.csv")
names<-c('TrialRep','Trial','AbsSuj','ID')
totaldf[,names] <- lapply(totaldf[,names],factor)

q <- quantile(totaldf$Inter_div_abs, probs = 0.99)
#mu <- quantile(totaldf$Inter_div_abs, probs = 0.5)
for(i in 1:nrow(totaldf)){#if outlier, p50=0
  tmp <- totaldf$Inter_div_abs[i]
  totaldf$Inter_div_abs[i] <- ifelse(totaldf$Inter_div_abs[i]>q,NA,tmp)
}

bprior2 <- get_prior(Inter_div_abs~ (1|TrialRep) + (1|ID/AbsSuj), data=totaldf)

bprior2$prior[1] <- "student_t(3, 2, 3)"
bprior2$prior[2] <- "student_t(3, 0, 3)"


total_bayesian_mlm2 <- brm(Inter_div_abs~ TrialRep + (1|ID/AbsSuj), data=totaldf, family=student,
                           prior=bprior2, chains=4, cores=4, iter=20000, warmup=2000, init='0', control=list(adapt_delta=0.9, max_treedepth = 10))

Thanks in advance, any insight or guidance will be highly appreciated.

2

I finally understood the family=student was not appropriate in this case. The distribution was closer to family=gaussian or even family="skew_normal" and this was specially important in the model with a single predictor. I finally decided to compare the models through pp_check() and other model comparators like loo_compare() for instance.

So takeaway message, probably obvious for most of you (but a learning for me), we need to always take account of different diagnostics to specify appropriately the model.

1 Like