How to cite the dataset in APA style

Van de Schoot, R. (2020). PhD-delay Dataset for Online Stats Training [Data set]. Zenodo. https://doi.org/10.5281/zenodo.3999424

1. Setting a goal: A research question

• Earning a Ph.D. degree is a long and enduring process. Many factors are involved in timely completion, abrupt termination, or some delay of the Ph.D. life. For the current tutorial, we examine how age is related to the Ph.D. delay.
• Let’s think about the relationship between age and Ph.D. delay. Do you think they are linearly related? Are there any different ideas such that the age is non-linearly associated with the delay? We expect the relationship between the completion time and age are non-linear. This might be because, at a certain point in your life (say mid-thirties), family life takes up more of your time compared to when you are in the twenties.
• To that end, we examine whether the age of the Ph.D. recipients predicts a delay in their Ph.D. projects. The regression model is consequently the one we should adopt to answer the research question. We use the difference variable (B3_difference_extra) as the dependent variable and the age (E22_Age) and the age-squared (E22_Age_Squared) as the independent variables for the regression model.

2. Understanding the data

3. Adding the JAGS button

• To follow the WAMBS Checklist in JASP, we call our helper, JAGS. We can invite JAGS using the plus button.
• Click the plus button on the top right of the JASP screen -> Check JAGS

• Can you see that the JAGS button has been added to the top bar? We are ready to go to the first stage of the WAMBS Checklist!

Point 1. Do you understand the priors?

• We first need to think about the prior distributions and their hyperparameters for our regression model. Hyperparameters are the parameters for the prior distributions.
• The regression model for the current example can be written as: \(Difference = b_{0} + beta_{Age} * X_{Age} + beta_{Age-squared} * X_{Age-squared}\)
• In our regression model, parameters are the intercept, two regression coefficients, and residual variance. Thus, we should specify the prior distributions for the intercept, two regression coefficients, and residual variance, respectively.

Point 2. Does the trace plot exhibit convergence?

• First, we only use 10 burn-in samples and take 50 samples. In other words, the number of iterations is 60.
• Click Advanced in the control panel -> Type 50 in No. samples (i.e., the number of samples) and 10 in No. burnin samples (i.e., the number of burn-in samples)
• Go again to the code that we have pasted. Press Ctrl and Enter on the keyboard together to apply the code. With this action, you are running the model that we specified.

• To see the trace plots for the parameters, move intercept, beta1, beta2, and tau from the Parameters in the model section to the Show results for these parameters section.

• In the control panel, check Trace in the Plots tab.

• With the trace plots, we can examine whether the parameter estimates from multiple chains are stable. If the trace plots are indistinguishable, this indicates that we do not need to increase the number of iterations.
• According to the trace plots of the four parameters, we see that the chains have not converged into one another yet. The ideal trace plot should exhibit one fat caterpillar. Thus, we need to increase the number of iterations.
• In addition to the trace plots, the convergence diagnostics tell us whether the chains have converged. One such convergence diagnostic is the Gelman and Rubin diagnostic. In the output panel, JASP provides the Gelman and Rubin diagnostic. Where? The Rhat column of the MCMC Summary table!

• If the multiple chains have converged neatly, the point estimate (Point est.) and the upper limit of the confidence interval (Upper CI; also referred to as the upper confidence limit) should be close to 1. According to the MCMC Summary table, the point estimate and the upper confidence limit of beta1, for example, is not close to 1. Therefore, we decide to increase the number of iterations.
• This time, we use 25000 as the number of burn-in samples and 50000 as the number of samples to take. In total, we use 75000 iterations.

• From the trace plots, we see the fat caterpillar. Let’s also check the Gelman and Rubin diagnostic from the Rhat column of the MCMC Summary table.

• The point estimate and the upper confidence limit are 1.000 for all parameters.
• Therefore, the trace plots and the Gelman and Rubin diagnostics exhibit convergence.

Point 3. Does convergence remain after doubling the number of iterations?

• Following the recommendation in the WAMBS Checklist, we double the number of iterations. This is because the convergence we saw might be a local convergence. Type 50000 as the number of burn-in samples and 100000 as the number of samples. In other words, we use 150000 iterations.

Point 4. Does the posterior distribution histogram have enough information?

• Before drawing the posterior distribution histogram, uncheck Trace in the Plots tab to save the processing time.
• Check Histogram to draw the posterior distribution histograms for parameters.

• The histograms look smooth and have no gaps or other abnormalities. This indicates the additional iterations are not necessary. Of course, you can increase the number of iterations if you are not satisfied.
  • Please note that the posterior distributions do not have to be symmetrical, but they seem to be in the current example.
• Let’s compare these histograms with those based on the first analysis with very few iterations (with 10 burn-in samples and 50 additionally drawn samples).
• Type 50 in No. samples and 10 in No. burnin samples to get the histograms with 60 iterations.

