The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Hypothesis-Testing

Setup

First, lets load the package:

library(bvarnet)
library(qgraph)

Data

Now, we can load the example data:

data(studentlife)

There is some missing data in the dataset. The models default options handle this by themselves. For a further elaboration on this, you can read vignette("Missing-Data").

Model Estimation

This time, well simply use the default priors, for the ordinal model they are

fit <- bvar(
  id_col = "id",
  time_col = "day",
  y_cols = c("anxious", "calm", "conventional", "critical", "dependable"),
  x_cols = c("sleep_hour"),
  re_cols = NULL,
  re_temporal = FALSE,
  K = 1,
  data = studentlife,
  family = c("ordinal"),
  priors = set_priors(),
  iter = 4000,
  warmup = 1000,
  chains = 4,
  cores = 4,
  seed = 1337
)

Model Output

Again, lets get the model overview and summary:

print(fit)
#> BVAR Network fit
#> ======================================== 
#> Family:      ordinal
#> Outcomes (p): 5 
#> Lags (K):     1 
#> Fixed eff.:   1 
#> Observations: 111 
#> Rhat max:    1.001
#> Divergences: 2  WARNING: check model/priors.
#> Priors:       beta ~ Normal(0, 1), phi ~ Normal(0, 0.5), kappa ~ Normal(0, 2) (all defaults)
#> Total time:  16.2 sec
#> ========================================
summary(fit)
#> BVAR Network Summary
#> ================================================== 
#> Family: ordinal | p=5 | K=1 | n=111
#> Rhat max: 1.001 | Divergences: 2
#>   WARNING: divergent transitions detected — check model/priors.
#> 
#> --- Fixed Effect ---
#>  predictor  outcome      mean   median q5     q95   rhat ess_bulk ess_tail
#>  sleep_hour anxious      -0.005 -0.005 -0.060 0.049 1    15771.43 12695.05
#>  sleep_hour calm         -0.032 -0.032 -0.088 0.022 1    17766.75 12604.63
#>  sleep_hour conventional -0.016 -0.015 -0.076 0.042 1    17485.60 12642.45
#>  sleep_hour critical     -0.008 -0.008 -0.056 0.038 1    14049.46 12084.18
#>  sleep_hour dependable   -0.003 -0.003 -0.060 0.053 1    18991.13 12742.34
#> 
#> 
#> --- Autoregressive ---
#>  predictor         outcome      mean  median q5     q95   rhat  ess_bulk ess_tail
#>  lag1_anxious      anxious      0.069 0.070  -0.260 0.395 1.000 16772.27 12193.66
#>  lag1_calm         calm         0.165 0.166  -0.145 0.478 1.000 12365.66 11825.79
#>  lag1_conventional conventional 0.391 0.390   0.036 0.750 1.001 15129.26 12680.89
#>  lag1_critical     critical     0.238 0.235  -0.020 0.503 1.000 18369.60 11582.80
#>  lag1_dependable   dependable   0.427 0.425   0.129 0.733 1.000 16335.00 11476.86
#> 
#> 
#> --- Cross-lagged ---
#>  predictor         outcome      mean   median q5     q95   rhat ess_bulk ess_tail
#>  lag1_calm         anxious      -0.249 -0.247 -0.574 0.070 1    12852.77 11468.22
#>  lag1_conventional anxious      -0.028 -0.027 -0.370 0.314 1    13704.43 12188.62
#>  lag1_critical     anxious       0.296  0.295  0.017 0.582 1    19347.83 13118.49
#>  lag1_dependable   anxious       0.085  0.084 -0.194 0.365 1    16984.34 12402.25
#>  lag1_anxious      calm         -0.028 -0.030 -0.343 0.290 1    16803.47 12540.98
#>  lag1_conventional calm         -0.199 -0.200 -0.542 0.146 1    13545.69 11926.93
#>  lag1_critical     calm         -0.173 -0.172 -0.453 0.099 1    19113.18 11669.83
#>  lag1_dependable   calm          0.167  0.164 -0.098 0.438 1    16138.68 11770.26
#>  lag1_anxious      conventional -0.313 -0.313 -0.663 0.032 1    16432.09 12229.12
#>  lag1_calm         conventional -0.063 -0.062 -0.396 0.266 1    13903.03 12044.49
#> 
#> ... 10 more rows. Use extract_temporal(fit, effect = "cl") for full output.
#> 
#> --- Threshold ---
#>  predictor               outcome mean   median q5     q95    rhat  ess_bulk ess_tail
#>  kappa(anxious, c1)      —       -1.392 -1.384 -2.106 -0.710 1.000 10840.42 10882.22
#>  kappa(calm, c1)         —       -1.787 -1.745 -2.665 -1.035 1.000 13231.26 11291.59
#>  kappa(conventional, c1) —       -1.553 -1.527 -2.347 -0.847 1.000 14172.12 10875.45
#>  kappa(critical, c1)     —        0.026  0.035 -0.497  0.524 1.001 12451.60 11045.22
#>  kappa(dependable, c1)   —       -2.043 -1.991 -3.108 -1.139 1.001 11713.93 10367.08
#>  kappa(anxious, c2)      —        0.405  0.406 -0.163  0.966 1.000 16633.22 13007.72
#>  kappa(calm, c2)         —       -1.215 -1.204 -1.855 -0.619 1.000 17831.13 12444.68
#>  kappa(conventional, c2) —       -1.079 -1.075 -1.709 -0.463 1.000 16960.92 12615.81
#>  kappa(critical, c2)     —        0.423  0.423 -0.034  0.889 1.000 16313.25 13806.83
#>  kappa(dependable, c2)   —       -1.082 -1.076 -1.759 -0.435 1.000 16347.72 12825.94
#> 
#> ... 10 more rows. Use extract_param(fit, type = "Threshold") for full output.
#> 
#> ==================================================
#> Use extract_param() or extract_param(fit, type = "...") for the full parameter table.
#> Use extract_network_matrix() for the temporal network matrix.

