Dyadic Analysis Workshop, SPR 2025

Author

Eran Bar-Kalifa, Ben-Gurion Univ. Israel

Intro to 3-level MLM

Libraries

pacman::p_load(dplyr,tidyr, rio,ggplot2,psych, nlme,GGally,sjPlot,osfr,vtable,emmeans,RSA,car)

Loading the data from OSF

file <- osf_retrieve_file("68467733c28a254a2fd6a394")
#Change to your local dowonloads folder
setwd("C:\\Users\\97254\\Downloads")
osf_download(file, progress = T,conflicts = "overwrigt") 
temp1 <- import("SimData3Level.csv")

names(temp1) <- c("Alliance","PosEmo","Client","Therapist")
temp1 <- temp1 %>% group_by(Therapist,Client) %>% 
  mutate(Session=row_number()) %>% 
  ungroup()
setwd("G:/My Drive/Research/Conferences & Workshops/2025- SPR/Workshop")

Descriptive stats

head(temp1)

# A tibble: 6 × 5
  Alliance PosEmo Client Therapist Session
     <dbl>  <dbl>  <int>     <int>   <int>
1    1.22    5.84      1         1       1
2    3.02    3.45      1         1       2
3    0.739   5.12      1         1       3
4    2.08    5.50      1         1       4
5    1.34    5.61      1         1       5
6    0.420   3.68      1         1       6

st(temp1)

Summary Statistics
Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 75	Max
Alliance	16000	6	3.1	-8	3.8	8	18
PosEmo	16000	4	1	-0.12	3.3	4.7	7.8
Client	16000	500	289	1	251	750	1000
Therapist	16000	50	29	1	26	75	100
Session	16000	8.5	4.6	1	4.8	12	16

Association between Clients’ Positive Emotions and Alliance

ggplot(temp1 %>% filter(Client < 500), aes(x = PosEmo, y = Alliance)) +
  geom_point(alpha = 0.6) +  # Scatter plot
  geom_smooth(method = "lm", se = FALSE, color = "red", size = 1.2) +  # Regression line
  labs(
    title = "Scatter Plot with Regression Line",
    x = "Positive Emotions",
    y = "Alliance"
  ) +
  theme_minimal()

Regression model

corr.test(temp1 %>% select(Alliance,PosEmo))

Call:corr.test(x = temp1 %>% select(Alliance, PosEmo))
Correlation matrix 
         Alliance PosEmo
Alliance     1.00   0.09
PosEmo       0.09   1.00
Sample Size 
[1] 16000
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
         Alliance PosEmo
Alliance        0      0
PosEmo          0      0

 To see confidence intervals of the correlations, print with the short=FALSE option

Model1 <- lm(Alliance~PosEmo,data=temp1)
tab_model(Model1,show.r2=F)

	Alliance
Predictors	Estimates	CI	p
(Intercept)	4.85	4.65 – 5.05	<0.001
PosEmo	0.28	0.23 – 0.32	<0.001
Observations	16000

Scatter plot for the first 10 clients

ggplot(temp1 %>% filter(Client <= 10), aes(x = PosEmo, y = Alliance)) +
  geom_point(alpha = 0.6,aes(color=as.factor(Client))) +  # Scatter plot
  geom_smooth(method = "lm", se = FALSE, color = "red", size = 1.2) +  # Regression line
  labs(
    title = "Scatter Plot with Regression Line",
    x = "Positive Emotions",
    y = "Alliance"
  ) +
  theme_minimal()

temp1 <- temp1 %>% mutate(residuals=Model1$residuals, 
                          fitted=Model1$fitted.values)
# Plot residuals vs fitted, color by therapist
ggplot(temp1 %>% filter(Client<10), aes(x = fitted, y = residuals, color = as.factor(Client))) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Residuals vs Fitted: OLS ignoring nesting",
       x = "Fitted values", y = "Residuals") +
  theme_minimal()

Spaghetti Plot Per Client

ggplot(temp1 %>% filter(Client<50), aes(x = PosEmo, y = Alliance, group = Client)) +
  #geom_line(alpha = 0.5) +
  geom_smooth(method = lm, se = FALSE, fullrange = TRUE, aes(group = interaction(Client)), color = "gray", size = 0.5) +  # Individual regression lines
  geom_smooth(aes(group = 1), method = "lm", se = FALSE, color = "red", size = 1.2) +
  labs(title = "Spaghetti Plot Per Client",
       x = "Positive Emotions", y = "Alliance") +
  theme_minimal()

Model2 <- lme(Alliance~PosEmo, 
              random = ~ 1 + PosEmo| Therapist/Client,
              data = temp1)

foo <- as.data.frame(Model2$coefficients$random$Therapist)
Model2 <- lme(Alliance~PosEmo, 
              random = ~ 1 + PosEmo| Therapist/Client,
              data = temp1)
