BJCP: Known-groups and convergent validities of the ToM-TB in children with ASD

---

S1: Deviations from protocol

1. In the manuscript:

• The title of the paper has been rephrased to describe more precisely the type of validity explored.
• To ensure consistency throughout the manuscript, the term ‘convergent validity’ was preferred over ‘concurrent validity’.

2. In the methods:

• The SCQ was not used as an inclusion criteria for the ASD sample. This deviation was carried out due to the COVID-19 pandemic as we needed to keep the time passed in the hospital by the participants as short as possible. Moreover, the scores of the SCQ were judged not particularly relevant for our study (for the ASD sample) given that each child with ASD had received a clinical diagnosis of ASD made by a multidisciplinary team and confirmed by standard instruments (ADI-R and/or ADOS-2).

3. In the statistical analyses

• It was initially planned to explore the association of the socialization subscale of the Vineland-II with the two ToM tests by combining the two groups (ASD vs. CONT) if the group did not interact with the Vineland scores. However, we judged that - given the large variations in ToM tests and in Vineland scores across groups - it would be more meaningful to present the results separately for each group.
• Various analyses were added

---

dat = read.delim("BJCP_dat.txt")
dat_long = dat %>%
  tidyr::pivot_longer(cols = c(TOMTB, SST), names_to = "ToM_test", values_to = "ToM")
dat_long$ToM_test = factor(dat_long$ToM_test, levels= c("ToM-TB" = "TOMTB", "SST"))
levels(dat_long$ToM_test)[1] <- "ToM-TB"

S2: Preliminary analyses

1. Table1

table1::table1(~ age + sex + 
                 IQ + 
                 Vineland + 
                 BRIEF +
                 EVIP | group,
                 render.missing = NULL, 
                 data = df.table, overall = "FALSE", 
               extra.col=list(`P-value` = pvalue,
                              `SMD/V` = es))

	ASD (N=34)	CONT (N=34)	P-value	SMD/V
Age
Mean (SD)	8.44 (1.28)	8.54 (1.29)	0.731	0.0838
Median [Min, Max]	8.54 [6.08, 10.8]	8.42 [6.25, 10.9]
Sex
Boys	23 (67.6%)	19 (55.9%)	0.318	0.0908
Girls	11 (32.4%)	15 (44.1%)
Full-scale IQ
Mean (SD)	112 (15.8)	111 (14.4)	0.792	0.0642
Median [Min, Max]	113 [82.0, 147]	113 [81.0, 137]
Vineland (socialization)
Mean (SD)	58.9 (24.3)	110 (18.9)	<0.001	2.37
Median [Min, Max]	56.5 [20.0, 119]	108 [50.0, 140]
BRIEF (total score)
Mean (SD)	69.3 (15.1)	53.4 (10.1)	<0.001	1.24
Median [Min, Max]	70.0 [44.0, 98.0]	50.0 [38.0, 76.0]
EVIP (total score)
Mean (SD)	119 (16.3)	118 (12.9)	0.902	0.0299
Median [Min, Max]	119 [86.0, 149]	121 [90.0, 144]

2. Internal reliability

KR20 for ToM-TB

DescTools::CronbachAlpha(
  dat[, c("TOMTB_A_1", "TOMTB_A_2", "TOMTB_A_3", "TOMTB_A_4", 
         "TOMTB_B_5", "TOMTB_C_6", "TOMTB_D_7", "TOMTB_D_8",
         "TOMTB_E_9", "TOMTB_F_10", "TOMTB_G_11", "TOMTB_G_12",
         "TOMTB_G_13", "TOMTB_H_14", "TOMTB_I_15")])

## [1] 0.6557157

Cronbach’s alpha for SST

DescTools::CronbachAlpha(
  dat[,c("sst1", "sst2", "sst3", "sst4", 
         "sst5", "sst6", "sst7", "sst8")])

## [1] 0.6576006

3. Correlations

ASD

ggstatsplot::ggcorrmat(dat[,
     c('IQ', 'BRIEF', 'EVIP')],
    p.adjust.method = "none",
    ggtheme = ggplot2::theme_bw(),
    colors = c("#0016FF", "white", "#FF0000"),
    pch = 13)

---

S3: Known-groups

1. Main analysis

A. Plot

ggplot(dat_long, aes(x = group, y = ToM)) + 
  geom_jitter(binaxis='y',
            stackdir = 'center', size = 1.5, alpha = 0.3,
            position = position_jitter(width = 0.05, height = 0.1)) +
  geom_boxplot(width = 0.5, alpha = 0)+
  facet_wrap(~ToM_test) +
  theme_bw() +
  labs(y = "Performance at ToM tests",
       x = "") +
  theme(text = element_text(size = 14),
        legend.position = "none",
        axis.title.y = element_text(size=14)) +
  scale_y_continuous(limits = c(0, 16), breaks = seq(0, 16, 4)) 

---

B. Multivariate effect

Known-groups validity for the two ToM tests combined

