# Choosing a sampling distribution for left skewed data

**URL:** <https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497>\
**Category:** brms\
**Created:** [March 15, 2024, 1:26am UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497 "2024-03-15T01:26:35Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 15, 2024, 1:26am UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/1 "2024-03-15T01:26:35Z")

</div>

I’m trying to decide on a sampling distribution for a response variable that is left skewed and slightly bimodal. The response variable is fatty acid concentration of tissues, so observed response values are continuous, and positive.

When I fit a simple intercept model with family = gaussian(link=“identity”), pp\_check() indicates negative posterior predictions for the response, which doesn’t make biological sense.

```no-highlight
condition_master %>%
  filter(lme == "EBS", 
         !vial_id %in% c("2019-65","2019-67","2019-68","2019-71","2019-66"), 
         maturity != 1) -> ebs.dat

model_normal <- brm(Total_FA_Conc_WWT ~ 1, family = gaussian(link="identity"), data = ebs.dat)

```

I tried to remedy this by truncating the response, but I’m not entirely sure I understand the implications of doing this, and the potential trade offs of instead using a model that allows for negative posterior predictions. I also looked at skew normal, gamma and lognormal distributions, and compared how well each model captures the range, mean, min and max of the data. I was unable to run a truncated skew normal model, as it appeared that predictions could not be evaluated properly at the lower end, leading to NA’s.

```no-highlight
# fitting a test brms model with a Gaussian likelihood - truncating response at 0, else models predict values < 0
  #ie. responses out of bounds are discarded
model_normal <- brm(Total_FA_Conc_WWT | trunc(lb = 0) ~ 1, family = gaussian(link="identity"), data = ebs.dat)
summary(model_normal)

# fitting a test brms model with a skew normal likelihood
model_skew <- brm(Total_FA_Conc_WWT ~ 1, family = skew_normal(), data = ebs.dat)
summary(model_skew)
#If bounded- predictions cannot be evaluated properly at the lower bound, leading to NAs?
#If not bounded, negative posterior predictions for response

# fitting a test brms model with a gamma likelihood
model_gamma <- brm(Total_FA_Conc_WWT ~ 1, family = "gamma", data = ebs.dat)
summary(model_gamma)

# fitting a test brms model with a lognormal likelihood
model_log <- brm(Total_FA_Conc_WWT ~ 1, family = lognormal(), data = ebs.dat)
summary(model_log)

```

Should have titled, but plots are from top to bottom, L to R: gaussian truncated, skew normal, gamma, and lognormal.

```no-highlight
# posterior predictive checking
  pp_check(model_normal, ndraws = 1e2) + pp_check(model_skew, ndraws = 1e2) +
 pp_check(model_gamma, ndraws = 1e2) + pp_check(model_log, ndraws = 1e2) 
#skew normal captures mean and variance best, none of the models picking up apparent bimodality of data 

```

 ![Rplot](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/a/a/aa0b32af9aceb35e5c43fb0e9ec522a5d8a4221d.png)

```no-highlight
# posterior predictive checking - boxplots
pp_check(model_normal, type = "boxplot", ndraws = 20) + pp_check(model_skew, type = "boxplot", ndraws = 20) +
  pp_check(model_gamma, type = "boxplot", ndraws = 20) + pp_check(model_log, type = "boxplot", ndraws = 20) 

```

 ![Rplot01](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/5/5/5576ed46a980f197fc88c769a71050009034399e.png)

```no-highlight
#and means
pp_check(model_normal, type = "stat", stat = "mean") + pp_check(model_skew, type = "stat", stat = "mean") +
  pp_check(model_gamma, type = "stat", stat = "mean") + pp_check(model_log, type = "stat", stat = "mean")

```

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/e/c/ecc5604e66aeacde73abe6f2be939bb764484b81.png)

loo\_compare(criterion=“waic”) indicates that predictive accuracy is highest for the gaussian truncated model, but I assume that the skew normal model isn’t top solely because I haven’t effectively truncated the response? Log transforming the raw response doesn’t appear to help much either so I’m uncertain which road to take. How “bad” is bad in regards to posterior predictive checks? And should I just be using family=“gaussian” with truncation given the loo results for my model runs that include covariates?