ggplot(temp1 %>% mutate(THERAPIST=as.factor(Therapist)) %>% 
         filter(Therapist%in%c(which.max(foo$PosEmo),
                               which.min(foo$PosEmo))), 
       aes(x = PosEmo, y = Alliance, group = Client)) +
  geom_smooth(
    method = lm, se = FALSE, fullrange = TRUE,
    aes(group = interaction(Client), color = THERAPIST),
    size = 0.5
  ) +  # Individual regression lines colored by therapist
  geom_smooth(
    aes(group = 1),
    method = "lm", se = FALSE, color = "red", size = 1.2
  ) +  # Overall regression line
  labs(
    title = "Spaghetti Plot Per Client and Therapist",
    x = "Positive Emotions", y = "Alliance"
  ) +
  theme_minimal()

3-level multi-level model

Model2 <- lme(Alliance~PosEmo, 
              random = ~ 1 + PosEmo| Therapist/Client,
              data = temp1)
summary(Model2)

Linear mixed-effects model fit by REML
  Data: temp1 
       AIC      BIC    logLik
  49516.24 49585.36 -24749.12

Random effects:
 Formula: ~1 + PosEmo | Therapist
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 1.0331539 (Intr)
PosEmo      0.3157697 0.075 

 Formula: ~1 + PosEmo | Client %in% Therapist
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 0.9765288 (Intr)
PosEmo      0.4621859 0.382 
Residual    0.9610456       

Fixed effects:  Alliance ~ PosEmo 
               Value  Std.Error    DF  t-value p-value
(Intercept) 4.903089 0.11278744 14999 43.47194       0
PosEmo      0.262182 0.03571402 14999  7.34115       0
 Correlation: 
       (Intr)
PosEmo 0.04  

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max 
-3.649017326 -0.647822488  0.001778542  0.646254809  3.701130771 

Number of Observations: 16000
Number of Groups: 
            Therapist Client %in% Therapist 
                  100                  1000

tab_model(Model2,show.r2=F)

	Alliance
Predictors	Estimates	CI	p
(Intercept)	4.90	4.68 – 5.12	<0.001
PosEmo	0.26	0.19 – 0.33	<0.001
Random Effects
σ²	0.92
τ₀₀ _Client	1.03
τ₀₀ _Therapist	0.98
τ₁₁ _{Client.PosEmo}	0.32
τ₁₁ _{Therapist.PosEmo}	0.46
ρ₀₁ _Client	0.07
ρ₀₁ _Therapist	0.38
ICC	0.95
N _Client	1000
N _Therapist	100
Observations	16000

T&B Model

Loading the data from OSF

file <- osf_retrieve_file("6846bf952b204a6ecd5396a6")
#Change to your local dowonloads folder
setwd("C:\\Users\\97254\\Downloads")
osf_download(file, progress = T,conflicts = "overwrigt") 
temp2 <- import("T_B_Data.csv")

names(temp2) <- c("J_PosEmo","C_PosEmo","T_PosEmo","Experience","Depression","Client","Therapist")
temp2 <- temp2 %>% group_by(Therapist,Client) %>% 
  mutate(Session=row_number()) %>% 
  ungroup() %>% relocate(Client:Session,.before = J_PosEmo)

Descriptive stats

head(temp2)

# A tibble: 6 × 8
  Client Therapist Session J_PosEmo C_PosEmo T_PosEmo Experience Depression
   <int>     <int>   <int>    <dbl>    <dbl>    <dbl>      <dbl>      <dbl>
1      1         1       1     14.1     5.91     6.09      -2.61     -0.548
2      1         1       2     12.2     5.19     5.20      -2.61     -0.548
3      1         1       3     10.4     4.59     5.37      -2.61     -0.548
4      1         1       4     14.2     6.68     5.96      -2.61     -0.548
5      1         1       5     13.7     5.94     6.16      -2.61     -0.548
6      1         1       6     13.6     6.33     5.35      -2.61     -0.548

st(temp2)

Summary Statistics
Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 75	Max
Client	16000	500	289	1	251	750	1000
Therapist	16000	50	29	1	26	75	100
Session	16000	8.5	4.6	1	4.8	12	16
J_PosEmo	16000	4.5	4.4	-14	1.5	7.3	20
C_PosEmo	16000	5	1	0.94	4.3	5.7	9
T_PosEmo	16000	6	1	1.8	5.3	6.7	9.8
Experience	16000	0.052	1.1	-2.6	-0.56	0.87	2.8
Depression	16000	0.023	0.97	-2.6	-0.62	0.68	3.2

Centering the Variables on the Truth Variable

temp2 <- temp2 %>% group_by(Therapist,Client) %>% 
  mutate(
    Avg_Truth=mean(C_PosEmo),
    Judg=J_PosEmo-mean(C_PosEmo),
         Truth=C_PosEmo-mean(C_PosEmo),
         Bias=T_PosEmo-mean(C_PosEmo)
  ) %>% 
  ungroup()
st(temp2)

