R Notebook

Predicting Math Anxiety

Outcome (Numerical): (Math anxiety) score_AMAS_total (Outcome Variable)

Dichotomous Predictor: sex (Covariate)

Multi-categorical Predictor: sample (Moderator)

Predictor 1: (General trait anxiety) score_GAD (Mediator)

Predictor 2: (Math self-concept) score_SDQ_M (Focal Predictor)

Predictor 3: (Neuroticism) score_BFI_N (Covariate)

# install.packages("readr")
# install.packages("psych")
# install.packages("dplyr")
# install.packages("tidyr")
# install.packages("car")
# install.packages("lm.beta")
# install.packages("interactions")
# install.packages("emmeans")
# install.packages("ggeffects")
# install.packages("mediation")

# load the following packages
library(readr)
library(psych)
library(dplyr)
library(tidyr)
library(car)
library(lm.beta)
#library(emmeans)
#library(ggeffects)
library(interactions)
library(mediation)

Descripitive Statisitics

# turn off scientific notation 
options(scipen = 999)

# examine the dataset
AMATUS_dataset_processed

# Data frame

Anxiety_data <- AMATUS_dataset_processed %>% 
  dplyr::select(
  score_AMAS_total,
  sex,
  sample,
  score_GAD,
  score_SDQ_M,
  score_BFI_N
)  %>%
  
  drop_na("score_AMAS_total",
                            "score_GAD",
                            "score_SDQ_M",
                            "score_BFI_N",
                            "sex",
                            "sample") %>%
  
  dplyr::mutate(
    sex_d = dplyr::recode(sex, 'f' = 0,
                        'm' =1),
    ger_stu = if_else(sample == "german_students", 1, 0),
    ger_tea = if_else(sample == "german_teachers", 1, 0),
    bel_tea = if_else(sample == "belgian_teachers", 1, 0),
  )
Anxiety_data$sample <- relevel(as.factor(Anxiety_data$sample), ref = "german_students")
Anxiety_data

describe(Anxiety_data[, c("score_AMAS_total",
                            "score_GAD",
                            "score_SDQ_M",
                            "score_BFI_N",
                            "sex_d",
                            "sample")])

describe(Anxiety_data)

table(Anxiety_data$sex)

  f   m 
820 286 

table(Anxiety_data$sample)

 german_students belgian_teachers  german_teachers 
             848              127              131 

#Bivariate cor. between num. variables.
cor(Anxiety_data$score_AMAS_total, Anxiety_data$score_GAD)

[1] 0.349942

cor(Anxiety_data$score_AMAS_total, Anxiety_data$score_SDQ_M)

[1] -0.6429761

cor(Anxiety_data$score_AMAS_total, Anxiety_data$score_BFI_N)

[1] 0.3432629

prop.table(table(Anxiety_data$sample))*100

 german_students belgian_teachers  german_teachers 
        76.67269         11.48282         11.84448 

prop.table(table(Anxiety_data$sex))*100

       f        m 
74.14105 25.85895 

Anxiety_data %>% 
  group_by(sex) %>%
  summarise(mean_anxiety = mean(score_AMAS_total))

Anxiety_data %>% 
  group_by(sample) %>%
  summarise(mean_anxiety = mean(score_AMAS_total))

Multiple Linear Regressions models

# Multi-categorical predictor model, German student as reference group.

Anxiety_lm_GAD_sex_sample <- lm(score_AMAS_total ~  score_GAD + 
                                      score_SDQ_M + 
                                      score_BFI_N +
                                      sex_d + 
                                      ger_tea + 
                                      bel_tea, 
                 data = Anxiety_data)

summary(Anxiety_lm_GAD_sex_sample)

Call:
lm(formula = score_AMAS_total ~ score_GAD + score_SDQ_M + score_BFI_N + 
    sex_d + ger_tea + bel_tea, data = Anxiety_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.2076  -3.4000  -0.3279   3.0469  23.3077 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept) 27.23036    0.93261  29.198 < 0.0000000000000002 ***
