Bayesian Response Surface Analysis (BRSA)

Overview

In the preceding tutorial, Foundations of Response Surface Analysis for Congruence Hypothesis, we introduced the use of response surface analysis (RSA) in testing the congruence hypothesis and achieved specifying, fitting, visualizing, and interpreting a polynomial regression model with the RSA package for both the broad and strict congruence hypothesis.

In the current tutorial, we will first move away from the RSA package to replicate specifying, fitting, visualizing, and interpreting a polynomial regression model with more flexible functions/packages in R and introduce the Bayesian Response Surface Analysis (BRSA). More specifically, we provide

  1. the frequentist RSA using lm(), glm(), and sem(), and

  2. introduce the Bayesian RSA to move from null hypothesis testing to a model comparison framework and equivalence tests of the posterior distributions.

It should be noted that this tutorial is closely connected to the first tutorial on implementing response surface analysis and examining congruence hypothesis. We assume that readers are familiar with the concepts introduced in the first tutorial including the first principal axis (PA1), Line of Congruence (LOC), and Line of InCongruence (LOIC) as many of the concepts and procedures would not be repeated in this tutorial.

Why other packages?

The existing RSA package offers a convenient way of implementing various common response surface models in R. However, as one became more interested in diverse data structures (e.g., multilevel data), more than two predictors and/or more than one outcome in dyadic data, and varying non-linear dynamics other than second- (current example) and third-order (i.e., cubic RSA, Humberg et al., 2022) polynomial regressions, some other functions and packages in R could provide more flexibility for those extensions, including lm() and glm() from the stats package (R Core Team, 2021), and sem() from the lavaan package (Rosseel et al., 2012).

More importantly, we move away from the frequentist approach in congruence hypothesis testing and demonstrate the Bayesian approach in fitting, visualizing, and interpreting the congruence hypothesis mainly with the brms package (Bürkner, 2017).

We propose the Bayesian Response Surface Analysis (BRSA) as a better alternative to the frequentist RSA for its superiority in obtaining the “true” parameter distributions for congruence hypothesis testing.

In addition, as one will see in the following section, parameter transformation can be unexciting and prone to errors (especially when calculating the partial derivatives) when we try to obtain the point estimate and its standard deviation for hypothesis testing.

We thus also demonstrate here how the Bayesian parameter estimation and hypothesis testing, instead of the frequentist approach, can make this process less painful. For more on the Bayesian parameter estimation, one could also consult the series of tutorials, Bayesian Methods, by Fischer (2024).

We will now fit the same model with the same data

  1. in all three functions (i.e., lm(), glm(), sem()) and estimate the parameters needed for interpreting the congruence hypothesis in the frequentist approach; and

  2. in Bayesian RSA using brm and demonstrate the congruence hypothesis testing in the Bayesian framework.

Bayesian Response Surface Analysis (BRSA)

As we presented in the interpreting RSA with the RSA package in the previous tutorial, the congruence hypothesis testing, as shown in the table below, involves a test that the first principle axis should not significantly differ from the Line of Congruence (LOC), using \(p_{10}=0\) and \(p_{11}=1\) as the first two tests for congruence effects.

Table 1

Results from Response Surface Analysis Examining How the Congruence of Perceived Self Agency and Perceived Other Agency Influences Happiness

Following the frequentist hypothesis testing, the RSA model estimation would suggest that \(p_{10}\) is not significant different than 0 and that \(p_{11}\) is not significant different than 1. However, this might not always be interpreted as “the first principle axis did not significantly differ from the LOC.” The current result would suggest that the standard error of the parameter \(p_{10}\) is over ten times of its point estimate. The wide confidence interval that includes 0 might only indicate the instability and uncertainty of the estimation of \(p_{10}\), instead of a more stable and certain estimation of \(p_{10}\approx0\).

To make up for such concerns, we here propose the Bayesian Response Surface Analysis where we could use model comparison and posterior distributions for congruence hypothesis testing. Unlike the frequentist approach that estimates how well sample-based estimates and statistical tests would perform under repeated sampling from the same population (with confidence intervals based on standard errors), Bayesian parameter estimation returns a probability distribution of the parameter itself (Kruschke, 2010).

Thus, Bayesian posterior distributions facilitate the direct hypothesis testing of how much of the distribution of a parameter of interest (e.g., \(p_{10}\)) overlaps with the region practical equivalence (ROPE; e.g., \([-0.05, 0.05]\) as practically equivalent to \(p_{10}\approx0\)) to accept or reject the congruence hypothesis. More specifically, we will calculate the the posterior distributions for the parameters of interest, as detailed in the table above, and quantify the proportions of the parameter within the region practical equivalence (ROPE) to make claims about the surface and examine the congruence hypothesis.

In addition, Bayesian model comparison can provide a direct comparison of how different models perform in predicting the outcome using leave-one-out cross validation. More specifically, we will compare a freely estimated polynomial regression model (i.e., the model estimated in the preceding tutorial) and another two polynomial regression models with constraints imposing the broad and the strict congruence hypothesis tests (i.e., \(p_{10} = 0\), \(p_{11} = 1\), and \(a_{3} = 0\) for broad congruence; \(a_{1} = 0\) and \(a_{2} = 0\) for strict congruence). If imposing such constraints would not render the model of significantly worse fit (or even render the model of better fit) than a freely estimated model, we could then accept the congruence hypothesis.

RSA with lm(), glm(), and sem()

Load Necessary Libraries

library(knitr)            # knitr tables
library(kableExtra)       # knitr table formatting
library(psych)            # data descriptives
library(RSA)              # Response Surface Analysis
library(plot3D)           # visualize response surfaces
library(plotly)           # visualize response surfaces
library(lavaan)           # structural equation modeling
library(brms)             # Bayesian parameter estimation
library(posterior)        # Bayesian posterior distribution
library(patchwork)        # add plots together
library(sandwich)         # calculate robust SEs
library(lmtest)           # linear regression tests
library(bayestestR)       # region of practical equivalence (ROPE)
library(see)              # visualize ROPE
library(parameters)       # parameter standardization 
library(tidyverse)        # data structures & visualization

Step 1: Specify Congruence Hypothesis

We will continue to use the iSAHIB dataset and explore the effect of the fit between the perceived agency of the self (SelfAgency_X, predictor 1) and the perceived agency of the interaction partner (OtherAgency_Y, predictor 2) in dyadic interpersonal interactions on happiness (Happy_Z, the outcome).

We again test the congruent hypothesis - the more congruent between the perceived self agency and the perceived other agency, the happier oneself would be.

To request the iSAHIB data, please have a look at the study materials and documentary - and send a note to Nilam Ram.

#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/refs/heads/main/iSAHIB_RSA/isahib_RSA.csv"

#read in the .csv file using the url() function
isahib <- read.csv(file=url(filepath), header=TRUE)

#look at the data
head(isahib, 10)
##    SelfAgency_X OtherAgency_Y Happy_Z
## 1            82            22      57
## 2            97            21      66
## 3            55            51      67
## 4            82            13      73
## 5            32            61      94
## 6            14            95      83
## 7            56            52      78
## 8            34            71      74
## 9            50            50      67
## 10           51            48      77
#describe the data
describe(isahib)
##               vars   n  mean    sd median trimmed   mad min max range  skew
## SelfAgency_X     1 866 51.51 15.76   53.0   51.85 14.83   2  97    95 -0.22
## OtherAgency_Y    2 866 55.93 17.07   57.5   56.45 17.05   4  99    95 -0.28
## Happy_Z          3 866 81.41 11.33   84.0   83.44  5.93   5  98    93 -2.91
##               kurtosis   se
## SelfAgency_X     -0.19 0.54
## OtherAgency_Y    -0.14 0.58
## Happy_Z          11.77 0.39
#plot the distribution of each variable and the correlation between them
psych::pairs.panels(isahib)

#plot raw data scatter plot in three dimensions
plot_ly(x = isahib$SelfAgency_X, y = isahib$OtherAgency_Y, z = isahib$Happy_Z, 
        type = "scatter3d", mode = "markers", color = isahib$Happy_Z) %>%
   layout(scene = list(xaxis=list(title = "Self Agency", nticks=10, range=c(0,100)),
                       yaxis=list(title = "Other Agency", nticks=10, range=c(0,100)),
                       zaxis=list(title = "Happiness", nticks=10, range=c(0,100))))
## Warning: Ignoring 3 observations

Step 2: Data Centering

We again center predictors on the scale midpoint (Barranti et al., 2017). In our case, this means centering at 50, the midpoint on a 0 to 100 scale.

#center predictors with the scale midpoint
isahib <- isahib %>%
  mutate(cSelfAgency_X = SelfAgency_X - 50,
         cOtherAgency_Y = OtherAgency_Y - 50)

Step 3: Fit the Polynomial Regression Model and Plot the Response Surface

glm() and lm()

#glm() and lm() function in the stats package
lm_happy_agency <- lm(formula = Happy_Z ~ cSelfAgency_X * cOtherAgency_Y + 
                        I(cSelfAgency_X^2) + I(cOtherAgency_Y^2),
                      data = isahib)

glm_happy_agency <- glm(formula = Happy_Z ~ cSelfAgency_X * cOtherAgency_Y + 
                          I(cSelfAgency_X^2) + I(cOtherAgency_Y^2),
                        data = isahib,
                        family = gaussian(link = "identity")) #conditional distribution and link function

# model output
summary(lm_happy_agency)
## 
## Call:
## lm(formula = Happy_Z ~ cSelfAgency_X * cOtherAgency_Y + I(cSelfAgency_X^2) + 
##     I(cOtherAgency_Y^2), data = isahib)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -72.850  -2.269   1.963   6.076  25.272 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  84.3008571  0.5348587 157.613  < 2e-16 ***
## cSelfAgency_X                -0.0787281  0.0242744  -3.243 0.001227 ** 
## cOtherAgency_Y               -0.0628345  0.0231707  -2.712 0.006825 ** 
## I(cSelfAgency_X^2)           -0.0042935  0.0012754  -3.366 0.000795 ***
## I(cOtherAgency_Y^2)          -0.0040760  0.0010747  -3.793 0.000159 ***
## cSelfAgency_X:cOtherAgency_Y  0.0002195  0.0011029   0.199 0.842292    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.91 on 860 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.07896,    Adjusted R-squared:  0.0736 
## F-statistic: 14.74 on 5 and 860 DF,  p-value: 6.796e-14
summary(glm_happy_agency)
## 
## Call:
## glm(formula = Happy_Z ~ cSelfAgency_X * cOtherAgency_Y + I(cSelfAgency_X^2) + 
##     I(cOtherAgency_Y^2), family = gaussian(link = "identity"), 
##     data = isahib)
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  84.3008571  0.5348587 157.613  < 2e-16 ***
## cSelfAgency_X                -0.0787281  0.0242744  -3.243 0.001227 ** 
## cOtherAgency_Y               -0.0628345  0.0231707  -2.712 0.006825 ** 
## I(cSelfAgency_X^2)           -0.0042935  0.0012754  -3.366 0.000795 ***
## I(cOtherAgency_Y^2)          -0.0040760  0.0010747  -3.793 0.000159 ***
## cSelfAgency_X:cOtherAgency_Y  0.0002195  0.0011029   0.199 0.842292    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 118.9638)
## 
##     Null deviance: 111079  on 865  degrees of freedom
## Residual deviance: 102309  on 860  degrees of freedom
##   (3 observations deleted due to missingness)
## AIC: 6604
## 
## Number of Fisher Scoring iterations: 2

