Differences between brms and GLS estimates with lagged DV and autoregressive error

I’m fitting some models both in Stan (using brms) and GLS. One would expect that both lead to (almost the) same results which is true for most of the cases, but I have trouble with replicating one model that includes an AR (1) error and a lagged dependent variable. I think the issue is due to the lagged dependent variable in combination with the AR (1) error.

I think the results make my problem much clearer. So here are some models comparing Stan and GLS:

  1. Model without lagged DV and without AR (1) error

Stan:

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: adawareness_logit ~ cv_tv_ttsd_level_log + cv_print_ttsd_level_log + cv_digital_ttsd_level_log + cv_radio_ttsd_level_log + cv_other_ttsd_level_log + buzz_rescaled_log_lag1 + cv_comp_share_log + cv_comp_share_log_lag1 + google_trend_log + q2 + q3 + q4 + time 
   Data: df_analysis[["mercedes-benz"]] (Number of observations: 199) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                    -1.61      0.53    -2.68    -0.61 1.00     3951     3196
cv_tv_ttsd_level_log          0.01      0.01    -0.00     0.02 1.00     3384     3082
cv_print_ttsd_level_log       0.01      0.01    -0.02     0.04 1.00     4029     2443
cv_digital_ttsd_level_log    -0.00      0.01    -0.03     0.02 1.00     3801     2889
cv_radio_ttsd_level_log       0.01      0.01    -0.00     0.02 1.00     4020     2744
cv_other_ttsd_level_log      -0.01      0.01    -0.03     0.01 1.00     4578     3013
buzz_rescaled_log_lag1        0.93      0.34     0.24     1.58 1.00     3836     3115
cv_comp_share_log            -0.25      0.48    -1.17     0.68 1.00     2941     2692
cv_comp_share_log_lag1        1.14      0.31     0.51     1.75 1.00     3555     2850
google_trend_log              0.03      0.12    -0.21     0.28 1.00     4055     2932
q2                           -0.01      0.03    -0.07     0.04 1.00     2735     2876
q3                           -0.02      0.03    -0.07     0.04 1.00     2650     2845
q4                            0.03      0.03    -0.04     0.09 1.00     2131     2807
time                         -0.00      0.00    -0.00    -0.00 1.00     3612     2902

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.12      0.01     0.11     0.14 1.00     4058     2859

Samples were drawn 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).

GLS:

Generalized least squares fit by REML
  Model: adawareness_logit ~ cv_tv_ttsd_level_log + cv_print_ttsd_level_log +      cv_digital_ttsd_level_log + cv_radio_ttsd_level_log + cv_other_ttsd_level_log +      buzz_rescaled_log_lag1 + cv_comp_share_log + cv_comp_share_log_lag1 +      google_trend_log + q2 + q3 + q4 + time 
  Data: data.frame(df_analysis[["mercedes-benz"]] %>% dplyr::select(adawareness_logit_lag1,      adawareness_logit, cv_tv_ttsd_level_log, cv_print_ttsd_level_log,      cv_digital_ttsd_level_log, cv_radio_ttsd_level_log, cv_other_ttsd_level_log,      buzz_rescaled_log_lag1, cv_comp_share_log, cv_comp_share_log_lag1,      google_trend_log, q2, q3, q4, time)) 
        AIC       BIC   logLik
  -167.9525 -119.6472 98.97627

Coefficients:
                               Value Std.Error   t-value p-value