Summary Statistics
Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 75	Max
Client	16000	500	289	1	251	750	1000
Therapist	16000	50	29	1	26	75	100
Session	16000	8.5	4.6	1	4.8	12	16
J_PosEmo	16000	4.5	4.4	-14	1.5	7.3	20
C_PosEmo	16000	5	1	0.94	4.3	5.7	9
T_PosEmo	16000	6	1	1.8	5.3	6.7	9.8
Experience	16000	0.052	1.1	-2.6	-0.56	0.87	2.8
Depression	16000	0.023	0.97	-2.6	-0.62	0.68	3.2
Avg_Truth	16000	5	0.25	4.2	4.8	5.2	5.9
Judg	16000	-0.55	4.4	-19	-3.5	2.3	15
Truth	16000	-0.0000000000000000029	0.97	-4	-0.66	0.65	3.9
Bias	16000	0.99	1	-3.3	0.29	1.7	4.7

head(temp2 %>% select(contains("PosEmo"),Avg_Truth,Judg,Truth,Bias),n=20)

# A tibble: 20 × 7
   J_PosEmo C_PosEmo T_PosEmo Avg_Truth   Judg   Truth   Bias
      <dbl>    <dbl>    <dbl>     <dbl>  <dbl>   <dbl>  <dbl>
 1    14.1      5.91     6.09      4.95  9.15   0.965   1.15 
 2    12.2      5.19     5.20      4.95  7.25   0.250   0.256
 3    10.4      4.59     5.37      4.95  5.49  -0.357   0.423
 4    14.2      6.68     5.96      4.95  9.25   1.74    1.02 
 5    13.7      5.94     6.16      4.95  8.72   0.995   1.21 
 6    13.6      6.33     5.35      4.95  8.69   1.39    0.408
 7    15.9      4.73     7.67      4.95 10.9   -0.210   2.73 
 8    13.8      5.00     5.77      4.95  8.86   0.0559  0.827
 9    15.0      4.88     7.14      4.95 10.1   -0.0615  2.19 
10    16.8      6.59     7.17      4.95 11.9    1.64    2.23 
11    13.9      3.48     5.95      4.95  8.96  -1.46    1.00 
12    12.0      4.53     6.02      4.95  7.03  -0.411   1.07 
13    13.4      4.32     6.23      4.95  8.42  -0.624   1.28 
14    12.3      4.20     5.83      4.95  7.37  -0.746   0.884
15    12.1      3.27     6.17      4.95  7.16  -1.67    1.23 
16    12.2      3.46     5.69      4.95  7.29  -1.48    0.740
17     6.97     3.27     6.05      5.14  1.82  -1.87    0.901
18     5.70     4.22     5.02      5.14  0.560 -0.919  -0.126
19     7.93     6.59     6.53      5.14  2.79   1.45    1.39 
20     6.08     6.10     7.43      5.14  0.934  0.961   2.29

Model 1: T&B

TB1 <- lme(Judg~Truth, 
              random = ~ 1 + Truth| Therapist/Client,
              
              data = temp2)
tab_model(TB1,show.r2=F)

	Judg
Predictors	Estimates	CI	p
(Intercept)	-0.55	-1.13 – 0.03	0.063
Truth	0.32	0.24 – 0.41	<0.001
Random Effects
σ²	1.26
τ₀₀ _Client	2.78
τ₀₀ _Therapist	3.20
τ₁₁ _Client.Truth	0.38
τ₁₁ _{Therapist.Truth}	0.50
ρ₀₁ _Client	0.54
ρ₀₁ _Therapist	0.66
ICC	0.84
N _Client	1000
N _Therapist	100
Observations	16000

Model 2: Moderated T&B, Level-2 Moderator (Client’s Depression)

TB2 <- lme(
  Judg ~ Truth * Depression,
  random = list(
    Therapist = ~ Depression+Truth,
    Client = ~ Truth
  ),

  data = temp2
)
summary(TB2)

Linear mixed-effects model fit by REML
  Data: temp2 
       AIC      BIC    logLik
  55228.91 55336.43 -27600.45

Random effects:
 Formula: ~Depression + Truth | Therapist
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr         
(Intercept) 2.7961278 (Intr) Dprssn
Depression  0.4302474 0.350        
Truth       0.3845697 0.545  0.184 

 Formula: ~Truth | Client %in% Therapist
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 3.0286570 (Intr)
Truth       0.4691316 0.634 
Residual    1.1214130       

Fixed effects:  Judg ~ Truth * Depression 
                      Value  Std.Error    DF   t-value p-value
(Intercept)      -0.5390417 0.29585330 14998 -1.821990  0.0685
Truth             0.3263256 0.04229916 14998  7.714707  0.0000
Depression       -0.9551760 0.11220278   899 -8.512944  0.0000
Truth:Depression -0.1603469 0.01863527 14998 -8.604483  0.0000
 Correlation: 
                 (Intr) Truth  Dprssn
Truth             0.540              
Depression        0.121  0.059       
Truth:Depression -0.004 -0.010  0.486

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-4.61138349 -0.59998947 -0.01230018  0.60297132  4.93291035 

Number of Observations: 16000
Number of Groups: 
            Therapist Client %in% Therapist 
                  100                  1000

