Descriptive statistics

• Let’s see the descriptive statistics of the variables to check whether all the data points make sense.
• Click Descriptives at the top bar to proceed with -> Move all the variables to space under the Variables section -> Check Frequency tables (nominal and ordinal variables) -> Click Statistics -> Check additional descriptive statistics (S. E. mean under Dispersion, Median under Central Tendency, and Skewness and Kurtosis under Distribution) for better understanding.

• For all the variables, the sample size is 333 without any missing values.
• For the difference variable (B3_difference_extra), the mean and the standard error of the mean are 9.97 and 0.79.
• For the age variable (E22_Age), the mean and the standard error of the mean are 31.68 and 0.38.
• For the age-squared variable (E22_Age_Squared), the mean and the standard error of the mean are 1050.22 and 35.97.
• For the dichotomous variables (E4_having_child and E21_sex), frequency tables are presented.
• Given the descriptive statistics, all the data points substantively make sense.

Data visualization

• Let’s plot the expected relationship between age and the delay for visual inspection. To do that, we need to remain only two variables, B3_difference_extra and E22_Age, under the Variables section.

• Click Plots and check Scatter Plots -> Under Scatter Plots, uncheck Show confidence interval 95.0%. We can see the scatter plot between the difference variable on the x-axis and the age variable on the y-axis.
• Although we can use this scatter plot for our visual inspection purposes, it might be better to locate the difference variable on the y-axis and the age variable on the x-axis since we regress the difference variable on the age variable. We can simply do that by changing the variable order under the Variables section: set E22_Age to be the first and B3_difference_extra the latter.

• With this plot, we can see the non-linear relationship between the difference variable and the age variable.
• It is also possible to examine the linear relationship by drawing the linear regression line. How? Uncheck Smooth and check Linear under Add regression line. Can you see the linear line?

Regression coefficients

• We follow the default prior specifications that the R package blavaan uses with the JAGS target (Merkle, & Rosseel, 2018). You can find more information in the tutorial Bayesian regression in blavaan (using JAGS).
• In blavaan, the normal prior \(\mathcal{N}(\mu, \ 1/\sigma^{2})\), where \(\mu\) is the mean and \(\ 1/\sigma^{2}\) is the precision of the distribution, is used for the regression coefficient.
• “Wait, what is the precision? Should it be the variance?” A great question! JAGS uses the precision instead of the variance. The precision is the inverse of the variance. For example, the precision of 0.01 is equal to the variance of 100 (the inverse of 0.01 is 100).
• Blavaan uses 0 and 1e-2 (1e-2 is equal to \(1 × 10^{-2}\), which is 0.01) as the mean and the precision for the prior distribution for regression coefficients. Since the mean 0 and the precision 1e-2 are the parameters for the prior distribution, they are hyperparameters. In summary, we set \(\mathcal{N}(0, 1/100)\) as our prior distributions for both regression coefficients.

Intercept

• Although we leave intercept and residual variance untouched, we have to specify the prior distributions for the intercept and the residual variance to run JAGS (otherwise, JAGS does not work!).
• We again refer to the default priors that blavaan uses for the intercept: the normal prior distribution \(\mathcal{N}(\mu_0, \ 1/\sigma^{2}_{0})\), where \(\mu_0\) is the mean and \(\ 1/\sigma^{2}_{0}\) is the precision of the distribution.
• Blavaan uses 0 and 1e-3 (1e-3 is equal to \(1 × 10^{-3}\), which is 0.001) as hyperparameters for the mean and the precision, respectively. That is, we set \(\mathcal{N}(0, 1/1000)\) as our prior distribution for the intercept.

Residual variance

• In a normal regression model, the inverse gamma distribution is used as a prior for the variance thanks to the conjugacy (Lynch, 2007). Keep in mind, however, that JAGS uses precision instead of the variance. Therefore, we have to consider prior distributions for the precision.
• To do that, we should use the gamma distribution as a prior for the precision. In other words, using a gamma distribution on the precision is the same as using an inverse gamma distribution on the variance (think about the inverse relationship!). Note that blavaan also places prior distributions on the precision instead of the variance for the residual variance.
• The gamma distribution looks like \(Gamma(\kappa_0,\theta_0)\), where \(\kappa_0\) is the shape parameter and \(\theta_0\) the rate parameter of the distribution.
• As hyperparameters, let’s use 0.5 for both shape and rate parameter. This setting is an uninformative prior for the precision, which has been found to perform well. In other words, we set \(Gamma(.5, .5)\) as our prior distribution for the precision.

JAGS grammar: A primer

• In writing the JAGS code, we have to follow the JAGS grammar. Although we list a couple of them, readers can find the extensive JAGS manual at the mcmc-jags project. We also refer readers to RJAGS: how to get started that illustrates how to run the JAGS but in R.
• The two main pillars that consist of JAGS model specification are the likelihood and the prior.
  • For the likelihood, we specify the model of interest: the regression model in our case.
  • For the prior, we specify the type of distribution (e.g., normal distribution, gamma distribution) with its hyperparameters.
• To that end, the JAGS model that we use can be specified as:

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(0.5,0.5)
}

