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

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)

belgian_teachers  german_students  german_teachers 
             127              848              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

belgian_teachers  german_students  german_teachers 
        11.48282         76.67269         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)                       30.24723    1.73205  17.463 < 0.0000000000000002 ***
score_SDQ_M                       -1.49600    0.13623 -10.982 < 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 ***
samplegerman_students             -3.57646    1.69805  -2.106               0.0354 *  
samplegerman_teachers             -0.48665    2.32821  -0.209               0.8345    
score_SDQ_M:samplegerman_students  0.27873    0.14605   1.908               0.0566 .  
score_SDQ_M:samplegerman_teachers  0.11962    0.20010   0.598               0.5501    
---
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 = german_teachers: 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  -1.38   0.15    -9.34   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

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  

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  

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  

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.0511897   -0.0113704 < 0.00000000000000022 ***
ADE            -1.2672256   -1.3663976   -1.1700950 < 0.00000000000000022 ***
Total Effect   -1.2957205   -1.3995183   -1.1936573 < 0.00000000000000022 ***
Prop. Mediated  0.0219915    0.0089175    0.0390627 < 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 <-AMATUS_dataset_processed %>%
  mutate(
     prey_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_age = influence(Anxiety_lm_GAD_sex_sample)$coef[, 2],
    DB_sex = influence(Anxiety_lm_GAD_sex_sample)$coef[, 3],
    DB_UK = influence(Anxiety_lm_GAD_sex_sample)$coef[, 4],
    DB_US = influence(Anxiety_lm_GAD_sex_sample)$coef[, 5],
    DB_NZ = influence(Anxiety_lm_GAD_sex_sample)$coef[, 6]
  )

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)