--- title: "Data Analysis of Predictive coding in ASD" author: "Z. Shi, L. Theisinger, F. Allenmark, R. Pistorius, H. Müller, C. Falter-Wagner" output: html_notebook --- # Data Structure 1. `/experiments`: Experimental codes and instructions This sub-folder contains Matlab codes and instructions for the duration reproduction task. The sequences of the duration reproductions are stored in the sub-folder `/experiments/seqs`. Those sequences were used for matched participants. 2. `/data`: raw data files - `rawdata.csv`: Raw reproduction trials - `outliers.csv`: those produced almost flat reproduction 3. `/figures`: store figures used in the paper. # Data Analysis ## 1. load raw data ```{r packages, message=FALSE, warning=FALSE, include=FALSE} # load packages library(tidyverse) library(ez) library(cowplot) library(BayesFactor) # in case need Bayes factor analysis library(bayestestR) library(rstatix) # using tidyverse friendly statistics library(GGally) # ---- read data and preparation ----- rawdata = read_csv('./data/rawdata.csv') rawdata$group = toupper(rawdata$group) # flag saveFig = FALSE ``` Show the raw trial structure: ```{r raw trial structure} glimpse(rawdata) ``` In the raw trials, there are several important columns which are relevant for further analyses. - `Duration`: The test durations generated by computer. 'pdur' is the actual presented durations by the computer. There were some fluctuations, but within 5 ms (within 1 refresh frame). - `sub`, `group`, `sequence` are anonymized subject number, group and the duration sequence used in the experiment. - `Reproduction`, `rep_err`: the reproduced durations and the reproduction errors compared to the given duration. - `itd`, `preDuration`: the inter-trial difference in a given trial, and the duration used in the previous trial. These were used in the sequential effect analysis. ### Participants ```{r participants} # Total participants print(length (unique(rawdata$sub))) # total sequences (pairs) print(length(unique(rawdata$sequence))) ``` ## The sampled durations and sequences Let's first illustrate the sequence and the distribution of the sample durations. ```{r sequence} # ---- illustrate one sequence ----- fig11 = ggplot(rawdata, aes(Duration)) + geom_histogram(binwidth = 0.1, fill = I("white"), col = I('black')) + theme_classic() + xlab('Duration (secs)') # a typical sequence sub1 = rawdata %>% filter(sub == 'ara27') fig12 = ggplot(sub1, aes(trlNo, Duration, color = Volatility)) + geom_line() + xlab("Trial Sequence") + theme_classic()+ theme(legend.position = 'top') # histogram of the typical sequence fig13 = ggplot(sub1, aes(y=Duration, color = Volatility)) + geom_density() + theme_classic() + theme(legend.position = 'top') fig13 fig1 = plot_grid(fig12,fig13, rel_widths =c(3,1)) fig1 if (saveFig){ ggsave("figures/fig_sequence.png", fig1, width=4, height=3.5) ggsave("figures/fig_sequence.pdf", fig1, width=4, height=3.5) } ``` ## Explorative data analysis and outlier detection We first estimate two key signatures - the central tendency index (ci) and the sequential dependence index (si). ```{r linear_models} # using error for regression, -slope is the central tendency index. ct_model <- function(df){ lm(rep_err ~ Duration, data = df) } # sequential effect using the duration from the previous trial seq_model <- function(df){ lm(rep_err ~ preDuration, data = df) } # exclude extreme trials (beyond [1/3D, 3D]) vdata = rawdata %>% filter(Reproduction > Duration/3, Reproduction < 3*Duration) # calculate slope for the central tendency as well as the sequential dependence. slopes_ct <- vdata %>% group_by(sub, sequence, Volatility, group, Order) %>% nest() %>% # nested data mutate(model = map(data, ct_model)) %>% # linear regression mutate(slope = map(model, broom::tidy)) %>% # get estimates out unnest(slope, .drop = TRUE) %>% # remove raw data select(-std.error,-statistic, -p.value) %>% # remove unnecessary columns spread(term, estimate) %>% # spread estimates rename(intercept_ct = `(Intercept)`, slope = Duration) %>% select(-data, -model) #sequential effect slopes_seq <- vdata %>% filter(!is.na(preDuration)) %>% group_by(sub, sequence, Volatility, group, Order) %>% nest() %>% # nested data mutate(model = map(data, seq_model)) %>% # linear regression mutate(slope = map(model, broom::tidy)) %>% # get estimates out unnest(slope, .drop = TRUE) %>% # remove raw data select(-std.error,-statistic, -p.value) %>% # remove unnecessary columns spread(term, estimate) %>% # spread estimates rename(intercept_seq = `(Intercept)`, seq_slope = preDuration) %>% select(-data, -model) # merge two tables together slopes_all = left_join(slopes_ct, slopes_seq, by = c("sub", "sequence", "Volatility","group","Order")) %>% mutate(ci = -slope, si = seq_slope) #central tendency index # estimated general biases (over-/under-estimates) # individual mean interval (middle point): most of them were around 1 by design mInterval = vdata %>% group_by(sub) %>% summarise(mDur = mean(Duration)) # join the mean interval, and estimate the general bias slopes = slopes_all %>% left_join(., mInterval, by = c('sub')) %>% mutate(gBias = (intercept_ct + slope * mDur)*1000) # change factor order for plotting slopes$Volatility = factor(slopes$Volatility) slopes$Order = factor(slopes$Order) slopes$group = factor(slopes$group) head(slopes) ``` Let's plot the histograms of the central tendency index and the sequential dependence index, which show some outliers. We visualize the outliers using 3-sigma rule. The 3-sigma rule only rule out 0.3% of the population if we assume the population is normal. ```{r} ci_3sig = mean(slopes$ci) + c(-1,1)*3*sd(slopes$ci) si_3sig = mean(slopes$si) + c(-1,1)*3*sd(slopes$si) hist1 = ggplot(slopes, aes(x = ci, y = ..density..)) + geom_histogram(colour = 1, fill = 'white', bins = 20) + geom_vline(xintercept = ci_3sig, linetype = 'dashed', color = 'red') + theme_classic() + xlab('Central Tendency Index') hist2 = ggplot(slopes, aes(x = si, y = ..density..)) + geom_histogram(colour = 1, fill = 'white', bins = 20) + geom_vline(xintercept = ci_3sig, linetype = 'dashed', color = 'red') + theme_classic() + xlab('Sequential Dependence Index') hist_fig = plot_grid(hist1,hist2, nrow = 2) hist_fig if (saveFig){ ggsave("figures/hist_fig.png", hist_fig, width=4, height=3.5) ggsave("figures/hist_fig.pdf", fig1, width=4, height=3.5) } ``` Note, the above outlier detection method is agnostic to the groups and conditions! Let's find out those outliers and exclude by sequences for further analyses (given that ASD and TD groups were paired). ```{r} slopes %>% ungroup() %>% filter(ci > ci_3sig[2] | ci < ci_3sig[1] | si > si_3sig[2] | si < si_3sig[1]) %>% select(sequence) -> outlier_seq outlier_seq$sequence ``` Now let's visualize the outliers together with their matched pairs. ```{r outliers} # ---- plot outliers ---- fig_outlier = rawdata %>% ungroup() %>% filter(sequence %in% outlier_seq$sequence) %>% filter(Reproduction > Duration/3, Reproduction < 3*Duration) %>% group_by(sub, group,sequence, Volatility, Duration) %>% summarise(mRep = mean(Reproduction), sdr=sd(Reproduction), n = n(), se = sd(Reproduction)/sqrt(n-1)) %>% filter(n>5) %>% #approx. for linear regression without averaging first ggplot(aes(Duration, mRep, color = group, group = interaction(group, Volatility), shape = Volatility, linetype = Volatility)) + geom_point() + geom_smooth(method = 'lm', se = FALSE) + facet_wrap(~sequence, ncol = 4) + theme_classic() + theme(legend.pos = 'bottom') + geom_abline(slope = 1, linetype = 3) + xlab('Duration (Secs)') + ylab(' Reproduction (Secs)') + theme(strip.background = element_blank(), strip.text.x = element_blank()) fig_outlier if (saveFig){ ggsave("figures/fig_outlier.png", fig_outlier, width=5, height=5) ggsave("figures/fig_outlier.pdf", fig_outlier, width=5, height=5) } ``` Mean slopes and related statistics: ```{r} slopes %>% filter(sequence %in% outlier_seq$sequence) %>% arrange(group, Volatility) vslopes = slopes %>% filter(!(sequence %in% outlier_seq$sequence)) # outliers oslopes = slopes %>% filter(sequence %in% outlier_seq$sequence) # ANOVA ezANOVA(data = oslopes, dv = ci, wid = sub, within = Volatility, between = group) ``` ### Outliers and sequential dependence Visualize the strong sequential dependence. ```{r} fig_sdep31 = rawdata %>% filter(sequence %in% c(31) ) %>% # the extreme sequential dependence pair mutate_at("preDuration", round, 1) %>% group_by(sequence, group, preDuration) %>% summarise(rep_err = mean(rep_err)) %>% ggplot(aes(preDuration, rep_err, color = group) ) + geom_point() + geom_smooth(method = 'lm') + theme_classic() + geom_abline(slope = 0, linetype = 2) + # facet_wrap(~sequence) + xlab('Duration trial n-1') + ylab('Repr. Error (Secs)') + theme(strip.background = element_blank(), strip.text.x = element_blank(), legend.position = 'bottom') fig_sdep31 fig_o = plot_grid(fig_outlier, fig_sdep31, nrow = 1, labels = c('a','b'), rel_widths = c(2.4,1)) if (saveFig){ ggsave("figures/fig_o.png", fig_o, width=7, height=3) ggsave("figures/fig_o.pdf", fig_o, width=7, height=3) } ``` ## Typical reproduction performance And here are the two typical participants produced errors from sequence 12: ```{r typical participants} # ---- reproduction figure - an example ----- # mean analysis mrep = rawdata %>% ungroup() %>% filter(!(sequence %in% outlier_seq$sequence)) %>% filter(Reproduction > Duration/3, Reproduction < 3*Duration) %>% group_by(sub, group,Order, sequence, Volatility, Duration) %>% summarise(mRep = mean(Reproduction), sdr=sd(Reproduction), n = n(), se = sd(Reproduction)/sqrt(n-1)) %>% filter(n>3) # individual examples # asd fig_asd = mrep%>% filter( sequence == '11', group == 'ASD') %>% ggplot(aes(Duration, mRep, color = Volatility, group = Volatility, shape = Volatility)) + geom_point(size = 2) + #geom_line(aes(linetype = Volatility)) + geom_smooth(method = 'lm', aes(fill = Volatility), se = FALSE) + geom_errorbar(aes(ymin = mRep - se, ymax = mRep +se), width = 0.05) + theme_classic() + theme(strip.background = element_blank()) + geom_abline(slope = 1, linetype = 2) + theme(legend.position = 'none', plot.margin = unit(c(0,0,0,0),'cm')) + ylab('') + xlab('') #td fig_td = mrep%>% filter( sequence == '11', group == 'TD') %>% ggplot(aes(Duration, mRep, color = Volatility, group = Volatility, shape = Volatility)) + geom_point(size = 2) + geom_smooth(method = 'lm', aes(fill = Volatility), se = FALSE) + geom_errorbar(aes(ymin = mRep - se, ymax = mRep +se), width = 0.05) + theme_classic() + theme(strip.background = element_blank()) + geom_abline(slope = 1, linetype = 2) + theme(legend.position = 'none', plot.margin = unit(c(0,0,0,0),'cm')) + ylab('') + xlab('') # group mean reproductions fig_repd_gr = mrep%>% ggplot(aes(Duration, mRep, color = Volatility, group = Volatility, shape = Volatility)) + geom_point(alpha = 0.2) + geom_smooth(method = 'lm', se = FALSE, aes(fill = Volatility)) + theme_classic() + theme(strip.background = element_blank()) + geom_abline(slope = 1, linetype = 2) + facet_wrap(~group) + theme(legend.position = c(0.1,0.85)) + xlab('Duration (Secs)') + ylab('Reproduction (Secs)') # plot individual example as inset. fig_mrep = ggdraw() + draw_plot(fig_repd_gr) + draw_plot(fig_asd, x = 0.3, y = 0.13, width = .2, height = .3) + draw_plot(fig_td, x = 0.75, y = 0.13, width = .2, height = .3) fig_mrep if(saveFig){ ggsave("figures/fig_reproduction.png", fig_mrep, width=7, height=3.5) ggsave("figures/fig_reproduction.pdf", fig_mrep, width=7, height=3.5) } ``` As we can see from the patterns above. The ASD participant produced relative flat errors, while the TD participant showed a strong central tendency effect (shorts being overestimated and longs being underestimated). ### Visulize CTE and sequential dependence Let's visualize the biases (central tendency and serial dependence) for two groups (ASD vs. TD) ```{r plot_biases} pd = position_dodge(width = 0.05) # plot CTI and SI together fig_biases = vslopes%>% group_by(group, Volatility, Order) %>% summarise(msi = mean(si), n = n(), se_si = sd(si)/sqrt(n), mci = mean(ci), se_ci = sd(ci)/sqrt(n)) %>% ggplot(aes(msi, mci, color = group, shape = Volatility, group=interaction(Order, group))) + geom_hline(yintercept = 0, color = 'gray', linetype = 'dashed') + geom_vline(xintercept = 0, color = 'gray', linetype = 'dashed') + geom_point(size = 2) + geom_line(aes(linetype = Order)) + geom_errorbar(aes(ymin = mci - se_ci, ymax = mci + se_ci), width = 0.01 ) + geom_errorbarh(aes(xmin = msi - se_si, xmax = msi + se_si), height = 0.01) + xlab('Serial Dependence') + ylab('Central Tendency') + scale_y_continuous(labels = scales::percent) + scale_x_continuous(labels = scales::percent) + guides(color = guide_legend(title = 'Group'), linetype = guide_legend(title = 'Volatility Order'), ) + theme_classic() + theme(legend.position = 'bottom') fig_biases if (saveFig){ ggsave("figures/fig_biases.png", fig_biases, width=5, height=4) ggsave("figures/fig_biases.pdf", fig_biases, width=5, height=4) } ``` By visual inspection, individuals with ASD exhibited less central tendency relative to their matched TD controls (the red lines below the cyan lines in the above figure), while the local serial dependence was relatively comparable between two groups except in the low volatility condition when that condition started first. Interestingly, the central tendency and serial dependence were similar for both groups when the high-volatility session started first (the solid lines), while they differed when the low-volatility session started first (the dashed lines). Given that the main difference was shown in CTE. We further plot the mean CTEs. ```{r} #pd = position_dodge(width = 0.5) # separate plots for appendix fig_cti = vslopes %>% group_by(group, Volatility, Order) %>% summarise(mci = mean(ci), n = n(), se = sd(ci)/sqrt(n)) %>% ggplot(aes(Volatility, mci, shape = Order, color = Order, group = Order)) + geom_line() + geom_point(size= 2) + #geom_bar(stat = 'identity', position = pd, width = 0.5) + facet_wrap(~group)+ geom_errorbar(aes(ymin = mci - se, ymax = mci + se), width = 0.2) + theme_classic() + theme(legend.position = 'bottom', strip.background = element_blank()) + xlab('Volatility') + ylab('CTI') + scale_y_continuous(labels = scales::percent) + guides(color = guide_legend(title = 'Order'), fill = guide_legend(title = 'Order')) fig_cti ``` Let's do the sequential dependence as well. ```{r} pd = position_dodge(width = 0.5) # separate plots for appendix fig_sdi = vslopes %>% group_by(group, Volatility, Order) %>% summarise(msi = mean(si), n = n(), se = sd(si)/sqrt(n)) %>% ggplot(aes(Volatility, msi, shape = Order, color = Order, group = Order)) + geom_line() + geom_point(size = 2) + #geom_bar(stat = 'identity', position = pd, width = 0.5) + facet_wrap(~group)+ geom_errorbar(aes(ymin = msi - se, ymax = msi + se), width = 0.2) + theme_classic() + theme(legend.position = 'bottom', strip.background = element_blank()) + xlab('Volatility') + ylab('SDI') + scale_y_continuous(labels = scales::percent) + guides(color = guide_legend(title = 'Order'), fill = guide_legend(title = 'Order')) fig_sdi fig3 = plot_grid(fig_cti, fig_sdi, nrow = 1, labels = c("a","b")) fig3 if (saveFig){ ggsave('figures/fig3.pdf',fig3, width = 7, height = 3.5) ggsave('figures/fig3.png',fig3, width = 7, height = 3.5) } ``` ### Statistics 1. Central tendency effect Average CTIs ```{r} vslopes %>% group_by(group) %>% summarise(mcti = mean(ci)) vslopes %>% group_by(Volatility) %>% summarise(mcti = mean(ci)) vslopes %>% group_by(Order) %>% summarise(mcti = mean(ci)) # calculate mean elevated CTIs between HV First vs. LV First vslopes %>% group_by(group, Volatility, Order) %>% summarise(mcti = mean(ci)) %>% pivot_wider(names_from = Order, values_from = mcti) %>% mutate(dCTI = `HV First` - `LV First`) %>% group_by(group) %>% summarise(md = mean(dCTI)) ``` ```{r ANOVAs_cti} vslopes = as.data.frame(vslopes) # required by rstatix # ---- central tendency index---- # repeated measures ANOVA on central tendency index anova1 = anova_test(data = vslopes, dv = ci, wid = sub, within = Volatility, between = c(group, Order)) anova1 ## separate for Volatility order anova1a = anova_test(data = vslopes %>% filter(group == 'ASD'), dv = ci, wid = sub, within = Volatility, between = Order) anova1a anova1b = anova_test(data = vslopes %>% filter(group == 'TD'), dv = ci, wid = sub, within = Volatility, between = Order) anova1b ``` 2. Bayes factor analysis for pair-wise comparison The above analysis showed that the main difference came from the central tendency. Here we further get the central tendency bias and Bayes factor analyses: ```{r cti_bayes} # --- Bayes t-tests bftest = function(df){ df = as.data.frame(df) # get the means rdf = df%>% summarise(mci = mean(ci), se_ci = sd(ci)/sqrt(n()), msi = mean(si), se_si = sd(si)/sqrt(n())) rdf$ci_bf = ttestBF(df$ci, mu = 0) %>% extractBF() %>% .$bf rdf$si_bf = ttestBF(df$si, mu = 0) %>% extractBF() %>% .$bf return(rdf) } vslopes %>% group_by(group, Volatility, Order) %>% nest() %>% mutate(bf = map(data, bftest)) %>% unnest(bf, .drop = TRUE) # group comparison for the low-vol-first, high-vol session vslopes %>% filter(Volatility == 'High Vola.', Order == 'LV First') %>% as.data.frame() ->v11 t_test(data = v11, formula = ci ~ group) ttestBF(data = v11, formula = ci ~ group) vslopes %>% filter(Volatility == 'Low Vola.', Order == 'LV First') %>% as.data.frame() ->v12 t_test(data = v12, formula = ci ~ group) ttestBF(data = v12, formula = ci ~ group) ``` 3. ANOVA analyses for the serial dependence indices: ```{r ANOVAs_sdi} # ---- Serial dependence index---- anova2 = ezANOVA(data = vslopes, dv = si, wid = sub, within = .(Volatility), between = .(group, Order)) anova2$ANOVA # Bayes factors bf = anovaBF(si ~ Volatility + group + Order, data = vslopes, whichRandom = "sub") bayesfactor_inclusion(bf) #bayes inclusion values ## separate for session order anova2a = ezANOVA(data = vslopes %>% filter(Order == 'LV First'), dv = si, wid = sub, within = .(Volatility), between = .(group)) anova2a$ANOVA anova2b = ezANOVA(data = vslopes %>% filter(Order == 'HV First'), dv = si, wid = sub, within = .(Volatility), between = .(group)) anova2b$ANOVA ``` 4. General biases Additionally we examined the general biases: ```{r gBias} # ---- descriptive of general bias ---- vslopes %>% group_by(group) %>% summarise(mg = mean(gBias), se = sd(gBias)/sqrt(n())) # ANOVA shows no differences anova2 = ezANOVA(data = vslopes, dv = gBias, wid = sub, within = .(Volatility), between = .(group, Order)) anova2$ANOVA bf = anovaBF(gBias ~ group + Volatility + Order, data = as.data.frame(vslopes), whichRandom = "sub") bayesfactor_inclusion(bf) #bayes inclusion values ``` We then visualize the general biases ```{r plotgBias} # --- general bias fig_bias = vslopes %>% group_by(group, Volatility, Order) %>% summarise(mbias = mean(gBias), n = n(), se = sd(gBias)/sqrt(n)) %>% ggplot(aes(Volatility, mbias, color = Order,linetype = group, group = interaction(Order, group), shape = group)) + geom_point(size = 3, position = pd) + geom_line(position = pd) + geom_errorbar(aes(ymin = mbias - se, ymax = mbias + se), width = 0.2, position = pd) + theme_classic() + theme(legend.position = 'bottom') + xlab('Volatility') + ylab('Mean overestimation (ms)') fig_bias ``` There was no difference in two groups in general biases, although both groups were positive overestimated. Next we examined the reproduced variability: ```{r rep_var} # ---- Reproduction variability ----- msds_r <- mrep %>% filter(n>5) %>% group_by(group, Order, Volatility, sub) %>% summarise(msd = mean(sdr), n=n(), msd_se = sd(sdr)/sqrt(n)) sd_ANOVA <- ezANOVA(msds_r, dv=msd, wid=sub, between=.(group,Order), within=Volatility) sd_ANOVA$ANOVA mmsds_r <- msds_r %>% summarize(mmsd=mean(msd*1000), n = n(), se = sd(msd*1000)/sqrt(n)) pd = position_dodge(width = 0.5) fig_sd <- ggplot(mmsds_r, aes(Volatility, mmsd, shape = Order, color = Order, group = Order)) + geom_line() + geom_point() + #geom_bar(stat = 'identity', position = pd, width = 0.5) + geom_errorbar(aes(ymin = mmsd - se, ymax = mmsd + se), width = 0.3) + facet_wrap(~group) + ylab('Mean Standard Deviation (ms)') + xlab('Volatility') + theme_classic() # coord_cartesian(ylim = c(140, 200)) fig_sd if(saveFig){ ggsave('figures/fig_sd.png', fig_sd, width=5, height=4) ggsave('figures/fig_sd.pdf', fig_sd, width=5, height=4) } ``` ## Correlation analysis Next we import individual participant ratings and do correlation analyses. ```{r parinfo} pinfo = read.csv('data/parinfo.csv') # join with ctis, sdis res = left_join(vslopes, pinfo, by = c('group','sequence')) res$group = as.factor(res$group) # scatter plot to see any relations fig_corr = ggduo(res, c('AQ','EQ','BDI'), c('ci','si'), mapping = aes(color = group, shape = Order), types = list(continuous = 'smooth_lm', se = FALSE), showStrips = FALSE, legend = 3, columnLabelsY = c('CTI','SDI')) + theme_classic() + theme(legend.position = 'bottom', strip.background = element_blank()) if(saveFig){ ggsave('figures/fig_corr.png', fig_corr, width=5, height=4) ggsave('figures/fig_corr.pdf', fig_corr, width=5, height=4) } fig_corr ``` Get linear regression out ```{r} ci_AQ <- function(df){ lm(ci ~ AQ, data = df) } ci_EQ <- function(df){ lm(ci ~ EQ, data = df) } ci_BDI <- function(df){ lm(ci ~ BDI, data = df) } si_AQ <- function(df){ lm(si ~ AQ, data = df) } si_EQ <- function(df){ lm(si ~ EQ, data = df) } si_BDI <- function(df){ lm(si ~ BDI, data = df) } # regression reg_res <- res %>% group_by(Volatility, group, Order) %>% nest() %>% # nested data mutate(ci_AQ = map(data, ci_AQ), ci_EQ = map(data, ci_EQ), ci_BDI = map(data, ci_BDI), si_AQ = map(data, si_AQ), si_EQ = map(data, si_EQ), si_BDI = map(data, si_BDI)) %>% mutate(mci_AQ = map(ci_AQ, broom::tidy), mci_EQ = map(ci_EQ, broom::tidy), mci_BDI = map(ci_BDI, broom::tidy), msi_AQ = map(si_AQ, broom::tidy), msi_EQ = map(si_EQ, broom::tidy), msi_BDI = map(si_BDI, broom::tidy) ) %>% unnest(cols = c(mci_AQ, mci_EQ, mci_BDI, msi_AQ, msi_EQ, msi_BDI), names_repair = 'unique' , .drop = TRUE) %>% select(-ci_AQ, -ci_BDI, -ci_EQ, -si_AQ, -si_EQ, -si_BDI, -data, -starts_with('st')) %>% filter(term == 'AQ') # remove intercept ``` Given there is no significant of slopes (after correction), but some difference in intercepts, we do general linear regression, without separating groups. ```{r} reg_res2 <- res %>% ungroup() %>% group_by(Order, Volatility) %>% nest() %>% # nested data mutate(ci_AQ = map(data, ci_AQ), ci_EQ = map(data, ci_EQ), ci_BDI = map(data, ci_BDI), si_AQ = map(data, si_AQ), si_EQ = map(data, si_EQ), si_BDI = map(data, si_BDI)) %>% mutate(mci_AQ = map(ci_AQ, broom::tidy), mci_EQ = map(ci_EQ, broom::tidy), mci_BDI = map(ci_BDI, broom::tidy), msi_AQ = map(si_AQ, broom::tidy), msi_EQ = map(si_EQ, broom::tidy), msi_BDI = map(si_BDI, broom::tidy) ) %>% unnest(cols = c(mci_AQ, mci_EQ, mci_BDI, msi_AQ, msi_EQ, msi_BDI), names_repair = 'unique' , .drop = TRUE) %>% select(-ci_AQ, -ci_BDI, -ci_EQ, -si_AQ, -si_EQ, -si_BDI, -data, -starts_with('st')) %>% filter(term == 'AQ') # remove intercept ``` Again, there is no significant of slopes.