(Intercept)               -1.6109716 0.5193790 -3.101727  0.0022
cv_tv_ttsd_level_log       0.0077812 0.0056776  1.370504  0.1722
cv_print_ttsd_level_log    0.0102662 0.0134993  0.760496  0.4479
cv_digital_ttsd_level_log -0.0025597 0.0134117 -0.190855  0.8488
cv_radio_ttsd_level_log    0.0077414 0.0063597  1.217244  0.2251
cv_other_ttsd_level_log   -0.0087907 0.0093418 -0.941005  0.3479
buzz_rescaled_log_lag1     0.9263150 0.3284313  2.820423  0.0053
cv_comp_share_log         -0.2482575 0.4788388 -0.518457  0.6048
cv_comp_share_log_lag1     1.1312420 0.2985923  3.788584  0.0002
google_trend_log           0.0313039 0.1210452  0.258613  0.7962
q2                        -0.0136007 0.0278711 -0.487986  0.6261
q3                        -0.0193100 0.0280559 -0.688270  0.4921
q4                         0.0264391 0.0316059  0.836526  0.4039
time                      -0.0015821 0.0001753 -9.026159  0.0000
  1. Model without lagged DV including AR (1) error

Stan:

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: adawareness_logit ~ cv_tv_ttsd_level_log + cv_print_ttsd_level_log + cv_digital_ttsd_level_log + cv_radio_ttsd_level_log + cv_other_ttsd_level_log + buzz_rescaled_log_lag1 + cv_comp_share_log + cv_comp_share_log_lag1 + google_trend_log + q2 + q3 + q4 + time 
         autocor ~ ar(p = 1)
   Data: df_analysis[["mercedes-benz"]] (Number of observations: 199) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Correlation Structures:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
ar[1]     0.32      0.08     0.17     0.47 1.00     2864     2840

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                    -1.66      0.60    -2.84    -0.49 1.00     3502     2978
cv_tv_ttsd_level_log          0.01      0.01    -0.01     0.02 1.00     3121     3024
cv_print_ttsd_level_log       0.01      0.01    -0.02     0.03 1.00     3152     2818
cv_digital_ttsd_level_log     0.00      0.01    -0.02     0.03 1.00     3201     2946
cv_radio_ttsd_level_log       0.01      0.01     0.00     0.03 1.00     3905     2551
cv_other_ttsd_level_log      -0.01      0.01    -0.02     0.01 1.00     4326     3043
buzz_rescaled_log_lag1        0.89      0.42     0.04     1.73 1.00     3046     2658
cv_comp_share_log            -0.19      0.49    -1.19     0.78 1.00     2254     2586
cv_comp_share_log_lag1        0.98      0.29     0.40     1.55 1.00     3800     3062
google_trend_log              0.04      0.14    -0.23     0.32 1.00     3363     2850
q2                           -0.01      0.04    -0.08     0.06 1.00     2289     2374
q3                           -0.01      0.04    -0.09     0.06 1.00     2140     2562
q4                            0.02      0.04    -0.06     0.09 1.00     2302     2444
time                         -0.00      0.00    -0.00    -0.00 1.00     4050     2960

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.12      0.01     0.11     0.13 1.00     3640     2825

Samples were drawn 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).

GLS:

Generalized least squares fit by REML
  Model: adawareness_logit ~ cv_tv_ttsd_level_log + cv_print_ttsd_level_log +      cv_digital_ttsd_level_log + cv_radio_ttsd_level_log + cv_other_ttsd_level_log +      buzz_rescaled_log_lag1 + cv_comp_share_log + cv_comp_share_log_lag1 +      google_trend_log + q2 + q3 + q4 + time 
  Data: data.frame(df_analysis[["mercedes-benz"]] %>% dplyr::select(adawareness_logit_lag1,      adawareness_logit, cv_tv_ttsd_level_log, cv_print_ttsd_level_log,      cv_digital_ttsd_level_log, cv_radio_ttsd_level_log, cv_other_ttsd_level_log,      buzz_rescaled_log_lag1, cv_comp_share_log, cv_comp_share_log_lag1,      google_trend_log, q2, q3, q4, time)) 
       AIC       BIC   logLik
  -183.683 -132.1573 107.8415

Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):
      Phi 
0.3190317 

Coefficients:
                               Value Std.Error   t-value p-value