KG_mult <- lm(cbind(TOMTB, SST) ~ group + IQ + EVIP + BRIEF, data = dat) 
car::Manova(KG_mult, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df    Pr(>F)    
## group  1  0.137156   4.9277      2     62 0.0103251 *  
## IQ     1  0.217124   8.5976      2     62 0.0005064 ***
## EVIP   1  0.002341   0.0727      2     62 0.9299162    
## BRIEF  1  0.011480   0.3600      2     62 0.6991065    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

---

C. Univariate effects

Known-groups validity for ToM-TB

KG_TOMTB <- lm(TOMTB ~ group + IQ + EVIP + BRIEF, data = dat)
broom::tidy(KG_TOMTB)

## # A tibble: 5 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  3.08       2.34      1.32   0.192   
## 2 groupCONT    1.74       0.556     3.12   0.00270 
## 3 IQ           0.0693     0.0178    3.90   0.000233
## 4 EVIP        -0.00105    0.0188   -0.0558 0.956   
## 5 BRIEF        0.00876    0.0191    0.459  0.648

emm.TOMTB <- emmeans(KG_TOMTB, "group")
smd(x = emm.TOMTB, level = 97.5, df = dat, test = "TOMTB")

## [1] "SMD = 0.81 [0.24, 1.37]"

Known-groups validity for SST

KG_SST <- lm(SST ~ group + IQ + EVIP + BRIEF, data = dat)
broom::tidy(KG_SST)

## # A tibble: 5 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept) -0.663      3.72      -0.178 0.859  
## 2 groupCONT    1.50       0.885      1.69  0.0961 
## 3 IQ           0.0821     0.0283     2.90  0.00512
## 4 EVIP         0.00984    0.0300     0.328 0.744  
## 5 BRIEF       -0.0148     0.0304    -0.485 0.629

emm.SST <- emmeans(KG_SST, "group")
smd(x = emm.SST, level = 97.5, df = dat, test = "SST")

## [1] "SMD = 0.46 [-0.09, 1.01]"

---

2. Comparison

Magnitude of group differences between two ToM tests

emm.KG <- emmeans(KG_mult, ~ group | ToMTest, mult.name = "ToMTest")
contrast(emm.KG, interaction = "consec", by = NULL)

##  group_consec ToMTest_consec estimate    SE df t.ratio p.value
##  CONT - ASD   SST - TOMTB       -0.24 0.842 63  -0.285  0.7765

---

3. Replication (unadjusted model)

Known-groups validity for the two ToM tests combined

KG_mult_crude <- lm(cbind(TOMTB, SST) ~ group, data = dat) 
car::Manova(KG_mult_crude, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df Pr(>F)  
## group  1   0.12104   4.4757      2     65 0.0151 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Known-groups validity for ToM-TB

KG_TOMTB_crude <- lm(TOMTB ~ group, data = dat)
broom::tidy(KG_TOMTB_crude)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)    11.4      0.369     30.8  7.23e-41
## 2 groupCONT       1.53     0.522      2.93 4.62e- 3

emm.TOMTB_crude <- emmeans(KG_TOMTB_crude, "group")
smd(x = emm.TOMTB_crude, level = 97.5, df = dat, test = "TOMTB")

## [1] "SMD = 0.71 [0.15, 1.27]"

Known-groups validity for SST

KG_SST_crude <- lm(SST ~ group, data = dat)
broom::tidy(KG_SST_crude)

## # A tibble: 2 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)     8.71     0.560     15.6  9.21e-24
## 2 groupCONT       1.65     0.791      2.08 4.13e- 2

emm.SST_crude <- emmeans(KG_SST_crude, "group")
smd(x = emm.SST_crude, level = 97.5, df = dat, test = "SST")

## [1] "SMD = 0.5 [-0.05, 1.06]"

---

4. Replication (BES)

Backward Elimination Selection of covariates

KG_BWD <- lm(cbind(TOMTB, SST) ~ group + IQ, data = dat) 
car::Manova(KG_BWD, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df    Pr(>F)    
## group  1   0.17000    6.554      2     64  0.002574 ** 
## IQ     1   0.26851   11.746      2     64 4.516e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

KG_TOMTB_BWD <- lm(TOMTB ~ group + IQ, data = dat)
broom::tidy(KG_TOMTB_BWD)

## # A tibble: 3 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)   3.57      1.76        2.02 0.0470   
## 2 groupCONT     1.60      0.459       3.48 0.000911 
## 3 IQ            0.0693    0.0154      4.50 0.0000291

emm.TOMTB.BWD <- emmeans(KG_TOMTB_BWD, "group")
smd(x = emm.TOMTB.BWD, level = 97.5, df = dat, test = "TOMTB")

## [1] "SMD = 0.74 [0.18, 1.3]"

KG_SST_BWD <- lm(SST ~ group + IQ, data = dat)
broom::tidy(KG_SST_BWD)

## # A tibble: 3 x 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)  -0.940     2.81      -0.335 0.739   
## 2 groupCONT     1.73      0.732      2.36  0.0211  
## 3 IQ            0.0858    0.0246     3.49  0.000860

emm.SST.BWD <- emmeans(KG_SST_BWD, "group")
smd(x = emm.SST.BWD, level = 97.5, df = dat, test = "SST")

## [1] "SMD = 0.53 [-0.02, 1.08]"

5. New multivariate method

library(MVLM)

## Warning: package 'MVLM' was built under R version 4.1.3

dat$group=factor(dat$group)

mvlm.res <- mvlm(as.matrix(dat$TOMTB, dat$SST) ~ group + IQ + EVIP + BRIEF, data = dat)

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in CompQuadForm::davies(0, lambda = gamma, h = nu, lim = lim, acc =
## acc): Consider playing with 'lim' or 'acc'.

## Warning in mvlm(as.matrix(dat$TOMTB, dat$SST) ~ group + IQ + EVIP + BRIEF, : Adjusted sample size = 63
## Asymptotic properties of the null distribution may not hold.
## This can result in overly conservative p-values.
## Increased sample size is recommended.

summary(mvlm.res)

##                 Statistic Numer.DF Pseudo.R.Square    p.value    
## Omnibus Effect   7.668810        4       3.275e-01 4.2335e-05 ***
## (Intercept)      3.034039        1                   0.086415 .  
## group            9.754514        1       1.041e-01  0.0027021 ** 
## IQ              15.237902        1       1.627e-01  0.0002334 ***
## EVIP             0.003108        1       3.318e-05    0.95572    
## BRIEF            0.210671        1       2.249e-03    0.64782    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6. PPS analysis

Model 1

library(MatchIt)

## Warning: package 'MatchIt' was built under R version 4.1.3

library(randomForest)

## Warning: package 'randomForest' was built under R version 4.1.3

## randomForest 4.7-1.1

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:dplyr':
## 
##     combine

## The following object is masked from 'package:ggplot2':
## 
##     margin

library(marginaleffects)

## Warning: package 'marginaleffects' was built under R version 4.1.3

dat$ASD = factor(ifelse(dat$group == "ASD", 1, 0))
# dat$ASD = factor(dat$group)

prop1 <- matchit(ASD ~ IQ + EVIP + BRIEF, data = dat,
                 method = "nearest",
                 distance = "glm",
                 caliper = 0.2)

plot(prop1, type = "jitter", interactive = FALSE)

plot(prop1, type = "density", interactive = FALSE,
     which.xs = ~IQ + EVIP + BRIEF)

summary(prop1)

## 
## Call:
## matchit(formula = ASD ~ IQ + EVIP + BRIEF, data = dat, method = "nearest", 
##     distance = "glm", caliper = 0.2)
## 
## Summary of Balance for All Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance        0.6501        0.3499          1.2030     1.2510    0.3131
## IQ            112.3529      111.3824          0.0612     1.2195    0.0279
## EVIP          118.7353      118.2941          0.0270     1.5954    0.0492
## BRIEF          69.3235       53.3824          1.0587     2.2056    0.2964
##          eCDF Max
## distance   0.5000
## IQ         0.0882
## EVIP       0.1176
## BRIEF      0.5294
## 
## Summary of Balance for Matched Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance        0.5002        0.4861          0.0564     1.1166    0.0261
## IQ            113.2222      113.2778         -0.0035     1.9385    0.0468
## EVIP          115.2222      114.2778          0.0578     1.6287    0.0711
## BRIEF          59.8333       58.8333          0.0664     1.4998    0.0629
##          eCDF Max Std. Pair Dist.
## distance   0.1111          0.0875
## IQ         0.1667          1.0832
## EVIP       0.2222          1.1252
## BRIEF      0.2222          0.2730
## 
## Sample Sizes:
##           Control Treated
## All            34      34
## Matched        18      18
## Unmatched      16      16
## Discarded       0       0

plot(summary(prop1))

dat.prop1 <- match.data(prop1)

fit.prop1.TOM <- lm(TOMTB ~ ASD * (IQ + EVIP + BRIEF), 
          data = dat.prop1, weights = weights)

comp.prop1.TOM <- comparisons(fit.prop1.TOM,
                   variables = "ASD",
                   vcov = ~subclass,
                   newdata = subset(dat.prop1, ASD == 1),
                   wts = "weights")
summary(comp.prop1.TOM)

##   Group Term Contrast Effect Std. Error z value  Pr(>|z|)  2.5 %  97.5 %
## 1   ASD  ASD    1 - 0 -1.875     0.5866  -3.197 0.0013898 -3.025 -0.7254
## 
## Model type:  lm 
## Prediction type:  response

fit.prop1.SST <- lm(SST ~ ASD * (IQ + EVIP + BRIEF), 
          data = dat.prop1, weights = weights)

comp.prop1.SST <- comparisons(fit.prop1.SST,
                   variables = "ASD",
                   vcov = ~subclass,
                   newdata = subset(dat.prop1, ASD == 1),
                   wts = "weights")
summary(comp.prop1.SST)

##   Group Term Contrast Effect Std. Error z value Pr(>|z|)  2.5 % 97.5 %
## 1   ASD  ASD    1 - 0 -1.657      1.069   -1.55  0.12109 -3.753 0.4381
## 
## Model type:  lm 
## Prediction type:  response

Model 2

prop2 <- matchit(ASD ~ IQ + EVIP + BRIEF, data = dat,
                 method = "nearest",
                 distance = "glm",
                 replace = TRUE)

plot(prop2, type = "jitter", interactive = FALSE)

plot(prop2, type = "density", interactive = FALSE,
     which.xs = ~IQ + EVIP + BRIEF)

summary(prop2)

## 
## Call:
## matchit(formula = ASD ~ IQ + EVIP + BRIEF, data = dat, method = "nearest", 
##     distance = "glm", replace = TRUE)
## 
## Summary of Balance for All Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance        0.6501        0.3499          1.2030     1.2510    0.3131
## IQ            112.3529      111.3824          0.0612     1.2195    0.0279
## EVIP          118.7353      118.2941          0.0270     1.5954    0.0492
## BRIEF          69.3235       53.3824          1.0587     2.2056    0.2964
##          eCDF Max
## distance   0.5000
## IQ         0.0882
## EVIP       0.1176
## BRIEF      0.5294
## 
## Summary of Balance for Matched Data:
##          Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance        0.6501        0.6279          0.0892     1.0786    0.0394
## IQ            112.3529      108.7059          0.2301     1.3691    0.0929
## EVIP          118.7353      106.7059          0.7360     1.1088    0.1984
## BRIEF          69.3235       63.9706          0.3555     2.5066    0.1107
##          eCDF Max Std. Pair Dist.
## distance   0.2941          0.1471
## IQ         0.2059          1.1284
## EVIP       0.3529          1.4955
## BRIEF      0.4118          0.5157
## 
## Sample Sizes:
##               Control Treated
## All             34.        34
## Matched (ESS)    6.35      34
## Matched         13.        34
## Unmatched       21.         0
## Discarded        0.         0

plot(summary(prop2))

dat.prop2 <- get_matches(prop2)
fit.prop2.TOM <- lm(TOMTB ~ ASD * (IQ + EVIP + BRIEF), 
          data = dat.prop2, weights = weights)
comp.prop2.TOM <- comparisons(fit.prop2.TOM,
                   variables = "ASD",
                   vcov = ~subclass + ID,
                   newdata = subset(dat.prop2, ASD == 1),
                   wts = "weights")
summary(comp.prop2.TOM)

##   Group Term Contrast Effect Std. Error z value  Pr(>|z|)  2.5 %  97.5 %
## 1   ASD  ASD    1 - 0 -1.695      0.534  -3.174 0.0015018 -2.742 -0.6485
## 
## Model type:  lm 
## Prediction type:  response

fit.prop2.SST <- lm(SST ~ ASD * (IQ + EVIP + BRIEF), 
          data = dat.prop2, weights = weights)

comp.prop2.SST <- comparisons(fit.prop2.SST,
                   variables = "ASD",
                   vcov = ~subclass + ID,
                   newdata = subset(dat.prop2, ASD == 1),
                   wts = "weights")
summary(comp.prop2.SST)

##   Group Term Contrast Effect Std. Error z value Pr(>|z|)  2.5 % 97.5 %
## 1   ASD  ASD    1 - 0 -2.292      1.286  -1.783 0.074661 -4.811 0.2281
## 
## Model type:  lm 
## Prediction type:  response

5. Replication (alternative scoring)

ToM-TB scored as SST

KG_mult_justif <- lm(cbind(TOMTB_justif, SST) ~ group + IQ + EVIP + BRIEF, data = dat) 
car::Manova(KG_mult_justif, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df   Pr(>F)   
## group  1  0.101032   3.4840      2     62 0.036819 * 
## IQ     1  0.148470   5.4050      2     62 0.006858 **
## EVIP   1  0.033325   1.0687      2     62 0.349702   
## BRIEF  1  0.005587   0.1742      2     62 0.840576   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

KG_TOMTB_justif <- lm(TOMTB_justif ~ group + IQ + EVIP + BRIEF, data = dat)
emm.TOMTB.justif <- emmeans(KG_TOMTB_justif, "group")
broom::tidy(KG_TOMTB_justif)

## # A tibble: 5 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept) -4.54       2.79      -1.63   0.109 
## 2 groupCONT    1.67       0.664      2.51   0.0146
## 3 IQ           0.0531     0.0212     2.50   0.0149
## 4 EVIP         0.0328     0.0225     1.46   0.149 
## 5 BRIEF        0.00356    0.0228     0.156  0.876

smd(x = emm.TOMTB.justif, level = 97.5, df = dat, test = "TOMTB_justif")

## [1] "SMD = 0.67 [0.11, 1.22]"

SST scored as ToM-TB

KG_mult_dich <- lm(cbind(TOMTB, SST_dich) ~ group + IQ + EVIP + BRIEF, data = dat) 
car::Manova(KG_mult_dich, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df    Pr(>F)    
## group  1  0.162720   6.0247      2     62 0.0040643 ** 
## IQ     1  0.213029   8.3915      2     62 0.0005953 ***
## EVIP   1  0.000891   0.0276      2     62 0.9727524    
## BRIEF  1  0.003362   0.1046      2     62 0.9008571    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

KG_SST_dich <- lm(SST_dich ~ group + IQ + EVIP + BRIEF, data = dat)
emm.SST.dich <- emmeans(KG_SST_dich, "group")
broom::tidy(KG_SST_dich)

## # A tibble: 5 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)  1.59       2.03      0.784   0.436 
## 2 groupCONT    1.18       0.482     2.44    0.0173
## 3 IQ           0.0379     0.0154    2.46    0.0167
## 4 EVIP        -0.00386    0.0163   -0.236   0.814 
## 5 BRIEF        0.00160    0.0166    0.0968  0.923

smd(x = emm.SST.dich, level = 97.5, df = dat, test = "SST_dich")

## [1] "SMD = 0.69 [0.13, 1.25]"

---

6. Diagnostics

Plots

TOMTB

par(mar = c(4, 4, .1, .1))
plot(KG_TOMTB)

SST

par(mar = c(4, 4, .1, .1))
plot(KG_SST)

Outliers

TOMTB_cooksd <- cooks.distance(KG_TOMTB)
SST_cooksd <- cooks.distance(KG_SST)
TOMTB_out = as.numeric(which(TOMTB_cooksd > 4/68))
SST_out = as.numeric(which(SST_cooksd > 4/68))
outliers_ID = unique(c(TOMTB_out, SST_out))
length(unique(outliers_ID))

## [1] 6

Analyses

KG_CD <- lm(cbind(TOMTB, SST) ~ group + IQ + EVIP + BRIEF, data = dat[-outliers_ID, ]) 

KG_TOMTB_CD <- lm(TOMTB ~ group + IQ + EVIP + BRIEF, 
                  data = dat[-outliers_ID, ])

KG_SST_CD <- lm(SST ~ group + IQ + EVIP + BRIEF, 
                data = dat[-outliers_ID, ])

Multivariate

car::Manova(KG_CD, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##       Df test stat approx F num Df den Df    Pr(>F)    
## group  1  0.090025   2.7701      2     56 0.0712558 .  
## IQ     1  0.261334   9.9062      2     56 0.0002073 ***
## EVIP   1  0.002439   0.0684      2     56 0.9339184    
## BRIEF  1  0.027918   0.8042      2     56 0.4525657    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ToM-TB

broom::tidy(KG_TOMTB_CD)

## # A tibble: 5 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)  2.75       2.48      1.11   0.273    
## 2 groupCONT    1.34       0.567     2.35   0.0221   
## 3 IQ           0.0814     0.0190    4.29   0.0000710
## 4 EVIP        -0.00322    0.0186   -0.173  0.863    
## 5 BRIEF        0.00131    0.0189    0.0692 0.945

emm.TOMTB_CD <- emmeans(KG_TOMTB_CD, "group")
smd(x = emm.TOMTB_CD, level = 97.5, df = dat[-outliers_ID, ], test = "TOMTB")

## [1] "SMD = 0.65 [0.06, 1.23]"

SST

broom::tidy(KG_SST_CD)

## # A tibble: 5 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)  0.300      3.95      0.0760 0.940  
## 2 groupCONT    1.14       0.902     1.26   0.212  
## 3 IQ           0.0904     0.0302    2.99   0.00409
## 4 EVIP         0.00678    0.0295    0.230  0.819  
## 5 BRIEF       -0.0341     0.0301   -1.14   0.261

emm.SST_CD <- emmeans(KG_SST_CD, "group")
smd(x = emm.SST_CD, level = 97.5, df = dat[-outliers_ID, ], test = "SST")

## [1] "SMD = 0.37 [-0.21, 0.95]"

Robust regression

bartlett.test(TOMTB ~ group, data = dat)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  TOMTB by group
## Bartlett's K-squared = 4.9033, df = 1, p-value = 0.02681

bartlett.test(SST ~ group, data = dat)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  SST by group
## Bartlett's K-squared = 0.14753, df = 1, p-value = 0.7009

summary(robustbase::lmrob(TOMTB ~ group + IQ + EVIP + BRIEF, data = dat, setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = TOMTB ~ group + IQ + EVIP + BRIEF, data = dat, 
##     setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4706 -1.1603 -0.1339  1.5180  3.8008 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.9068420  2.5831742   1.125  0.26473    
## groupCONT    1.7616575  0.6096676   2.890  0.00528 ** 
## IQ           0.0695451  0.0194563   3.574  0.00068 ***
## EVIP        -0.0002205  0.0205648  -0.011  0.99148    
## BRIEF        0.0095707  0.0210017   0.456  0.65017    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 2.024 
## Multiple R-squared:  0.3062, Adjusted R-squared:  0.2622 
## Convergence in 7 IRWLS iterations
## 
## Robustness weights: 
##  46 weights are ~= 1. The remaining 22 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.7825  0.8941  0.9604  0.9356  0.9774  0.9953 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         1.471e-03         2.710e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

summary(robustbase::lmrob(SST ~ group + IQ + EVIP + BRIEF, data = dat, setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = SST ~ group + IQ + EVIP + BRIEF, data = dat, 
##     setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7635 -2.4918  0.5717  2.6026  5.2147 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -0.524582   3.967648  -0.132  0.89524   
## groupCONT    1.345280   0.948358   1.419  0.16096   
## IQ           0.088128   0.030361   2.903  0.00509 **
## EVIP         0.007585   0.031913   0.238  0.81291   
## BRIEF       -0.021297   0.032497  -0.655  0.51463   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 3.214 
## Multiple R-squared:  0.2084, Adjusted R-squared:  0.1582 
## Convergence in 8 IRWLS iterations
## 
## Robustness weights: 
##  49 weights are ~= 1. The remaining 19 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.6954  0.8768  0.9742  0.9317  0.9846  0.9977 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         1.471e-03         2.710e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

---

S4: Convergent (Vineland)

1. Main analysis

A. Plot

ggplot(dat_long, aes(x = ToM, y = VINELAND_SOC, color = group)) + 
  geom_point(size = 3, alpha = 0.6) + 
  geom_smooth(method = "lm", fullrange = TRUE, alpha = 0.3, aes(fill = group)) + 
  facet_grid(~ToM_test, scales="free_x") + 
  theme_bw() + 
  geom_smooth(fullrange = TRUE, linetype = "dashed", alpha = 0.6, se = FALSE) + 
  scale_color_manual(values=c("#8BC2FF", "#F1C12C")) +
  scale_fill_manual(values=c("#8BC2FF", "#F1C12C")) + 
  labs(y = "Social compentence",
       x = "Performance on ToM tests") +
  theme(legend.position = "top")

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

par(mar = c(2, 1, .1, .1))
ggstatsplot::ggcorrmat(dat[which(dat$group == "ASD"),
     c('VINELAND_SOC', 'TOMTB', 'SST')],
    p.adjust.method = "none",
    ggtheme = ggplot2::theme_bw(),
    colors = c("#0016FF", "white", "#FF0000"),
    title = "ASD group",
    pch = 13)

ggstatsplot::ggcorrmat(dat[which(dat$group == "CONT"),
     c('VINELAND_SOC', 'TOMTB', 'SST')],
    p.adjust.method = "none",
    ggtheme = ggplot2::theme_bw(),
    colors = c("#0016FF", "white", "#FF0000"),
    title = "CONT group",
    pch = 13)

---

B. Multivariate effect

CONV_mult <- lm(cbind(TOMTB, SST) ~ VINELAND_SOC * group, data = dat)
car::Manova(CONV_mult, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##                    Df test stat approx F num Df den Df Pr(>F)
## VINELAND_SOC        1  0.018381  0.58984      2     63 0.5574
## group               1  0.022420  0.72243      2     63 0.4896
## VINELAND_SOC:group  1  0.017846  0.57237      2     63 0.5671

---

C. Univariate effects

Post-hoc tests : groups and tests taken separately

emt.CONV.comp <- emtrends(CONV_mult, ~ ToMTest1 | group, 
                          var = "VINELAND_SOC", 
                          mult.name = "ToMTest1")

test(emt.CONV.comp, adjust = "none")

## group = ASD:
##  ToMTest1 VINELAND_SOC.trend     SE df t.ratio p.value
##  TOMTB               0.00619 0.0154 64   0.401  0.6897
##  SST                 0.01100 0.0237 64   0.464  0.6440
## 
## group = CONT:
##  ToMTest1 VINELAND_SOC.trend     SE df t.ratio p.value
##  TOMTB               0.02435 0.0199 64   1.225  0.2250
##  SST                -0.00149 0.0305 64  -0.049  0.9611

---

2. Comparison

contrast(emt.CONV.comp, interaction = "consec", by = "group")

## group = ASD:
##  ToMTest1_consec estimate     SE df t.ratio p.value
##  SST - TOMTB      0.00481 0.0207 64   0.232  0.8169
## 
## group = CONT:
##  ToMTest1_consec estimate     SE df t.ratio p.value
##  SST - TOMTB     -0.02585 0.0266 64  -0.971  0.3354

---

3. Diagnostics

Outliers

CONV_TOMTB <- lm(TOMTB ~ VINELAND_SOC * group, data = dat)
CONV_SST <- lm(SST ~ VINELAND_SOC * group, data = dat)

CONV_TOMTB_cooksd <- cooks.distance(CONV_TOMTB)
CONV_SST_cooksd <- cooks.distance(CONV_SST)
CONV_TOMTB_out = as.numeric(which(CONV_TOMTB_cooksd > 4/68))
CONV_SST_out = as.numeric(which(CONV_SST_cooksd > 4/68))
outliers_ID_CONV = unique(c(CONV_TOMTB_out, CONV_SST_out))
length(unique(outliers_ID))

## [1] 6

CONV_mult_CD <- lm(cbind(TOMTB, SST) ~ VINELAND_SOC * group, 
                   data = dat[-outliers_ID_CONV, ])
car::Manova(CONV_mult_CD, test.statistic = "Pillai")

## 
## Type II MANOVA Tests: Pillai test statistic
##                    Df test stat approx F num Df den Df Pr(>F)
## VINELAND_SOC        1  0.006088  0.17151      2     56 0.8428
## group               1  0.029839  0.86118      2     56 0.4282
## VINELAND_SOC:group  1  0.043443  1.27166      2     56 0.2883

emt.CONV.comp.CD <- emtrends(CONV_mult_CD, ~ ToMTest1 | group, 
                          var = "VINELAND_SOC", 
                          mult.name = "ToMTest1")

test(emt.CONV.comp.CD, adjust = "none")

## group = ASD:
##  ToMTest1 VINELAND_SOC.trend     SE df t.ratio p.value
##  TOMTB              -0.00889 0.0192 57  -0.463  0.6449
##  SST                 0.00275 0.0270 57   0.102  0.9191
## 
## group = CONT:
##  ToMTest1 VINELAND_SOC.trend     SE df t.ratio p.value
##  TOMTB               0.04129 0.0260 57   1.590  0.1173
##  SST                 0.01947 0.0365 57   0.533  0.5960

Robust regression

summary(robustbase::lmrob(TOMTB ~ VINELAND_SOC, 
                          data = dat[which(dat$group == "ASD"),], setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = TOMTB ~ VINELAND_SOC, data = dat[which(dat$group == 
##     "ASD"), ], setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3915 -1.5096  0.4365  2.4697  3.5671 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.192593   1.241215   9.017 2.67e-10 ***
## VINELAND_SOC  0.004144   0.019445   0.213    0.833    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 2.687 
## Multiple R-squared:  0.001468,   Adjusted R-squared:  -0.02974 
## Convergence in 6 IRWLS iterations
## 
## Robustness weights: 
##  26 weights are ~= 1. The remaining 8 ones are
##     11     18     23     42     53     56     61     64 
## 0.7340 0.9675 0.8721 0.8741 0.9613 0.9735 0.9617 0.9543 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         2.941e-03         2.165e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

summary(robustbase::lmrob(SST ~ VINELAND_SOC, 
                          data = dat[which(dat$group == "ASD"),], setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = SST ~ VINELAND_SOC, data = dat[which(dat$group == 
##     "ASD"), ], setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2847 -2.8740 -0.5691  2.8884  5.8895 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.00285    1.68328   4.754 4.05e-05 ***
## VINELAND_SOC  0.01245    0.02664   0.467    0.644    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 3.616 
## Multiple R-squared:  0.007106,   Adjusted R-squared:  -0.02392 
## Convergence in 7 IRWLS iterations
## 
## Robustness weights: 
##  22 weights are ~= 1. The remaining 12 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8327  0.9473  0.9784  0.9535  0.9934  0.9985 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         2.941e-03         2.165e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

summary(robustbase::lmrob(TOMTB ~ VINELAND_SOC, 
                          data = dat[which(dat$group == "CONT"),], setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = TOMTB ~ VINELAND_SOC, data = dat[which(dat$group == 
##     "CONT"), ], setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6134 -1.1293  0.2089  1.3860  2.5272 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.00593    1.81258   5.520 4.37e-06 ***
## VINELAND_SOC  0.02653    0.01621   1.636    0.112    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 1.739 
## Multiple R-squared:  0.07907,    Adjusted R-squared:  0.0503 
## Convergence in 5 IRWLS iterations
## 
## Robustness weights: 
##  25 weights are ~= 1. The remaining 9 ones are
##     19     25     30     31     46     57     58     65     66 
## 0.9223 0.8623 0.9949 0.9893 0.9636 0.7061 0.9981 0.9226 0.8722 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         2.941e-03         2.547e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

summary(robustbase::lmrob(SST ~ VINELAND_SOC, 
                          data = dat[which(dat$group == "CONT"),], setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = SST ~ VINELAND_SOC, data = dat[which(dat$group == 
##     "CONT"), ], setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6729 -1.6214  0.5454  1.6016  4.6203 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  9.799670   3.559516   2.753  0.00965 **
## VINELAND_SOC 0.006237   0.031915   0.195  0.84630   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 3.297 
## Multiple R-squared:  0.001295,   Adjusted R-squared:  -0.02991 
## Convergence in 8 IRWLS iterations
## 
## Robustness weights: 
##  23 weights are ~= 1. The remaining 11 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.7271  0.8604  0.9453  0.9052  0.9755  0.9976 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         2.941e-03         2.547e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

---

4. Equivalence tests

res_TOMTB_SST = data.frame(res = c("COR", "TOSTER"))
for (groups in c("ASD", "CONT")) {
  for (tests in c("TOMTB", "SST")) {
    res_TOMTB_SST[1 ,c(paste0(groups, "_", tests))] <- 
      # cor
      cor.test(~dat[which(dat$group == groups), "VINELAND_SOC"] + 
                 dat[which(dat$group == groups), tests])$estimate
      # TOST
        res_TOMTB_SST[2 ,c(paste0(groups, "_", tests))] <- 
          max(
            abs(
            TOSTER::TOSTr(n = 34, r = res_TOMTB_SST[1 ,c(paste0(groups, "_", tests))], 
              low_eqbound_r = -.40,
              high_eqbound_r = .40, plot = FALSE, alpha = 0.05, verbose = FALSE)$LL_CI_TOST
            ), 
            abs(
            TOSTER::TOSTr(n = 34, r = res_TOMTB_SST[1 ,c(paste0(groups, "_", tests))], 
              low_eqbound_r = -.40,
              high_eqbound_r = .40, plot = FALSE, alpha = 0.05, verbose = FALSE)$UL_CI_TOST
          ) 
        ) 
  }
}
res_TOMTB_SST

##      res  ASD_TOMTB    ASD_SST CONT_TOMTB    CONT_SST
## 1    COR 0.05970783 0.07932218  0.2702624 -0.00895625
## 2 TOSTER 0.34098206 0.35828210  0.5172452  0.29531662

---

S5: Convergent (ToM tests)

1. Main analysis

A. Plot

ggplot(dat, aes(x = TOMTB, y = SST, color = group)) + 
  geom_point(size = 3, alpha = 0.6) + 
  geom_smooth(aes(fill = group), method = "lm", fullrange = TRUE) + 
  geom_smooth(aes(fill = group), se = FALSE, linetype = "dashed") + 
  facet_wrap(~group, scale = "free_x") +
  theme_bw() + 
  scale_y_continuous(name = "Performance at the ToM tests") +
  scale_color_manual(values=c("#8BC2FF", "#F1C12C")) +
  scale_fill_manual(values=c("#8BC2FF", "#F1C12C")) + 
  labs(y = "Performance at the SST",
       x = "Performance at the ToM-TB") 

## `geom_smooth()` using formula 'y ~ x'

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

B. Correlation

CONC_ASD <- cor.test(~ TOMTB + SST, 
                     data = subset(dat, group == "ASD"),
                     conf.level = 1-(0.05/(2)))
CONC_CONT <- cor.test(~ TOMTB + SST, 
                      data = subset(dat, group != "ASD"),
                     conf.level = 1-(0.05/(2)))
broom::tidy(CONC_ASD)

## # A tibble: 1 x 8
##   estimate statistic p.value parameter conf.low conf.high method         alter~1
##      <dbl>     <dbl>   <dbl>     <int>    <dbl>     <dbl> <chr>          <chr>  
## 1    0.515      3.40 0.00182        32    0.166     0.750 Pearson's pro~ two.si~
## # ... with abbreviated variable name 1: alternative

broom::tidy(CONC_CONT)

## # A tibble: 1 x 8
##   estimate statistic p.value parameter conf.low conf.high method         alter~1
##      <dbl>     <dbl>   <dbl>     <int>    <dbl>     <dbl> <chr>          <chr>  
## 1    0.500      3.26 0.00262        32    0.145     0.740 Pearson's pro~ two.si~
## # ... with abbreviated variable name 1: alternative

2. Replication (advanced score)

CONC_ASD_adv <- cor.test(~ TOMTB_advanced + SST, 
                         data = subset(dat, group == "ASD"),
                     conf.level = 1-(0.05/(2)))
CONC_CONT_adv <- cor.test(~ TOMTB_advanced + SST, 
                          data = subset(dat, group != "ASD"),
                     conf.level = 1-(0.05/(2)))
broom::tidy(CONC_ASD_adv)

## # A tibble: 1 x 8
##   estimate statistic p.value parameter conf.low conf.high method         alter~1
##      <dbl>     <dbl>   <dbl>     <int>    <dbl>     <dbl> <chr>          <chr>  
## 1    0.316      1.88  0.0690        32  -0.0757     0.623 Pearson's pro~ two.si~
## # ... with abbreviated variable name 1: alternative

broom::tidy(CONC_CONT_adv)

## # A tibble: 1 x 8
##   estimate statistic p.value parameter conf.low conf.high method         alter~1
##      <dbl>     <dbl>   <dbl>     <int>    <dbl>     <dbl> <chr>          <chr>  
## 1    0.340      2.05  0.0488        32  -0.0479     0.639 Pearson's pro~ two.si~
## # ... with abbreviated variable name 1: alternative

3. Replication (regression)

CONC_INT <- lm(TOMTB ~ SST * group, data = dat)
test(emtrends(CONC_INT, ~group, var = "SST"), adjust = "bonferroni")

##  group SST.trend    SE df t.ratio p.value
##  ASD       0.385 0.097 64   3.974  0.0004
##  CONT      0.270 0.104 64   2.600  0.0231
## 
## P value adjustment: bonferroni method for 2 tests

4. Diagnostics

Outliers

SST_TOMTB <- lm(TOMTB ~ SST * group, data = dat)

SST_TOMTB_cooksd <- cooks.distance(SST_TOMTB)
SST_TOMTB_out = as.numeric(which(SST_TOMTB_cooksd > 4/68))
outliers_ID_SST_TOMTB = as.numeric(which(SST_TOMTB_out > 4/68))
length(unique(outliers_ID))

## [1] 6

CONV_SST_TOMTB_CD <- CONC_INT <- lm(TOMTB ~ SST * group, data = dat[-outliers_ID_SST_TOMTB, ])
test(emtrends(CONV_SST_TOMTB_CD, ~group, var = "SST"), adjust = "bonferroni")

##  group SST.trend     SE df t.ratio p.value
##  ASD       0.390 0.0961 60   4.053  0.0003
##  CONT      0.276 0.1020 60   2.706  0.0177
## 
## P value adjustment: bonferroni method for 2 tests

CONV_SST_TOMTB_CD <- CONC_INT <- lm(TOMTB ~ SST * group, data = dat)
test(emtrends(CONV_SST_TOMTB_CD, ~group, var = "SST"), adjust = "bonferroni")

##  group SST.trend    SE df t.ratio p.value
##  ASD       0.385 0.097 64   3.974  0.0004
##  CONT      0.270 0.104 64   2.600  0.0231
## 
## P value adjustment: bonferroni method for 2 tests

Robust regression

summary(robustbase::lmrob(TOMTB ~ VINELAND_SOC, 
                          data = dat[which(dat$group == "ASD"),],
                          setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = TOMTB ~ VINELAND_SOC, data = dat[which(dat$group == 
##     "ASD"), ], setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3915 -1.5096  0.4365  2.4697  3.5671 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  11.192593   1.241215   9.017 2.67e-10 ***
## VINELAND_SOC  0.004144   0.019445   0.213    0.833    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 2.687 
## Multiple R-squared:  0.001468,   Adjusted R-squared:  -0.02974 
## Convergence in 6 IRWLS iterations
## 
## Robustness weights: 
##  26 weights are ~= 1. The remaining 8 ones are
##     11     18     23     42     53     56     61     64 
## 0.7340 0.9675 0.8721 0.8741 0.9613 0.9735 0.9617 0.9543 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         2.941e-03         2.165e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

summary(robustbase::lmrob(SST ~ VINELAND_SOC, 
                          data = dat[which(dat$group == "ASD"),],
                          setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = SST ~ VINELAND_SOC, data = dat[which(dat$group == 
##     "ASD"), ], setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2847 -2.8740 -0.5691  2.8884  5.8895 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   8.00285    1.68328   4.754 4.05e-05 ***
## VINELAND_SOC  0.01245    0.02664   0.467    0.644    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 3.616 
## Multiple R-squared:  0.007106,   Adjusted R-squared:  -0.02392 
## Convergence in 7 IRWLS iterations
## 
## Robustness weights: 
##  22 weights are ~= 1. The remaining 12 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8327  0.9473  0.9784  0.9535  0.9934  0.9985 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         2.941e-03         2.165e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

summary(robustbase::lmrob(TOMTB ~ VINELAND_SOC, 
                          data = dat[which(dat$group == "CONT"),],
                          setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = TOMTB ~ VINELAND_SOC, data = dat[which(dat$group == 
##     "CONT"), ], setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6134 -1.1293  0.2089  1.3860  2.5272 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  10.00593    1.81258   5.520 4.37e-06 ***
## VINELAND_SOC  0.02653    0.01621   1.636    0.112    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 1.739 
## Multiple R-squared:  0.07907,    Adjusted R-squared:  0.0503 
## Convergence in 5 IRWLS iterations
## 
## Robustness weights: 
##  25 weights are ~= 1. The remaining 9 ones are
##     19     25     30     31     46     57     58     65     66 
## 0.9223 0.8623 0.9949 0.9893 0.9636 0.7061 0.9981 0.9226 0.8722 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         2.941e-03         2.547e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)

summary(robustbase::lmrob(SST ~ VINELAND_SOC, 
                          data = dat[which(dat$group == "CONT"),],
                          setting = "KS2014"))

## 
## Call:
## robustbase::lmrob(formula = SST ~ VINELAND_SOC, data = dat[which(dat$group == 
##     "CONT"), ], setting = "KS2014")
##  \--> method = "SMDM"
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.6729 -1.6214  0.5454  1.6016  4.6203 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  9.799670   3.559516   2.753  0.00965 **
## VINELAND_SOC 0.006237   0.031915   0.195  0.84630   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Robust residual standard error: 3.297 
## Multiple R-squared:  0.001295,   Adjusted R-squared:  -0.02991 
## Convergence in 8 IRWLS iterations
## 
## Robustness weights: 
##  23 weights are ~= 1. The remaining 11 ones are summarized as
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.7271  0.8604  0.9453  0.9052  0.9755  0.9976 
## Algorithmic parameters: 
##       tuning.chi1       tuning.chi2       tuning.chi3       tuning.chi4 
##        -5.000e-01         1.500e+00                NA         5.000e-01 
##                bb       tuning.psi1       tuning.psi2       tuning.psi3 
##         5.000e-01        -5.000e-01         1.500e+00         9.500e-01 
##       tuning.psi4        refine.tol           rel.tol         scale.tol 
##                NA         1.000e-07         1.000e-07         1.000e-10 
##         solve.tol       eps.outlier             eps.x warn.limit.reject 
##         1.000e-07         2.941e-03         2.547e-10         5.000e-01 
## warn.limit.meanrw 
##         5.000e-01 
##      nResample         max.it       best.r.s       k.fast.s          k.max 
##           1000            500             20              2           2000 
##    maxit.scale      trace.lev            mts     compute.rd      numpoints 
##            200              0           1000              0             10 
## fast.s.large.n 
##           2000 
##               setting                   psi           subsampling 
##              "KS2014"                 "lqq"         "nonsingular" 
##                   cov compute.outlier.stats 
##             ".vcov.w"                "SMDM" 
## seed : int(0)