Windows 10  
brms\_2.19.0

---

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 15, 2024, 2:12am UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/2 "2024-03-15T02:12:25Z")

</div>

Also including distribution of minimum and maximum values for posterior distributions

```no-highlight
pp_check(model_normal, type = "stat", stat = "min") + pp_check(model_skew, type = "stat", stat = "min") +
  pp_check(model_gamma, type = "stat", stat = "min") + pp_check(model_log, type = "stat", stat = "min")
#skew normal is predicting negative values 

```

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

```no-highlight
#now maximum values
pp_check(model_normal, type = "stat", stat = "max") + pp_check(model_skew, type = "stat", stat = "max") +
    pp_check(model_gamma, type = "stat", stat = "max") + pp_check(model_log, type = "stat", stat = "max")
#lognormal is way overshooting max 

```

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/3/2/3287ecfcc2e64a7b8847d14d2ac4ef87e9394b90.png)

---

<div class="post-metadata">

**Author:** ![AWoodward](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/awoodward/32/13015_2.png) [@AWoodward](https://discourse.mc-stan.org/u/AWoodward)\
**Post date:** [March 15, 2024, 3:41am UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/3 "2024-03-15T03:41:58Z")

</div>

Hello @Erin_Fedewa_NOAA_Federal, I suspect that for a positively constrained variable like concentration, your model will need to include some structural constraint to prevent negative predictions, because as you say the normal response model will predict negative concentrations, so the normal and skewnormal response distributions are probably unreasonable at face-value. If the concentration is in fact left-skewed, then I suppose that the major weakness of lognormal or gamma response models will be the right-skew of the distribution, which seems apparent in the PPC. Left-skewed distributions for biologic concentrations are fairly unusual. What does this variable actually look like (say on the log scale)? Might a flexible distribution like the skewnormal on the log-tranformed responses be an option? Personally I would probably take this approach rather than truncation.

A side point is that I presume you will want to work with some covariates here; is it feasible that the left skew of the distribution is explainable by a major covariate?

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 15, 2024, 12:53pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/4 "2024-03-15T12:53:32Z")

</div>

- The common terminology would say that your data are right skewed (see [Skewness - Wikipedia](https://en.wikipedia.org/wiki/Skewness))
- Have you used `loo()` and `loo_compare()` to check which is providing the best predictions?
- I recommend to use also `pit_ecdf` plots.
- Can you post a reproducible example of failure with truncated skew normal? I think this is something that should work, but may require some changes in the code to make it more robust
- Eye balling the current plots, the skew normal looks the best, so getting that truncated skew normal would be useful. You could also test skew normal for log(y)

---

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 15, 2024, 3:54pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/5 "2024-03-15T15:54:42Z")

</div>

Thank you for the advice!

This is the response variable on log scale:  
`ebs.dat %>% ggplot(aes(log(Total_FA_Conc_WWT))) + geom_density()`  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/5/9/595cee2d371d7bac93a05338ab4ab2899dc920bb.png)  
I went ahead and tried skew normal on the log-transformed response and the PPC plots definately look much better, so maybe I’ll go that route

```no-highlight
#skew normal with log response 
model_skew_log <- brm(log(Total_FA_Conc_WWT) ~ 1, family = skew_normal(), data = ebs.dat)

```

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

I had also had the same thought about a covariate explaining the skewed distribution, and imagine that would be very probable. I’ll play around with some model runs including covariates as well, and reassess the PPC plots

---

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 15, 2024, 4:14pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/6 "2024-03-15T16:14:34Z")

</div>

The truncated skew normal model runs fine, I just get a warning for NAs being excluded from the model. I also can’t run pp\_check() on this model, I assume due to the NAs?

```no-highlight
#Bounded skew normal 
model_skew_bound <- brm(Total_FA_Conc_WWT | trunc(lb = 0) ~ 1, family = skew_normal(), data = ebs.dat)
#If bounded- predictions cannot be evaluated properly at the lower bound, leading to NAs?

Warning message:
Rows containing NAs were excluded from the model.

> pp_check(model_skew_bound, ndraws = 1e2)
Error in `validate_predictions()`:
! NAs not allowed in predictions.
Run `rlang::last_trace()` to see where the error occurred.
There were 50 or more warnings (use warnings() to see the first 50)
> 

```

