A tutorial to select participants and
create two balanced groups
using propensity scores
Janet Y. Bang, Megha Sharda, and Aparna Nadig
last update: 31 May 2020
A template of this workflow can be downloaded here.
All data and code for this workflow can be downloaded here.
This tutorial follows the workflow below to select participants for matched groups (Bang, Sharda, & Nadig, in preparation). Please refer to Code buttons to show/hide code as desired.
include_graphics("diagram_4steps.png")
Load libraries and set theme
# load libraries
library(tidyverse)
library(MatchIt)
library(GGally)
library(compute.es)
library(compareGroups)
library(gridExtra)
# set theme to black and white
theme_set(theme_bw())Read in data
# read in data
demo <- read_csv("./data/demographics.csv")
celf <- read_csv("./data/celf.csv")
leiter <- read_csv("./data/leiter.csv")
scq <- read_csv("./data/scq.csv")
vineland <- read_csv("./data/vineland.csv")
Step 1: Assess data
The final possible sample size includes 43 TD children and 25 children with ASD.
Five TD children and three children from the ASD group were tested but excluded from the study.
Select all possible participants to include for matching and combine data into one dataframe (df)
# create new df
df <- demo %>%
filter(keep == 1) %>% # final children to include
left_join(celf, by = "participant") %>%
left_join(leiter, by = "participant") %>%
left_join(scq, by = "participant") %>%
left_join(vineland, by = "participant")
# sample size count
df %>% group_by(group) %>% count()## # A tibble: 2 x 2
## # Groups: group [2]
## group n
## <chr> <int>
## 1 ASD 25
## 2 TD 43
Step 2: Select covariates
Prior to data collection, we considered multiple covariates that could influence children’s performance on the experimental tasks. Ultimately, we could not include all possible covariates given our sample size, thus we describe how we decided on our final two covariates of age and IQ.
The covariates we considered were based on their known relationships with referential gaze following and word learning in experimental tasks. These covariates included: age, nonverbal IQ, language ability, sex, and parental education. We then decided to exclude covariates of sex and parental education. Though both variables are known to be related to language abilities in the general population, we prioritized age, nonverbal IQ, and language ability due to their known relations with referential gaze following and word learning in children with ASD.
We examined the distributions and interrelations between age, nonverbal IQ, and language abilities seen in the figure below, comparing similarities and differences in their distributions between groups. From these plots we concluded that while distributions were similar on age and IQ, wide heterogeneity was seen in the ASD group on language variables that were not seen to the same extent in the TD group. Whereas the ASD group ranged widely from 2 SD below the mean to 2 SD above the mean on the normed language measures, there were no children in the TD group 2 SD below the mean. For Recalling Sentences, the distribution of scores in the TD group was skewed with 73% at or above the mean, whereas this was the case for only 40% children with ASD. In contrast to the different distributions between groups on Word Classes and Recalling Sentences, the performance of children with ASD and TD children were more similarly distributed on Word Associations.
Thus on measures of semantic language abilities, children with ASD versus TD children demonstrated both similarities (as measured by the Word Associations subtest) and weaknesses (as measured by the Word Classes subtest), whereas a large proportion of children with ASD had weaker structural language abilities (as measured by the Recalling Sentences subtest).
Select participant variables
# create df with relevant variables
df2 <- df %>%
select(participant, keep, group, age_yrs, language, sex, block1, parental_ed,
leiter_composite, celf4_RS_scaled, celf4_WC_rec_scaled, celf4_WC_exp_scaled, celf4_WC_total_scaled,
celf4_WA_total, scq_total, social_standScore)Create matrix of variables to consider for propensity scores
# create df with variables for the matrix
mtx <- df2 %>%
select("age_yrs", "leiter_composite", "celf4_RS_scaled", "celf4_WC_total_scaled", "celf4_WA_total", "group") %>%
rename("Age (yr)" = age_yrs, "Leiter (IQ)" = leiter_composite, "CELF4-RS" = celf4_RS_scaled,
"CELF4-WC" = celf4_WC_total_scaled, "CELF4-WA" = celf4_WA_total, "Group" = group)
# age: years
# Leiter: Nonverbal IQ, composite score (normed: M = 100, SD = 15)
# CELF4 - RS: CELF-4 Recalling Sentences, scaled score (normed: M = 10, SD = 3)
# CELF4 - WC: CELF-4 Word Classes Total (Receptive and Expressive), scaled score (normed: M = 10, SD = 3)
# CELF4 - WA: CELF-4 Word Associations, total score
# create matrix plot
corr_mtx <- ggpairs(data = mtx,
columns = 1:6,
title = "Correlation Matrix of Potential Covariates",
mapping = aes(colour = Group, alpha = 1),
axisLabels = "show") +
theme(text= element_text(size = 15))
# print and save plot
corr_mtx
CELF-4-RS = CELF-4 Recalling Sentences (scaled scores); CELF-4-WC = CELF-4 Word Classes (scaled scores) and CELF-4-WA = CELF-4 Word Associations (raw scores, because no normative data is provided). The density plots arranged diagonally and the histograms at the bottom can be used to compare the distribution of scores in ASD and TD groups. The scatterplots and correlations on either side of the density plots can be used to examine the strength of the relation between variables both within and across groups. The boxplots and bar graphs on the right-hand side provide another way to compare distributions between groups.
ggsave("./plots/corr_mtx.pdf", height = 8, width = 11, units = "in")Taken together, the above matrix indicated that matching on all variables of age, nonverbal IQ, and the three language measures would likely result in less than half of our sample of children with ASD to be matched to a group with TD children. To retain as many children with ASD as possible, we matched on age and nonverbal IQ because both children with ASD and TD children were similarly distributed on these measures and both covariates have demonstrated relations with our primary experimental manipulation: how children learn using referential gaze. The choice of two covariates followed the guidelines by Blackford (2007), who suggests an approximate ratio of 5-10 participants per covariate.
Though one of our outcomes measures of interest was word learning, and children with stronger language abilities have been shown to have better performance on word learning tasks, we concluded that the overarching goal was to examine how children use referential gaze to learn in two different contexts (word learning and action understanding).
Additionally, our experimental manipulation tested how children with ASD treated referential gaze in contrast to an arrow cue, and there was no theoretical or evidence-based rationale suggesting that having stronger language abilities meant that children would learn better with one cue vs. another (or the contrary, that weaker language abilities meant that children would learn worse with one cue versus another). Thus, we reasoned that regardless of language ability, it was still possible to test whether children learned new words differently with a referential gaze cue vs. an arrow cue.
Another bonus of including all children with ASD was that the sample would reflect part of the heterogeneity seen in language abilities (keeping in mind that our sample was selected to have nonverbal IQ in the normal range, and therefore does not represent a full range of language abilities). We decided to further investigate the role of language ability on our experimental measures as a part of our exploratory analyses.
Step 3: Conduct matching
We used the MatchIt package (Ho et al., 2011) in R version 3.6.1 (R Core Team, 2019) to apply both the nearest neighbor and optimal matching methods.
Both methods resulted in the same 25 TD children chosen as matches to the 25 children with ASD.
As seen in Table 2, nonverbal IQ is significantly higher in the TD group (n = 43) than in the ASD group (n = 25; p = .039), although on age both groups share similar means and standard deviations (p = .570). Thus, our next step was to select participants such that groups were balanced on both age and nonverbal IQ.
Examine descriptive statistics of age and IQ before matching
# subset only age and IQ
df2_age_IQ <- df2 %>%
select(group, age_yrs, leiter_composite)
# compare groups
df2_age_IQ_compare <- compareGroups(group ~ age_yrs + leiter_composite , data = df2_age_IQ)
df2_age_IQ_compare_tbl <- createTable(df2_age_IQ_compare, digits = 2, type = 1, show.n = TRUE)
print(df2_age_IQ_compare_tbl, header.labels = c(p.overall = "p-value"))##
## --------Summary descriptives table by 'group'---------
##
## _________________________________________________________
## ASD TD p-value N
## N=25 N=43
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## age_yrs 8.93 (1.34) 8.75 (1.14) 0.570 68
## leiter_composite 107.16 (13.61) 114.53 (14.18) 0.039 68
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Create dataframe with covariates to be used for propensity scores
# create df
df_ps_n68 <- df2 %>%
select(participant, group, age_yrs, leiter_composite) %>%
mutate(group.int = ifelse(group == "ASD", 1, 0)) # group must be an integer for MatchItConduct Nearest neighbor algorithm: n = 68
# apply matching
m_near_n68 <- matchit(group.int ~ age_yrs + leiter_composite,
data = df_ps_n68, method = "nearest")
# get df of matched participants and associated covariates
m_near_n68_df <- match.data(m_near_n68)
head(m_near_n68_df)## participant group age_yrs leiter_composite group.int distance weights
## 1 302 ASD 8.83 123 1 0.2575106 1
## 2 304 ASD 8.75 123 1 0.2539244 1
## 3 306 ASD 10.00 104 1 0.5088229 1
## 5 310 ASD 9.50 119 1 0.3254607 1
## 7 312 ASD 10.75 109 1 0.4991269 1
## 9 314 ASD 10.17 92 1 0.6439131 1Conduct Optimal matching algorithm: n = 68
# apply matching
m_opt_n68 <- matchit(group.int ~ age_yrs + leiter_composite,
data = df_ps_n68, method = "optimal")
# get df of matched participants and associated covariates
m_opt_n68_df <- match.data(m_opt_n68)
head(m_opt_n68_df)## participant group age_yrs leiter_composite group.int distance weights
## 1 302 ASD 8.83 123 1 0.2575106 1
## 2 304 ASD 8.75 123 1 0.2539244 1
## 3 306 ASD 10.00 104 1 0.5088229 1
## 5 310 ASD 9.50 119 1 0.3254607 1
## 7 312 ASD 10.75 109 1 0.4991269 1
## 9 314 ASD 10.17 92 1 0.6439131 1
## subclass
## 1 1
## 2 22
## 3 25
## 5 10
## 7 17
## 9 18Compare participants selected with nearest neighbor vs. optimal matching: n = 68
This side-by-side comparison demonstrates that with nearest neighbor and optimal matching methods, the same set of 25 TD children were chosen as matches to the 25 children with ASD.
# create dfs of ASD and TD participants selected from nearest neighbor and optimal matching methods
# English = 300 series
# French = 400 series
# ASD participants = even numbers
# TD participants = odd numbers
m_near_n68_part <- m_near_n68_df %>%
dplyr::select(c(participant, group)) %>%
rename(m_near_participants = participant,
m_near_group = group) %>%
arrange(m_near_group)
m_opt_n68_part <- m_opt_n68_df %>%
dplyr::select(c(participant, group)) %>%
rename(m_opt_participants = participant,
m_opt_group = group) %>%
arrange(m_opt_group)
# compare participants selected from nearest neighbor vs. optimal matching methods
check_n68_df <- cbind(m_near_n68_part, m_opt_n68_part)
check_n68_df## m_near_participants m_near_group m_opt_participants m_opt_group
## 1 302 ASD 302 ASD
## 2 304 ASD 304 ASD
## 3 306 ASD 306 ASD
## 4 310 ASD 310 ASD
## 5 312 ASD 312 ASD
## 6 314 ASD 314 ASD
## 7 318 ASD 318 ASD
## 8 322 ASD 322 ASD
## 9 324 ASD 324 ASD
## 10 326 ASD 326 ASD
## 11 330 ASD 330 ASD
## 12 332 ASD 332 ASD
## 13 402 ASD 402 ASD
## 14 404 ASD 404 ASD
## 15 406 ASD 406 ASD
## 16 408 ASD 408 ASD
## 17 410 ASD 410 ASD
## 18 412 ASD 412 ASD
## 19 414 ASD 414 ASD
## 20 416 ASD 416 ASD
## 21 418 ASD 418 ASD
## 22 420 ASD 420 ASD
## 23 424 ASD 424 ASD
## 24 428 ASD 428 ASD
## 25 430 ASD 430 ASD
## 26 315 TD 315 TD
## 27 317 TD 317 TD
## 28 319 TD 319 TD
## 29 323 TD 323 TD
## 30 325 TD 325 TD
## 31 329 TD 329 TD
## 32 331 TD 331 TD
## 33 333 TD 333 TD
## 34 341 TD 341 TD
## 35 347 TD 347 TD
## 36 403 TD 403 TD
## 37 407 TD 407 TD
## 38 409 TD 409 TD
## 39 413 TD 413 TD
## 40 415 TD 415 TD
## 41 417 TD 417 TD
## 42 421 TD 421 TD
## 43 425 TD 425 TD
## 44 427 TD 427 TD
## 45 429 TD 429 TD
## 46 433 TD 433 TD
## 47 435 TD 435 TD
## 48 439 TD 439 TD
## 49 443 TD 443 TD
## 50 445 TD 445 TD
Step 4: Diagnose matching
We used visual inspection and statistics to examine how well groups were matched on their propensity scores. As seen in the figure below, visual inspection of propensity score plots depicted the same ASD participant with a high propensity score (Matched Treatment Units; this participant has a propensity score of approximately 0.8) without a close match among the selected matches in the TD group (Matched Control Units).
We next examined matching of propensity scores based on cut-off values proposed in the literature: a maximum standardized mean difference (d) of approximately .25, and variance ratios (vr) within the range of .5 to 2. The d value was close to the maximum of .25, although the vr was within the acceptable range (d = .24, vr = 1.46). Given the high standardized mean difference and the outlier seen in the propensity score plot, we removed the outlier ASD participant and conducted nearest neighbor and optimal matching methods with a revised sample of 24 children with ASD and 43 TD children.
Propensity scores distribution
# nearest neighbor (the same participant is identified using the optimal matching method)
plot(m_near_n68, type = "jitter", interactive = F)
Matched Treatment Units = Children with ASD; Matched Control Units = selected matches of TD children; Unmatched Control Units = remaining unmatched TD children. Propensity scores calculated using the nearest neighbor and optimal matching methods resulted in the same values. This plot demonstrates the distribution of propensity scores when including covariates of age and IQ for all 25 children with ASD and 43 TD children. We see a similar distribution of propensity scores for Matched Treatment Units and Matched Control units, ranging from scores of 0.2 to above 0.6. Among the Matched Treatment Units, there appears to be one outlier where a child with ASD was assigned a propensity score of approximately 0.8.
# participant with high propensity score in ASD group
m_near_n68_df %>% filter(distance > .7)## participant group age_yrs leiter_composite group.int distance weights
## 1 330 ASD 11.42 80 1 0.8027958 1Calculate standardized mean differences and variance ratios
# NOTE: both nearest neighbor and optimal matching methods have the same participants, so either method will provide the same scores. Scores from both methods are presented here for clarity.
# nearest neighbor
# a) create vectors with data in asd and td group
near_asd_n68 <- m_near_n68_df %>%
filter(group == "ASD")
near_td_n68 <- m_near_n68_df %>%
filter(group == "TD")
# b) standardized mean difference of propensity score
# (Mtreatment - Mcontrol) / SDtreatment; (Stuart, 2010; Rubin, 2001; Rosenbaum & Rubin, 1985 - The American Statistician)
(mean(near_asd_n68$distance) - mean(near_td_n68$distance)) / sd(near_asd_n68$distance)## [1] 0.2410058# c) variance ratio
(sd(near_asd_n68$distance)^2) / (sd(near_td_n68$distance)^2)## [1] 1.456683# optimal matching
# a) create vectors with data in asd and td group
opt_asd_n68 <- m_opt_n68_df %>%
filter(group == "ASD")
opt_td_n68 <- m_opt_n68_df %>%
filter(group == "TD")
# b) standardized mean difference of propensity score
(mean(opt_asd_n68$distance) - mean(opt_td_n68$distance)) / sd(opt_asd_n68$distance)## [1] 0.2410058# c) variance ratio
(sd(opt_asd_n68$distance)^2) / (sd(opt_td_n68$distance)^2)## [1] 1.456683
Step 3 with revised sample (24 ASD, 42 TD):
Conduct matching
The second iteration using the revised sample resulted in 24 children with ASD matched to 24 TD children. In this iteration, there was a difference between the two methods in the set of TD children selected for the matched group, which was a difference in one child.
Remove participant in ASD group with high propensity score
# create df that removes ASD participant
df_ps_n67 <- df2 %>%
select(participant, group, age_yrs, leiter_composite) %>%
mutate(group.int = ifelse(group == "ASD", 1, 0)) %>%
filter(participant != 330)
# check count of participants ASD = 24 and TD = 43
df_ps_n67 %>% count(group)## # A tibble: 2 x 2
## group n
## <chr> <int>
## 1 ASD 24
## 2 TD 43Conduct Nearest neighbor algorithm: n = 67
# apply matching
m_near_n67 <- matchit(group.int ~ age_yrs + leiter_composite,
data = df_ps_n67, method = "nearest")
# get df of matched participants and associated covariates
m_near_n67_df <- match.data(m_near_n67)Conduct Optimal matching algorithm: n = 67
# apply matching
m_opt_n67 <- matchit(group.int ~ age_yrs + leiter_composite,
data = df_ps_n67, method = "optimal")
# get df of matched participants and associated covariates
m_opt_n67_df <- match.data(m_opt_n67)Compare participants selected with nearest neighbor vs. optimal matching: n = 67
This side-by-side comparison demonstrates that nearest neighbor and optimal matching methods resulted in two different sets of 24 TD children selected as matches to the 24 children with ASD. Nearest neighbor matching includes participant 351, whereas optimal matching includes participant 323. This single difference in participants resulted in a stronger matched set of TD children to children with ASD using the optimal matching method versus the nearest neighbor method.
# extract list of ASD and TD participants selected from nearest neighbor and optimal matching methods; the list is presented in ascending order by participant id
# English = 300 series
# French = 400 series
# ASD participants = even numbers
# TD participants = odd numbers
m_near_n67_part <- m_near_n67_df %>%
dplyr::select(c(participant, group)) %>%
rename(m_near_participants = participant,
m_near_group = group) %>%
arrange(m_near_group)
m_opt_n67_part <- m_opt_n67_df %>%
dplyr::select(c(participant, group)) %>%
rename(m_opt_participants = participant,
m_opt_group = group) %>%
arrange(m_opt_group)
# compare participants selected from nearest neighbor vs. optimal matching methods
check_n67_df <- cbind(m_near_n67_part, m_opt_n67_part)
check_n67_df## m_near_participants m_near_group m_opt_participants m_opt_group
## 1 302 ASD 302 ASD
## 2 304 ASD 304 ASD
## 3 306 ASD 306 ASD
## 4 310 ASD 310 ASD
## 5 312 ASD 312 ASD
## 6 314 ASD 314 ASD
## 7 318 ASD 318 ASD
## 8 322 ASD 322 ASD
## 9 324 ASD 324 ASD
## 10 326 ASD 326 ASD
## 11 332 ASD 332 ASD
## 12 402 ASD 402 ASD
## 13 404 ASD 404 ASD
## 14 406 ASD 406 ASD
## 15 408 ASD 408 ASD
## 16 410 ASD 410 ASD
## 17 412 ASD 412 ASD
## 18 414 ASD 414 ASD
## 19 416 ASD 416 ASD
## 20 418 ASD 418 ASD
## 21 420 ASD 420 ASD
## 22 424 ASD 424 ASD
## 23 428 ASD 428 ASD
## 24 430 ASD 430 ASD
## 25 315 TD 315 TD
## 26 321 TD 321 TD
## 27 325 TD 323 TD
## 28 329 TD 325 TD
## 29 331 TD 329 TD
## 30 333 TD 331 TD
## 31 341 TD 333 TD
## 32 347 TD 341 TD
## 33 349 TD 347 TD
## 34 351 TD 349 TD
## 35 403 TD 403 TD
## 36 407 TD 407 TD
## 37 409 TD 409 TD
## 38 413 TD 413 TD
## 39 415 TD 415 TD
## 40 421 TD 421 TD
## 41 425 TD 425 TD
## 42 427 TD 427 TD
## 43 429 TD 429 TD
## 44 433 TD 433 TD
## 45 435 TD 435 TD
## 46 439 TD 439 TD
## 47 443 TD 443 TD
## 48 445 TD 445 TD
Step 4 with revised sample (24 ASD, 42 TD): Diagnose matching
The optimal method with 24 children ASD and 24 TD children established balanced groups on the desired covariates of age and IQ.
The final step in diagnosing groups is to determine how well groups are matched on each covariate included in the propensity score, as well as any other pre-existing variables that may be of interest in the study.
When examining propensity score distributions, there were no clear outliers with either method. However, an examination of standardized mean differences and variance ratios indicated that the optimal matching method was better than the nearest neighbor results (optimal: d = .14, vr = 1.05; nearest neighbor: d = .24, vr = 1.32). These findings suggest that the optimal matching method with 24 children per group resulted in better balanced groups versus the nearest neighbor method, as well as the matched groups in the first iteration with 25 ASD and 25 TD children.
After finalizing our matched samples using propensity scores, we evaluated our variables. Guidelines to evaluate well matched groups on each variable included examination of boxplots (where one would observe significant overlap when groups are well-balanced), p values > .5, Cohen’s d close to 0, and variance ratios close to 1 (5,13,14). Cohen’s d was calculated using the compute.es package with formulas in line with Kover and Atwood (2013). The use of Cohen’s d and variance ratios are recommended as alternatives to inferential statistics such as p values, due to difficulties with establishing equivalence with inferential statistics.
Propensity scores distribution
# nearest neighbor
plot(m_near_n67, type = "jitter", interactive = F)
# optimal matching
plot(m_opt_n67, type = "jitter", interactive = F)
Calculate standardized mean differences and variance ratios
# nearest neighbor
# a) create vectors with data in asd and td group
near_asd_n67 <- m_near_n67_df %>%
filter(group == "ASD")
near_td_n67 <- m_near_n67_df %>%
filter(group == "TD")
# b) standardized mean difference of propensity score
(mean(near_asd_n67$distance) - mean(near_td_n67$distance)) / sd(near_asd_n67$distance)## [1] 0.2417592# c) variance ratio
(sd(near_asd_n67$distance)^2) / (sd(near_td_n67$distance)^2)## [1] 1.324467# optimal matching
# a) create vectors with data in asd and td group
opt_asd_n67 <- m_opt_n67_df %>%
filter(group == "ASD")
opt_td_n67 <- m_opt_n67_df %>%
filter(group == "TD")
# b) standardized mean difference of propensity score
(mean(opt_asd_n67$distance) - mean(opt_td_n67$distance)) / sd(opt_asd_n67$distance)## [1] 0.1443307# c) variance ratio
(sd(opt_asd_n67$distance)^2) / (sd(opt_td_n67$distance)^2)## [1] 1.04522
Get list of matched participants and all covariates to report in study
As seen in the table and figure below, our final revised sample with the optimal matching method resulted in two successfully balanced groups according to criteria listed above on our covariates of interest, age and IQ. We next examined other variables not included in our propensity scores, but may be related to group diagnosis and/or performance on outcome measures.
The matched groups met the cut-off for p values > .5 on the ratio of English- to French-speaking children. On measures of sex, parental education level, and CELF-4 Word Association, groups were not significantly different (p’s between .136 and .461), but these values did not meet recommended matching cutoffs of p > .5.
As expected prior to selecting participants, groups were significantly different in their distribution on language measures of Recalling Sentences and Word Classes. Additionally, as expected due to diagnoses, groups were significantly different on social skills measures of the SCQ and VABS-II Socialization domain.
Lastly, we verified the distribution of children for randomized experimental factors (i.e., block order). As seen in the table below including descriptive statistics, the same proportion of children had both block orders.
# get list of participants
matchit_list <- m_opt_n67_df %>%
select(participant)
# create df with participants in matched groups
df_matchit <- df2 %>%
right_join(matchit_list, by = "participant") %>%
select(-c(keep, celf4_WC_rec_scaled, celf4_WC_exp_scaled)) %>%
mutate(group = factor(group), language = factor(language),
sex = factor(sex), parental_ed = factor(parental_ed),
block1 = factor(block1))
# check number of rows of participants
nrow(df_matchit)## [1] 48Calculate inferential statistics
# remove participant column to run compareGroups
df_matchit2 <- df_matchit %>%
select(-participant)
# compareGroups uses Fisher's exact test as a default when cell values < 5
compare_group <- compareGroups(group ~ . , data = df_matchit2)
compare_group_tbl <- createTable(compare_group, digits = 2, type = 1, show.n = TRUE)
print(compare_group_tbl, header.labels = c(p.overall = "p-value"))##
## --------Summary descriptives table by 'group'---------
##
## ______________________________________________________________
## ASD TD p-value N
## N=24 N=24
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## age_yrs 8.83 (1.26) 8.70 (1.12) 0.713 48
## language: 1.000 48
## English 45.83% 41.67%
## French 54.17% 58.33%
## sex: 0.461 48
## F 12.50% 25.00%
## M 87.50% 75.00%
## block1: 1.000 48
## order1 50.00% 54.17%
## order2 50.00% 45.83%
## parental_ed: 0.136 48
## 0 50.00% 25.00%
## 1 50.00% 75.00%
## leiter_composite 108.29 (12.65) 109.50 (13.24) 0.748 48
## celf4_RS_scaled 8.08 (4.16) 11.17 (2.18) 0.003 48
## celf4_WC_total_scaled 9.74 (3.74) 12.08 (3.06) 0.024 47
## celf4_WA_total 29.92 (15.01) 33.29 (11.17) 0.382 48
## scq_total 20.88 (5.83) 4.42 (2.62) <0.001 48
## social_standScore 76.83 (11.64) 110.00 (11.88) <0.001 48
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Calculate Cohen’s d and variance ratios
# create separate vectors for groups
asd <- df_matchit %>%
filter(group == "ASD")
td <- df_matchit %>%
filter(group == "TD")
# age_yrs - cohen's d and variance ratios
age_yrs_es <- mes(mean(asd$age_yrs), mean(td$age_yrs),
sd(asd$age_yrs), sd(td$age_yrs),
nrow(asd), nrow(td), verbose = F)
age_yrs <- age_yrs_es[, 4]
age_yrs_es_ci <- age_yrs_es[, 6:7]
age_yrs_es_vr <- (sd(asd$age_yrs)^2)/(sd(td$age_yrs)^2)
# celf4_RS_scaled - cohen's d and variance ratios
celf4_rs_es <- mes(mean(asd$celf4_RS_scaled), mean(td$celf4_RS_scaled),
sd(asd$celf4_RS_scaled), sd(td$celf4_RS_scaled),
nrow(asd), nrow(td), verbose = F)
celf4_rs <- celf4_rs_es[, 4]
celf4_rs_es_ci <- celf4_rs_es[, 6:7]
celf4_rs_es_vr <- (sd(asd$celf4_RS_scaled)^2)/(sd(td$celf4_RS_scaled)^2)
# celf4_WC_total_scaled - cohen's d and variance ratios
# removing child with no score for WC
asd2 <- asd %>%
filter(participant != 420)
celf4_wc_es <- mes(mean(asd2$celf4_WC_total_scaled), mean(td$celf4_WC_total_scaled),
sd(asd2$celf4_WC_total_scaled), sd(td$celf4_WC_total_scaled),
nrow(asd2), nrow(td), verbose = F)
celf4_wc <- celf4_wc_es[, 4]
celf4_wc_es_ci <- celf4_wc_es[, 6:7]
celf4_wc_es_vr <- (sd(asd2$celf4_WC_total_scaled)^2)/(sd(td$celf4_WC_total_scaled)^2)
# celf4_WA_total - cohen's d and variance ratios
celf4_wa_es <- mes(mean(asd$celf4_WA_total), mean(td$celf4_WA_total),
sd(asd$celf4_WA_total), sd(td$celf4_WA_total),
nrow(asd), nrow(td), verbose = F)
celf4_wa <- celf4_wa_es[, 4]
celf4_wa_es_ci <- celf4_wa_es[, 6:7]
celf4_wa_es_vr <- (sd(asd$celf4_WA_total)^2)/(sd(td$celf4_WA_total)^2)
# leiter_composite - cohen's d and variance ratios
leiter_es <- mes(mean(asd$leiter_composite), mean(td$leiter_composite),
sd(asd$leiter_composite), sd(td$leiter_composite),
nrow(asd), nrow(td), verbose = F)
leiter <- leiter_es[, 4]
leiter_es_ci <- leiter_es[, 6:7]
leiter_es_vr <- (sd(asd$leiter_composite)^2)/(sd(td$leiter_composite)^2)
# scq_total - cohen's d and variance ratios
scq_es <- mes(mean(asd$scq_total), mean(td$scq_total),
sd(asd$scq_total), sd(td$scq_total),
nrow(asd), nrow(td), verbose = F)
scq <- scq_es[, 4]
scq_es_ci <- scq_es[, 6:7]
scq_es_vr <- (sd(asd$scq_total)^2)/(sd(td$scq_total)^2)
# social_standScore - cohen's d and variance ratios
vineland_es <- mes(mean(asd$social_standScore), mean(td$social_standScore),
sd(asd$social_standScore), sd(td$social_standScore),
nrow(asd), nrow(td), verbose = F)
vineland <- vineland_es[, 4]
vineland_es_ci <- vineland_es[, 6:7]
vineland_es_vr <- (sd(asd$social_standScore)^2)/(sd(td$social_standScore)^2)
es_d <- rbind(age_yrs, celf4_rs, celf4_wc, celf4_wa, leiter, scq, vineland)
es_ci <- rbind(age_yrs_es_ci, celf4_rs_es_ci, celf4_wc_es_ci, celf4_wa_es_ci, leiter_es_ci, scq_es_ci, vineland_es_ci)
es_vr <- rbind(age_yrs_es_vr, celf4_rs_es_vr, celf4_wc_es_vr, celf4_wa_es_vr, leiter_es_vr, scq_es_vr, vineland_es_vr)
effect_sizes <- cbind(es_d, es_ci, es_vr)
effect_sizes %>%
rename( d = 1, d_CI_lower_bound = 2, d_CI_upper_bound = 3, vr = 4)## d d_CI_lower_bound d_CI_upper_bound vr
## age_yrs 0.11 -0.46 0.67 1.2742083
## celf4_rs -0.93 -1.52 -0.33 3.6387195
## celf4_wc -0.69 -1.28 -0.10 1.4939979
## celf4_wa -0.26 -0.82 0.31 1.8049142
## leiter -0.09 -0.66 0.47 0.9119877
## scq 3.64 2.72 4.57 4.9458817
## vineland -2.82 -3.62 -2.02 0.9591544Create violin plots to examine continuous variables
# create violin plots
p1 <- ggplot(df_matchit, aes(group, age_yrs, fill = group)) +
geom_violin(alpha = .7) +
geom_jitter(alpha = .3) +
theme(legend.position = "none", axis.text=element_text(size=15)) +
labs(y = "Age (yrs)", x = "")
p2 <- ggplot(df_matchit, aes(group, leiter_composite, fill = group)) +
geom_violin(alpha = .7) +
geom_jitter(alpha = .3) +
theme(legend.position = "none", axis.text=element_text(size=15)) +
labs(y = "Nonverbal IQ (Leiter)", x = "")
p3 <- ggplot(df_matchit, aes(group, celf4_WA_total, fill = group)) +
geom_violin(alpha = .7) +
geom_jitter(alpha = .3) +
theme(legend.position = "none", axis.text=element_text(size=15)) +
labs(y = "CELF-4 Word Associations", x = "")
p4 <- ggplot(df_matchit, aes(group, celf4_RS_scaled, fill = group)) +
geom_violin(alpha = .7) +
geom_jitter(alpha = .3) +
theme(legend.position = "none", axis.text=element_text(size=15)) +
labs(y = "CELF-4 Recalling Sentences", x = "")
p5 <- ggplot(df_matchit, aes(group, celf4_WC_total_scaled, fill = group)) +
geom_violin(alpha = .7) +
geom_jitter(alpha = .3) +
theme(legend.position = "none", axis.text=element_text(size=15)) +
labs(y = "CELF-4 Word Classes", x = "")
p6 <- ggplot(df_matchit, aes(group, scq_total, fill = group)) +
geom_violin(alpha = .7) +
geom_jitter(alpha = .3) +
theme(legend.position = "none", axis.text=element_text(size=15)) +
labs(y = "Social Communication Questionnaire", x = "")
p7 <- ggplot(df_matchit, aes(group, social_standScore, fill = group)) +
geom_violin(alpha = .7) +
geom_jitter(alpha = .3) +
theme(legend.position = "none", axis.text=element_text(size=15)) +
labs(y = "Vineland Socialization Domain", x = "")
# print violin plots
grid.arrange(p1, p2, p3, p4, p5, p6, p7, nrow = 2)
Points represent observations per participant. For age and nonverbal IQ, matching was achieved according to criteria of p > .5, Cohen’s d close to 0, and variance ratios close to 1. CELF-4 Word Associations did not meet criteria of p > .5, but between groups the distribution on this variable appears similar between groups. Groups are significantly different on other language measures of CELF-4 Recalling Sentences and Word Classes, as well as the Social Communication Questionnaire (SCQ) and the Vineland Socialization Domain.
Examine matched groups when including age, nonverbal IQ, and language (Recalling Sentences)
To examine the consequences of matching on age, nonverbal IQ and language, we conducted the nearest neighbor and optimal matching methods with all three variables. We chose the Recalling Sentences subtest to represent our language variable, because on Word Classes one child with ASD was unable to complete the measure and on Word Associations the distribution between groups in the full sample appeared similar (see above matrix).
With either the nearest neighbor or optimal matching method, including all three variables of age, nonverbal IQ and language resulted in different sets of TD children as potential matches for the 25 children with ASD. Matched groups with both matching methods that used three covariates resulted in less well-balanced matches than with just two covariates of age and IQ.
For example, as seen in the figure below, an examination of the distribution of propensity scores with the optimal matching method demonstrated that 8 children with ASD were outside the range of propensity scores relative to the rest of the children with ASD and TD children. The table below describes the propensity scores and each covariate when including age, nonverbal IQ, and CELF-4 Recalling Sentences. Only the variable of age meets the desired cutoff of p > .5, and the distribution of CELF-4 Recalling Sentences still appears substantially different between both groups (propensity score Cohen’s d > .5, propensity score variance ratios > 3, and ps < .5 for two of three variables).
Due to the poor balancing when including all three proposed covariates, this evidence supports balancing on two covariates of age and nonverbal IQ to retain as many children in the sample as possible. Additionally, because language abilities of children with ASD were not categorically poorer across all three measures relative to TD children, it is unclear on which language measure to match when intercorrelations between language measures ranged widely across the full sample (rs = .28 - .72), within ASD (rs = .34 - .84) and within TD children (rs = .09 - .47).
Therefore, matching groups on two covariates of age and nonverbal IQ appear to be both theoretically supported based on prior studies and empirically supported by the current evidence with our sample.
Select covariates
df_ps_n68_3cov <- df2 %>%
select(participant, group, age_yrs, leiter_composite, celf4_RS_scaled) %>%
mutate(group.int = ifelse(group == "ASD", 1, 0)) # group must be an integer for MatchItConduct Nearest neighbor algorithm: n = 68 including age, IQ, and language (Recalling Sentences)
# apply matching
m_near_n68_3cov <- matchit(group.int ~ age_yrs + leiter_composite + celf4_RS_scaled,
data = df_ps_n68_3cov, method = "nearest")
# obtain list of matched participants
m_near_n68_3cov_df <- match.data(m_near_n68_3cov)
m_near_n68_3cov_df$participant## [1] 302 304 306 307 310 311 312 313 314 315 317 318 319 321 322 323 324 326 330
## [20] 331 332 333 339 341 347 351 401 402 403 404 406 408 410 411 412 414 416 417
## [39] 418 420 421 423 424 425 428 430 431 437 439 445Diagnose matching - View propensity scores
# view propensity scores
plot(m_near_n68_3cov, type = "jitter")
Matched Treatment Units = Children with ASD; Matched Control Units = selected matches of TD children; Unmatched Control Units = remaining unmatched TD children. This plot depicts the distribution of propensity scores when including covariates of age, IQ, and CELF-4 Recalling Sentences for all 25 children with ASD and 43 TD children. Propensity scores were calculated using the optimal matching method. There are 8 children with ASD who appear to be outliers relative to the propensity scores for TD children.
## [1] "To identify the units, use first mouse button; to stop, use second."## integer(0)Diagnose matching - Standardized mean difference and Variance ratio
# standardized mean difference
# a) create vectors with data in asd and td group
near_asd_n68_3cov <- m_near_n68_3cov_df %>%
filter(group == "ASD")
near_td_n68_3cov <- m_near_n68_3cov_df %>%
filter(group == "TD")
# b) standardized mean difference of propensity score
# (Mtreatment - Mcontrol) / SDtreatment; (Stuart, 2010; Rubin, 2001; Rosenbaum & Rubin, 1985 - The American Statistician)
(mean(near_asd_n68_3cov$distance) - mean(near_td_n68_3cov$distance)) / sd(near_asd_n68_3cov$distance)## [1] 0.5498751# variance ratio
(sd(near_asd_n68_3cov$distance)^2) / (sd(near_td_n68_3cov$distance)^2)## [1] 3.572499Diagnose matching - Descriptive and inferential statistics on age, IQ, and language
# subset only age and IQ
age_IQ_language_3cov <- m_near_n68_3cov_df %>%
select(group, age_yrs, leiter_composite, celf4_RS_scaled)
# compare groups
age_IQ_language_compare <- compareGroups(group ~ age_yrs + leiter_composite + celf4_RS_scaled, data = age_IQ_language_3cov)
age_IQ_language_compare_tbl <- createTable(age_IQ_language_compare, digits = 2, type = 1, show.n = TRUE)
print(age_IQ_language_compare_tbl, header.labels = c(p.overall = "p-value"))##
## --------Summary descriptives table by 'group'---------
##
## _________________________________________________________
## ASD TD p-value N
## N=25 N=25
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## age_yrs 8.93 (1.34) 8.73 (1.24) 0.587 50
## leiter_composite 107.16 (13.61) 113.00 (14.89) 0.154 50
## celf4_RS_scaled 7.84 (4.25) 10.00 (2.10) 0.029 50
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Conduct Optimal matching algorithm: n = 68 including age, IQ, and language (Recalling Sentences)
# apply matching
m_opt_n68_3cov <- matchit(group.int ~ age_yrs + leiter_composite + celf4_RS_scaled,
data = df_ps_n68_3cov, method = "optimal")
# obtain list of matched participants
m_opt_n68_3cov_df <- match.data(m_opt_n68_3cov)
m_opt_n68_3cov_df$participant## [1] 302 304 306 310 311 312 313 314 315 317 318 319 321 322 323 324 326 330 331
## [20] 332 333 339 341 347 349 351 401 402 403 404 406 408 410 411 412 414 416 418
## [39] 420 421 423 424 425 428 430 431 433 437 439 445Diagnose matching - View propensity scores
# view propensity scores
plot(m_opt_n68_3cov, type = "jitter")
## [1] "To identify the units, use first mouse button; to stop, use second."## integer(0)Diagnose matching - Standardized mean difference and Variance ratio
# standardized mean difference
# a) create vectors with data in asd and td group
opt_asd_n68_3cov <- m_opt_n68_3cov_df %>%
filter(group == "ASD")
opt_td_n68_3cov <- m_opt_n68_3cov_df %>%
filter(group == "TD")
# b) standardized mean difference of propensity score
# (Mtreatment - Mcontrol) / SDtreatment; (Stuart, 2010; Rubin, 2001; Rosenbaum & Rubin, 1985 - The American Statistician)
(mean(opt_asd_n68_3cov$distance) - mean(opt_td_n68_3cov$distance)) / sd(opt_asd_n68_3cov$distance)## [1] 0.5427627# variance ratio
(sd(opt_asd_n68_3cov$distance)^2) / (sd(opt_td_n68_3cov$distance)^2)## [1] 3.717604Diagnose matching - Descriptive and inferential statistics on age, IQ, and language
# subset only age and IQ
age_IQ_language_3cov <- m_opt_n68_3cov_df %>%
select(group, age_yrs, leiter_composite, celf4_RS_scaled)
# compare groups
age_IQ_language_compare <- compareGroups(group ~ age_yrs + leiter_composite + celf4_RS_scaled, data = age_IQ_language_3cov)
age_IQ_language_compare_tbl <- createTable(age_IQ_language_compare, digits = 2, type = 1, show.n = TRUE)
print(age_IQ_language_compare_tbl, header.labels = c(p.overall = "p-value"))##
## --------Summary descriptives table by 'group'---------
##
## _________________________________________________________
## ASD TD p-value N
## N=25 N=25
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## age_yrs 8.93 (1.34) 8.69 (1.31) 0.520 50
## leiter_composite 107.16 (13.61) 113.80 (15.10) 0.109 50
## celf4_RS_scaled 7.84 (4.25) 9.88 (1.94) 0.036 50
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Diagnose matching - Descriptive and inferential statistics on age, IQ, and language - removing 8 ASD outliers
# subset only age and IQ
age_IQ_language_3cov_rm8asd <- m_opt_n68_3cov_df %>%
filter(distance < .7) %>%
select(group, age_yrs, leiter_composite, celf4_RS_scaled)
# compare groups
age_IQ_language_compare_rm8asd <- compareGroups(group ~ age_yrs + leiter_composite + celf4_RS_scaled, data = age_IQ_language_3cov_rm8asd)
age_IQ_language_compare_rm8asd_tbl <- createTable(age_IQ_language_compare_rm8asd, digits = 2, type = 1, show.n = TRUE)
print(age_IQ_language_compare_rm8asd_tbl, header.labels = c(p.overall = "p-value"))##
## --------Summary descriptives table by 'group'---------
##
## _________________________________________________________
## ASD TD p-value N
## N=17 N=25
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
## age_yrs 8.65 (1.02) 8.69 (1.31) 0.903 42
## leiter_composite 112.53 (11.62) 113.80 (15.10) 0.760 42
## celf4_RS_scaled 10.12 (2.96) 9.88 (1.94) 0.773 42
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
Session Info
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] gridExtra_2.3 compareGroups_4.2.0 SNPassoc_1.9-2
## [4] mvtnorm_1.0-11 survival_3.1-8 haplo.stats_1.7.9
## [7] compute.es_0.2-5 GGally_1.4.0 MatchIt_3.0.2
## [10] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3
## [13] purrr_0.3.3 readr_1.3.1 tidyr_1.0.0
## [16] tibble_2.1.3 ggplot2_3.2.1 tidyverse_1.3.0
## [19] knitr_1.26
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.0-10 minqa_1.2.4 colorspace_1.4-1
## [4] flextable_0.5.7 htmlTable_1.13.2 base64enc_0.1-3
## [7] fs_1.3.1 rstudioapi_0.11 mice_3.7.0
## [10] farver_2.0.1 MatrixModels_0.4-1 fansi_0.4.0
## [13] lubridate_1.7.4 xml2_1.2.2 codetools_0.2-16
## [16] splines_3.6.1 zeallot_0.1.0 Formula_1.2-3
## [19] jsonlite_1.6.1 nloptr_1.2.1 broom_0.5.2
## [22] cluster_2.1.0 dbplyr_1.4.2 compiler_3.6.1
## [25] httr_1.4.1 backports_1.1.5 assertthat_0.2.1
## [28] Matrix_1.2-17 lazyeval_0.2.2 cli_2.0.2
## [31] acepack_1.4.1 htmltools_0.4.0 quantreg_5.54
## [34] tools_3.6.1 gtable_0.3.0 glue_1.4.0
## [37] reshape2_1.4.3 optmatch_0.9-13 Rcpp_1.0.3
## [40] cellranger_1.1.0 vctrs_0.2.0 writexl_1.2
## [43] nlme_3.1-140 xfun_0.11 lme4_1.1-21
## [46] rvest_0.3.5 lifecycle_0.1.0 polspline_1.1.17
## [49] pan_1.6 MASS_7.3-51.4 zoo_1.8-6
## [52] scales_1.1.0 hms_0.5.2 sandwich_2.5-1
## [55] SparseM_1.78 RColorBrewer_1.1-2 HardyWeinberg_1.6.3
## [58] yaml_2.2.0 gdtools_0.2.1 epitools_0.5-10
## [61] rms_5.1-4 rpart_4.1-15 reshape_0.8.8
## [64] latticeExtra_0.6-28 stringi_1.4.3 highr_0.8
## [67] checkmate_1.9.4 zip_2.0.4 boot_1.3-22
## [70] truncnorm_1.0-8 chron_2.3-55 systemfonts_0.1.1
## [73] rlang_0.4.5 pkgconfig_2.0.3 Rsolnp_1.16
## [76] evaluate_0.14 lattice_0.20-38 labeling_0.3
## [79] htmlwidgets_1.5.1 tidyselect_0.2.5 plyr_1.8.4
## [82] magrittr_1.5 R6_2.4.1 RItools_0.1-17
## [85] generics_0.0.2 Hmisc_4.3-0 mitml_0.3-7
## [88] multcomp_1.4-12 DBI_1.0.0 pillar_1.4.2
## [91] haven_2.2.0 foreign_0.8-71 withr_2.1.2
## [94] abind_1.4-5 nnet_7.3-12 modelr_0.1.5
## [97] crayon_1.3.4 jomo_2.6-10 utf8_1.1.4
## [100] uuid_0.1-2 officer_0.3.6 rmarkdown_1.17
## [103] grid_3.6.1 readxl_1.3.1 data.table_1.12.6
## [106] svd_0.5 webshot_0.5.2 reprex_0.3.0
## [109] digest_0.6.25 xtable_1.8-4 munsell_0.5.0
## [112] viridisLite_0.3.0 kableExtra_1.1.0