score_GAD    0.24500    0.04541   5.396         0.0000000836 ***
score_SDQ_M -1.26723    0.04843 -26.168 < 0.0000000000000002 ***
score_BFI_N  0.14165    0.03144   4.505         0.0000073394 ***
sex_d       -1.69373    0.35113  -4.824         0.0000016079 ***
ger_tea      1.28745    0.46584   2.764              0.00581 ** 
bel_tea      0.46216    0.47624   0.970              0.33204    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.903 on 1099 degrees of freedom
Multiple R-squared:  0.4943,    Adjusted R-squared:  0.4916 
F-statistic: 179.1 on 6 and 1099 DF,  p-value: < 0.00000000000000022

anova(Anxiety_lm_GAD_sex_sample)

Analysis of Variance Table

Response: score_AMAS_total
              Df  Sum Sq Mean Sq  F value                Pr(>F)    
score_GAD      1  6399.3  6399.3 266.1468 < 0.00000000000000022 ***
score_SDQ_M    1 17865.7 17865.7 743.0358 < 0.00000000000000022 ***
score_BFI_N    1   683.7   683.7  28.4362          0.0000001175 ***
sex_d          1   691.7   691.7  28.7686          0.0000000994 ***
ger_tea        1   168.7   168.7   7.0182              0.008184 ** 
bel_tea        1    22.6    22.6   0.9417              0.332044    
Residuals   1099 26424.6    24.0                                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

lm.beta(Anxiety_lm_GAD_sex_sample)

Call:
lm(formula = score_AMAS_total ~ score_GAD + score_SDQ_M + score_BFI_N + 
    sex_d + ger_tea + bel_tea, data = Anxiety_data)

Standardized Coefficients::
(Intercept)   score_GAD score_SDQ_M score_BFI_N       sex_d     ger_tea     bel_tea 
         NA  0.14747323 -0.57697605  0.12464582 -0.10789144  0.06052326  0.02143577 

# Model without focal predictor

Anxiety_lm_sex_sample <- lm(score_AMAS_total ~ score_GAD + 
                              score_BFI_N+ 
                              sex_d + 
                              ger_tea + 
                              bel_tea, 
                 data = Anxiety_data)

summary(Anxiety_lm_sex_sample)

Call:
lm(formula = score_AMAS_total ~ score_GAD + score_BFI_N + sex_d + 
    ger_tea + bel_tea, data = Anxiety_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.4700  -4.6780  -0.7561   4.1943  24.3665 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)  9.90479    0.83640  11.842 < 0.0000000000000002 ***
score_GAD    0.37458    0.05748   6.517       0.000000000109 ***
score_BFI_N  0.20359    0.03992   5.099       0.000000401321 ***
sex_d       -2.37219    0.44592  -5.320       0.000000125792 ***
ger_tea      1.32295    0.59321   2.230               0.0259 *  
bel_tea      0.52973    0.60645   0.874               0.3826    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.244 on 1100 degrees of freedom
Multiple R-squared:  0.1793,    Adjusted R-squared:  0.1755 
F-statistic: 48.05 on 5 and 1100 DF,  p-value: < 0.00000000000000022

summary(Anxiety_lm_GAD_sex_sample)$r.squared - summary(Anxiety_lm_sex_sample)$r.squared 

[1] 0.3150722

anova(Anxiety_lm_GAD_sex_sample, Anxiety_lm_sex_sample)

Analysis of Variance Table

Model 1: score_AMAS_total ~ score_GAD + score_SDQ_M + score_BFI_N + sex_d + 
    ger_tea + bel_tea
Model 2: score_AMAS_total ~ score_GAD + score_BFI_N + sex_d + ger_tea + 
    bel_tea
  Res.Df   RSS Df Sum of Sq      F                Pr(>F)    
1   1099 26425                                              
2   1100 42889 -1    -16465 684.76 < 0.00000000000000022 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Nonlinear

# Quadratic Model with squared score
# Math self-concept

Anxiety_lm_GAD_sex_sample_quad <- lm(score_AMAS_total ~ score_SDQ_M + I(score_SDQ_M^2) +
                                      score_GAD + 
                                      score_BFI_N + 
                                      sex_d + 
                                      ger_tea + 
                                      bel_tea, 
                 data = Anxiety_data)



summary(Anxiety_lm_GAD_sex_sample_quad)