sem()

#sem() function in the lavaan package
model_sem <- '
Happy_Z ~ b1 * cSelfAgency_X + b2 * cOtherAgency_Y + b3 * cSelfAgency_X:cSelfAgency_X + b4 * cSelfAgency_X:cOtherAgency_Y + b5 * cOtherAgency_Y:cOtherAgency_Y
          
Happy_Z ~ 1
'

sem_happy_agency <- sem(model_sem,
                        estimator = "ML",
                        se = "standard",
                        data = isahib)

summary(sem_happy_agency)
## lavaan 0.6-19 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         7
## 
##                                                   Used       Total
##   Number of observations                           866         869
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate    Std.Err  z-value  P(>|z|)
##   Happy_Z ~                                             
##     cSlfAgn_X (b1)     -0.079    0.024   -3.255    0.001
##     cOthrAg_Y (b2)     -0.063    0.023   -2.721    0.007
##     cSA_X:SA_ (b3)     -0.004    0.001   -3.378    0.001
##     cSA_X:OA_ (b4)      0.000    0.001    0.200    0.842
##     cOA_Y:OA_ (b5)     -0.004    0.001   -3.806    0.000
## 
## Intercepts:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Happy_Z            84.301    0.533  158.162    0.000
## 
## Variances:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Happy_Z           118.140    5.677   20.809    0.000

Compare model output from glm()/lm() with sem()

The model outputs for the regression coefficients’ estimations are the “same” across all four functions with slightly different corresponding standard errors (SEs) and two different forms of significance testing. The differences come from the two different estimation processes used. Functions lm() and glm() use ordinary least squares (OLS) while functions sem(), which is the background function for RSA(), uses maximum likelihood estimation (MLE).

In the current case of a polynomial regression with one outcome and no latent variables, both methods would produce the same outcome for regression coefficients. However, those two methods differ in their SE estimation. The SE estimation in sem() under the normal assumption would always be smaller than its equivalent OLS estimation results from lm() and glm(). One could also specify likelihood = "wishart" instead of se="standard" to obtain estimation from the unbiased covariance matrix for even smaller discrepancies. Those discrepancies would also be smaller with larger sample size (for a more detailed discussion, see: Satorra & Bentler, 1994).

Those differences in SE estimation pertain to all the hypothesis testing but often not influential enough to alter the interpretation of the significance level with a 5% threshold. The two different tests parallel to OLS and MLE distinction are t-test and z-test, each correspondent of their own statistical assumptions about estimating the sample parameter or the population parameter.

After fitting the polynomial regressions in all four packages/functions, we turn to visualize the 3D surface. Since the brm() is used for parameter estimation, we demonstrate here how to plot the surface with the lm() and glm() object

lm() plot

#create the data grid
x_grid <- seq(-50, 50, by = 1) #use the range of the centered predictor
y_grid <- seq(-50, 50, by = 1)
newdat <- expand.grid(x_grid, y_grid)
colnames(newdat) <- c("cSelfAgency_X", "cOtherAgency_Y")

#predict the z value from the lm() model
pred <- predict(lm_happy_agency, 
                newdata = newdat, 
                se.fit = TRUE)

#calculate the confidence intervals for predictions
alpha <- 0.05 # 95% confidence interval
ymin <- pred$fit - qnorm(1 - alpha/2) * pred$se.fit
ymax <- pred$fit + qnorm(1 - alpha/2) * pred$se.fit
fit  <- pred$fit # the predicted values

#output for visualization
z_matrix <- matrix(fit, nrow = length(x_grid), byrow = TRUE)
ci_low_matrix <- matrix(ymin, nrow = length(x_grid), byrow = TRUE)
ci_up_matrix  <- matrix(ymax, nrow = length(x_grid), byrow = TRUE)

#plot
plot_ly(x = x_grid, y = y_grid) %>% 
  add_surface(z = z_matrix,
              colorscale = list(c(0,1),c("red","blue"))) %>% 
  add_surface(z = ci_low_matrix, opacity = 0.5, showscale = FALSE, 
              colorscale = list(c(0,1),c("grey","grey"))) %>% 
  add_surface(z = ci_up_matrix, opacity = 0.5, showscale = FALSE, 
              colorscale = list(c(0,1),c("grey","grey")))

glm() plot

#create the data grid
x_grid <- seq(-52, 52, by = 1)
y_grid <- seq(-52, 52, by = 1)
newdat <- expand.grid(x_grid, y_grid)
colnames(newdat) <- c("cSelfAgency_X", "cOtherAgency_Y")

#predict the z value from the glm() model
pred <- predict(glm_happy_agency, 
                newdata = newdat, 
                type = "link", # standard errors for confidence interval
                se.fit = TRUE)

#calculate the confidence intervals for predictions
alpha <- 0.05 # 95% confidence interval
ymin <- glm_happy_agency$family$linkinv(pred$fit - qnorm(1 - alpha/2) * pred$se.fit)
ymax <- glm_happy_agency$family$linkinv(pred$fit + qnorm(1 - alpha/2) * pred$se.fit)
fit  <- glm_happy_agency$family$linkinv(pred$fit) 

#output for visualization
z_matrix <- matrix(fit, length(x_grid), byrow = TRUE)
ci_low_matrix <- matrix(ymin, nrow = length(x_grid), byrow = TRUE)
ci_up_matrix  <- matrix(ymax, nrow = length(x_grid), byrow = TRUE)

plot_ly(x = x_grid, y = y_grid) %>% 
  add_surface(z = z_matrix,
              colorscale = list(c(0,1),c("red","blue"))) %>% 
  add_surface(z = ci_low_matrix, opacity = 0.5, showscale = FALSE, 
              colorscale = list(c(0,1),c("grey","grey"))) %>% 
  add_surface(z = ci_up_matrix, opacity = 0.5, showscale = FALSE, 
              colorscale = list(c(0,1),c("grey","grey")))

sem() plot

#create the data grid
x_grid <- seq(-52, 52, by = 1)
y_grid <- seq(-52, 52, by = 1)
newdat <- expand.grid(x_grid, y_grid)
colnames(newdat) <- c("cSelfAgency_X", "cOtherAgency_Y")

#for lavPredictY(), include the outcome variable column although their values will be ignored). 
newdat$Happy_Z <- NA 

#predict the z value from the sem() model
newdat$Happy_Z <- lavPredictY(sem_happy_agency, 
                              newdata = newdat)

#output for visualization
z_matrix <- matrix(newdat$Happy_Z, nrow = length(x_grid), byrow = TRUE)
plot_ly(x = x_grid, y = y_grid) %>% 
  add_surface(z_matrix,
              colorscale = list(c(0,1),c("red","blue")))

Surface constructed using lm() is equivalent or glm(), both are the interactive version of the surface constructed using the rsa() object with some differences in the functions for calculating the confidence intervals as glm() requires specification of the link function.

Additionally, we added the estimated 95% confidence interval to demonstrate the possibilities of playing around with visualizing response surfaces. One difference is that there are currently no built-in functions in lavaan to obtain the standard errors (SEs) of predicted values. They can be obtained with more computation (e.g., Bootstrap). In addition, brm() here are used primarily for parameter estimation. We thus did not include the surface visualization using model outputs from brm() or the procedure of calculating the 95% confidence interval from the sem() model output for brevity.

Step 4: Parameter Transformation, Model Constraints, & Posterior Distribution

In this step, we provide more flexibility in examining the statistical significance for congruence hypothesis testing. Parameter transformation aims at calculating the point estimate and confidence interval and is applicable to lm(), glm(), and sem(). We here demonstrate additional two: model constrains for sem() and posterior distributions for brm().

Parameter Transformation

Three parameters (i.e., \(a_{3}\), \(a_{4}\), \(a_{5}\)) are needed for the broad sense of congruence hypothesis and additional two parameters (i.i., \(a_{1}\), \(a_{2}\)) can also be included for the strict sense of congruence hypothesis (i.e., without the main effects of the predictors). Those parameters are transformed using the model parameters from the line of congruence (LOC), first principal axis, and the line of incongruence (LOIC). One could refer to the preceding tutorial on basic RSA for the transformation process and motivation behind it.

We here present how we can obtain the point estimation along with confidence intervals using lm() and sem(). We did not repeat the procedure for glm() because (1) glm() is equivalent to lm() in the present case and (2) it is relatively straightforward to extract the model coefficients and the variance-covariance matrices from lm() output to the glm() output while the other steps of calculation and estimation remain the same.

\[ \begin{split} a_{1} & = b_{1} + b_{2} \\ \\ a_{2} & = b_{3}+b_{4}+b_{5} \\ \\ a_{3} & = b_{1} - b_{2} \\ \\ a_{4} & = b_{3}-b_{4}+b_{5} \\ \\ a_{5} & = b_{3}-b_{5} \\ \end{split} \]

Parameter Transformation with lm()/glm()

a1

# using lm() output for a1

#extract the regression coefficient for computing the estimate of a1
lm_coefs <- coef(lm_happy_agency)
b1 <- lm_coefs["cSelfAgency_X"]
b2 <- lm_coefs["cOtherAgency_Y"]

a1 <- b1 + b2
a1
## cSelfAgency_X 
##    -0.1415627
#calculate the partial derivatives of a1 with respect to each coefficient
grad_a1 <- c(1, 1)

#using standard covariance matrix for delta method
standard_cov <- vcov(lm_happy_agency)
relevant_stand_cov <- standard_cov[c("cSelfAgency_X", 
                                     "cOtherAgency_Y"), 
                                   c("cSelfAgency_X", 
                                     "cOtherAgency_Y")]

# Variance of a1
var_stand_a1 <- t(grad_a1) %*% relevant_stand_cov %*% grad_a1

# Standard error of a1
se_stand_a1 <- sqrt(var_stand_a1)
se_stand_a1
##            [,1]
## [1,] 0.03207802
# confidence intervals
ci_lower_a1 <- unname(a1) - qnorm(0.975) * unname(se_stand_a1)
ci_lower_a1
##            [,1]
## [1,] -0.2044344
ci_upper_a1 <- unname(a1) + qnorm(0.975) * unname(se_stand_a1)
ci_upper_a1
##             [,1]
## [1,] -0.07869091
# t-test 
t_a1 <- unname(a1) / unname(se_stand_a1)
p_value_a1 <- 2 * (1 - pt(abs(t_a1), df = df.residual(lm_happy_agency)))
p_value_a1
##              [,1]
## [1,] 1.148623e-05

a2

# using lm() output for a2

#extract the regression coefficient for computing the estimate of a2
lm_coefs <- coef(lm_happy_agency)
b3 <- lm_coefs["I(cSelfAgency_X^2)"]
b4 <- lm_coefs["cSelfAgency_X:cOtherAgency_Y"]
b5 <- lm_coefs["I(cOtherAgency_Y^2)"]