# Simple slopes of PosEmo at different levels of Depression
emtrends(TB2, ~ Depression, var = "Truth", at = list(Depression = c(-1, 0, 1)))

 Depression Truth.trend     SE    df lower.CL upper.CL
         -1       0.487 0.0464 14998   0.3957    0.578
          0       0.326 0.0423 14998   0.2434    0.409
          1       0.166 0.0461 14998   0.0757    0.256

Degrees-of-freedom method: containment 
Confidence level used: 0.95

emmip(TB2, Depression ~ Truth,
      at = list(Truth = seq(min(temp2$Truth), max(temp2$Truth), length.out = 20),
                Depression = c(-1, 0, 1)),
      CIs = TRUE)

Model 3: Moderated T&B, Level-3 Moderator (Therapist’s Experience)

TB3 <- lme(
  Judg ~ Truth * Depression+Truth*Experience,
  random = list(
    Therapist = ~ Depression+Truth,
    Client = ~ Truth
  ),
  
  data = temp2
)
summary(TB3)

Linear mixed-effects model fit by REML
  Data: temp2 
       AIC      BIC    logLik
  55215.31 55338.19 -27591.65

Random effects:
 Formula: ~Depression + Truth | Therapist
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr         
(Intercept) 2.6351611 (Intr) Dprssn
Depression  0.4311979 0.361        
Truth       0.3365549 0.457  0.195 

 Formula: ~Truth | Client %in% Therapist
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev   Corr  
(Intercept) 3.028508 (Intr)
Truth       0.469176 0.634 
Residual    1.121415       

Fixed effects:  Judg ~ Truth * Depression + Truth * Experience 
                      Value  Std.Error    DF   t-value p-value
(Intercept)      -0.5839883 0.28098461 14997 -2.078364  0.0377
Truth             0.3175634 0.03802809 14997  8.350760  0.0000
Depression       -0.9539401 0.11216979   899 -8.504429  0.0000
Experience        0.8567644 0.24645606    98  3.476337  0.0008
Truth:Depression -0.1602282 0.01860836 14997 -8.610549  0.0000
Truth:Experience  0.1674730 0.03362128 14997  4.981161  0.0000
 Correlation: 
                 (Intr) Truth  Dprssn Exprnc Trth:D
Truth             0.465                            
Depression        0.124  0.061                     
Experience       -0.046 -0.021  0.004              
Truth:Depression -0.004 -0.011  0.485 -0.001       
Truth:Experience -0.021 -0.047  0.002  0.460  0.000

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-4.6130337 -0.5995688 -0.0118541  0.6012764  4.9300565 

Number of Observations: 16000
Number of Groups: 
            Therapist Client %in% Therapist 
                  100                  1000

# Simple slopes of PosEmo at different levels of Experience
emtrends(TB3, ~ Experience, var = "Truth", at = list(Experience = c(-1, 0, 1)))

 Experience Truth.trend     SE    df lower.CL upper.CL
         -1       0.146 0.0519 14997   0.0446    0.248
          0       0.314 0.0380 14997   0.2393    0.388
          1       0.481 0.0496 14997   0.3842    0.579

Degrees-of-freedom method: containment 
Confidence level used: 0.95

emmip(TB3, Experience ~ Truth,
      at = list(Truth = seq(min(temp2$Truth), max(temp2$Truth), length.out = 20),
                Experience = c(-1, 0, 1)),
      CIs = TRUE)

Model 4: Adding the Bias Force (therapist’s emotions)

TB4 <- lme(Judg~Truth+Bias, 
              random = ~ 1 + Truth+Bias| Therapist/Client,
              
              data = temp2)
tab_model(TB4,show.r2=F)

	Judg
Predictors	Estimates	CI	p
(Intercept)	-1.02	-1.55 – -0.48	<0.001
Truth	0.33	0.24 – 0.41	<0.001
Bias	0.48	0.40 – 0.55	<0.001
Random Effects
σ²	0.72
τ₀₀ _Client	2.56
τ₀₀ _Therapist	2.92
τ₁₁ _Client.Truth	0.37
τ₁₁ _Client.Bias	0.34
τ₁₁ _{Therapist.Truth}	0.58
τ₁₁ _{Therapist.Bias}	0.44
ρ₀₁	0.63
	0.63
	0.73
	0.62
ICC	0.93
N _Client	1000
N _Therapist	100
Observations	16000

Dyadic Data

file <- osf_retrieve_file("68493f15c941fea4dfa4db8c")
#Change to your local dowonloads folder
setwd("C:\\Users\\97254\\Downloads")
osf_download(file, progress = T,conflicts = "overwrigt") 
temp3 <- import("T_B_Dyadic_Data.csv")

names(temp3) <- c("M_J_Vul","W_J_Vul","W_Vul","M_Vul","Couple")
temp3 <- temp3 %>% group_by(Couple) %>%
  mutate(Session=row_number()) %>%
  ungroup() %>% relocate(Couple:Session,.before = M_J_Vul)

Descriptive Stat

head(temp3)

# A tibble: 6 × 6
  Couple Session M_J_Vul W_J_Vul W_Vul M_Vul
   <int>   <int>   <dbl>   <dbl> <dbl> <dbl>
