Fitting Gamma model with brms

Interfaces brms
1

Hi to all,

I have problems fitting the following gamma model:

# Library imports

library(data.table)
library(simstudy)
library(brms)
#> Lade nötiges Paket: Rcpp
#> Loading 'brms' package (version 2.9.0). Useful instructions
#> can be found by typing help('brms'). A more detailed introduction
#> to the package is available through vignette('brms_overview').

# Define the data

def <- defData(varname = "male", formula = .5, dist = "binary", id = "eye_id")
def <- defData(def, varname = "oct", formula = .5, dist = "binary")

defC <- defDataAdd(varname = "thickness", 
                   formula = "390 + 8 * male - 10 * oct", 
                   variance = 66^2
)

defC <- defDataAdd(defC, varname = "pressure", 
                   formula = 2.86,
                   variance = .563^2,
                   dist = "gamma",
                   link = "identity")

# Generate the data

set.seed(41948)
dd <- genData(1000, def)
dd <- addCorFlex(dd, defC, rho = .6, corstr = "cs")

dd
#>       eye_id male oct thickness  pressure
#>    1:      1    1   0  461.6095 2.3230065
#>    2:      2    0   1  382.8064 2.9249466
#>    3:      3    1   1  450.2778 1.9641713
#>    4:      4    0   0  368.7375 1.2235662
#>    5:      5    0   0  431.0898 1.8071625
#>   ---                                    
#>  996:    996    1   1  314.1276 1.3455995
#>  997:    997    0   1  413.1119 2.3405988
#>  998:    998    0   0  203.4344 0.7652045
#>  999:    999    1   1  280.7003 0.9592640
#> 1000:   1000    0   0  419.2663 4.4517026

# Check means, sd, and correlation

dd[, .(mean(male), mean(oct))]
#>       V1    V2
#> 1: 0.493 0.502

dd[, .(mean(thickness), sd(thickness)), keyby = .(male, oct)]
#>    male oct       V1       V2
#> 1:    0   0 389.3978 68.68572
#> 2:    0   1 382.5514 69.09854
#> 3:    1   0 403.4099 66.39633
#> 4:    1   1 388.8900 67.21697

dd[, .(mean(pressure), sd(pressure)), keyby = .(male, oct)]
#>    male oct       V1       V2
#> 1:    0   0 2.985553 1.689790
#> 2:    0   1 2.813059 1.590758
#> 3:    1   0 2.990881 1.659774
#> 4:    1   1 2.954744 1.616318

dd[, .(mean(pressure), sd(pressure))]
#>          V1       V2
#> 1: 2.936086 1.638918

dd[, cor(thickness, pressure)]
#> [1] 0.5869474

# Fit Gamma model

fit <- brm(pressure ~ thickness + male + oct, family = Gamma(link = "identity"), data = dd)
#> Compiling the C++ model
#> Start sampling
#> 
#> SAMPLING FOR MODEL '0b1a6aa24224242a80a515bd48d206ca' NOW (CHAIN 1).
#> Chain 1: Rejecting initial value:
#> Chain 1:   Error evaluating the log probability at the initial value.
#> Chain 1: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.071384, but must be > 0!  (in 'modeldeea13a01997_0b1a6aa24224242a80a515bd48d206ca' at line 37)
#> 
#> Chain 1: Rejecting initial value:
#> Chain 1:   Error evaluating the log probability at the initial value.
#> Chain 1: Exception: gamma_lpdf: Inverse scale parameter[2] is -0.164441, but must be > 0!  (in 'modeldeea13a01997_0b1a6aa24224242a80a515bd48d206ca' at line 37)
#> 
#> Chain 1: Rejecting initial value:
#> Chain 1:   Error evaluating the log probability at the initial value.
#> Chain 1: Exception: gamma_lpdf: Inverse scale parameter[2] is -0.0536331, but must be > 0!  (in 'modeldeea13a01997_0b1a6aa24224242a80a515bd48d206ca' at line 37)
#> 
#> Chain 1: Rejecting initial value:
#> Chain 1:   Error evaluating the log probability at the initial value.
#> Chain 1: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.0102431, but must be > 0!  (in 'modeldeea13a01997_0b1a6aa24224242a80a515bd48d206ca' at line 37)
#> 
#> Chain 1: Rejecting initial value:
#> Chain 1:   Error evaluating the log probability at the initial value.
#> Chain 1: Exception: gamma_lpdf: Inverse scale parameter[2] is -0.981855, but must be > 0!  (in 'modeldeea13a01997_0b1a6aa24224242a80a515bd48d206ca' at line 37)
#> 
#> Chain 1: Rejecting initial value:
#> Chain 1:   Error evaluating the log probability at the initial value.
#> Chain 1: Exception: gamma_lpdf: Inverse scale parameter[1] is -0.0306889, but must be > 0!  (in 'modeldeea13a01997_0b1a6aa24224242a80a515bd48d206ca' at line 37)
.
.
.
#> 
#> Chain 1: 
#> Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
#> Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
#> Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.
#> character(0)
#> error occurred during calling the sampler; sampling not done

Created on 2019-09-02 by the reprex package (v0.2.1)

I tried setting the initial values to “0”, but this does not help understandably (according to the error message).

How can I start the sampler?

Thanks and kind regards,
Pascal

Please also provide the following information in addition to your question:

  • Operating System: macOS Sierra 10.12.6
  • brms Version: 2.9.0
2

Hi, I suspect the use of the identity link is the problem here. Combined with that fact that a linear predictor can take on negative values (including at initialisation), you’re allowing the mean (or scale, not sure how it as implemented in brms) of the gamma distribution to be negative. Use a log link function to transform the linear predictor, or define a strictly positive non-linear predictor.

2 Likes
3

Ah, I see! Thank you for the clarification!