(Intercept)               -1.6790615 0.5940652 -2.826392  0.0052
cv_tv_ttsd_level_log       0.0070462 0.0061614  1.143603  0.2543
cv_print_ttsd_level_log    0.0073586 0.0133033  0.553142  0.5808
cv_digital_ttsd_level_log  0.0010626 0.0134279  0.079131  0.9370
cv_radio_ttsd_level_log    0.0132868 0.0060888  2.182167  0.0304
cv_other_ttsd_level_log   -0.0052724 0.0091015 -0.579292  0.5631
buzz_rescaled_log_lag1     0.8941425 0.4137180  2.161237  0.0320
cv_comp_share_log         -0.1846266 0.4870345 -0.379083  0.7051
cv_comp_share_log_lag1     0.9765684 0.2881332  3.389295  0.0009
google_trend_log           0.0481074 0.1369736  0.351216  0.7258
q2                        -0.0093892 0.0354115 -0.265145  0.7912
q3                        -0.0118310 0.0361723 -0.327074  0.7440
q4                         0.0173099 0.0389121  0.444846  0.6570
time                      -0.0016105 0.0002395 -6.724811  0.0000
  1. Model with lagged DV and without AR (1) error

Stan:

Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: adawareness_logit ~ adawareness_logit_lag1 + cv_tv_ttsd_level_log + cv_print_ttsd_level_log + cv_digital_ttsd_level_log + cv_radio_ttsd_level_log + cv_other_ttsd_level_log + buzz_rescaled_log_lag1 + cv_comp_share_log + cv_comp_share_log_lag1 + google_trend_log + q2 + q3 + q4 + time 
   Data: df_analysis[["mercedes-benz"]] (Number of observations: 199) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                    -1.20      0.51    -2.20    -0.21 1.00     3931     3002
adawareness_logit_lag1        0.29      0.07     0.15     0.42 1.00     4169     2593
cv_tv_ttsd_level_log          0.01      0.01    -0.00     0.02 1.00     3358     3326
cv_print_ttsd_level_log       0.01      0.01    -0.02     0.04 1.00     4147     3229
cv_digital_ttsd_level_log    -0.00      0.01    -0.03     0.02 1.00     3659     3319
cv_radio_ttsd_level_log       0.01      0.01    -0.00     0.02 1.00     4587     3117
cv_other_ttsd_level_log      -0.01      0.01    -0.03     0.01 1.00     4589     3129
buzz_rescaled_log_lag1        0.66      0.33     0.03     1.30 1.00     3363     2691
cv_comp_share_log            -0.20      0.47    -1.10     0.71 1.00     2774     2692
cv_comp_share_log_lag1        0.87      0.29     0.30     1.45 1.00     3977     2679
google_trend_log              0.03      0.12    -0.20     0.26 1.00     3943     2809
q2                           -0.01      0.03    -0.07     0.04 1.00     2361     2877
q3                           -0.01      0.03    -0.06     0.04 1.00     2700     3042
q4                            0.02      0.03    -0.04     0.08 1.00     2459     2795
time                         -0.00      0.00    -0.00    -0.00 1.00     4315     3407

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.12      0.01     0.11     0.13 1.00     3906     2621

Samples were drawn 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).

GLS:

Generalized least squares fit by REML
  Model: adawareness_logit ~ adawareness_logit_lag1 + cv_tv_ttsd_level_log +      cv_print_ttsd_level_log + cv_digital_ttsd_level_log + cv_radio_ttsd_level_log +      cv_other_ttsd_level_log + buzz_rescaled_log_lag1 + cv_comp_share_log +      cv_comp_share_log_lag1 + google_trend_log + q2 + q3 + q4 +      time 
  Data: data.frame(df_analysis[["mercedes-benz"]] %>% dplyr::select(adawareness_logit_lag1,      adawareness_logit, cv_tv_ttsd_level_log, cv_print_ttsd_level_log,      cv_digital_ttsd_level_log, cv_radio_ttsd_level_log, cv_other_ttsd_level_log,      buzz_rescaled_log_lag1, cv_comp_share_log, cv_comp_share_log_lag1,      google_trend_log, q2, q3, q4, time)) 
        AIC       BIC   logLik
  -180.3339 -128.8949 106.1669