1      1       1   -4.20   17.1   6.84  7.35
2      1       2   -4.32   14.9   6.12  6.65
3      1       3   -3.05   12.3   6.61  4.73
4      1       4   -3.21   11.2   5.09  4.30
5      1       5   -3.71   10.2   5.19  4.30
6      1       6   -2.26    9.24  2.83  3.96

st(temp3)

Summary Statistics
Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 75	Max
Couple	1600	50	29	1	26	75	100
Session	1600	8.5	4.6	1	4.8	12	16
M_J_Vul	1600	3.6	2.6	-5.7	2.1	5.2	12
W_J_Vul	1600	3.7	2.5	-3.3	2.1	4.9	18
W_Vul	1600	5.1	0.98	1.8	4.4	5.7	8.6
M_Vul	1600	5	1	1.1	4.3	5.7	8.1

corr.test(temp3)

Call:corr.test(x = temp3)
Correlation matrix 
        Couple Session M_J_Vul W_J_Vul W_Vul M_Vul
Couple    1.00    0.00   -0.03   -0.01  0.01  0.01
Session   0.00    1.00    0.01    0.00 -0.03 -0.03
M_J_Vul  -0.03    0.01    1.00   -0.58  0.12  0.04
W_J_Vul  -0.01    0.00   -0.58    1.00  0.07  0.15
W_Vul     0.01   -0.03    0.12    0.07  1.00  0.35
M_Vul     0.01   -0.03    0.04    0.15  0.35  1.00
Sample Size 
[1] 1600
Probability values (Entries above the diagonal are adjusted for multiple tests.) 
        Couple Session M_J_Vul W_J_Vul W_Vul M_Vul
Couple    0.00    1.00    1.00       1  1.00     1
Session   1.00    0.00    1.00       1  1.00     1
M_J_Vul   0.24    0.81    0.00       0  0.00     1
W_J_Vul   0.56    0.97    0.00       0  0.04     0
W_Vul     0.65    0.29    0.00       0  0.00     0
M_Vul     0.63    0.24    0.15       0  0.00     0

 To see confidence intervals of the correlations, print with the short=FALSE option

Data stacking

temp3_long <- temp3 %>%
  mutate(M_T_Vul=W_Vul,
         W_T_Vul=M_Vul) %>% 
  pivot_longer(
    cols = c(M_J_Vul, W_J_Vul, M_T_Vul, W_T_Vul),
    names_to = c("Person", ".value"),
    names_pattern = "([MW])_(.*)"
  ) %>% mutate(Female=ifelse(Person=="W",1,0),
               Male=ifelse(Person=="M",1,0)
               ) %>% 
  #centering
  group_by(Couple,Person) %>% 
  mutate(avg_T_Vul=mean(T_Vul),
         J_Vulmc=J_Vul-mean(T_Vul),
         T_Vulmc=T_Vul-mean(T_Vul)
         ) %>% ungroup()
head(temp3_long,30)

# A tibble: 30 × 12
   Couple Session W_Vul M_Vul Person J_Vul T_Vul Female  Male avg_T_Vul J_Vulmc
    <int>   <int> <dbl> <dbl> <chr>  <dbl> <dbl>  <dbl> <dbl>     <dbl>   <dbl>
 1      1       1  6.84  7.35 M      -4.20  6.84      0     1      5.66   -9.86
 2      1       1  6.84  7.35 W      17.1   7.35      1     0      5.28   11.8 
 3      1       2  6.12  6.65 M      -4.32  6.12      0     1      5.66   -9.98
 4      1       2  6.12  6.65 W      14.9   6.65      1     0      5.28    9.66
 5      1       3  6.61  4.73 M      -3.05  6.61      0     1      5.66   -8.71
 6      1       3  6.61  4.73 W      12.3   4.73      1     0      5.28    6.99
 7      1       4  5.09  4.30 M      -3.21  5.09      0     1      5.66   -8.86
 8      1       4  5.09  4.30 W      11.2   4.30      1     0      5.28    5.94
 9      1       5  5.19  4.30 M      -3.71  5.19      0     1      5.66   -9.37
10      1       5  5.19  4.30 W      10.2   4.30      1     0      5.28    4.89
# ℹ 20 more rows
# ℹ 1 more variable: T_Vulmc <dbl>

Model 5: Dyadic T&B

TB4<- lme(J_Vulmc ~ -1 + 
                       Female+ Female:T_Vulmc+
                       Male+ Male:T_Vulmc
                       ,
           random = ~ -1 + 
                      Female+Female:T_Vulmc+
                      Male+Male:T_Vulmc| Couple,
           weights = varIdent(form=~1|Person),
           # corr=corAR1(form = ~1 | Couple/DiaryDay),
           correlation = corCompSymm(form = ~1|Couple/Session),
           data = temp3_long,na.action = na.exclude)
summary(TB4)

Linear mixed-effects model fit by REML
  Data: temp3_long 
       AIC      BIC    logLik
  9058.325 9161.509 -4512.162