a2 <- b3 + b4 + b5
a2
## I(cSelfAgency_X^2) 
##        -0.00814994
#calculate the partial derivatives of a2 with respect to each coefficient
grad_a2 <- c(1, 1, 1)

#using standard covariance matrix for delta method
standard_cov <- vcov(lm_happy_agency)
relevant_stand_cov <- standard_cov[c("I(cSelfAgency_X^2)", 
                                     "cSelfAgency_X:cOtherAgency_Y", 
                                     "I(cOtherAgency_Y^2)"), 
                                   c("I(cSelfAgency_X^2)", 
                                     "cSelfAgency_X:cOtherAgency_Y",
                                     "I(cOtherAgency_Y^2)")]

# Variance of a2
var_stand_a2 <- t(grad_a2) %*% relevant_stand_cov %*% grad_a2

# Standard error of a2
se_stand_a2 <- sqrt(var_stand_a2)
se_stand_a2
##             [,1]
## [1,] 0.001992289
# confidence intervals
ci_lower_a2 <- unname(a2) - qnorm(0.975) * unname(se_stand_a2)
ci_lower_a2
##             [,1]
## [1,] -0.01205475
ci_upper_a2 <- unname(a2) + qnorm(0.975) * unname(se_stand_a2)
ci_upper_a2
##              [,1]
## [1,] -0.004245126
# t-test 
t_a2 <- unname(a2) / unname(se_stand_a2)
p_value_a2 <- 2 * (1 - pt(abs(t_a2), df = df.residual(lm_happy_agency)))
p_value_a2
##              [,1]
## [1,] 4.703967e-05

a3

# using lm() output for a3

#extract the regression coefficient for computing the estimate of a3
lm_coefs <- coef(lm_happy_agency)
b1 <- lm_coefs["cSelfAgency_X"]
b2 <- lm_coefs["cOtherAgency_Y"]

a3 <- b1 - b2
a3
## cSelfAgency_X 
##   -0.01589362
#calculate the partial derivatives of a3 with respect to each coefficient
grad_a3 <- c(1, -1)

#using standard covariance matrix for delta method
standard_cov <- vcov(lm_happy_agency)
relevant_stand_cov <- standard_cov[c("cSelfAgency_X", 
                                     "cOtherAgency_Y"), 
                                   c("cSelfAgency_X", 
                                     "cOtherAgency_Y")]

# Variance of a3
var_stand_a3 <- t(grad_a3) %*% relevant_stand_cov %*% grad_a3

# Standard error of a3
se_stand_a3 <- sqrt(var_stand_a3)
se_stand_a3
##           [,1]
## [1,] 0.0349751
# confidence intervals
ci_lower_a3 <- unname(a3) - qnorm(0.975) * unname(se_stand_a3)
ci_lower_a3
##             [,1]
## [1,] -0.08444356
ci_upper_a3 <- unname(a3) + qnorm(0.975) * unname(se_stand_a3)
ci_upper_a3
##            [,1]
## [1,] 0.05265632
# t-test 
t_a3 <- unname(a3) / unname(se_stand_a3)
p_value_a3 <- 2 * (1 - pt(abs(t_a3), df = df.residual(lm_happy_agency)))
p_value_a3
##           [,1]
## [1,] 0.6496363

a4

# using lm() output for a4

#extract the regression coefficient for computing the estimate of a4
lm_coefs <- coef(lm_happy_agency)
b3 <- lm_coefs["I(cSelfAgency_X^2)"]
b4 <- lm_coefs["cSelfAgency_X:cOtherAgency_Y"]
b5 <- lm_coefs["I(cOtherAgency_Y^2)"]

a4 <- b3 - b4 + b5
a4
## I(cSelfAgency_X^2) 
##       -0.008588948
#calculate the partial derivatives of a4 with respect to each coefficient
grad_a4 <- c(1, -1, 1)

#using standard covariance matrix for delta method
standard_cov <- vcov(lm_happy_agency)
relevant_stand_cov <- standard_cov[c("I(cSelfAgency_X^2)", 
                                     "cSelfAgency_X:cOtherAgency_Y", 
                                     "I(cOtherAgency_Y^2)"), 
                                   c("I(cSelfAgency_X^2)", 
                                     "cSelfAgency_X:cOtherAgency_Y",
                                     "I(cOtherAgency_Y^2)")]

# Variance of a4
var_stand_a4 <- t(grad_a4) %*% relevant_stand_cov %*% grad_a4

# Standard error of a4
se_stand_a4 <- sqrt(var_stand_a4)
se_stand_a4
##             [,1]
## [1,] 0.001398085
# confidence intervals
ci_lower_a4 <- unname(a4) - qnorm(0.975) * unname(se_stand_a4)
ci_lower_a4
##             [,1]
## [1,] -0.01132914
ci_upper_a4 <- unname(a4) + qnorm(0.975) * unname(se_stand_a4)
ci_upper_a4
##              [,1]
## [1,] -0.005848751
# t-test 
t_a4 <- unname(a4) / unname(se_stand_a4)
p_value_a4 <- 2 * (1 - pt(abs(t_a4), df = df.residual(lm_happy_agency)))
p_value_a4
##              [,1]
## [1,] 1.234016e-09

a5

#extract the regression coefficient for computing the estimate of a5
lm_coefs <- coef(lm_happy_agency)
b3 <- lm_coefs["I(cSelfAgency_X^2)"]
b5 <- lm_coefs["I(cOtherAgency_Y^2)"]

a5 <- b3 - b5
a5
## I(cSelfAgency_X^2) 
##      -0.0002174958
#calculate the partial derivatives of a4 with respect to each coefficient
grad_a5 <- c(1, -1)

#using standard covariance matrix for delta method
standard_cov <- vcov(lm_happy_agency)
relevant_stand_cov <- standard_cov[c("I(cSelfAgency_X^2)", 
                                     "I(cOtherAgency_Y^2)"), 
                                   c("I(cSelfAgency_X^2)", 
                                     "I(cOtherAgency_Y^2)")]

# Variance of a5
var_stand_a5 <- t(grad_a5) %*% relevant_stand_cov %*% grad_a5

# Standard error of a4
se_stand_a5 <- sqrt(var_stand_a5)
se_stand_a5
##            [,1]
## [1,] 0.00195383
# confidence intervals
ci_lower_a5 <- unname(a5) - qnorm(0.975) * unname(se_stand_a5)
ci_lower_a5
##              [,1]
## [1,] -0.004046933
ci_upper_a5 <- unname(a5) + qnorm(0.975) * unname(se_stand_a5)
ci_upper_a5
##             [,1]
## [1,] 0.003611941
# t-test 
t_a5 <- unname(a5) / unname(se_stand_a5)
p_value_a5 <- 2 * (1 - pt(abs(t_a5), df = df.residual(lm_happy_agency)))
p_value_a5
##           [,1]
## [1,] 0.9113904

Summary Table

#create the table
table_lm <- data.frame(
  Parameter = c("a1", "a2", "a3", "a4", "a5"),
  Point_Estimation = c(sprintf("%.3f", a1), sprintf("%.3f", a2), sprintf("%.3f", a3), 
                       sprintf("%.3f", a4), sprintf("%.3f", a5)),  # Keeping no column name
  SE = c(sprintf("%.3f", se_stand_a1), sprintf("%.3f", se_stand_a2), sprintf("%.3f", se_stand_a3), 
         sprintf("%.3f", se_stand_a4), sprintf("%.3f", se_stand_a5)),
  CI = c(paste0("[", sprintf("%.3f", ci_lower_a1), ", ", sprintf("%.3f", ci_upper_a1), "]"),
         paste0("[", sprintf("%.3f", ci_lower_a2), ", ", sprintf("%.3f", ci_upper_a2), "]"),
         paste0("[", sprintf("%.3f", ci_lower_a3), ", ", sprintf("%.3f", ci_upper_a3), "]"),
         paste0("[", sprintf("%.3f", ci_lower_a4), ", ", sprintf("%.3f", ci_upper_a4), "]"),
         paste0("[", sprintf("%.3f", ci_lower_a5), ", ", sprintf("%.3f", ci_upper_a5), "]")),
  p_value = c(sprintf("%.3f", p_value_a1), sprintf("%.3f", p_value_a2), 
              sprintf("%.3f", p_value_a3), sprintf("%.3f", p_value_a4), 
              sprintf("%.3f", p_value_a5))
)

#change column names
colnames(table_lm)[2] <- "Point Estimate"
colnames(table_lm)[3] <- "$SE$"
colnames(table_lm)[5] <- "$p$ value"

#print table using knitr for better formatting
table_lm %>% 
  kable(caption = "Summary Table for Parameters Transformed from lm()") %>%
  column_spec(1, width = "2.5cm") %>%   # Adjust width for the "Parameter" column
  column_spec(2, width = "4cm") %>%     # Adjust width for "Point Estimate"
  column_spec(3, width = "3cm") %>%     # Adjust width for "SE"
  column_spec(4, width = "4cm") %>%     # Adjust width for "CI"
  column_spec(5, width = "3cm") %>%     # Adjust width for "*p* value"
  footnote(general = "$Note$: SE = Standard Error, CI = Confidence Interval.",
           general_title = "") %>%
  kable_styling(full_width = FALSE, html_font = "Times New Roman")
Summary Table for Parameters Transformed from lm()
Parameter Point Estimate \(SE\) CI \(p\) value
a1 -0.142 0.032 [-0.204, -0.079] 0.000
a2 -0.008 0.002 [-0.012, -0.004] 0.000
a3 -0.016 0.035 [-0.084, 0.053] 0.650
a4 -0.009 0.001 [-0.011, -0.006] 0.000
a5 -0.000 0.002 [-0.004, 0.004] 0.911
\(Note\): SE = Standard Error, CI = Confidence Interval.

Parameter Transformation with sem()

While we need to manually calculate the point estimate and confidence intervals for the lm() and glm() model outputs, sem() offers the flexibility to specify any parameters within the model and obtain the corresponding estimation directly.

#sem() function in the lavaan package
model_sem <- '
Happy_Z ~ b1 * cSelfAgency_X + b2 * cOtherAgency_Y + b3 * cSelfAgency_X:cSelfAgency_X + b4 * cSelfAgency_X:cOtherAgency_Y + b5 * cOtherAgency_Y:cOtherAgency_Y
          
Happy_Z ~ 1

a1 := b1 + b2
a2 := b3 + b4 + b5
a3 := b1 - b2
a4 := b3 - b4 + b5
a5 := b3 - b5

'

sem_happy_agency <- sem(model_sem,
                        estimator = "ML",
                        se = "standard",
                        data = isahib)

summary(sem_happy_agency)
## lavaan 0.6-19 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         7
## 
##                                                   Used       Total
##   Number of observations                           866         869
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate    Std.Err  z-value  P(>|z|)
##   Happy_Z ~                                             
##     cSlfAgn_X (b1)     -0.079    0.024   -3.255    0.001
##     cOthrAg_Y (b2)     -0.063    0.023   -2.721    0.007
##     cSA_X:SA_ (b3)     -0.004    0.001   -3.378    0.001
##     cSA_X:OA_ (b4)      0.000    0.001    0.200    0.842
##     cOA_Y:OA_ (b5)     -0.004    0.001   -3.806    0.000
## 
## Intercepts:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Happy_Z            84.301    0.533  158.162    0.000
## 
## Variances:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Happy_Z           118.140    5.677   20.809    0.000
## 
## Defined Parameters:
##                    Estimate    Std.Err  z-value  P(>|z|)
##     a1                 -0.142    0.032   -4.428    0.000
##     a2                 -0.008    0.002   -4.105    0.000
##     a3                 -0.016    0.035   -0.456    0.648
##     a4                 -0.009    0.001   -6.165    0.000
##     a5                 -0.000    0.002   -0.112    0.911

