0
$\begingroup$

First of all I apologize as some might see this a stupid question, but I cannot manage to find the answer.

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 plot composite

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 (only 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

Given that the datasets are based on a repetition of the same experiment I would expect them to be at the same scale as well, but this is not the case. I understand what I am modeling with this specification is Trial repetition 1 as Intercept and then model calculates the estimated differences in the next 2 levels, but I cannot figure out how to specify this simple model in order to have uncentered results and therefore, have a straightforward reading of the results.

Even if they come from different samples, in a sanity check calculating means from these two samples I obtain these results, which are clearly inside the scale I would expect, that is captured by the interaction model or model2.

#Simple Model: Group*TrialRep

> mean(totaldf1[totaldf1$TrialRep =="1",]$Inter_div_abs, na.rm=TRUE)
[1] 5.866142
> mean(totaldf1[totaldf1$TrialRep =="2",]$Inter_div_abs, na.rm=TRUE)
[1] 2.410974
> mean(totaldf1[totaldf1$TrialRep =="3",]$Inter_div_abs, na.rm=TRUE)
[1] 1.306225

#Interaction Model: Group*TrialRep

> mean(totaldf2[totaldf2$TrialRep =="1",]$Inter_div_abs, na.rm=TRUE)
[1] 5.607245
> mean(totaldf2[totaldf2$TrialRep =="2",]$Inter_div_abs, na.rm=TRUE)
[1] 3.41219
> mean(totaldf2[totaldf2$TrialRep =="3",]$Inter_div_abs, na.rm=TRUE)
[1] 2.906711

Final test, doing it in a frequentist approach, I get the fixed effects estimators close to what I would be expecting according to previous sanity check.

> m <- lmer(Inter_div_abs~ TrialRep + (1|ID/AbsSuj), data=totaldf)
boundary (singular) fit: see help('isSingular')
> 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  

I am posting the complete code for modeling I am using and I also upload the dataset and can be accessed from this link.

totaldf <- read.csv("totaldata2reduced.csv")

#Make sure this columns are factors
names<-c('TrialRep','Trial','AbsSuj','ID')
totaldf[,names] <- lapply(totaldf[,names],factor)

#Control of outliers
q <- quantile(totaldf$Inter_div_abs, probs = 0.99)
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)
}

#Specify prior
bprior2 <- get_prior(Inter_div_abs~ (1|TrialRep) + (1|ID/AbsSuj), data=totaldf)
#Edit prior
bprior2$prior[1] <- "student_t(3, 2, 3)"
bprior2$prior[2] <- "student_t(3, 0, 3)"

#Specify the model
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!

$\endgroup$
7
  • $\begingroup$ This question seems hard to answer without the data. Can you provide us with the data? Or a dataset which reproduces the issue and is okay to share even if it's not the "actual" data? $\endgroup$
    – dipetkov
    Feb 11, 2023 at 10:09
  • $\begingroup$ You can also consider asking the question on the Stan forum: discourse.mc-stan.org. $\endgroup$
    – dipetkov
    Feb 11, 2023 at 10:19
  • $\begingroup$ Thank you! I thought it could be a specification problem, not a programming problem, that's why I asked here. $\endgroup$ Feb 13, 2023 at 9:52
  • $\begingroup$ That may be the case. But then you provide few details about your data, variables, statistical problem formulation and no minimal reproducible example. $\endgroup$
    – dipetkov
    Feb 13, 2023 at 9:58
  • $\begingroup$ I understand, I also could have edited that, but the question was closed. Anyhow I am taking your advice and asking the question in Stan forum. Thanks again. $\endgroup$ Feb 13, 2023 at 10:11

1 Answer 1

1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.