Random effects:
 Formula: ~-1 + Female + Female:T_Vulmc + Male + Male:T_Vulmc | Couple
 Structure: General positive-definite, Log-Cholesky parametrization
               StdDev    Corr                
Female         2.2173090 Female Male   Fm:T_V
Male           2.3921426 -0.732              
Female:T_Vulmc 0.5049870  0.960 -0.770       
T_Vulmc:Male   0.4931943 -0.701  0.930 -0.688
Residual       0.8711177                     

Correlation Structure: Compound symmetry
 Formula: ~1 | Couple/Session 
 Parameter estimate(s):
      Rho 
0.2828992 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | Person 
 Parameter estimates:
        M         W 
1.0000000 0.9950838 
Fixed effects:  J_Vulmc ~ -1 + Female + Female:T_Vulmc + Male + Male:T_Vulmc 
                    Value  Std.Error   DF   t-value p-value
Female         -1.3629909 0.22278739 3097 -6.117900       0
Male           -1.4186450 0.24020354 3097 -5.906012       0
Female:T_Vulmc  0.3156361 0.05469998 3097  5.770314       0
T_Vulmc:Male    0.3544787 0.05423213 3097  6.536322       0
 Correlation: 
               Female Male   Fm:T_V
Male           -0.723              
Female:T_Vulmc  0.882 -0.708       
T_Vulmc:Male   -0.634  0.842 -0.562

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-3.16772997 -0.63694388 -0.01781197  0.63908689  3.42718090 

Number of Observations: 3200
Number of Groups: 100

tab_model(TB4,show.r2=F)

	J Vulmc
Predictors	Estimates	CI	p
Female	-1.36	-1.80 – -0.93	<0.001
Male	-1.42	-1.89 – -0.95	<0.001
Female × T Vulmc	0.32	0.21 – 0.42	<0.001
T Vulmc × Male	0.35	0.25 – 0.46	<0.001
Random Effects
σ²	0.76
τ₀₀
τ₀₀
τ₁₁ _Couple.Male	5.72
τ₁₁ _{Couple.Female:T_Vulmc}	0.26
τ₁₁ _{Couple.T_Vulmc:Male}	0.24
ρ₀₁	-0.73
	0.96
	-0.70
ICC	0.88
N _Couple	100
Observations	3200

RSA Model

Loading the data from OSF

file <- osf_retrieve_file("6857a9c8fc6aa8ceafa1c717")
#Change to your local dowonloads folder
setwd("C:\\Users\\97254\\Downloads")
osf_download(file, progress = T,conflicts = "overwrigt") 
temp4 <- import("RSA.csv")

setwd("G:/My Drive/Research/Conferences & Workshops/2025- SPR/Workshop")

Descriptive stats

head(temp4)

    id Alliance       Pos Perceived_Pos n_rows
1 1231     6.44  0.822973      0.822973     37
2 1231     3.38 -0.177027     -0.677027     37
3 1231     6.86  0.822973     -4.177027     37
4 1231     5.74 -4.077027     -1.077027     37
5 1231     2.70  0.822973      0.822973     37
6 1231     6.26  0.822973      0.822973     37

st(temp4)

Summary Statistics
Variable	N	Mean	Std. Dev.	Min	Pctl. 25	Pctl. 75	Max
id	6211	109702	237828	1231	4247	9708	988945
Alliance	6186	5.1	2.1	0	3.7	6.6	10
Pos	6202	0.049	1.6	-9.4	-0.62	0.99	5
Perceived_Pos	6202	-0.049	1.7	-9.2	-0.78	1	5
n_rows	6211	33	11	16	26	37	79

Person-mean centering and computation of RSA variables

temp4 <- temp4 %>% 
  group_by(id) %>% 
  mutate(x=Pos-mean(Pos,na.rm=T),
         y=Perceived_Pos-mean(Perceived_Pos,na.rm=T)
         
         ) %>% 
  ungroup() %>% 
  mutate(x2=x*x,xy=x*y,y2=y*y)

head(temp4,20)

# A tibble: 20 × 10
      id Alliance    Pos Perceived_Pos n_rows      x      y      x2     xy
   <int>    <dbl>  <dbl>         <dbl>  <int>  <dbl>  <dbl>   <dbl>  <dbl>
 1  1231     6.44  0.823         0.823     37  0.697  0.949  0.486   0.661
 2  1231     3.38 -0.177        -0.677     37 -0.303 -0.551  0.0916  0.167
 3  1231     6.86  0.823        -4.18      37  0.697 -4.05   0.486  -2.82 
 4  1231     5.74 -4.08         -1.08      37 -4.20  -0.951 17.7     4.00 
 5  1231     2.7   0.823         0.823     37  0.697  0.949  0.486   0.661
 6  1231     6.26  0.823         0.823     37  0.697  0.949  0.486   0.661
 7  1231     0.1  -5.68         -6.88      37 -5.80  -6.75  33.7    39.2  
 8  1231     3.76  0.823         0.823     37  0.697  0.949  0.486   0.661
 9  1231     5.24  0.823         0.823     37  0.697  0.949  0.486   0.661