• Did you notice that there are two pillars in the JAGS model specification, which are the likelihood and the prior? They are guided with the number sign (#).
• In the code, intercept, beta1, beta2, and tau are the parameters.
  • The regression coefficients for the age and the age-squared variable are beta1 and beta2, respectively.
  • Tau is the precision parameter, which is the inverse of the variance.
• In the likelihood, 333 in ‘for (i in 1:333)’ is equal to the sample size. The alphabet i corresponds to each row in the dataset.
• The tilde (~) sign means ‘follow’.
  • beta1 ~ dnorm(0,1e-2) means that the regression coefficient of the age variable follows the normal distribution with the mean of 0 and the precision of 0.01.
• Such terms as dnorm and dgamma refer to the density of the normal distribution and the density of the gamma distribution.
• Copy and paste the code into the Enter JAGS model below section.

Set a seed value and the number of iterations

• Before running this code, we set the seed and the number of iterations.
• Regarding a seed, let’s use a seed value of 123 for replicable results. To do that, click the Advanced option in the control panel. Below Repeatability, type 123 for the seed.
• The initial number of burnin samples and that of samples to take are 500 and 2000, respectively. For our example, we use 5000 burnin samples and take additional 20000 samples. Under MCMC parameters of the Advanced option, type 20000 in No. samples (i.e., the number of samples) and 5000 in No. burnin samples (i.e., the number of burin samples). In total, we do 25000 iterations.
• We keep using the seed value of 123 and the number of iterations for subsequent analyses.
• If you have followed what is illustrated, your JASP screen should look like the snapshot below.

• Hope we are looking at the same screen! Go again to the code that we have pasted (i.e., the Enter JAGS model below section). Press Ctrl and Enter on the keyboard together to apply the code. With this action, you are running the model that we specified!
• Shortly, we can see a list of parameters under the Parameters in model in the control panel. This means that the posterior distributions for regression parameters are obtained by JAGS.

• As we are interested in the posterior distributions of two regression coefficients (i.e., the age and the age-squared variable), move beta1 and beta2 into the Show results for these parameters section on the right. After a few seconds, JASP presents the MCMC Summary table on the output panel.

Interpretation of the output

• We can find a summary of regression coefficients in the MCMC Summary table. The table presents the posterior distribution of regression coefficients after considering the default priors for age and age-squared variable and the likelihood of the data.
  • The posterior mean of the regression coefficient of age (beta1) is 2.356. We can interpret this value such that a one-year increase in age adds about 2.356 delays in Ph.D. projects on average. The 95% credible interval of age is [1.283, 3.420], which means that there is a 95% probability that the regression coefficient of age lies in the population with the corresponding credible interval. We also notice that a 95% credible interval does not include 0, which shows 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.023. This can be interpreted such that one unit increase of the age-squared variable leads to a decrease of 0.023 unit in the Ph.D. delay on average. The 95% credible interval of [-0.034, -0.012] indicates that we are 95% sure that the regression coefficient of age-squared lies within the corresponding interval in the population. Did you notice that 0 is not included in the interval? Therefore, it is fairly sure that there is an effect of the age-squared variable in predicting the delay.
• If you are curious about examining convergence diagnostics, please go to WAMBS Checklist in JASP (using JAGS).

beta1 ~ N(3, 1/0.4) and beta2 ~ N(0, 1/0.1)

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(3,1/0.4)
    beta2 ~ dnorm(0,1/0.1)
    tau ~ dgamma(0.5,0.5)
}

beta1 ~ N(3, 1/1000) and beta2 ~ N(0, 1/1000)

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(3,1/1000)
    beta2 ~ dnorm(0,1/1000)
    tau ~ dgamma(0.5,0.5)
}

beta1 ~ N(20, 1/0.4) and beta2 ~ N(20, 1/0.1)

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(20,1/0.4)
    beta2 ~ dnorm(20,1/0.1)
    tau ~ dgamma(0.5,0.5)
}

beta1 ~ N(20, 1/1000) and beta2 ~ N(20, 1/1000)

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(20,1/1000)
    beta2 ~ dnorm(20,1/1000)
    tau ~ dgamma(0.5,0.5)
}

Summary of the posterior mean and the posterior standard deviation

Relative bias

• Relative bias can be used to express the difference between the default prior and the informative prior.
• The formula for the relative bias follows: \(bias = 100*\frac{(model \; with \; specified \; priors\; – \; model \; with \; default \;priors)}{model \; with \; default \; priors}\)
• To preserve clarity, we will just calculate the bias of the two regression coefficients and only compare the model with the default prior and the model that uses \(\mathcal{N}(20, 1/0.4)\) for the age variable and \(\mathcal{N}(20, 1/0.1)\) for the age-squared variable as informative priors.
  • For the age variable, the relative bias is \(100*\frac{(11.215\; – \; 2.356)}{2.356}\). This equals 376.02%.
  • For the age-squared variable, the relative bias is \(100*\frac{(-0.114\; – \; (-0.023))}{-0.023}\), which equals 395.65%.
• The influence of the informative prior is around 376.02% for the age variable and 395.65% for the age-squared variable. This signals a tremendous change in the results with different prior specifications. For the concrete examination, let’s take a look at the parameter estimates.

Parameter estimates

• See the posterior mean, the posterior standard deviation, and the 95% credible interval from the MCMC Summary table.

• For both variables, there are changes in the parameter estimates. For instance, a one-year increase in age adds about 2.356 delays in Ph.D. projects on average under the default prior. When the informative prior is used, a one-year increase in age adds about 11.215 delays on average!