Re. my response to @AWoodward, maybe I should pursue log transformation of the response, and the skew normal distribution instead of truncation? I used loo\_compare() on the truncated Gaussian, skew normal, gamma and lognormal, and the gaussian truncated had the highest predictive accuracy, but I assume that’s because the skew normal model was predicting negative response values. I assume I can’t use loo to compare a truncated vrs a log-response model because the response variables differ? (see script below). I will look into pit\_ecdf plots, that’s a new one for me but thank you for the suggestion!

```no-highlight
> loo_compare(model_normal, model_skew, model_skew_log, model_gamma, model_log, criterion = "waic")
                          elpd_diff se_diff
model_skew_log 0.0 0.0
model_normal -2142.4 20.3
model_skew -2153.5 16.7
model_gamma -2153.8 17.7
model_log -2217.2 18.0
Warning message:
Not all models have the same y variable. ('yhash' attributes do not match) 
> 

```

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 15, 2024, 4:34pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/7 "2024-03-15T16:34:40Z")

</div>

> [@Erin\_Fedewa\_NOAA\_Federal](#):
>
> The truncated skew normal model runs fine, I just get a warning for NAs being excluded from the model. I also can’t run pp\_check() on this model, I assume due to the NAs?

Strange.

> [@Erin\_Fedewa\_NOAA\_Federal](#):
>
> I assume I can’t use loo to compare a truncated vrs a log-response model because the response variables differ?

You can if you take into account the Jacobian transformation. There is an example of that in Regression and Other Stories chapter 12, with code at [Regression and Other Stories: Mesquite](https://avehtari.github.io/ROS-Examples/Mesquite/mesquite.html#Predict_log(weight)_given_log_transformed_predictors).

---

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 15, 2024, 5:32pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/8 "2024-03-15T17:32:23Z")

</div>

Thanks for the helpful example @avehtari. I used the Jacobian transformation for Gaussian log(y) and skew normal log(y) models, and then ran loo\_compare to assess the Gaussian log(y), skew normal log(y), Gaussian truncated, and skew normal truncated. I’m still trying to sort out what’s going on with the truncated skew normal model, and am not even able to run loo() for it

```no-highlight
> loo_skew_trunc <- loo(model_skew_bound)
Error in if (varx == 0) { : missing value where TRUE/FALSE needed
In addition: There were 50 or more warnings (use warnings() to see the first 50)

```

I dropped the skew normal truncated model from loo\_compare() because of this, and the skew normal log(y) model comes out on top.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/3/7/37627c7f2f005d10f2a2fce8242116c1424d2b91.png)  
If I use the skew normal distribution with a log transformed response, should I be worried that my PPC plot (below) still doesn’t look great? I’ll continue to assess once I run more complex models with covariates, and maybe these added covariates will help explain the underlying distribution of the data.  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/e/a/eaf552737d4a9af81652f1650832b754acc1ac78.png)

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 15, 2024, 6:41pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/9 "2024-03-15T18:41:57Z")

</div>

I think the dens\_overlay plot for skew normal looks better than for skew normal log. It would be good to get the truncated model to work, and I could have a look, but a reproducible example would help.

I now realized that your formula had `~ 1`, so that you didn’t have any covariates. If you do have covariates, they can change completely which distribution is best describing the conditional distribution.

---

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 15, 2024, 8:05pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/10 "2024-03-15T20:05:25Z")

</div>