Call:
lm(formula = score_AMAS_total ~ score_SDQ_M + I(score_SDQ_M^2) + 
    score_GAD + score_BFI_N + sex_d + ger_tea + bel_tea, data = Anxiety_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.8534  -3.4056  -0.3035   3.0707  23.1467 

Coefficients:
                 Estimate Std. Error t value             Pr(>|t|)    
(Intercept)      25.51715    1.80617  14.128 < 0.0000000000000002 ***
score_SDQ_M      -0.92319    0.31438  -2.937              0.00339 ** 
I(score_SDQ_M^2) -0.01585    0.01431  -1.108              0.26829    
score_GAD         0.24765    0.04547   5.447         0.0000000632 ***
score_BFI_N       0.14071    0.03145   4.474         0.0000084766 ***
sex_d            -1.70638    0.35128  -4.858         0.0000013604 ***
ger_tea           1.26188    0.46636   2.706              0.00692 ** 
bel_tea           0.45788    0.47621   0.962              0.33651    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.903 on 1098 degrees of freedom
Multiple R-squared:  0.4949,    Adjusted R-squared:  0.4917 
F-statistic: 153.7 on 7 and 1098 DF,  p-value: < 0.00000000000000022

summary(Anxiety_lm_GAD_sex_sample_quad)$r.squared - summary(Anxiety_lm_GAD_sex_sample)$r.squared

[1] 0.0005643132

anova(Anxiety_lm_GAD_sex_sample ,Anxiety_lm_GAD_sex_sample_quad)

Analysis of Variance Table

Model 1: score_AMAS_total ~ score_GAD + score_SDQ_M + score_BFI_N + sex_d + 
    ger_tea + bel_tea
Model 2: score_AMAS_total ~ score_SDQ_M + I(score_SDQ_M^2) + score_GAD + 
    score_BFI_N + sex_d + ger_tea + bel_tea
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1   1099 26425                           
2   1098 26395  1    29.489 1.2267 0.2683

Moderation

#
#

Anxiety_lm_GAD_sex_sample_mod <- lm(score_AMAS_total ~ score_SDQ_M +
                                      score_GAD + 
                                      score_BFI_N + 
                                      sex_d + 
                                      sample +
                                      score_SDQ_M*sample, 
                 data = Anxiety_data)

summary(Anxiety_lm_GAD_sex_sample_mod)

Call:
lm(formula = score_AMAS_total ~ score_SDQ_M + score_GAD + score_BFI_N + 
    sex_d + sample + score_SDQ_M * sample, data = Anxiety_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.2146  -3.3878  -0.3429   3.1310  23.3755 

Coefficients:
                                   Estimate Std. Error t value             Pr(>|t|)    
(Intercept)                        26.67077    0.97209  27.436 < 0.0000000000000002 ***
score_SDQ_M                        -1.21727    0.05455 -22.316 < 0.0000000000000002 ***
score_GAD                           0.24611    0.04539   5.422         0.0000000726 ***
score_BFI_N                         0.14078    0.03141   4.482         0.0000081837 ***
sex_d                              -1.70208    0.35087  -4.851         0.0000014059 ***
samplebelgian_teachers              3.57646    1.69805   2.106               0.0354 *  
samplegerman_teachers               3.08981    1.82727   1.691               0.0911 .  
score_SDQ_M:samplebelgian_teachers -0.27873    0.14605  -1.908               0.0566 .  
score_SDQ_M:samplegerman_teachers  -0.15911    0.15616  -1.019               0.3085    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.898 on 1097 degrees of freedom
Multiple R-squared:  0.4963,    Adjusted R-squared:  0.4926 
F-statistic: 135.1 on 8 and 1097 DF,  p-value: < 0.00000000000000022

summary(Anxiety_lm_GAD_sex_sample_mod)$r.squared - summary(Anxiety_lm_GAD_sex_sample)$r.squared

[1] 0.001957103

anova(Anxiety_lm_GAD_sex_sample ,Anxiety_lm_GAD_sex_sample_mod)

Analysis of Variance Table

Model 1: score_AMAS_total ~ score_GAD + score_SDQ_M + score_BFI_N + sex_d + 
    ger_tea + bel_tea
Model 2: score_AMAS_total ~ score_SDQ_M + score_GAD + score_BFI_N + sex_d + 
    sample + score_SDQ_M * sample
  Res.Df   RSS Df Sum of Sq      F Pr(>F)
1   1099 26425                           
2   1097 26322  2    102.27 2.1311 0.1192

sim_slopes(Anxiety_lm_GAD_sex_sample_mod, 
           pred = score_SDQ_M, 
           modx = sample)

SIMPLE SLOPES ANALYSIS

Slope of score_SDQ_M when sample = german_students: 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  -1.22   0.05   -22.32   0.00

Slope of score_SDQ_M when sample = belgian_teachers: 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  -1.50   0.14   -10.98   0.00

Slope of score_SDQ_M when sample = german_teachers: 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  -1.38   0.15    -9.34   0.00

interact_plot(Anxiety_lm_GAD_sex_sample_mod, 
           pred = score_SDQ_M, 
           modx = sample,
           interval = TRUE)

Mediation

# Outcome Variable:score_AMAS_total (Math anxiety) 
# Mediator: score_GAD (General trait anxiety)
# Focal Predictor: score_SDQ_M (Math self-concept)
# Moderator: sample
# Covariate: score_BFI_N (Neuroticism) 
# Covariate: sex

Anxiety_lm_TE <- lm(score_AMAS_total ~ score_SDQ_M + 
                                      score_BFI_N +
                                      sex_d + 
                                      ger_tea + 
                                      bel_tea, 
                 data = Anxiety_data)

Anxiety_lm_TE

Call:
lm(formula = score_AMAS_total ~ score_SDQ_M + score_BFI_N + sex_d + 
    ger_tea + bel_tea, data = Anxiety_data)

Coefficients:
(Intercept)  score_SDQ_M  score_BFI_N        sex_d      ger_tea      bel_tea  
    28.2663      -1.2957       0.2414      -1.5992       1.3309       0.8816  

summary(Anxiety_lm_TE)

Call:
lm(formula = score_AMAS_total ~ score_SDQ_M + score_BFI_N + sex_d + 
    ger_tea + bel_tea, data = Anxiety_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.5668  -3.4082  -0.3589   3.2188  25.0138 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept) 28.26634    0.92422  30.584 < 0.0000000000000002 ***
score_SDQ_M -1.29572    0.04875 -26.579 < 0.0000000000000002 ***
score_BFI_N  0.24140    0.02576   9.373 < 0.0000000000000002 ***
sex_d       -1.59925    0.35515  -4.503           0.00000741 ***
ger_tea      1.33094    0.47169   2.822              0.00486 ** 
bel_tea      0.88159    0.47582   1.853              0.06418 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.966 on 1100 degrees of freedom
Multiple R-squared:  0.4809,    Adjusted R-squared:  0.4786 
F-statistic: 203.8 on 5 and 1100 DF,  p-value: < 0.00000000000000022

Anxiety_lm_DE <- lm(score_AMAS_total ~ score_GAD+ 
                                      score_SDQ_M + 
                                      score_BFI_N +
                                      sex_d + 
                                      ger_tea + 
                                      bel_tea, 
                 data = Anxiety_data)

Anxiety_lm_DE

Call:
lm(formula = score_AMAS_total ~ score_GAD + score_SDQ_M + score_BFI_N + 
    sex_d + ger_tea + bel_tea, data = Anxiety_data)

Coefficients:
(Intercept)    score_GAD  score_SDQ_M  score_BFI_N        sex_d      ger_tea      bel_tea  
    27.2304       0.2450      -1.2672       0.1417      -1.6937       1.2875       0.4622  

summary(Anxiety_lm_DE)

Call:
lm(formula = score_AMAS_total ~ score_GAD + score_SDQ_M + score_BFI_N + 
    sex_d + ger_tea + bel_tea, data = Anxiety_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-17.2076  -3.4000  -0.3279   3.0469  23.3077 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept) 27.23036    0.93261  29.198 < 0.0000000000000002 ***
score_GAD    0.24500    0.04541   5.396         0.0000000836 ***
score_SDQ_M -1.26723    0.04843 -26.168 < 0.0000000000000002 ***
score_BFI_N  0.14165    0.03144   4.505         0.0000073394 ***
sex_d       -1.69373    0.35113  -4.824         0.0000016079 ***
ger_tea      1.28745    0.46584   2.764              0.00581 ** 
bel_tea      0.46216    0.47624   0.970              0.33204    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.903 on 1099 degrees of freedom
Multiple R-squared:  0.4943,    Adjusted R-squared:  0.4916 
F-statistic: 179.1 on 6 and 1099 DF,  p-value: < 0.00000000000000022

GAD_lm_SDQ <- lm(score_GAD ~ +  score_SDQ_M + 
                                      score_BFI_N +
                                      sex_d + 
                                      ger_tea + 
                                      bel_tea, 
                 data = Anxiety_data)

GAD_lm_SDQ

Call:
lm(formula = score_GAD ~ +score_SDQ_M + score_BFI_N + sex_d + 
    ger_tea + bel_tea, data = Anxiety_data)

Coefficients:
(Intercept)  score_SDQ_M  score_BFI_N        sex_d      ger_tea      bel_tea  
     4.2285      -0.1163       0.4071       0.3856       0.1775       1.7120  

summary(GAD_lm_SDQ)

Call:
lm(formula = score_GAD ~ +score_SDQ_M + score_BFI_N + sex_d + 
    ger_tea + bel_tea, data = Anxiety_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.4540 -2.1361 -0.4019  1.6516 13.0460 

Coefficients:
            Estimate Std. Error t value             Pr(>|t|)    
(Intercept)  4.22848    0.60601   6.978     0.00000000000518 ***
score_SDQ_M -0.11631    0.03196  -3.639             0.000287 ***
score_BFI_N  0.40713    0.01689  24.108 < 0.0000000000000002 ***
sex_d        0.38565    0.23287   1.656             0.097995 .  
ger_tea      0.17748    0.30928   0.574             0.566200    
bel_tea      1.71196    0.31200   5.487     0.00000005068829 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.256 on 1100 degrees of freedom
Multiple R-squared:  0.3841,    Adjusted R-squared:  0.3813 
F-statistic: 137.2 on 5 and 1100 DF,  p-value: < 0.00000000000000022

Anxiety_mediation <- mediate(model.m = GAD_lm_SDQ,
                          model.y = Anxiety_lm_DE,
                          treat = "score_SDQ_M",
                          mediator = "score_GAD",
                          boot = TRUE,
                          sims = 1000)

summary(Anxiety_mediation) 

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

                 Estimate 95% CI Lower 95% CI Upper               p-value    
ACME           -0.0284949   -0.0528058   -0.0103261 < 0.00000000000000022 ***
ADE            -1.2672256   -1.3690575   -1.1687130 < 0.00000000000000022 ***
Total Effect   -1.2957205   -1.3957466   -1.1960217 < 0.00000000000000022 ***
Prop. Mediated  0.0219915    0.0078397    0.0412136 < 0.00000000000000022 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 1106 


Simulations: 1000 

Regression Diagnostics

Anxiety_lm_diag <-Anxiety_data %>%
  mutate(
    Anxiety_pred = fitted(Anxiety_lm_GAD_sex_sample),
    res = residuals(Anxiety_lm_GAD_sex_sample),
    h = hatvalues(Anxiety_lm_GAD_sex_sample),
    tres = rstudent(Anxiety_lm_GAD_sex_sample),
    Cook = cooks.distance(Anxiety_lm_GAD_sex_sample),
    DB_b0 = influence(Anxiety_lm_GAD_sex_sample)$coef[, 1],
    DB_GAD = influence(Anxiety_lm_GAD_sex_sample)$coef[, 2],
    DB_SDQ = influence(Anxiety_lm_GAD_sex_sample)$coef[, 3],
    DB_BFI = influence(Anxiety_lm_GAD_sex_sample)$coef[, 4],
    DB_sex = influence(Anxiety_lm_GAD_sex_sample)$coef[, 5],
    DB_ger_tea = influence(Anxiety_lm_GAD_sex_sample)$coef[, 6],
    DB_bel_tea = influence(Anxiety_lm_GAD_sex_sample)$coef[, 7]
  )


Anxiety_lm_diag

residualPlots(Anxiety_lm_GAD_sex_sample)

            Test stat Pr(>|Test stat|)
score_GAD     -0.7358           0.4620
score_SDQ_M   -1.1076           0.2683
score_BFI_N   -0.8373           0.4026
sex_d         -0.9117           0.3621
ger_tea        0.3200           0.7490
bel_tea       -1.0502           0.2939
Tukey test     0.3533           0.7239

qqPlot(Anxiety_lm_GAD_sex_sample)

[1]  78 736

outlierTest(Anxiety_lm_GAD_sex_sample)