• The difference becomes clear. For example, take a look at the histogram of beta1. The histogram is not smooth and has a gap. The additional iterations are thus necessary in the 60-iteration example.
• Please type 25000 as the number of burn-in samples and 50000 as the number of samples again to proceed with the next point. Uncheck Histogram as well to save the processing time.

Point 5. Do the chains exhibit a strong degree of autocorrelation?

• One evidence of convergence is no autocorrelation within chains. We can inspect if the autocorrelations exist with the autocorrelation plot.
• To draw the autocorrelation plot, check Autocorrelation in the Plots tab.

• The autocorrelation plots show nearly zero autocorrelations within the initial lags.
• If you detect the strong autocorrelation even after a few lags, you have to run the analysis with a lot of samples to identify the whole parameter space (see Link & Eaton, 2012 for more information on autocorrelation check).

Point 6. Do the posterior distributions make substantive sense?

• By plotting the posterior distributions, we can examine whether they are unimodal (i.e., one peak), whether they are centered around one value, whether they give a realistic estimate, and whether they make substantive sense compared to our prior beliefs (i.e., prior distributions).
• In the Plots tab, uncheck Autocorrelation and check Density.

• For the intercept, the highest peak is around -35. This makes substantive sense because the age of 0 for Ph.D. recipients is impossible.
• For both regression coefficients (i.e., beta1 and beta2), the posterior distributions show plausible ranges for the value of regression coefficients.
• For the precision (i.e., tau), the ranges cover positive values. This makes sense because the precision is the inverse of the variance, which cannot be negative.
• Before going to the next stage, uncheck Density to save the processing time.

Point 7. Do different specifications of the variance priors influence the results?

• To understand the effects of different prior specifications of precision on the results, let’s use different hyperparameter values for precision.
• The previous setting used \(Gamma(.5, .5)\) as the prior for the precision. This time, we use \(Gamma(.01, .01)\) as a prior.

model{
    # likelihood
    for (i in 1:333){
    mu[i] <- intercept + beta1 * E22_Age[i] + beta2 * E22_Age_Squared[i]
    B3_difference_extra[i] ~ dnorm(mu[i],tau)
    }

    # priors
    intercept ~ dnorm(-35,1/20)
    beta1 ~ dnorm(0.8,1/5)
    beta2 ~ dnorm(0,1/10)
    tau ~ dgamma(0.01,0.01)
}

• Don’t forget to press Ctrl and Enter together to apply the changes!

• The following table presents the comparison of posterior means between different prior specifications alongside the relative bias.

• The relative bias for the intercept is almost 0 whereas those for other parameters are 0. Therefore, the results are robust.

Point 8. Is there a notable effect of the prior when compared with non-informative priors?

• To compare the effect of the non-informative priors with that of the informative priors, we need to run the analysis with the non-informative priors. We then compare the output with what we have done so far with the informative priors.
• For the non-informative priors, we borrow the default priors that the R package blavaan (Merkle & Rosseel, 2018) uses with the JAGS target.
  • For more information, we refer readers to JASP for Bayesian analyses with informative priors (using JAGS) and Bayesian Regression in Blavaan (using JAGS).
• Following the default prior settings that blavaan employs, we are going to specify the noninformative priors as:
  • \(intercept \sim \mathcal{N}(0, 1e-3)\)
  • \(\beta_1 \sim \mathcal{N}(0, 1e-2)\)
  • \(\beta_2 \sim \mathcal{N}(0, 1e-2)\)
  • \(tau \sim Gamma(1, 0.5)\)
• For your convenience, we provide the code that you can copy and paste!

model{
    # likelihood
    for (i in 1:333){
    mu[i] <- intercept + beta1 * E22_Age[i] + beta2 * E22_Age_Squared[i]
    B3_difference_extra[i] ~ dnorm(mu[i],tau)
    }

    # priors
    intercept ~ dnorm(0,1e-3)
    beta1 ~ dnorm(0,1e-2)
    beta2 ~ dnorm(0,1e-2)
    tau ~ dgamma(1,0.5)
}

• Press Ctrl and Enter together to apply the changes.

• The subsequent table summarizes the comparison of posterior means between non-informative and informative priors alongside the relative bias.

• Except for the precision, the informative priors have influences on the posterior means ranging from 9.52% to 12.52% compared to the non-informative priors.
• Importantly, you have to keep in mind where the priors come from. If your informative priors come from accurate and reliable sources, you should use them. Simultaneously, we recommend you to have good arguments behind the usage of the informative priors. If you do not have trustful sources for the prior distributions, however, you should use the non-informative priors.

Point 9. Are the results stable from a sensitivity analysis?