10  1231     4.26 -0.977        -1.48      37 -1.10  -1.35   1.22    1.49 
11  1231     3.68  0.823         0.823     37  0.697  0.949  0.486   0.661
12  1231     7     0.823         0.823     37  0.697  0.949  0.486   0.661
13  1231     6.52  0.823         0.723     37  0.697  0.849  0.486   0.592
14  1231     5.74  0.823         0.823     37  0.697  0.949  0.486   0.661
15  1231     5.44  0.623         0.823     37  0.497  0.949  0.247   0.472
16  1231     6.46  0.523         0.823     37  0.397  0.949  0.158   0.377
17  1231     6.1   0.823         0.823     37  0.697  0.949  0.486   0.661
18  1231     3.68 -2.38         -2.78      37 -2.50  -2.65   6.26    6.64 
19  1231     6.48  0.823         0.823     37  0.697  0.949  0.486   0.661
20  1231     7.32  0.823         0.823     37  0.697  0.949  0.486   0.661
# ℹ 1 more variable: y2 <dbl>

Plotting Terciles of difference scores

# Step 1: Create the difference between the two variables
temp <- temp4%>%
  mutate(difference = x  - y )

# Step 2: Calculate quantiles (terciles)
tercile_cutoffs <- quantile(temp$difference, probs = c(1/3, 2/3), na.rm = TRUE)

# Step 3: Assign tercile based on the cutoffs
temp <- temp %>%
  mutate(tercile = case_when(
    difference <= tercile_cutoffs[1] ~ "Lower",
    difference > tercile_cutoffs[1] & difference <= tercile_cutoffs[2] ~ "Medium",
    difference > tercile_cutoffs[2] ~ "High"
  ))

# Step 4: Plot histogram of the difference scores by tercile
ggplot(temp, aes(x = difference, fill = tercile)) +
  geom_histogram(bins = 30, position = "identity", alpha = 0.7) +
  labs(title = "Histogram of Difference Scores by Tercile",
       x = "Difference (Pos vs. Perceived)",
       y = "Count",
       fill = "Tercile") +
  theme_minimal()

Multilevel Response Surface Analysis (MLRSA)

MLRSA<- lme(Alliance ~ 1 + x+y+x2+xy+y2,
               
               random = ~ 1 + x+y| id,
               data = temp4,na.action = na.exclude,
               )
summary(MLRSA)

Linear mixed-effects model fit by REML
  Data: temp4 
      AIC      BIC    logLik
  23427.1 23514.57 -11700.55

Random effects:
 Formula: ~1 + x + y | id
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr         
(Intercept) 1.1826257 (Intr) x     
x           0.2177434  0.054       
y           0.1597083  0.060 -0.646
Residual    1.5011529              

Fixed effects:  Alliance ~ 1 + x + y + x2 + xy + y2 
                Value  Std.Error   DF  t-value p-value
(Intercept)  5.130118 0.08579764 5966 59.79323  0.0000
x            0.397927 0.02687157 5966 14.80847  0.0000
y            0.212022 0.02321558 5966  9.13276  0.0000
x2          -0.012335 0.00753353 5966 -1.63738  0.1016
xy           0.041036 0.00885452 5966  4.63442  0.0000
y2          -0.021120 0.00613937 5966 -3.44012  0.0006
 Correlation: 
   (Intr) x      y      x2     xy    
x   0.006                            
y  -0.010 -0.678                     
x2 -0.089  0.356 -0.230              
xy  0.047 -0.019 -0.031 -0.599       
y2 -0.102 -0.246  0.442 -0.091 -0.487

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-4.64901893 -0.55237973  0.04123639  0.61576559  3.87027444 

Number of Observations: 6179
Number of Groups: 208

tab_model(MLRSA,show.r2=F)

	Alliance
Predictors	Estimates	CI	p
(Intercept)	5.13	4.96 – 5.30	<0.001
x	0.40	0.35 – 0.45	<0.001
y	0.21	0.17 – 0.26	<0.001
x2	-0.01	-0.03 – 0.00	0.102
xy	0.04	0.02 – 0.06	<0.001
y2	-0.02	-0.03 – -0.01	0.001
Random Effects
σ²	2.25
τ₀₀ _id	1.40
τ₁₁ _id.x	0.05
τ₁₁ _id.y	0.03
ρ₀₁	0.05
	0.06
ICC	0.40
N _id	208
Observations	6179

Plotting surface

# Plotting 
(fixed_effects <- fixef(MLRSA))

(Intercept)           x           y          x2          xy          y2 
 5.13011816  0.39792673  0.21202230 -0.01233522  0.04103556 -0.02112013

plotRSA(x=fixed_effects[2],
         y=fixed_effects[3],
         x2=fixed_effects[4],
         xy=fixed_effects[5],
         y2=fixed_effects[6]
         
         )

Getting estimates and their SE using the deltaMethod

