# Brms model linear vs non linear

**URL:** <https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033>\
**Category:** Modeling\
**Tags:** rstan, fitting-issues, brms\
**Created:** [May 6, 2024, 1:03am UTC](https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033 "2024-05-06T01:03:46Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Felix\_Requena](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/felix_requena/32/18454_2.png) [@Felix\_Requena](https://discourse.mc-stan.org/u/Felix_Requena)\
**Post date:** [May 6, 2024, 1:03am UTC](https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033/1 "2024-05-06T01:03:46Z")

</div>

Good morning,

I am running both linear and non-linear models in brms, and my results are quite similar. First, I am running my model using my response variable, along with two predictors and one random effect. Below are the results I obtained:

Family: lognormal  
Links: mu = identity; sigma = identity  
Formula: Social\_dist ~ 1 + X50 + Num\_Unique\_Healthy\_Ids + (1 | p | Name)  
Data: final\_data (Number of observations: 30538)  
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 4;  
total post-warmup draws = 4000

Group-Level Effects:  
~Name (Number of levels: 292)  
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk\_ESS Tail\_ESS  
sd(Intercept) 0.30 0.01 0.27 0.33 1.00 2133 3041

Population-Level Effects:  
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk\_ESS  
Intercept 1.45 0.02 1.41 1.49 1.00 1027  
X50 0.10 0.01 0.08 0.12 1.00 3443  
Num\_Unique\_Healthy\_Ids -0.11 0.01 -0.12 -0.09 1.00 3532  
Tail\_ESS  
Intercept 2006  
X50 3890  
Num\_Unique\_Healthy\_Ids 3737

Family Specific Parameters:  
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk\_ESS Tail\_ESS  
sigma 0.83 0.00 0.82 0.83 1.00 4041 3974

Draws were sampled using sampling(NUTS). For each parameter, Bulk\_ESS  
and Tail\_ESS are effective sample size measures, and Rhat is the potential  
scale reduction factor on split chains (at convergence, Rhat = 1).

I have also created a non-linear model. After extensive research, I’m not entirely certain if it’s the correct approach for a non-linear model, but I have included quadratic effects for both predictors.

Family: lognormal  
Links: mu = identity; sigma = identity  
Formula: Social\_dist ~ 1 + X50^2 + Num\_Unique\_Healthy\_Ids^2 + (1 | p | Name)  
Data: final\_data (Number of observations: 30538)  
Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 4;  
total post-warmup draws = 4000

Group-Level Effects:  
~Name (Number of levels: 292)  
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk\_ESS Tail\_ESS  
sd(Intercept) 0.30 0.01 0.27 0.33 1.00 2033 3073

Population-Level Effects:  
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk\_ESS Tail\_ESS  
Intercept 1.45 0.02 1.41 1.49 1.00 1065 2007  
X50 0.10 0.01 0.08 0.12 1.00 3478 3663  
Num\_Unique\_Healthy\_Ids -0.11 0.01 -0.13 -0.09 1.00 3618 3818

Family Specific Parameters:  
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk\_ESS Tail\_ESS  
sigma 0.83 0.00 0.82 0.83 1.00 3740 3356

Draws were sampled using sampling(NUTS). For each parameter, Bulk\_ESS  
and Tail\_ESS are effective sample size measures, and Rhat is the potential  
scale reduction factor on split chains (at convergence, Rhat = 1).

This is the LOO comparison.

> loo\_compare(loo\_linear, loo\_nonlinear)  
> elpd\_diff se\_diff  
> fit\_log 0.0 0.0  
> fit\_nonlinear -1.0 0.4

After doing some reading, I think the models have been executed well, but I’m not entirely sure. Could you please provide me with some assistance? Thank you very much!

Cheers,  
Felix

---

<div class="post-metadata">

**Author:** ![wpetry](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wpetry/32/4852_2.png) [@wpetry](https://discourse.mc-stan.org/u/wpetry)\
**Post date:** [May 6, 2024, 3:03am UTC](https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033/2 "2024-05-06T03:03:43Z")

</div>

Welcome to the Stan community, @Felix_Requena!

It looks like you’re fitting two candidate models: one with untransformed predictors and one where the predictors have been squared. Importantly, both models still fall within the linear model framework. (`brms` has a separate syntax for fitting a wider range of non-linear functional forms, but simple polynomial terms are better fit using the standard linear syntax as you do here.) Squaring the predictors will throw away information on the sign of the predictors, which is ok if only magnitude matters. Also, it should make your fit more sensitive to the data points with large values of the predictor variables. For some problems, it often is useful to include linear and polynomial terms for the predictors in the model. It’s difficult to say what’s reasonable in your case without knowing more about your specific motivations and having domain expertise in your field.

On a minor note, your group-level intercept for `Name` doesn’t need the `(... | p | ...)`. It’s not hurting anything, but isn’t being used either. This syntax is useful for setting up correlations between the group-level effects for models with multivariate responses. `(1 | Name)` will do just fine for your models.

Both of your candidate models seem to fit without any obvious problems based on the R-hat statistics you’ve provided. Any warning messages about divergences wouldn’t show up in these summary tables, but would be reason to investigate before proceeding further. A [posterior predictive check](http://paul-buerkner.github.io/brms/reference/pp_check.brmsfit.html) is another good step to add to your post-fitting workflow.

Lastly, you use leave-one-out cross validation to measure the expected predictive accuracy of the models and to compare the two fits. By this metric, `fit_log` is the better model.

You’re doing a great job thinking carefully, asking questions, and sharing details for others to help you navigate your issue. Please feel welcome to add more to the thread if this didn’t answer your question and to open new topics as challenges arise.

---

<div class="post-metadata">

**Author:** ![Felix\_Requena](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/felix_requena/32/18454_2.png) [@Felix\_Requena](https://discourse.mc-stan.org/u/Felix_Requena)\
**Post date:** [May 6, 2024, 4:55am UTC](https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033/3 "2024-05-06T04:55:48Z")

</div>

Thank you very much for the warm welcome and your assistance. I’ll strive to clarify the advertisement further.

To delve into my models, I’m exploring how my dependent variable changes based on mixed effects, represented by the number of healthy individuals within a home range, denoted as X50. I’ve developed a same model for disease individuals.

I’m endeavoring to construct both nonlinear and linear models. In animals, it has been previously observed that when there are more healthy individuals, animals tend to be closer, suggesting a linear relationship. However, with the increasing number of individuals contracting diseases each year, I’m investigating the possibility of a nonlinear relationship. Additionally, my findings deviate somewhat from prior research; specifically, I’m noticing that animals closer to diseased animals tend to draw closer, whereas conventional wisdom suggests that more diseased animals should lead to greater separation, and healthy animals should cluster together.

As the results were pretty similar between linear and non linear model (see first comment), I did more reading and I found the possibility to use poli term to do more flexibility to the non model, and I found some non linear effects. When I compared (LOO) with the original model (see first comment) this model is better fit.

This are outcomes that I have obtained in my new models  
1- Healthy individuals

```no-highlight
H_fit_nonlinear_improved1 <-bf(Social_dist ~ 1 + poly(X50, 2) + poly(Num_Unique_Healthy_Ids, 2) + (1|Name)) + lognormal()

```

Outcome

```no-highlight
Group-Level Effects: 
~Name (Number of levels: 292) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.27 0.01 0.25 0.30 1.00 2163 3246

Population-Level Effects: 
                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.44 0.02 1.41 1.48 1.00 1245 2306
polyX5021 17.16 1.74 13.75 20.64 1.00 3496 3724
polyX5022 -7.08 1.40 -9.87 -4.39 1.00 3624 3679
polyNum_Unique_Healthy_Ids21 -18.88 1.55 -21.91 -15.86 1.00 2995 3430
polyNum_Unique_Healthy_Ids22 6.52 1.31 3.96 9.06 1.00 3995 3932

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.83 0.00 0.82 0.83 1.00 3709 3518

```

Posterior predictive check

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/d/5/d5a7d558ce26fd990634f5df4e53772a6f490c5b.png)

2- Diseased individuals

```no-highlight
fit_nonlinear_improved1 <- bf(Social_dist ~ 1 + poly(X50, 2) + poly(Num_Unique_Diseased_Ids, 2) + (1|Name)) + lognormal()

```

Outcome

```no-highlight
Group-Level Effects: 
~Name (Number of levels: 292) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.33 0.02 0.30 0.37 1.00 1840 2805

Population-Level Effects: 
                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.45 0.02 1.41 1.49 1.00 835 1586
polyX5021 13.31 1.63 10.09 16.47 1.00 3128 3739
polyX5022 -4.61 1.41 -7.42 -1.92 1.00 3321 3420
polyNum_Unique_Diseased_Ids21 -13.32 1.65 -16.61 -10.07 1.00 3003 3739
polyNum_Unique_Diseased_Ids22 2.29 1.29 -0.25 4.76 1.00 3700 3749

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.83 0.00 0.82 0.83 1.00 3755 3667

```

Posterior predictive check

 ![b87e6e4f-a591-407b-a0a9-af832f612d27](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/9/d/9d0e76f843590bc13237d34a07af54d6237fc184.png)

The posterior checks are not perfect, but the data fits well with the line, and from what I see, the models are well-built. Could you please advise on which model would be better to use in my case? I have been working on these models for four months, and I am unsure if the results I am obtaining are accurate.

Thank you very much.  
Cheers,  
F

---

<div class="post-metadata">

**Author:** ![wpetry](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wpetry/32/4852_2.png) [@wpetry](https://discourse.mc-stan.org/u/wpetry)\
**Post date:** [May 6, 2024, 2:08pm UTC](https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033/4 "2024-05-06T14:08:49Z")

</div>

Thanks for the additional information. It sounds like you have good rationale for using squared predictors in your model, and I agree that `poly(..., 2)` is a good way to express these in your formula. (Note that it’s equivalent to adding a linear/untransformed term for the predictor as I was asking about above, i.e.,`y ~ 1 + x + I(x^2)`.)

> [@Felix\_Requena](#):
>
> Could you please advise on which model would be better to use in my case?

If I understand you correctly, you’ve now fit four candidate models: `fit_log`, `fit_nonlinear`, and `H_fit_nonlinear_improved1` include the same predictor variables just with increasing complexity. The fourth model, `fit_nonlinear_improved1`, mirrors the third model but substitutes the diseased individuals for healthy individuals as a predictor.

There isn’t a universal meaning of “better” when it comes to modeling. Better for what purpose? If your goal is prediction, you can compare all four models with LOO to find the model with the best predictive performance. If you want to evaluate the evidence in favor/against a given predictor, you should simply include it in the model and evaluate the posterior distribution of its coefficient. For example, it sounds like there is an interesting hypothesis about the effects of healthy vs. diseased individuals on the response. For this, you’d need both predictors in a single model. Note that this approach is very different than the model selection route you’ve been taking. Doing model selection prior to inference on the same data set will bias your inferences ([Treddenick et al. 2021](https://doi.org/10.1002/ecy.3336)), although there are some regularization approaches that will downweight predictors that appear to be less important (lasso, ridge, horseshoe priors).

All of your models seem to be fitting without obvious technical problems, which is where the Stan forums are most useful. Which one _you_ should use is a scientific question and will rely on your goals and expertise to move ahead.

---

<div class="post-metadata">

**Author:** ![Felix\_Requena](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/felix_requena/32/18454_2.png) [@Felix\_Requena](https://discourse.mc-stan.org/u/Felix_Requena)\
**Post date:** [May 8, 2024, 10:46pm UTC](https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033/5 "2024-05-08T22:46:45Z")

</div>

Sorry for the delay in replying. Thank you very much for your help. I tried to follow one of your pieces of advice—I decided to include both predictors in a model with an interaction effect. After discussing it with my supervisor, it seems better this way. I am now using this model.

```no-highlight
fit_log_all_int <- bf(Social_dist ~ 1 + X50 +
                    Num_Unique_Diseased_Ids*Num_Unique_Healthy_Ids +  
                      (1|Name)) + lognormal()

```

```no-highlight
Population-Level Effects: 
                                               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.37 0.03 1.32 1.42 1.00 1895 2815
X50 0.13 0.01 0.11 0.15 1.00 3658 3797
Num_Unique_Diseased_Ids -0.05 0.01 -0.07 -0.03 1.00 3437 3629
Num_Unique_Healthy_Ids -0.17 0.01 -0.19 -0.15 1.00 3392 3809
Num_Unique_Diseased_Ids:Num_Unique_Healthy_Ids 0.03 0.01 0.02 0.04 1.00 3805 3892

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.82 0.00 0.82 0.83 1.00 3625 3895

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

```

I think up until now the outcome is giving me what I need. The posterior checks are exactly the same as my previous comments. Thank you so much again. One last piece of advice?  
Best,  
F

---

<div class="post-metadata">

**Author:** ![wpetry](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wpetry/32/4852_2.png) [@wpetry](https://discourse.mc-stan.org/u/wpetry)\
**Post date:** [May 8, 2024, 11:24pm UTC](https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033/6 "2024-05-08T23:24:33Z")

</div>

> [@Felix\_Requena](#):
>
> One last piece of advice?

I don’t understand—what part are you looking for feedback on?

---

<div class="post-metadata">

**Author:** ![Felix\_Requena](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/felix_requena/32/18454_2.png) [@Felix\_Requena](https://discourse.mc-stan.org/u/Felix_Requena)\
**Post date:** [May 9, 2024, 1:53am UTC](https://discourse.mc-stan.org/t/brms-model-linear-vs-non-linear/35033/7 "2024-05-09T01:53:38Z")

</div>

Sorry, I think I sent a previous version of the draft. I believe I am getting to the point, but I will also spend some time thinking deeply about it. Thank you very much for all your help.

Cheers,  
Felix