• The idea of the sensitivity analysis is that it tells us how robust the estimates are depending on the specifications of the prior distribution.
• To investigate the stability of your results with the sensitivity analysis, you can adjust the hyperparameters of the priors upward and downward, thus re-estimating the model with varied priors to check for robustness (see Van Erp, Mulder, & Oberski, 2018 for more information on the sensitivity analysis).

“If informative or weakly-informative priors are used, then we suggest running a sensitivity analysis of these priors. When subjective priors are in place, then there might be a discrepancy between results using different subjective prior settings. A sensitivity analysis for priors would entail adjusting the entire prior distribution (i.e., using a completely different prior distribution than before) or adjusting hyperparameters upward and downward and re-estimating the model with these varied priors. Several different hyperparameter specifications can be made in a sensitivity analysis, and results obtained will point toward the impact of small fluctuations in hyperparameter values. […] The purpose of this sensitivity analysis is to assess how much of an impact the location of the mean hyperparameter for the prior has on the posterior. […] Upon receiving results from the sensitivity analysis, assess the impact that fluctuations in the hyperparameter values have on the substantive conclusions. Results may be stable across the sensitivity analysis, or they may be highly instable based on substantive conclusions. Whatever the finding, this information is important to report in the results and discussion sections of a paper. We should also reiterate here that original priors should not be modified, despite the results obtained.”

Point 10. Is the Bayesian way of interpreting and reporting model results used?

• Finally, it is time to interpret and report models that we constructed with the informative priors. You can find the JAGS model specification below (it is the same code that we used at the beginning of Stage 2). Copy and paste the code into the Enter JAGS model below section. Don’t forget to press Ctrl and Enter!

model{
    # likelihood
    for (i in 1:333){
    mu[i] <- intercept + beta1 * E22_Age[i] + beta2 * E22_Age_Squared[i]
    B3_difference_extra[i] ~ dnorm(mu[i],tau)
    }

    # priors
    intercept ~ dnorm(-35,1/20)
    beta1 ~ dnorm(0.8,1/5)
    beta2 ~ dnorm(0,1/10)
    tau ~ dgamma(0.5,0.5)
}

• In the output panel, the results are presented in the MCMC Summary table.

• Although we present how to interpret the parameter estimates, please refer to Van de Schoot and Depaoli (2014) for a summary of how to interpret and report models thoroughly.
  • The posterior mean of the regression coefficient of age (beta1) is 2.141. We interpret this value such that a one-year increase in age adds about 2.141 delays in the Ph.D. projects on average. The 95% credible interval of age is [1.731, 2.550], which means that there is a 95% probability that the regression coefficient of age lies in the population with the corresponding credible interval. Since the 95% credible interval does not include 0, there is the evidence of the effect of age in predicting the Ph.D. delay.
  • The posterior mean of the regression coefficient of age-squared (beta2) is -0.021. This can be interpreted such that one unit increase of the age-squared variable leads to a decrease of 0.021 unit in the Ph.D. delay on average. The 95% credible interval of age-squared is [-0.026, -0.015]. This indicates that we are 95% sure that the regression coefficient of the age-squared variable lies within the corresponding interval in the population. Therefore, it is fairly sure that there is an effect of the age-squared variable in predicting the delay.

References

Depaoli, S., & Van de Schoot, R. (2017). Improving transparency and replication in Bayesian statistics: The WAMBS-Checklist. Psychological methods, 22(2), 240-261. https://doi.org/10.1037/met0000065

JASP Team (2020). JASP (Version 0.13.1)[Computer software].

Merkle, E., & Rosseel, Y. (2018). blavaan: Bayesian Structural Equation Models via Parameter Expansion. Journal of Statistical Software, 85(4), 1–30. https://doi.org/10.18637/jss.v085.i04

Lynch, S. M. (2007). Introduction to applied Bayesian statistics and estimation for social scientists. Springer Science & Business Media.

Link, W. A., & Eaton, M. J. (2012). On thinning of chains in MCMC. Methods in ecology and evolution, 3(1), 112-115. https://doi.org/10.1111/j.2041-210X.2011.00131.x

Van de Schoot, R. (2020). PhD-delay Dataset for Online Stats Training [Data set]. Zenodo. https://doi.org/10.5281/zenodo.3999424

Van de Schoot, R., & Depaoli, S. (2014). Bayesian analyses: Where to start and what to report. The European Health Psychologist, 16(2), 75-84.

Van de Schoot, R., Yerkes, M. A., Mouw, J. M., & Sonneveld, H. (2013). What took them so long? Explaining PhD delays among doctoral candidates. PloS one, 8(7), e68839. https://doi.org/10.1371/journal.pone.0068839

Van Erp, S., Mulder, J., & Oberski, D. L. (2018). Prior sensitivity analysis in default Bayesian structural equation modeling. Psychological Methods, 23(2), 363-388. https://doi.org/10.1037/met0000162