Constrained Model Specification for sem()

Instead of the a freely estimated model, we also estimate two models corresponding to the broad sense and the strict sense of congruence for model comparisons in testing congruence.

Broad sense

model_sem_broad <- '
# a3 = b1 - b2 = 0
# b1 = b2

# a5 = b3 - b5 = 0
# b3 = b5

Happy_Z ~ b1 * cSelfAgency_X + b1 * cOtherAgency_Y + b3 * cSelfAgency_X:cSelfAgency_X + b4 * cSelfAgency_X:cOtherAgency_Y + b3 * cOtherAgency_Y:cOtherAgency_Y

Happy_Z ~ 1

a4 := b3 - b4 + b3
'

sem_happy_agency_broad <- sem(model_sem_broad,
                              estimator = "ML",
                              se = "standard",
                              data = isahib)

summary(sem_happy_agency_broad)
## lavaan 0.6-19 ended normally after 13 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         7
##   Number of equality constraints                     2
## 
##                                                   Used       Total
##   Number of observations                           866         869
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.265
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.876
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate    Std.Err  z-value  P(>|z|)
##   Happy_Z ~                                             
##     cSlfAgn_X (b1)     -0.070    0.016   -4.431    0.000
##     cOthrAg_Y (b1)     -0.070    0.016   -4.431    0.000
##     cSA_X:SA_ (b3)     -0.004    0.001   -6.430    0.000
##     cSA_X:OA_ (b4)      0.000    0.001    0.138    0.890
##     cOA_Y:OA_ (b3)     -0.004    0.001   -6.430    0.000
## 
## Intercepts:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Happy_Z            84.309    0.531  158.633    0.000
## 
## Variances:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Happy_Z           118.176    5.679   20.809    0.000
## 
## Defined Parameters:
##                    Estimate    Std.Err  z-value  P(>|z|)
##     a4                 -0.008    0.001   -6.243    0.000

Strict sense

model_sem_strict <- '
# a1 = b1 + b2 = 0 
# a3 = b1 - b2 = 0
# b1 = 0
# b2 = 0

Happy_Z ~ 0 * cSelfAgency_X + 0 * cOtherAgency_Y + b3 * cSelfAgency_X:cSelfAgency_X + b4 * cSelfAgency_X:cOtherAgency_Y + b3 * cOtherAgency_Y:cOtherAgency_Y

Happy_Z ~ 1

# a2 = b3 + b4 + b5 = 0
# a5 = b3 - b5 = 0
# b3 = b5
b4 == -2 * b3

a4 := b3 - b4 + b3
'

sem_happy_agency_strict <- sem(model_sem_strict,
                              estimator = "ML",
                              se = "standard",
                              data = isahib)

summary(sem_happy_agency_strict)
## lavaan 0.6-19 ended normally after 12 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         5
##   Number of equality constraints                     2
## 
##                                                   Used       Total
##   Number of observations                           866         869
## 
## Model Test User Model:
##                                                       
##   Test statistic                                43.717
##   Degrees of freedom                                 4
##   P-value (Chi-square)                           0.000
## 
## Parameter Estimates:
## 
##   Standard errors                             Standard
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate    Std.Err  z-value  P(>|z|)
##   Happy_Z ~                                             
##     cSlfAgn_X           0.000                           
##     cOthrAg_Y           0.000                           
##     cSA_X:SA_ (b3)     -0.002    0.000   -5.287    0.000
##     cSA_X:OA_ (b4)      0.004    0.001    5.287    0.000
##     cOA_Y:OA_ (b3)     -0.002    0.000   -5.287    0.000
## 
## Intercepts:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Happy_Z            82.363    0.420  196.333    0.000
## 
## Variances:
##                    Estimate    Std.Err  z-value  P(>|z|)
##    .Happy_Z           124.257    5.971   20.809    0.000
## 
## Defined Parameters:
##                    Estimate    Std.Err  z-value  P(>|z|)
##     a4                 -0.007    0.001   -5.287    0.000
## 
## Constraints:
##                                                |Slack|
##     b4 - (-2*b3)                                 0.000

Step 5: Congruence Hypothesis Testing and Interpretation

We here present three different approaches in testing the statistical significance of our congruence hypothesis: the congruence between perceived self agency (\(X\), SelfAgency_X) and other agency (\(Y\),OtherAgency_Y) is related to higher happiness (\(Z\), Happy_Z).

The null hypothesis testing approach is identical to the first tutorial where we used the RSA package. The model comparison approach compares a freely estimated model, a model with constraints pertaining to the broad sense of congruence, and a model with constraints pertaining to the strict sense of congruence, which we demonstrate using sem(). The posterior distribution approach relies on the true distributions of the parameters of interests from brm() to make statistical judgments on the congruence hypothesis.

Null Hypothesis Testing

#create the table
table_lm <- data.frame(
  Parameter = c("a1", "a2", "a3", "a4", "a5"),
  Null_Hypothesis = c(
    "H<sub>0</sub>: <i>a</i><sub>1</sub> = 0", 
    "H<sub>0</sub>: <i>a</i><sub>2</sub> = 0", 
    "H<sub>0</sub>: <i>a</i><sub>3</sub> = 0", 
    "H<sub>0</sub>: <i>a</i><sub>4</sub> &lt; 0",
    "H<sub>0</sub>: <i>a</i><sub>5</sub> = 0"
  ),
  Alternative = c(
    "H<sub>1</sub>: <i>a</i><sub>1</sub> ≠ 0", 
    "H<sub>1</sub>: <i>a</i><sub>2</sub> ≠ 0", 
    "H<sub>1</sub>: <i>a</i><sub>3</sub> ≠ 0", 
    "H<sub>1</sub>: <i>a</i><sub>4</sub> ≠ 0",
    "H<sub>1</sub>: <i>a</i><sub>5</sub> ≠ 0"
  ),
  Point_Estimation = c(sprintf("%.3f", a1), sprintf("%.3f", a2), sprintf("%.3f", a3), 
                       sprintf("%.3f", a4), sprintf("%.3f", a5)),  # Keeping no column name
  SE = c(sprintf("%.3f", se_stand_a1), sprintf("%.3f", se_stand_a2), sprintf("%.3f", se_stand_a3), 
         sprintf("%.3f", se_stand_a4), sprintf("%.3f", se_stand_a5)),
  CI = c(paste0("[", sprintf("%.3f", ci_lower_a1), ", ", sprintf("%.3f", ci_upper_a1), "]"),
         paste0("[", sprintf("%.3f", ci_lower_a2), ", ", sprintf("%.3f", ci_upper_a2), "]"),
         paste0("[", sprintf("%.3f", ci_lower_a3), ", ", sprintf("%.3f", ci_upper_a3), "]"),
         paste0("[", sprintf("%.3f", ci_lower_a4), ", ", sprintf("%.3f", ci_upper_a4), "]"),
         paste0("[", sprintf("%.3f", ci_lower_a5), ", ", sprintf("%.3f", ci_upper_a5), "]")),
  p_value = c(sprintf("%.3f", p_value_a1), sprintf("%.3f", p_value_a2), 
              sprintf("%.3f", p_value_a3), sprintf("%.3f", p_value_a4), 
              sprintf("%.3f", p_value_a5))
)

#change column names
colnames(table_lm)[2] <- "Null Hypothesis"
colnames(table_lm)[3] <- "Alternative"
colnames(table_lm)[4] <- "Point Estimate"
colnames(table_lm)[5] <- "$SE$"
colnames(table_lm)[7] <- "$p$ value"

#print table using knitr for better formatting
table_lm %>% 
  kable(caption = "Summary Table for Parameters Transformed from lm()",
        format = "html", escape = FALSE) %>%
  column_spec(1, width = "2.5cm") %>%   # Adjust width for the "Parameter" column
  column_spec(2, width = "4cm") %>%    # Adjust width for the "Null Hypothesis" column
  column_spec(3, width = "4cm") %>%    # Adjust width for the "Alternative" column
  column_spec(4, width = "4cm") %>%     # Adjust width for "Point Estimate"
  column_spec(5, width = "3cm") %>%     # Adjust width for "SE"
  column_spec(6, width = "4cm") %>%     # Adjust width for "CI"
  column_spec(7, width = "3cm") %>%     # Adjust width for "*p* value"
  footnote(general = "$Note$: SE = Standard Error, CI = Confidence Interval.",
           general_title = "") %>%
  kable_styling(full_width = FALSE, html_font = "Times New Roman")
Summary Table for Parameters Transformed from lm()
Parameter Null Hypothesis Alternative Point Estimate \(SE\) CI \(p\) value
a1 H0: a1 = 0 H1: a1 ≠ 0 -0.142 0.032 [-0.204, -0.079] 0.000
a2 H0: a2 = 0 H1: a2 ≠ 0 -0.008 0.002 [-0.012, -0.004] 0.000
a3 H0: a3 = 0 H1: a3 ≠ 0 -0.016 0.035 [-0.084, 0.053] 0.650
a4 H0: a4 < 0 H1: a4 ≠ 0 -0.009 0.001 [-0.011, -0.006] 0.000
a5 H0: a5 = 0 H1: a5 ≠ 0 -0.000 0.002 [-0.004, 0.004] 0.911
\(Note\): SE = Standard Error, CI = Confidence Interval.

As we can see, the current data confirmed the broad sense of congruence (\(a_{3} = -0.016\), \(p = 0.650\); \(a_{4} = -0.009\), \(p < 0.001\); \(a_{5} = -0.000\), \(p = 0.911\)) with a linear and curvilinear main effects of the predictors (\(a_{1} = -0.142\), \(p < 0.001\); \(a_{2} = -0.008\), \(p < 0.001\)).

More specifically, this participant felt happier when the perceived other agency is more congruent with a given self agency. However, there is an optimum level of congruence - even when the perceived other agency is congruent with (or equal to) a given self agency, this participant would feel less happier if the perceived self agency (or the perceived other agency) is lower or higher congruent than the optimum level (i.e., stationary point \((X_{0}, Y_{0})\)).

Model Comparison for sem()

We have already specified all three models for comparison (i.e., freely estimated, broad congruence, and strict congruence). However, one still has to specify \(a_{4}\) and test the null hypothesis of \(a_{4} = 0\) as the alternative hypothesis of \(a_{4} < 0\) cannot be added as a model constraint for sem() (One could add such a constraint in the Bayesian SEM using bsem() from blavaan). Thus, we will first determine whether \(a_{4}\) is significantly less than zero following the null hypothesis testing approach. Then, we can use anova() to perform the likelihood ratio test (LRT) to determine whether the additional constraints in the model significantly worsen model fit.

