Fitting autoregressive (AR1) model for multivariate binary response

Interfaces brms
1

I have three binary response, y1bin , y2bin and y3bin and one predictor, wealthindex . I have two grouping variables, years and hhid . I want to fit a mixed multivariate model and also fit an AR1 one model. See my formula specification below

`model <- brm(mvbind(y1bin, y2bin, y3bin) ~ 0 + intercept + wealthindex + (0 + 1|g|years) + (0 + 1|p|hhid) 
, autocor = list(cor_ar(form = ~years, p = 1, cov = TRUE), cor_ar(form = ~years, p = 1, cov = TRUE), cor_ar(form = ~years, p = 1, cov = TRUE))
 , data = df , family = list(bernoulli(link = "logit"), bernoulli(link = "logit"), bernoulli(link = "logit")) , prior = priors )`

I tried this but I getting the following error. I will appreciate any help on how to resolve this.
Thanks,

Duplicate declaration of variable, name=sderr; attempt to redeclare as real in parameter; previously declared as real in parameter
error in ‘model5c8a4b5f584c_file5c8a8ab4743’ at line 334, column 21

332: vector<lower=-1,upper=1>[Kar_y2bin] ar_y2bin; // autoregressive effects
333: vector[N_y2bin] zerr_y2bin; // unscaled residuals
334: real<lower=0> sderr; // SD of residuals
^
335: vector[K_y3bin] b_y3bin; // population-level effects
-------------------------------------------------

  • Operating System: ubuntu 18.10
  • brms Version: 2.10.0
2

Thanks! Should be fixed now on github.

3

Thanks for your help. However, I installed the github version of brms but now I am getting a different error, like name=chol_cor . See bellow

Duplicate declaration of variable, name=chol_cor; attempt to redeclare as matrix in transformed parameter; previously declared as matrix in transformed parameter
error in ‘model75bd3cdea9ab_file75bd7792b1cb’ at line 353, column 57

351: vector[N_y1bin] err_y1bin; // actual residuals

352: // cholesky factor of the autocorrelation matrix
353: matrix[max(nobs_tg_y2bin), max(nobs_tg_y2bin)] chol_cor;
^
354: vector[N_y2bin] err_y2bin; // actual residuals

Error in stanc(model_code = paste(program, collapse = “\n”), model_name = model_cppname, :

failed to parse Stan model ‘file75bd7792b1cb’ due to the above error.
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Duplicate declaration of variable, name=chol_cor; attempt to redeclare as matrix in transformed parameter; previously declared as matrix in transformed parameter
error in ‘model75bd54ecf89e_file75bd7792b1cb’ at line 353, column 57

351: vector[N_y1bin] err_y1bin; // actual residuals

352: // cholesky factor of the autocorrelation matrix
353: matrix[max(nobs_tg_y2bin), max(nobs_tg_y2bin)] chol_cor;
^
354: vector[N_y2bin] err_y2bin; // actual residuals

Error in stanc(model_code = paste(program, collapse = “\n”), model_name = model_cppname, :

4

My bad sorry. Can you try it again?

5

Thanks for the fix.
Seems to work now though model still running. I will let you know if it works fine or otherwise.

Thanks

6

Hi Paul,

Although the cor_ar model runs, I am still having troubles in getting simulation based true values. When I use multiple continuous outcomes on the same model (with gausian family since the outcomes are continuous), I get the true values. True values I mean simulation based pre-set values. However, in the binary case, I am not able to capture true values. in fact the ar values shrink and the intercept values (and sd) shoot up. On the other hand if I apply binary responses (the same data) and ignore the cor_ar, I get the true values. I am not sure where the problem is but please could check if everything seems fine from your side or if there something wrong with my code. Below is my model for the binary case

(model <- brm(
     mvbind(y1bin, y2bin, y3bin) ~ 0 + intercept + wealthindex + (1|g|years) + (1|q|hhid)
        , autocor = list(cor_ar(form = ~years|hhid, p = 1, cov = TRUE)
           , cor_ar(form = ~years|hhid, p = 1, cov = TRUE)
           , cor_ar(form = ~years|hhid, p = 1, cov = TRUE)
        )
        , data = df
        , family = list(bernoulli(link = "logit")
           , bernoulli(link = "logit")
           , bernoulli(link = "logit")
        ) 
        , warmup = 1e3
        , iter = 1e4
        , chains = 4
        , cores = parallel::detectCores()
        , control = list(adapt_delta = 0.95)
        , seed = 7777
        , prior = priors
  )

Thank you so much

7

I think the problem is that through AR for bernoulli models we have to add an error term to the linear predictor. However, an error term is not identified for a binary response and hence the model runs into a whole lot of problems. Basically, I don’t think AR structures are currently sensibly possible with bernoulli outcomes. There may be some way to fix the residual SD parameter but this is not yet implemented in brms.