Hypothesis Testing with Bayes Factors

To estimate Bayes factors for a single parameter, or a set of parameters, we can use the bf_table() function. Bayes factors are computed via the Savage–Dickey density ratio (SDDR), which, for a point null \(H_0 : \theta = \theta_0\) nested in \(H_1 : \theta \neq \theta_0\) with shared nuisance parameters, gives

\[ \mathrm{BF}_{01} \;=\; \frac{p(\theta = \theta_0 \mid y)}{p(\theta = \theta_0)}, \]

i.e. the ratio of the posterior to the prior density evaluated at the null value. bf_table() reports \(\mathrm{BF}_{10} = 1 / \mathrm{BF}_{01}\) by approximating the posterior density at \(\theta_0 = 0\) from the MCMC draws.

To perform a hypothesis test on effects involving a specific variable, we can call the bf_table() function, specifying which variable we want to look at using the variable argument. Lets investigate if the hours of sleep the individual got in the previous night (sleep_hour) has an influence on the network structure:

bf_tab <- bf_table(fit, variable = "sleep_hour")
bf_tab
#>                   type  predictor      outcome    BF10
#> 1         Fixed Effect sleep_hour      anxious 0.03386
#> 2         Fixed Effect sleep_hour         calm 0.05187
#> 3         Fixed Effect sleep_hour conventional 0.03760
#> 4         Fixed Effect sleep_hour     critical 0.02951
#> 5         Fixed Effect sleep_hour   dependable 0.03403
#> 6 Fixed Effect (joint) sleep_hour            — 0.00000

The output now shows 6 rows. The five “Fixed Effect” rows show if the the covariate sleep_hour has an influence on each node separately. The “Fixed Effect (join)” row displays a single Bayes factor that in this case tests the hypothesis

H: Is there a difference in the baseline of the network structure that is explained by the variable sleep_hour

As we can see, the BF_{10} = bf_tab[]. Therefore we can conlude that the variable sleep_hour does not influence the baseline level of the network.

Joint Bayes Factors

For more elaborate analyses, or to report more results at once, we can get a table with all single and joint Bayes factors. For this we can call the bf_table(fit) function

bf_tab <- bf_table(fit)
bf_tab
#>                      type         predictor      outcome     BF10
#> 1          Autoregressive      lag1_anxious      anxious  0.42563
#> 2          Autoregressive         lag1_calm         calm  0.54824
#> 3          Autoregressive lag1_conventional conventional  2.23181
#> 4          Autoregressive     lag1_critical     critical  1.01615
#> 5          Autoregressive   lag1_dependable   dependable  6.07838
#> 6  Autoregressive (joint)            all_ar            —  2.49597
#> 7            Cross-lagged         lag1_calm      anxious  0.85207
#> 8            Cross-lagged lag1_conventional      anxious  0.41647
#> 9            Cross-lagged     lag1_critical      anxious  1.52971
#> 10           Cross-lagged   lag1_dependable      anxious  0.38114
#> 11           Cross-lagged      lag1_anxious         calm  0.38559
#> 12           Cross-lagged lag1_conventional         calm  0.66942
#> 13           Cross-lagged     lag1_critical         calm  0.55497
#> 14           Cross-lagged   lag1_dependable         calm  0.54173
#> 15           Cross-lagged      lag1_anxious conventional  1.28950
#> 16           Cross-lagged         lag1_calm conventional  0.40758
#> 17           Cross-lagged     lag1_critical conventional  0.39803
#> 18           Cross-lagged   lag1_dependable conventional  0.36941
#> 19           Cross-lagged      lag1_anxious     critical  0.38522
#> 20           Cross-lagged         lag1_calm     critical  0.33946
#> 21           Cross-lagged lag1_conventional     critical  0.37277
#> 22           Cross-lagged   lag1_dependable     critical  0.49510
#> 23           Cross-lagged      lag1_anxious   dependable  2.31460
#> 24           Cross-lagged         lag1_calm   dependable  0.46823
#> 25           Cross-lagged lag1_conventional   dependable 24.44204
#> 26           Cross-lagged     lag1_critical   dependable  0.73579
#> 27   Cross-lagged (joint)            all_cl            —  0.00333
#> 28       Temporal (joint)           all_phi            —  0.00000
#> 29           Fixed Effect        sleep_hour      anxious  0.03386
#> 30           Fixed Effect        sleep_hour         calm  0.05187
#> 31           Fixed Effect        sleep_hour conventional  0.03760
#> 32           Fixed Effect        sleep_hour     critical  0.02951
#> 33           Fixed Effect        sleep_hour   dependable  0.03403
#> 34   Fixed Effect (joint)        sleep_hour            —  0.00000

The table can be used to easily get a set of Bayes factors. If we have hypotheses concerning the temporal structure, we can call the bf_table(fit) function using the argument type = "temporal".

bf_temp <- bf_table(fit, type = "temporal")
bf_temp
#>               type predictor outcome BF10
#> 1 Temporal (joint)   all_phi       —    0

The following arguments are available as filters for the type argument: “ar”, “cl”, “intercepts”, “fe”, “lag_fe”, and “temporal”.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.