# estimation of a4 for the freely estimated model
free_a4 <- parameterEstimates(sem_happy_agency)[
  parameterEstimates(sem_happy_agency)$label == "a4", ]

# estimation of a4 for broad congruence
broad_a4 <- parameterEstimates(sem_happy_agency_broad)[
  parameterEstimates(sem_happy_agency_broad)$label == "a4", ]

# estimation of a4 for strict congruence
strict_a4 <-parameterEstimates(sem_happy_agency_strict)[
  parameterEstimates(sem_happy_agency_strict)$label == "a4", ]

# Combine results into a single data frame with a new "Model" column
results_table <- data.frame(
  Model = c("Freely estimated", "Broad congruence", "Strict congruence"),
  Estimate = c(sprintf("%.3f", free_a4$est), 
               sprintf("%.3f", broad_a4$est), 
               sprintf("%.3f", strict_a4$est)),
  SE = c(sprintf("%.3f", free_a4$se), 
         sprintf("%.3f", broad_a4$se), 
         sprintf("%.3f", strict_a4$se)),
  p_value = c(sprintf("%.3f", free_a4$pvalue), 
                sprintf("%.3f", broad_a4$pvalue), 
                sprintf("%.3f", strict_a4$pvalue))
) 

colnames(results_table)[3] <- "$SE$"
colnames(results_table)[4] <- "$p$ value"

# print as a formatted table
kable(results_table, caption = "Summary of a4 across Models",
      format = "html", escape = FALSE) %>%
  column_spec(1, width = "4cm") %>%
  column_spec(2, width = "2.8cm") %>%
  column_spec(3, width = "2.8cm") %>%
  column_spec(4, width = "2.8cm") %>%
  footnote(general = "$Note$: SE = Standard Error.",
           general_title = "") %>%
  kable_styling(full_width = FALSE, html_font = "Times New Roman")
Summary of a4 across Models
Model Estimate \(SE\) \(p\) value
Freely estimated -0.009 0.001 0.000
Broad congruence -0.008 0.001 0.000
Strict congruence -0.007 0.001 0.000
\(Note\): SE = Standard Error. 
#model comparison
anova(sem_happy_agency, sem_happy_agency_broad, sem_happy_agency_strict) %>%
  mutate(across(where(is.numeric), ~sprintf("%.3f", .)),  # format to 3 decimals
         across(everything(), ~ifelse(is.na(.), "", .))) %>%  # replace NAs with blank spaces
  kable(caption = "Comparisons of Freely Estimated, Broad Congruence Constraints, 
        and Strict Sense Constraints Models",
        format = "html", escape = FALSE) %>%
  kable_styling(full_width = FALSE, html_font = "Times New Roman")
Comparisons of Freely Estimated, Broad Congruence Constraints, and Strict Sense Constraints Models
Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
sem_happy_agency 0.000 6604.038 6637.385 0.000 NA NA NA NA
sem_happy_agency_broad 2.000 6600.303 6624.123 0.265 0.265 0.000 2.000 0.876
sem_happy_agency_strict 4.000 6639.755 6654.047 43.717 43.452 0.155 2.000 0.000

As we can see, the estimation of \(a_{4}\) is significantly less than zero across all three models (\(ps < .001\)). Furthermore, adding constraints of broad congruence did not significantly worsen model fit (\(p = .876\)), while the model with constraints indicating the strict sense of congruence has significantly lower model fit (\(p < .001\)). This result aligns with what we saw from the null hypothesis approach. One could again print the estimation of \(a_{1}\) and \(a_{2}\) to interpret the main effects on the congruence to obtain a similar account of what we described using the results from null hypotheses.

Bayesian RSA (BRSA)

Load Necessary Libraries

Most libraries are already loaded for the demonstration in the previous packages, we still included them here for the use of this section as a stand-alone component for BRSA.

#new packages
library(shiny)            # interactive data visualization
library(loo)              # Leave-One-Out Cross-Validation and WAIC for Bayesian Models
#previously loaded packes
library(knitr)            # knitr tables
library(kableExtra)       # knitr table formatting
library(psych)            # data descriptives
library(RSA)              # Response Surface Analysis
library(plot3D)           # visualize response surfaces
library(plotly)           # visualize response surfaces
library(lavaan)           # structural equation modeling
library(brms)             # Bayesian parameter estimation
library(posterior)        # Bayesian posterior distribution
library(patchwork)        # add plots together
library(sandwich)         # calculate robust SEs
library(lmtest)           # linear regression tests
library(bayestestR)       # region of practical equivalence (ROPE)
library(see)              # visualize ROPE
library(parameters)       # parameter standardization 
library(tidyverse)        # data structures & visualization

Step 1: Specify Congruence Hypothesis

Similarly, we continue to use the iSAHIB dataset and explore the effect of the fit between the perceived agency of the self (predictor 1) and the perceived agency of the interaction partner (predictor 2) in dyadic interpersonal interactions on happiness (the outcome). We again test the congruent hypothesis - the more congruent between the perceived self agency and the perceived other agency, the happier oneself would be.

#set filepath for data file
filepath <- "https://raw.githubusercontent.com/The-Change-Lab/collaborations/refs/heads/main/iSAHIB_RSA/isahib_RSA.csv"

#read in the .csv file using the url() function
isahib <- read.csv(file=url(filepath), header=TRUE)

#look at the data
head(isahib, 10)
##    SelfAgency_X OtherAgency_Y Happy_Z
## 1            82            22      57
## 2            97            21      66
## 3            55            51      67
## 4            82            13      73
## 5            32            61      94
## 6            14            95      83
## 7            56            52      78
## 8            34            71      74
## 9            50            50      67
## 10           51            48      77
#describe the data
describe(isahib)
##               vars   n  mean    sd median trimmed   mad min max range  skew
## SelfAgency_X     1 866 51.51 15.76   53.0   51.85 14.83   2  97    95 -0.22
## OtherAgency_Y    2 866 55.93 17.07   57.5   56.45 17.05   4  99    95 -0.28
## Happy_Z          3 866 81.41 11.33   84.0   83.44  5.93   5  98    93 -2.91
##               kurtosis   se
## SelfAgency_X     -0.19 0.54
## OtherAgency_Y    -0.14 0.58
## Happy_Z          11.77 0.39
#plot the distribution of each variable and the correlation between them
psych::pairs.panels(isahib)

#plot raw data scatter plot in three dimensions
plot_ly(x = isahib$SelfAgency_X, y = isahib$OtherAgency_Y, z = isahib$Happy_Z, 
        type = "scatter3d", mode = "markers", color = isahib$Happy_Z) %>%
   layout(scene = list(xaxis=list(title = "Self Agency", nticks=10, range=c(0,100)),
                       yaxis=list(title = "Other Agency", nticks=10, range=c(0,100)),
                       zaxis=list(title = "Happiness", nticks=10, range=c(0,100))))
## Warning: Ignoring 3 observations

Step 2: Data Centering

We again center predictors on the scale midpoint (Barranti et al., 2017). In our case, this means centering at 50, the midpoint on a 0 to 100 scale.

#center predictors with the scale midpoint
isahib <- isahib %>%
  mutate(cSelfAgency_X = SelfAgency_X - 50,
         cOtherAgency_Y = OtherAgency_Y - 50)

Step 3: Fit the Polynomial Regression Model and Plot the Response Surface

#brm() function in the brms package
brm_happy_agency <- brm(formula = Happy_Z ~ cSelfAgency_X * cOtherAgency_Y + 
                          I(cSelfAgency_X^2) + I(cOtherAgency_Y^2),
                        data = isahib, family = gaussian(),
                        chains = 4, cores = 4, iter = 2000)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Start sampling

The specification of the regression model is similar to what we did in glm() with some noticeable differences.

First, we specified both chains = and cores = to 4. The Markov chain Monte Carlo (MCMC) chain arguments here refers to the number of the sequence of samples we wish to draw from the posterior distribution independently in parallel. It is helpful for us to assess the convergence to the posterior distribution and reliability of the results.

The core arguments indicates the number of CPU cores we intend to use for such computationally intensive MCMC simulations. It is suggest to set “as many processors as the hardware and RAM allow (up to the number of chains; Bürkner, 2017).”

Second, we have the iter = arguments and set it to the default number 2000. It indicates the number of total iterations per chain. To determine an appropriate number of interactions, we could start with the default 2000 iterations and then look at the output to assess convergence and determine the adequacy of the interactions.

Here, we did not set any prior distributions and used the default and uninformative priors in the brms package to estimate the model parameters. This is decision is motivated by the facts that we have relatively little knowledge about the distribution of the model parameters and that we have no subjective beliefs about the results. This decision should allow a more “data-driven” approach as we allow the likelihood to dominate the posterior.

# model output
summary(brm_happy_agency)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Happy_Z ~ cSelfAgency_X * cOtherAgency_Y + I(cSelfAgency_X^2) + I(cOtherAgency_Y^2) 
##    Data: isahib (Number of observations: 866) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Regression Coefficients:
##                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                       84.31      0.53    83.26    85.39 1.00     3470
## cSelfAgency_X                   -0.08      0.02    -0.13    -0.03 1.00     3416
## cOtherAgency_Y                  -0.06      0.02    -0.11    -0.02 1.00     3162
## IcSelfAgency_XE2                -0.00      0.00    -0.01    -0.00 1.00     5468
## IcOtherAgency_YE2               -0.00      0.00    -0.01    -0.00 1.00     4390
## cSelfAgency_X:cOtherAgency_Y     0.00      0.00    -0.00     0.00 1.00     5087
##                              Tail_ESS
## Intercept                        3119
## cSelfAgency_X                    2346
## cOtherAgency_Y                   2612
## IcSelfAgency_XE2                 3068
## IcOtherAgency_YE2                3209
## cSelfAgency_X:cOtherAgency_Y     2976
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    10.90      0.26    10.43    11.44 1.00     3182     2327
## 
## Draws were sampled 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).
# plot posterior distribution and trace of the chains
plot(brm_happy_agency)

We could see from the model output that the Gelman-Rubin statistic (R-hat; Gelman & Rubin, 1992) for all parameter estimates are \(\le\) 1.01, indicating good convergence (Vehtari et al., 2021; however, also see Vats & Knudson, 2018).

Again with the trace plot of the chains, we could see that they mixed well with each other and are relatively stable as they fluctuate within a similar range across the iterations after warm-up (which is usually half of the specified iterations).

# more decimals for comparison with model output of other functions
summary(brm_happy_agency)$fixed %>%
  as.data.frame() %>%
  mutate_at(vars(Estimate), round, digits = 5) %>%
  print()
##                              Estimate   Est.Error     l-95% CI     u-95% CI
## Intercept                    84.30957 0.534175922 83.264269380 85.386227705
## cSelfAgency_X                -0.07845 0.024837036 -0.128724619 -0.030342231
## cOtherAgency_Y               -0.06359 0.023244345 -0.108878072 -0.016710124
## IcSelfAgency_XE2             -0.00427 0.001267358 -0.006838591 -0.001763284
## IcOtherAgency_YE2            -0.00408 0.001096436 -0.006150086 -0.001859723
## cSelfAgency_X:cOtherAgency_Y  0.00020 0.001069532 -0.001917102  0.002283184
##                                   Rhat Bulk_ESS Tail_ESS
## Intercept                    1.0010006 3470.467 3118.946
## cSelfAgency_X                1.0017406 3415.911 2346.498
## cOtherAgency_Y               1.0002062 3161.539 2612.091
## IcSelfAgency_XE2             0.9995979 5468.017 3067.990
## IcOtherAgency_YE2            1.0010571 4389.805 3209.255
## cSelfAgency_X:cOtherAgency_Y 1.0009456 5087.019 2975.679

