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:
- 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
- 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
- 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
- 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!