Coefficients:
                               Value Std.Error   t-value p-value
(Intercept)               -1.2147028 0.5044561 -2.407946  0.0170
adawareness_logit_lag1     0.2873558 0.0663221  4.332733  0.0000
cv_tv_ttsd_level_log       0.0093671 0.0054354  1.723345  0.0865
cv_print_ttsd_level_log    0.0105311 0.0128943  0.816723  0.4151
cv_digital_ttsd_level_log -0.0047696 0.0128206 -0.372024  0.7103
cv_radio_ttsd_level_log    0.0094874 0.0060880  1.558380  0.1209
cv_other_ttsd_level_log   -0.0080605 0.0089246 -0.903172  0.3676
buzz_rescaled_log_lag1     0.6529229 0.3199911  2.040441  0.0427
cv_comp_share_log         -0.2030397 0.4574921 -0.443810  0.6577
cv_comp_share_log_lag1     0.8657742 0.2917139  2.967888  0.0034
google_trend_log           0.0331197 0.1156196  0.286454  0.7749
q2                        -0.0132785 0.0266218 -0.498783  0.6185
q3                        -0.0098434 0.0268871 -0.366102  0.7147
q4                         0.0195172 0.0302313  0.645596  0.5193
time                      -0.0010964 0.0002015 -5.441858  0.0000
  1. Model with lagged DV and with AR (1) error

Stan:

Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: adawareness_logit ~ adawareness_logit_lag1 + cv_tv_ttsd_level_log + cv_print_ttsd_level_log + cv_digital_ttsd_level_log + cv_radio_ttsd_level_log + cv_other_ttsd_level_log + buzz_rescaled_log_lag1 + cv_comp_share_log + cv_comp_share_log_lag1 + google_trend_log + q2 + q3 + q4 + time 
         autocor ~ ar(p = 1)
   Data: df_analysis[["mercedes-benz"]] (Number of observations: 199) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Correlation Structures:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
ar[1]     0.14      0.25    -0.28     0.57 1.01      903     2390

Population-Level Effects: 
                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                    -1.43      0.65    -2.78    -0.28 1.00     1799     2221
adawareness_logit_lag1        0.16      0.22    -0.25     0.51 1.01      918     2087
cv_tv_ttsd_level_log          0.01      0.01    -0.01     0.02 1.00     2526     2175
cv_print_ttsd_level_log       0.01      0.01    -0.02     0.03 1.00     2782     2846
cv_digital_ttsd_level_log    -0.00      0.01    -0.03     0.03 1.00     2237     2737
cv_radio_ttsd_level_log       0.01      0.01    -0.00     0.02 1.00     2662     2851
cv_other_ttsd_level_log      -0.01      0.01    -0.02     0.01 1.00     3172     2920
buzz_rescaled_log_lag1        0.76      0.41     0.03     1.66 1.00     2163     1946
cv_comp_share_log            -0.16      0.48    -1.09     0.76 1.00     2749     2775
cv_comp_share_log_lag1        0.89      0.31     0.28     1.49 1.00     2863     2726
google_trend_log              0.04      0.13    -0.20     0.29 1.00     3302     2633
q2                           -0.01      0.03    -0.08     0.05 1.00     2078     2297
q3                           -0.01      0.03    -0.08     0.05 1.00     2123     2350
q4                            0.01      0.04    -0.06     0.08 1.00     2078     2068
time                         -0.00      0.00    -0.00    -0.00 1.01     1122     2485

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.12      0.01     0.11     0.13 1.00     2510     2769

Samples were drawn 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).

GLS:

Generalized least squares fit by REML
  Model: adawareness_logit ~ adawareness_logit_lag1 + cv_tv_ttsd_level_log +      cv_print_ttsd_level_log + cv_digital_ttsd_level_log + cv_radio_ttsd_level_log +      cv_other_ttsd_level_log + buzz_rescaled_log_lag1 + cv_comp_share_log +      cv_comp_share_log_lag1 + google_trend_log + q2 + q3 + q4 +      time 
  Data: data.frame(df_analysis[["mercedes-benz"]] %>% dplyr::select(adawareness_logit_lag1,      adawareness_logit, cv_tv_ttsd_level_log, cv_print_ttsd_level_log,      cv_digital_ttsd_level_log, cv_radio_ttsd_level_log, cv_other_ttsd_level_log,      buzz_rescaled_log_lag1, cv_comp_share_log, cv_comp_share_log_lag1,      google_trend_log, q2, q3, q4, time)) 
        AIC       BIC   logLik
  -178.5605 -123.9066 106.2802

Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):
        Phi 
-0.09342653 

Coefficients:
                               Value Std.Error   t-value p-value
(Intercept)               -1.0993852 0.4752762 -2.313150  0.0218
adawareness_logit_lag1     0.3615251 0.0641852  5.632527  0.0000
cv_tv_ttsd_level_log       0.0099741 0.0051739  1.927761  0.0554
cv_print_ttsd_level_log    0.0110537 0.0125873  0.878163  0.3810
cv_digital_ttsd_level_log -0.0055290 0.0124520 -0.444025  0.6575
cv_radio_ttsd_level_log    0.0079168 0.0059922  1.321183  0.1881
cv_other_ttsd_level_log   -0.0082435 0.0087033 -0.947160  0.3448
buzz_rescaled_log_lag1     0.5908379 0.2967428  1.991078  0.0480
cv_comp_share_log         -0.1862249 0.4459772 -0.417566  0.6768
cv_comp_share_log_lag1     0.8070543 0.2919039  2.764794  0.0063
google_trend_log           0.0303427 0.1089460  0.278511  0.7809
q2                        -0.0147119 0.0246686 -0.596380  0.5517
q3                        -0.0095221 0.0249191 -0.382120  0.7028
q4                         0.0172693 0.0282832  0.610586  0.5422
time                      -0.0009659 0.0001879 -5.141090  0.0000

While for model 1 to 3 the estimates are pretty similar, there is a significant difference for model 4 including both autoregressive effects of the dependent variable and the error. Is there anything I’m missing here in the brms model? I’m using the following syntax:

model.stan <- brm(bf(adawareness_logit ~adawareness_logit_lag1+
                                       cv_tv_ttsd_level_log+
                                       cv_print_ttsd_level_log+
                                       cv_digital_ttsd_level_log+
                                       cv_radio_ttsd_level_log+
                                       cv_other_ttsd_level_log+
                                       buzz_rescaled_log_lag1+
                                       cv_comp_share_log+
                                       cv_comp_share_log_lag1+
                                       google_trend_log+
                                       q2+
                                       q3+
                                       q4+
                                       time, autocor= ~ar(p=1)),  
                                  data=df_analysis, 
                                  cores=4,
                                  control=list(adapt_delta=0.95))

Thanks!

One thing I would suspect is that you end up having too little data to constrain all of the parameters, or that you’ve introduced a collinearity in your predictors (using lagged DV could IMHO have a somewhat similar effect to an ar(1) term, making the coefficients somewhat interchangeable).

This seems to be supported by the quite wide posterior intervals for both the ar and adawareness_logit_lag1 terms. You might want to look into pairs plot for those two parameters (see pairs.brmsfit).

If that’s indeed the case than the difference between GLS and and brms is to be expected - the maximum likelihood estimate becomes unstable (sensitive to small changes in the data), while the brms is able to capture that the parameters are not well informed and gives you those wide intervals. Note that the GLS results are well within the posterior 95% credible intervals reported by brms.

Does that make sense?