Instead of relying on standard errors, confidence intervals, or p-values, we often use credible interval to determine statistical significance for a Bayesian model.

One popular evaluation is to examine whether the 95% credible interval includes 0.

Later in the tutorial, we will also introduce another test that does not only involves the value of 0, but a region that is considered to be practically equivalent to 0. If 0 falls within the 95% credible interval, we often interpret the results as not statistically significant. As we can easily compare, the parameter estimates relatively align with what we obtained using OLS and ML.

brm() plot

#create data grid
x_grid <- seq(-52, 52, by = 1)
y_grid <- seq(-52, 52, by = 1)
newdat <- expand.grid(cSelfAgency_X = x_grid, cOtherAgency_Y = y_grid)

#posterior expected predictions on response scale
#posterior_epred returns predictions already on response scale
pred_epred <- posterior_epred(brm_happy_agency, newdata = newdat)

#calculate mean and 95% CI for each point
fit <- apply(pred_epred, 2, mean)
ci_low <- apply(pred_epred, 2, quantile, probs = 0.025)
ci_up  <- apply(pred_epred, 2, quantile, probs = 0.975)

#turn into matrices for plotly
z_matrix       <- matrix(fit, nrow = length(x_grid), byrow = TRUE)
ci_low_matrix  <- matrix(ci_low, nrow = length(x_grid), byrow = TRUE)
ci_up_matrix   <- matrix(ci_up, nrow = length(x_grid), byrow = TRUE)

#plot
plot_ly(x = x_grid, y = y_grid) %>% 
  add_surface(z = z_matrix,
              colorscale = list(c(0,1),c("red","blue"))) %>% 
  add_surface(z = ci_low_matrix, opacity = 0.5, showscale = FALSE, 
              colorscale = list(c(0,1),c("grey","grey"))) %>% 
  add_surface(z = ci_up_matrix, opacity = 0.5, showscale = FALSE, 
              colorscale = list(c(0,1),c("grey","grey")))

Step 4: Model Constraints & Posterior Distribution

We here specify the two different approaches for examining the congruence hypothesis - model comparisons and posterior distributions. We will first specify how to extract posterior distributions of the parameters of interest and then specify two additional models representing the broad congruence and the strict congruence.

Posterior Distributions

We here extract the posterior distributions of all five parameters of interests from brm(). Compared to hand calculation of matrices to obtain the standard error as well as the significance tests in the frequentist approach, the \(posterior\) offered an elegant way to obtaining the posterior draws of the transformed parameters directly.

a1

# extract posterior samples
posterior_draws <- posterior::as_draws_df(brm_happy_agency)

# calculate a1 <- b1 + b2
posterior_draws$a1 <- posterior_draws$'b_cSelfAgency_X' + 
  posterior_draws$'b_cOtherAgency_Y'

# mean of a1
mean(posterior_draws$a1)
## [1] -0.1420359
# 95% credible intervals for a1
quantiles_a1 <- quantile(posterior_draws$a1, probs = c(0.025, 0.975))
quantiles_a1
##        2.5%       97.5% 
## -0.20384763 -0.08031294
# plot the 95% credible intervals
#a1_95 <- posterior_draws$a1[posterior_draws$a1 >= quantiles_a1[1] & 
 #                             posterior_draws$a1 <= quantiles_a1[2]]

# plot the posterior distribution of a1
brm_a1 <- ggplot(data = as.data.frame(posterior_draws), aes(x = a1)) +
  geom_histogram(aes(y = ..density..), bins = 100, fill = "#7096ae", color = "black") +
  #geom_density(alpha = 0.3, fill = "dodgerblue") +
  geom_vline(xintercept = 0) + 
  geom_vline(xintercept = quantile(posterior_draws$a1, probs = 0.025)) +    
  annotate("text", x = quantile(posterior_draws$a1, probs = 0.001), 
           y = 10, label = "2.5%") + 
  geom_vline(xintercept = quantile(posterior_draws$a1, probs = 0.975)) +   
  annotate("text", x = quantile(posterior_draws$a1, probs = 0.999), 
           y = 10, label = "97.5%") + 
  labs(title = "Posterior Distribution of a1",
       x = "a1",
       y = "Density") +
  theme_minimal()

brm_a1
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

a2

# calculate a2 <- b3 + b4 + b5
posterior_draws$a2 <- posterior_draws$'b_IcSelfAgency_XE2' +
  posterior_draws$'b_cSelfAgency_X:cOtherAgency_Y' + 
  posterior_draws$'b_IcOtherAgency_YE2'

# mean and sd of a2
mean(posterior_draws$a2)
## [1] -0.008148231
sd(posterior_draws$a2)
## [1] 0.001927042
# 95% credible intervals for a2
quantiles_a2 <- quantile(posterior_draws$a2, probs = c(0.025, 0.975))
quantiles_a2
##         2.5%        97.5% 
## -0.012012504 -0.004351602
# plot the 95% credible intervals
#a2_95 <- posterior_draws$a2[posterior_draws$a2 >= quantiles_a2[1] & 
 #                             posterior_draws$a2 <= quantiles_a2[2]]

# plot the posterior distribution of a2
brm_a2 <- ggplot(data = as.data.frame(posterior_draws), aes(x = a2)) +
  geom_histogram(aes(y = ..density..), bins = 100, fill = "#7096ae", color = "black") +
  #geom_density(alpha = 0.3, fill = "dodgerblue") +
  geom_vline(xintercept = 0) + 
  geom_vline(xintercept = quantile(posterior_draws$a2, probs = 0.025)) +    
  annotate("text", x = quantile(posterior_draws$a2, probs = 0.001), 
           y = 150, label = "2.5%") + 
  geom_vline(xintercept = quantile(posterior_draws$a2, probs = 0.975)) +   
  annotate("text", x = quantile(posterior_draws$a2, probs = 0.999), 
           y = 150, label = "97.5%") + 
  labs(title = "Posterior Distribution of a2",
       x = "a2",
       y = "Density") +
  theme_minimal()

brm_a2

a3

# calculate a3 <- b1 - b2
posterior_draws$a3 <- posterior_draws$b_cSelfAgency_X - posterior_draws$b_cOtherAgency_Y 

# mean and sd of a3
mean(posterior_draws$a3)
## [1] -0.01486413
sd(posterior_draws$a3)
## [1] 0.03614357
# 95% credible intervals for a3
quantiles_a3 <- quantile(posterior_draws$a3, probs = c(0.025, 0.975))
quantiles_a3
##        2.5%       97.5% 
## -0.08369855  0.05618696
# plot the 95% credible intervals
#a3_95 <- posterior_draws$a3[posterior_draws$a3 >= quantiles_a3[1] & 
 #                             posterior_draws$a3 <= quantiles_a3[2]]

# plot the posterior distribution of a3
brm_a3 <- ggplot(data = as.data.frame(posterior_draws), aes(x = a3)) +
  geom_histogram(aes(y = ..density..), bins = 100, fill = "#7096ae", color = "black") +
  #geom_density(alpha = 0.3, fill = "dodgerblue") +
  geom_vline(xintercept = 0) + 
  geom_vline(xintercept = quantile(posterior_draws$a3, probs = 0.025)) +    
  annotate("text", x = quantile(posterior_draws$a3, probs = 0.001), 
           y = 8, label = "2.5%") + 
  geom_vline(xintercept = quantile(posterior_draws$a3, probs = 0.975)) +   
  annotate("text", x = quantile(posterior_draws$a3, probs = 0.999), 
           y = 8, label = "97.5%") + 
  labs(title = "Posterior Distribution of a3",
       x = "a3",
       y = "Density") +
  theme_minimal()

brm_a3

a4

# calculate a4 = b3 - b4 + b5
posterior_draws$a4 <- posterior_draws$'b_IcSelfAgency_XE2' -
  posterior_draws$'b_cSelfAgency_X:cOtherAgency_Y' + 
  posterior_draws$'b_IcOtherAgency_YE2'

# mean and sd of a4
mean(posterior_draws$a4)
## [1] -0.008554372
sd(posterior_draws$a4)
## [1] 0.001404694
# 95% credible intervals for a4
quantiles_a4 <- quantile(posterior_draws$a4, probs = c(0.025, 0.975))
quantiles_a4
##         2.5%        97.5% 
## -0.011270290 -0.005758328
# plot the 95% credible intervals
#a4_95 <- posterior_draws$a4[posterior_draws$a4 >= quantiles_a4[1] & 
 #                             posterior_draws$a4 <= quantiles_a4[2]]

# plot the posterior distribution of a4
brm_a4 <- ggplot(data = as.data.frame(posterior_draws), aes(x = a4)) +
  geom_histogram(aes(y = ..density..), bins = 100, fill = "#7096ae", color = "black") +
  #geom_density(alpha = 0.3, fill = "dodgerblue") +
  geom_vline(xintercept = 0) + 
  geom_vline(xintercept = quantile(posterior_draws$a4, probs = 0.025)) +    
  annotate("text", x = quantile(posterior_draws$a4, probs = 0.0001), 
           y = 250, label = "2.5%") + 
  geom_vline(xintercept = quantile(posterior_draws$a4, probs = 0.975)) +   
  annotate("text", x = quantile(posterior_draws$a4, probs = 0.999), 
           y = 250, label = "97.5%") + 
  labs(title = "Posterior Distribution of a4",
       x = "a4",
       y = "Density") +
  theme_minimal()

brm_a4

a5

# calculate a5 = b3 - b5
posterior_draws$a5 <- posterior_draws$'b_IcSelfAgency_XE2' -
  posterior_draws$'b_IcOtherAgency_YE2'

# mean and sd of a5
mean(posterior_draws$a5)
## [1] -0.0001821705
sd(posterior_draws$a5)
## [1] 0.00197922
# 95% credible intervals for a5
quantiles_a5 <- quantile(posterior_draws$a5, probs = c(0.025, 0.975))
quantiles_a5
##         2.5%        97.5% 
## -0.004152232  0.003706022
# plot the 95% credible intervals
#a5_95 <- posterior_draws$a5[posterior_draws$a5 >= quantiles_a5[1] & 
 #                             posterior_draws$a5 <= quantiles_a5[2]]

# plot the posterior distribution of a5
brm_a5 <- ggplot(data = as.data.frame(posterior_draws), aes(x = a5)) +
  geom_histogram(aes(y = ..density..), bins = 100, fill = "#7096ae", color = "black") +
  #geom_density(alpha = 0.3, fill = "dodgerblue") +
  geom_vline(xintercept = 0) + 
  geom_vline(xintercept = quantile(posterior_draws$a5, probs = 0.025)) +    
  annotate("text", x = quantile(posterior_draws$a5, probs = 0.001), 
           y = 200, label = "2.5%") + 
  geom_vline(xintercept = quantile(posterior_draws$a5, probs = 0.975)) +   
  annotate("text", x = quantile(posterior_draws$a5, probs = 0.999), 
           y = 200, label = "97.5%") + 
  labs(title = "Posterior Distribution of a5",
       x = "a5",
       y = "Density") +
  theme_minimal()

