Issue with initial values / constraints

[gowerc]

This happens for most models but taking the Claret Bruno as an example.

When trying to run a model I quite often get something like this:

Chain 3 Rejecting initial value:
Chain 3   Log probability evaluates to log(0), i.e. negative infinity.
Chain 3   Stan can't start sampling from this initial value.
Chain 3 User-specified initialization failed.
Chain 3  Try specifying new initial values, using partially specialized initialization, reducing the range of constrained values, or reparameterizing the model.
Chain 3 Initialization failed.

However when manually inspecting the initial values they all appear to be completely fine ...

> initialValues(jm, n_chains = 1)
[[1]]
[[1]]$lm_clbr_mu_b
[1] 4.310798

[[1]]$lm_clbr_mu_g
[1] -0.1743525

[[1]]$lm_clbr_mu_c
[1] -1.024193

[[1]]$lm_clbr_mu_p
[1] 0.7080146

[[1]]$lm_clbr_omega_b
[1] 0.1069801

[[1]]$lm_clbr_omega_g
[1] 0.1352902

[[1]]$lm_clbr_omega_c
[1] 0.2290127

[[1]]$lm_clbr_omega_p
[1] 0.3337157

[[1]]$lm_clbr_sigma
[1] 0.6920734

[[1]]$lm_clbr_ind_b
[1] 60.55209

[[1]]$lm_clbr_ind_g
[1] 1.006471

[[1]]$lm_clbr_ind_c
[1] 0.3466599

[[1]]$lm_clbr_ind_p
[1] 2.778457

I double checked all the constraints and parameterisations and they all appear to look fine to me...

Weirdly if I set the minimum value to be .Machine$double.eps^(1/3) (I was still getting the error with sqrt) instead of 0 I no longer get this error though I do get alternative warnings instead:

Chain 1 Iteration:    1 / 4200 [  0%]  (Warmup) 
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: lognormal_lpdf: Scale parameter[1] is inf, but must be positive finite! (in '/tmp/RtmpgVZ2BH/model-fabc8850f.stan', line 332, column 4 to column 103)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1 
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: lognormal_lpdf: Scale parameter[1] is inf, but must be positive finite! (in '/tmp/RtmpgVZ2BH/model-fabc8850f.stan', line 332, column 4 to column 103)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1 

Most Stan programs I look at just have the minimum value set at 0 so I'm not sure why this is causing such problems for us :(

Full code for example:

    sim_params <- list(
        sigma = 1,
        mu_b = log(60),
        mu_g = log(c(0.9, 1.1)),
        mu_c = log(c(0.45, 0.35)),
        mu_p = log(c(2.4, 1.8)),
        omega_b = 0.1,
        omega_g = c(0.1, 0.15),
        omega_c = c(0.2, 0.25),
        omega_p = c(0.3, 0.35),
        link_ttg = 0,
        link_dsld = 0,
        link_identity = 0,
        link_growth = 0,
        lambda = 0.5,
        lambda_cen = 1 / 9000,
        beta_cat_b = -0.1,
        beta_cat_c = 0.5,
        beta_cont = 0.3
    )

    set.seed(6128)
    ## Generate Test data with known parameters
    jlist <- SimJointData(
        design = list(
            SimGroup(170, "Arm-A", "Study-X"),
            SimGroup(170, "Arm-B", "Study-X")
        ),
        longitudinal = SimLongitudinalClaretBruno(
            times = c(1, 50, 100, 150, 200, 250, 300, 400, 500, 600, 700, 800, 900, 1000) / 365,
            sigma = sim_params$sigma,
            mu_b = sim_params$mu_b,
            mu_g = sim_params$mu_g,
            mu_c = sim_params$mu_c,
            mu_p = sim_params$mu_p,
            omega_b = sim_params$omega_b,
            omega_g = sim_params$omega_g,
            omega_c = sim_params$omega_c,
            omega_p = sim_params$omega_p,
            link_ttg = sim_params$link_ttg,
            link_dsld = sim_params$link_dsld,
            link_identity = sim_params$link_identity,
            link_growth = sim_params$link_growth,
            scaled_variance = FALSE
        ),
        survival = SimSurvivalExponential(
            time_max = 5,
            time_step = 1 / 365,
            lambda = sim_params$lambda,
            lambda_cen = 1 / 9000,
            beta_cat = c(
                "A" = 0,
                "B" = sim_params$beta_cat_b,
                "C" = sim_params$beta_cat_c
            ),
            beta_cont = sim_params$beta_cont
        ),
        .silent = TRUE
    )

    jm <- JointModel(
        longitudinal = LongitudinalClaretBruno(
            mu_b = prior_normal(sim_params$mu_b, 0.4),
            mu_g = prior_normal(mean(sim_params$mu_g), 0.4),
            mu_c = prior_normal(mean(sim_params$mu_c), 0.4),
            mu_p = prior_normal(mean(sim_params$mu_p), 0.4),
            omega_b = prior_lognormal(log(mean(sim_params$omega_b)), 0.4),
            omega_g = prior_lognormal(log(mean(sim_params$omega_g)), 0.4),
            omega_c = prior_lognormal(log(mean(sim_params$omega_c)), 0.4),
            omega_p = prior_lognormal(log(mean(sim_params$omega_p)), 0.4),
            sigma = prior_lognormal(log(mean(sim_params$sigma)), 0.4),
            centred = TRUE,
            scaled_variance = FALSE
        )
    )

    jdat <- DataJoint(
        subject = DataSubject(
            data = jlist@survival,
            subject = "subject",
            arm = "arm",
            study = "study"
        ),
        longitudinal = DataLongitudinal(
            data = jlist@longitudinal,
            formula = sld ~ time,
            threshold = 2
        )
    )

    ## Sample from JointModel
    set.seed(213)
    mp <- sampleStanModel(
                jm,
                data = jdat,
                iter_sampling = 3000,
                iter_warmup = 1200,
                chains = 5,
                parallel_chains = 5
            )