Sorry, forgot to attach the data on my last reply for a reproducible example!  
[total\_FA\_master.csv](https://discourse.mc-stan.org/uploads/short-url/w4cF1PhXjpQKLoIEO3EVrsLmGvo.csv) (286.9 KB)

```no-highlight
#data wrangling 
condition_master <- read.csv("./data/total_FA_master.csv")

condition_master %>%
 mutate(julian=yday(parse_date_time(start_date, "mdy", "US/Alaska"))) %>% #add julian date 
 filter(lme == "EBS", 
        !vial_id %in% c("2019-65","2019-67","2019-68","2019-71","2019-66"), 
        maturity != 1) %>%
 mutate(year = as.factor(year),
        sex = as.factor(sex),
        region = as.factor(sample_region),
        station = as.factor(gis_station),
        temperature = as.numeric(gear_temperature),
        fourth.root.cpue = as.numeric(cpue^0.25),
        fourth.root.invert = as.numeric(total_benthic_cpue^0.25)) -> ebs.dat

```

I was just starting to explore covariate structure and more complex model runs, so as a first pass, I ran a skew normal model that includes two population-level covariates and a group-level effect to explain the spatial structure of the sampling design.

```no-highlight
mod5_formula <- bf(Total_FA_Conc_WWT ~ cw + s(temperature, k=4) + (1 | region)) 

mod5 <- brm(mod5_formula,
            data = ebs.dat,
            family = skew_normal(),
            cores = 4, chains = 4, iter = 2500,
            save_pars = save_pars(all = TRUE),
            control = list(adapt_delta = 0.999, max_treedepth = 14))

```

I had one divergent transition after warmup, and the PPC plot actually looks pretty good. But still predicting negative response values, so maybe a skew normal truncated model is the best way to go moving forward.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/2/a/2aad787fae1abcf04f1ec9a06f75e51a8592a7d3.png)  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/3/c/3c5bf09c4289f48b39a9e360c0b7e5418411df20.png)

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 19, 2024, 7:56am UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/11 "2024-03-19T07:56:47Z")

</div>

Thanks for the example. I can replicate the error. I’ll investigate this.

---

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 19, 2024, 4:15pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/12 "2024-03-19T16:15:17Z")

</div>

Thank you!!

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [March 19, 2024, 4:29pm UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/13 "2024-03-19T16:29:58Z")

</div>

should be fixed in the github version of brms.

---

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 20, 2024, 5:13am UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/14 "2024-03-20T05:13:45Z")

</div>

@paul.buerkner @avehtari I installed the github version of brms and the truncated skew normal distribution model looks good. pp\_check() and loo\_compare() are running, and are pretty convincing that truncated skew normal is probably the best option for the dataset moving forward, though I will reassess as I add covariates. Thank you for all the assistance, greatly appreciated!  
 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/3X/8/b/8b32ba004441cc7e18f78dc6e54e4e22744dc8a4.png)

---

<div class="post-metadata">

**Author:** ![Erin\_Fedewa\_NOAA\_Federal](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/erin_fedewa_noaa_federal/32/18095_2.png) [@Erin\_Fedewa\_NOAA\_Federal](https://discourse.mc-stan.org/u/Erin_Fedewa_NOAA_Federal)\
**Post date:** [March 20, 2024, 7:08am UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/15 "2024-03-20T07:08:31Z")

</div>

@paul.buerkner @avehtari but as a quick follow up to this, it seems that you cannot yet compute posterior draws of the posterior predictive distribution for truncated skew normal models?

```no-highlight
> preds <- posterior_epred(mod5) #posterior draws
Error: posterior_epred values on the respone scale not yet implemented for truncated 'skew_normal' models.

```

---

<div class="post-metadata">

**Author:** ![avehtari](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/avehtari/32/5935_2.png) [@avehtari](https://discourse.mc-stan.org/u/avehtari)\
**Post date:** [March 20, 2024, 7:53am UTC](https://discourse.mc-stan.org/t/choosing-a-sampling-distribution-for-left-skewed-data/34497/16 "2024-03-20T07:53:37Z")

</div>

> [@Erin\_Fedewa\_NOAA\_Federal](#):
>
> it seems that you cannot yet compute posterior draws of the posterior predictive distribution for truncated skew normal models?

Yes, that is what the error message says. There seem to be equations for moments of (at least some version of) skew-normal, so it would be possible to implement this, but I guess @paul.buerkner just didn’t have time and no-one else has asked before.

Meanwhile, maybe `posterior_linpred(., transform = TRUE)` would be useful? It’s not giving the draws for expectation, but the draws of location parameter. If you need `epred`, you could use `posterior_linpred(., transform = TRUE)` and posterior draws of other skew-normal parameters to compute the expectation for each draw yourself using the equations you can find via search engines. If you are not in a hurry, I guess this will get in `brms` eventually