brm_a5

Constrained Model Specification

We here estimate two models with different constraints corresponding to the broad sense and the strict sense of congruence for model comparisons.

Broad sense

# a3 = b1 - b2 = 0
# b1 = b2

# a5 = b3 - b5 = 0
# b3 = b5

# a4 = b3 - b4 + b3 < 0

brm_broad <- bf(
  Happy_Z ~ b0 * 1 + 
    b1 * cSelfAgency_X + b1 * cOtherAgency_Y +
    b3 * cSelfAgency_X^2 + 
    b4 * cSelfAgency_X * cOtherAgency_Y + 
    b3 * cOtherAgency_Y^2,
  b0 ~ 1,
  b1 ~ 1,
  b3 ~ 1,
  b4 ~ 1,
  nl = TRUE
)

brm_happy_agency_broad <- brm(
  formula = brm_broad,
  data = isahib,
  family = gaussian(),
  chains = 4,
  cores = 4,
  iter = 2000)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Start sampling
summary(brm_happy_agency_broad)

Strict sense

# a1 = b1 + b2 = 0 
# a3 = b1 - b2 = 0
# b1 = 0
# b2 = 0
# b3 = b5
# b4 == -2 * b3

# a4 = b3 - b4 + b3 < 0

brm_strict <- bf(
  Happy_Z ~ b0 * 1 + 
    0 * cSelfAgency_X + 0 * cOtherAgency_Y +
    b3 * cSelfAgency_X^2 + 
    -2 * b3 * cSelfAgency_X * cOtherAgency_Y + 
    b3 * cOtherAgency_Y^2,
  b0 ~ 1,
  b3 ~ 1,
  nl = TRUE
)

brm_happy_agency_strict <- brm(
  formula = brm_strict,
  data = isahib,
  family = gaussian(),
  chains = 4,
  cores = 4,
  iter = 2000)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Start sampling
summary(brm_happy_agency_strict)

Step 5: Congruence Hypothesis Testing and Interpretation

Posterior Distributions for brm()

We here introduce two ways to interpret the results of the transformed parameter estimations using posterior distributions. Both of which involve testing the 95% credible intervals (often referred to as the high density region, Kruschke, 2018) but against different criterion: a single value of 0 or a region that is considered practically equivalent to being 0.

The first approach is rather similar to the frequentist approach in testing parameter significance. We can evaluate whether the 95% credible interval includes the value of 0 to make a decision of whether the estimated parameter is significantly different from 0 or not. The following plot summarizes the distributions of the five parameters of interest with their 95% credible intervals as well as a third black vertical line indicating the value of 0.

brm_a1 + brm_a2 +
  brm_a3 + brm_a4 + 
  brm_a5 + 
  plot_layout(ncol = 2)

According to the plot, we can easily tell that \(a_{1}\), \(a_{2}\), and \(a_{4}\) are less than 0 as their 95% credible intervals are all located on the left side of zero while \(a_{3}\) and \(a_{5}\) are not significantly different from 0. Those results are the same as the results obtained earlier from the frequentist approach and does not make direct claims of whether a certain parameter is approximating zero.

We thus demonstrate the second approach where we can make use of the Bayesian parameter estimation with posterior distributions to make claims of whether a parameter is practically equivalent to zero.

Unlike the congruence hypothesis testing in the frequentist approach, Bayesian parameter estimation offers the probabilistic distribution of the parameters. Instead of testing again one value (often zero) in the frequentist approach, the Bayesian framework also allowed us to assess the uncertainty on the estimation of parameters and test the posterior distributions against a specific range, where we consider the range represents practically no effect (or a negligible magnitude). This range is the region of practical equivalence (ROPE; Makowski et al., 2019).

When we test the congruence hypothesis against the null hypothesis (e.g., a certain parameter equals zero), the Bayesian approach tests again the null region (e.g., whether the 95% credible interval is included in ROPE).

We followed the recommendation of previous studies and used 95% credible intervals along with a “-0.1 to 0.1” range (Cohen, 1988; Kruschke & Liddell, 2018). In this case, ROPE of a parameter is calculated as \([-0.1*SD_{outcome}, 0.1*SD_{outcome}]\), where SD represents the standard deviation. To make the the interpretation of congruence more meaningful, researchers should consider the unit of change that is of interest, or that can be important practically.

To summarize the posterior distributions and ROPE for \(a_{1}\) to \(a_{5}\), we plot all five distributions in one panel with 95% credible intervals (black lines) and ROPE (red rectangular regions). We also summarized their means, standard deviations, 95% credible intervals, and ROPE in one table to facilitate interpretation.

brm_a1 + 
  annotate("rect", 
           xmin = -0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           xmax =  0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           ymin = 0, ymax = Inf, 
           fill = "#F27781", alpha = 0.7) + 
  brm_a2 +
  annotate("rect", 
           xmin = -0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           xmax =  0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           ymin = 0, ymax = Inf, 
           fill = "#F27781", alpha = 0.7) + 
  brm_a3 + 
  annotate("rect", 
           xmin = -0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           xmax =  0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           ymin = 0, ymax = Inf, 
           fill = "#F27781", alpha = 0.7) + 
  brm_a4 + 
  annotate("rect", 
           xmin = -0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           xmax =  0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           ymin = 0, ymax = Inf, 
           fill = "#F27781", alpha = 0.7) + 
brm_a5 + 
  annotate("rect", 
           xmin = -0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           xmax =  0.1 * sd(isahib$Happy_Z, na.rm = TRUE), 
           ymin = 0, ymax = Inf, 
           fill = "#F27781", alpha = 0.7) + 
  plot_layout(ncol = 2)

#parameter names
parameters <- c("a1", "a2", "a3", "a4", "a5")

#mean and standard deviation
mean_values <- c(mean(posterior_draws$a1), mean(posterior_draws$a2), 
                 mean(posterior_draws$a3), mean(posterior_draws$a4), 
                 mean(posterior_draws$a5))

sd_values <- c(sd(isahib$Happy_Z, na.rm = TRUE), sd(isahib$Happy_Z, na.rm = TRUE), 
               sd(isahib$Happy_Z, na.rm = TRUE), sd(isahib$Happy_Z, na.rm = TRUE), 
               sd(isahib$Happy_Z, na.rm = TRUE))

#95% credible intervals
quantiles_list <- list(
  quantiles_a1, quantiles_a2, quantiles_a3, 
  quantiles_a4, quantiles_a5)

#compute ROPE as [-0.1*SD, 0.1*SD]
rope_intervals <- lapply(sd_values, function(sd) c(-0.1 * sd, 0.1 * sd))

#table
result_table <- data.frame(
  Parameter = parameters,
  Mean = sprintf("%.4f", mean_values),
  SD = sprintf("%.4f", sd_values),
  `95% CIs` = sapply(quantiles_list, function(q) {
    paste0("[", sprintf("%.4f", q["2.5%"]), ", ", 
           sprintf("%.4f", q["97.5%"]), "]")}),
  `ROPE` = sapply(rope_intervals, function(r) {
    paste0("[", sprintf("%.4f", r[1]), ", ", 
           sprintf("%.4f", r[2]), "]")
})) 

colnames(result_table) <- c("Parameter", "Mean", "SD", "95% CIs", "ROPE")

result_table %>%
  kable(caption = "Summary of Bayesian posterior distributions and the region of practical equivalence",
      format = "html", escape = FALSE) %>%
  column_spec(1, width = "2.8cm") %>%
  column_spec(2, width = "2.8cm") %>%
  column_spec(3, width = "2.8cm") %>%
  column_spec(4, width = "4cm") %>%
  column_spec(5, width = "4cm") %>%
  footnote(general = "$Note$: SD = Standard Deviation; CIs = Credible Intervals.",
           general_title = "") %>%
  kable_styling(full_width = FALSE, html_font = "Times New Roman")
Summary of Bayesian posterior distributions and the region of practical equivalence
Parameter Mean SD 95% CIs ROPE
a1 -0.1420 11.3321 [-0.2038, -0.0803] [-1.1332, 1.1332]
a2 -0.0081 11.3321 [-0.0120, -0.0044] [-1.1332, 1.1332]
a3 -0.0149 11.3321 [-0.0837, 0.0562] [-1.1332, 1.1332]
a4 -0.0086 11.3321 [-0.0113, -0.0058] [-1.1332, 1.1332]
a5 -0.0002 11.3321 [-0.0042, 0.0037] [-1.1332, 1.1332]
\(Note\): SD = Standard Deviation; CIs = Credible Intervals.

As illustrated in the Figure and Table, the 95% credible intervals for all five parameters are located within the predetermined ROPE, meaning all five parameters are equivalent to 0 if we are to assume that only effect sizes larger than 0.2 are significant (Cohen’s definition of a small effect).

Contrary to the previous results, we now do not have enough evidence to accept the congruence hypothesis.

More interestingly, the ROPE approach would suggest that perceived agency of the self (\(X\), SelfAgency_X) and the other (\(Y\), SelfAgency_Y) are not practically meaningful predictors for happiness (\(Z\), Happy_Z) for this given data. Using the ROPE approach in the Bayesian framework would urge us to predefine and consider the practical meaning of the parameters and relationship between variables - a goal that aligns with the purpose of RSA in its original development to maximize an outcome.

Decide the ROPE