# Define a function that receives var_names and builds the expression list
build_and_evaluate_expressions <- function(model, vars) {
  
  # Get the fixed effects and variance-covariance matrix from the model
  fixed_effects <- fixef(model)
  vcov_matrix <- vcov(model)
  
  # Build the expressions dynamically based on the provided variable names
  
  
  expressions_list <- list(
    "a1" = paste(vars[1], "+", vars[2]),  
    "a2" = paste(vars[3], "+", vars[4], "+", vars[5]),  
    "a3" = paste(vars[1], "-", vars[2]),  
    "a4" = paste(vars[3], "-", vars[4], "+", vars[5]),  
    "X0" = paste("(", vars[2], "*", vars[4], "- 2 *", vars[1], "*", vars[5], ") / (4 *", 
                 vars[3], "*", vars[5], "-", vars[4], "^2)"),
    "Y0" = paste("(", vars[1], "*", vars[4], "- 2 *", vars[2], "*", vars[3], ") / (4 *", 
                 vars[3], "*", vars[5], "-", vars[4], "^2)"),
    "p10" = paste("((", vars[1], "*", vars[4], "- 2 *", vars[2], "*", vars[3], ") / (4 *", 
                  vars[3], "*", vars[5], "-", vars[4], "^2)) - ((", vars[5], "-", 
                  vars[3], "+ sqrt((", vars[3], "-", vars[5], ")^2 +", vars[4], "^2))/", 
                  vars[4], ") * ((", vars[2], "*", vars[4], "- 2 *", vars[1], "*", 
                  vars[5], ") / (4 *", vars[3], "*", vars[5], "-", vars[4], "^2))"),
    "p11" = paste("(", vars[5], "-", vars[3], "+ sqrt((", vars[3], "-", vars[5], ")^2 +",
                  vars[4], "^2)) /", vars[4]),
    "p20" = paste("((", vars[1], "*", vars[4], "- 2 *", vars[2], "*", vars[3], ") / (4 *",
                  vars[3], "*", vars[5], "-", vars[4], "^2)) - ((", vars[5], "-",
                  vars[3], "- sqrt((", vars[3], "-", vars[5], ")^2 +", vars[4], "^2))/",
                  vars[4], ") * ((", vars[2], "*", vars[4], "- 2 *", vars[1], "*",
                  vars[5], ") / (4 *", vars[3], "*", vars[5], "-", vars[4], "^2))"),
    "p21" = paste("(", vars[5], "-", vars[3], "- sqrt((", vars[3], "-", vars[5], ")^2 +",
                  vars[4], "^2)) /", vars[4])
  )
  
  # Initialize an empty dataframe to store results
  results_df <- data.frame(Expression = character(), Estimate = numeric(), SE = numeric(),CI_2.5= numeric(),CI_97.5= numeric(), Z = numeric(), p_value = numeric(), stringsAsFactors = FALSE)
  
  # Iterate over the expressions and apply deltaMethod for each
  for (expr_name in names(expressions_list)) {
    expr <- expressions_list[[expr_name]]
    
    # Apply deltaMethod using the current expression
    result <- deltaMethod(object = fixed_effects, 
                          g = expr, 
                          vcov. = vcov_matrix)
    
    # Calculate Z-score and p-value
    z_score <- result$Estimate / result$SE
    
    
    p_value <- round(2 * (1 - pnorm(abs(z_score))),4)
    
    # Append result to the dataframe
    results_df <- results_df %>%
      add_row(Expression = expr_name, Estimate = result$Estimate, SE = result$SE,CI_2.5=result$`2.5 %`,CI_97.5=result$`97.5 %`, Z = z_score, p_value = p_value)
  }
  
  # Return the results dataframe
  return(results_df)
}

# Define the male and female variable names 
vars <- c("x", "y", "x2", "xy", "y2")

# Apply the function to the model (in this case Model1.1)
results <- build_and_evaluate_expressions(MLRSA, vars)

# View the results
print(results)

   Expression      Estimate           SE        CI_2.5     CI_97.5         Z
1          a1   0.609949035  0.020385818  5.699936e-01  0.64990450 29.920263
2          a2   0.007580209  0.005627002 -3.448512e-03  0.01860893  1.347113
3          a3   0.185904428  0.045896704  9.594854e-02  0.27586032  4.050496
4          a4  -0.074490917  0.017241104 -1.082829e-01 -0.04069897 -4.320542
5          X0 -39.744061765 27.177601473 -9.301118e+01 13.52305831 -1.462383
6          Y0 -33.591116632 25.437462253 -8.344763e+01 16.26539324 -1.320537
7         p10  -1.454938813  0.927085202 -3.271992e+00  0.36211479 -1.569369
8         p11   0.808578097  0.204101582  4.085463e-01  1.20860985  3.961645
9         p20 -82.744144552 54.915981174 -1.903775e+02 24.88920073 -1.506741
10        p21  -1.236738917  0.312178095 -1.848597e+00 -0.62488109 -3.961645
   p_value
1   0.0000
2   0.1779
3   0.0001
4   0.0000
5   0.1436
6   0.1867
7   0.1166
8   0.0001
9   0.1319
10  0.0001