A detailed discussion on how to decide the “right” ROPE is outside the scope of the current tutorial. One could refer to Kruschke (2018) for a relatively more detailed discussion of defining ROPE limits. Though we do not expand on the discussion of deciding the ROPE limits, we present the visualization of how the change in the ROPE limits could alter the hypothesis testing using the Shiny app (https://zhenchaohu.shinyapps.io/rsa2_shinyapp_may12/).

The decision of “what is practically equivalent to zero” differs a lot depending on the outcome and the research question (e.g., randing from 1.25 as a FDA bioequivalence recommendation to a 0.2 industry conventions, see Kruschke, 2018).

For testing congruence hypothesis, we strongly recommend researchers to use predefined ROPE limits based on relevant literature (e.g., what is the expected effect size) and real-world outcome (e.g., how much of a change in the outcome is meaningful/risky) instead of adopting a post-hoc approach that risks unethical research practice such as p-hacking.

Model Comparisons for brm()

To compare all three Bayesian models (i.e., freely estimated, broad congruence, and strict congruence), we adopted the leave-one-out cross validation (LOO) proposed by Vehtari and colleagues (2017) using the loo package (Vehtari et al., 2024). Essentially, we are comparing the out-of-sample prediction accuracy of the fitted BRSA models - which model can do a better job in predicting the outcome.

loo_compare(
  loo(brm_happy_agency),
  loo(brm_happy_agency_broad),
  loo(brm_happy_agency_strict)
)
## Warning: Found 1 observations with a pareto_k > 0.7 in model
## 'brm_happy_agency'. We recommend to set 'moment_match = TRUE' in order to
## perform moment matching for problematic observations.
##                         elpd_diff se_diff
## brm_happy_agency_broad    0.0       0.0  
## brm_happy_agency         -2.8       0.8  
## brm_happy_agency_strict -18.5      10.6

As shown in the table, the results of the LOO for Bayesian model comparison compute two columns: the estimated difference in predictive accuracy (difference of expected log predictive density, elpd_diff) compared to the top rank model and the standard error of the difference in expected log predictive density (i.e., se_diff).

For our models here, the best performing model is the broad congruence model as negative elpd_diff indicates the model has less predictive accuracy compared to the broad congruence model.

However, following the guidelines by Sivula and colleagues (2023) and Vehtari (2024), differences smaller than 4 in the elpd scale for the freely estimated model (i.e., brm_happy_agency) are usually negligible and large difference in standard error (e.g., se_diff \(>\) \(1.96*|\)elpd_diff\(|\)) represent that the model (i.e., brm_happy_agency_strict where \(1.96*10.7=20.972>|-18.3|\)) is not meaningful different.

This result is in line with what we obtained with ROPE where we do not have enough evidence to accept the hypothesis and that the current models are not ideal for meaningful predicting the outcome.

Conclusion

We have achieved fitting, visualizing, and interpreting an RSA model with more flexibility in model specification (i.e., lm(), glm(), and sem()) using null hypothesis and model comparison by hand.

More importantly, we proposed the Bayesian RSA framework and accomplished using posterior distributions and predictive accuracy to examine the congruence hypothesis with the brms package. Further, as one follow along the tutorials, they are basis for specifying the polynomial regressions in more complex data structures (e.g., within-person or within-dyad dependencies) and including more nonlinear dynamics in parameter estimation. Have fun!

References

Box, G. E. P., & Wilson, K. B. (1951). On the Experimental Attainment of Optimum Conditions. Journal of the Royal Statistical Society: Series B (Methodological), 13(1), 1–38. https://doi.org/10.1111/j.2517-6161.1951.tb00067.x

Bürkner, P.-C. (2017). brms: An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software, 80, 1–28. https://doi.org/10.18637/jss.v080.i01

Bürkner, P.-C., Gabry, J., Kay, M., & Vehtari, A. (2025). posterior: Tools for Working with Posterior Distributions (Version 1.6.1). https://mc-stan.org/posterior/

Caldwell, D. F., & O’Reilly III, C. A. (1990). Measuring person-job fit with a profile-comparison process. Journal of Applied Psychology, 75(6), 648–657. https://doi.org/10.1037/0021-9010.75.6.648

C. Feng, G. (2017). Do difference scores make a difference on the third-person effect? Communications in Statistics - Simulation and Computation, 46(7), 5085–5104. https://doi.org/10.1080/03610918.2016.1143104

Fischer, J. R. (2024). Bayesian methods tutorials for psychological research [Undergraduate honors thesis, Stanford University]. https://thechangelab.stanford.edu/tutorials/bayesian-methods/

Hu, Z. (2025). Dynamic ROPE Visualization for Posterior Parameters (a1–a5) [R Shiny Application]. Dynamic ROPE Visualization for Posterior Parameters (A1–A5). Retrieved May 16, 2025, from https://zhenchaohu.shinyapps.io/rsa2_shinyapp_may12/

Kenny, D. A., Kashy, D. A., & Cook, W. L. (2006). Dyadic Data Analysis. Guilford Press.

Khuri, A. I., & Cornell, J. A. (1996). Response Surfaces: Designs and Analyses. Taylor & Francis Group.

Edwards, J. R. (2002). Alternatives to difference scores: Polynomial regression analysis and response surface methodology. In F. Drasgow & N. Schmitt (Eds.), Measuring and analyzing behavior in organizations: Advances in measurement and data analysis (pp. 350–400). Jossey-Bass/Wiley.

Edwards, J. R., & Parry, M. E. (1993). On the use of polynomial regression equations as an alternative to difference scores in organizational research. Academy of Management Journal, 36(6), 1577–1613. https://www.jstor.org/stable/256822

Gelman, A., & Rubin, D. B. (1992). Inference from Iterative Simulation Using Multiple Sequences. Statistical Science, 7(4), 457–472. https://doi.org/10.1214/ss/1177011136

Hoffman, J. A. (1984). Psychological separation of late adolescents from their parents. Journal of Counseling Psychology, 31(2), 170–178. https://doi.org/10.1037/0022-0167.31.2.170

Humberg, S., Dufner, M., Schönbrodt, F. D., Geukes, K., Hutteman, R., Küfner, A. C. P., van Zalk, M. H. W., Denissen, J. J. A., Nestler, S., & Back, M. D. (2019). Is accurate, positive, or inflated self-perception most advantageous for psychological adjustment? A competitive test of key hypotheses. Journal of Personality and Social Psychology, 116(5), 835–859. https://doi.org/10.1037/pspp0000204

Humberg, S., Nestler, S., & Back, M. D. (2019). Response Surface Analysis in Personality and Social Psychology: Checklist and Clarifications for the Case of Congruence Hypotheses. Social Psychological and Personality Science, 10(3), 409–419. https://doi.org/10.1177/1948550618757600

Humberg, S., Schönbrodt, F. D., Back, M. D., & Nestler, S. (2022). Cubic response surface analysis: Investigating asymmetric and level-dependent congruence effects with third-order polynomial models. Psychological Methods, 27(4), 622–649. https://doi.org/10.1037/met0000352

Kruschke, J. K. (2010). Bayesian data analysis. WIREs Cognitive Science, 1(5), 658–676. https://doi.org/10.1002/wcs.72

Kruschke, J. K. (2018). Rejecting or Accepting Parameter Values in Bayesian Estimation. Advances in Methods and Practices in Psychological Science, 1(2), 270–280. https://doi.org/10.1177/2515245918771304

Kruschke, J. K., & Liddell, T. M. (2018). The Bayesian New Statistics: Hypothesis testing, estimation, meta-analysis, and power analysis from a Bayesian perspective. Psychonomic Bulletin & Review, 25(1), 178–206. https://doi.org/10.3758/s13423-016-1221-4

Lüdecke, D., Ben-Shachar, M. S., Patil, I., & Makowski, D. (2020). Extracting, Computing and Exploring the Parameters of Statistical Models using R. Journal of Open Source Software, 5(53), 2445. https://doi.org/10.21105/joss.02445

Lüdecke, D., Patil, I., Ben-Shachar, M. S., Wiernik, B. M., Waggoner, P., & Makowski, D. (2021). see: An R Package for Visualizing Statistical Models. Journal of Open Source Software, 6(64), 3393. https://doi.org/10.21105/joss.03393

Makowski, D., Ben-Shachar, M. S., & Lüdecke, D. (2019). bayestestR: Describing Effects and their Uncertainty, Existence and Significance within the Bayesian Framework. Journal of Open Source Software, 4(40), 1541. https://doi.org/10.21105/joss.01541

Moskowitz, D. S., & Zuroff, D. C. (2005). Assessing Interpersonal Perceptions Using the Interpersonal Grid. Psychological Assessment, 17(2), 218–230. https://doi.org/10.1037/1040-3590.17.2.218

Pedersen, T. (2024). patchwork: The Composer of Plots (Version 1.3.0). https://CRAN.R-project.org/package=patchwork

R Core Team., (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

Revelle, W. (2024). psych: Procedures for Psychological, Psychometric, and Personality Research. Northwestern University. https://CRAN.R-project.org/package=psych

Satorra, A., & Bentler, P. M. (1994). Corrections to test statistics and standard errors in covariance structure analysis. In A. von Eye & C. C. Clogg (Eds.), Latent variables analysis: Applications for developmental research (pp. 399–419). Sage Publications, Inc.

Schönbrodt, F. D., & Humberg, S. (2023). RSA: An R package for response surface analysis (Version 0.10.6). https://cran.r-project.org/package=RSA

Sievert, C. (2020). Interactive Web-Based Data Visualization with R, plotly, and shiny. Chapman and Hall/CRC. https://doi.org/10.1201/9780429447273

Sivula, T., Magnusson, M., Matamoros, A. A., & Vehtari, A. (2023). Uncertainty in Bayesian Leave-One-Out Cross-Validation Based Model Comparison (arXiv:2008.10296). arXiv. https://doi.org/10.48550/arXiv.2008.10296

Soetaert, K. (2024). plot3D: Plotting Multi-Dimensional Data (Version 1.4.1). https://CRAN.R-project.org/package=plot3D

Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48, 1–36. https://doi.org/10.18637/jss.v048.i02

Trinh, T. K., & Kang, L. S. (2010). Application of Response Surface Method as an Experimental Design to Optimize Coagulation Tests. Environmental Engineering Research, 15(2), 63–70.

Vats, D., & Knudson, C. (2020). Revisiting the Gelman-Rubin Diagnostic (arXiv:1812.09384). arXiv. https://doi.org/10.48550/arXiv.1812.09384

Vehtari, A. (2023, December). Model checking: Comparison using loo vs. Loo_compare [Online post]. Stan Forums. https://discourse.mc-stan.org/t/model-checking-comparison-using-loo-vs-loo-compare/33694

Vehtari, A., Gabry, J., Magnusson, M., Yao, Y., Bürkner, P.-C., Paananen, T., & Gelman, A. (2024). loo: Efficient leave-one-out cross validation and WAIC for Bayesian models (Version 2.8.0). https://mc-stan.org/loo/

Vehtari, A., Gelman, A., & Gabry, J. (2017). Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Statistics and Computing, 27(5), 1413–1432. https://doi.org/10.1007/s11222-016-9696-4

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2021). Rank-Normalization, Folding, and Localization: An Improved Rˆ for Assessing Convergence of MCMC (with Discussion). Bayesian Analysis, 16(2), 667–718. https://doi.org/10.1214/20-BA1221

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., Grolemund, G., Hayes, A., Henry, L., Hester, J., Kuhn, M., Pedersen, T. L., Miller, E., Bache, S. M., Müller, K., Ooms, J., Robinson, D., Seidel, D. P., Spinu, V., … Yutani, H. (2019). Welcome to the Tidyverse. Journal of Open Source Software, 4(43), 1686. https://doi.org/10.21105/joss.01686

Xie, Y. (2025). knitr: A General-Purpose Package for Dynamic Report Generation in R (Version 1.50). https://yihui.org/knitr/

Zeileis, A. (2004). Econometric Computing with HC and HAC Covariance Matrix Estimators. Journal of Statistical Software, 11, 1–17. https://doi.org/10.18637/jss.v011.i10

Zeileis, A. (2006). Object-oriented Computation of Sandwich Estimators. Journal of Statistical Software, 16, 1–16. https://doi.org/10.18637/jss.v016.i09

Zeileis, A., & Hothorn, T. (2002). Diagnostic Checking in Regression Relationships. R News, 2(3), 7–10.

Zeileis, A., Köll, S., & Graham, N. (2020). Various Versatile Variances: An Object-Oriented Implementation of Clustered Covariances in R. Journal of Statistical Software, 95, 1–36. https://doi.org/10.18637/jss.v095.i01

Zhu, H. (2024). kableExtra: Construct Complex Table with “kable” and Pipe Syntax (Version 1.4.0). https://CRAN.R-project.org/package=kableExtra