From 6448dccd8213a228e06435d94d941c0544355766 Mon Sep 17 00:00:00 2001
From: Anthony Dick <31479113+anthonystevendick@users.noreply.github.com>
Date: Tue, 16 Apr 2019 09:41:11 -0400
Subject: [PATCH] Update bilingual_analysis.r

This file was updated to reflect the analysis for the accepted paper in Nature Human Behavior.
---
 bilingual_analysis.r | 546 ++++++++++++++++++++++++++++---------------
 1 file changed, 352 insertions(+), 194 deletions(-)

diff --git a/bilingual_analysis.r b/bilingual_analysis.r
index cb636da..3cf1e5b 100644
--- a/bilingual_analysis.r
+++ b/bilingual_analysis.r
@@ -1,6 +1,7 @@
+#Updated 4/16/19 by Anthony Dick
 
 rm(list=ls())
-setwd("<path to data>")
+setwd("<path to data")
 list.files()
 
 ####################
@@ -91,6 +92,62 @@ accult_q2_y[which(accult_q3_other_y=="British")]<-"0"
 accult_q2_y[which(accult_q3_other_y=="English")]<-"0"
 SSRTr<-(beh_sst_ssrt_mean_total) * -1 #Stop-Signal reaction time, recode: recode/reverse score SSRT (make lower SSRT mean poorer performance)
 
+#Missing Data Summary
+
+sum(is.na(nda17$age))
+sum(is.na(nda17$female))
+sum(is.na(nda17$household.income))
+sum(is.na(nda17$high.educ))
+sum(is.na(nda17$married))
+sum(is.na(nda17$race.ethnicity))
+sum(is.na(nda17$nihtbx_reading_uncorrected))
+sum(is.na(nda17$nihtbx_picture_uncorrected))
+sum(is.na(nda17$nihtbx_list_uncorrected))
+sum(is.na(nda17$nihtbx_pattern_uncorrected))
+sum(is.na(nda17$nihtbx_picvocab_uncorrected))
+sum(is.na(nda17$nihtbx_flanker_uncorrected))
+sum(is.na(nda17$nihtbx_cardsort_uncorrected))
+sum(is.na(nda17$beh_sst_ssrt_mean_total))
+sum(is.na(nda17$abcd_site))
+sum(is.na(nda17$rel_family_id))
+
+#predictor independent variables (i.e., the ultimate outcome variables of interest)
+
+predictor_iv = c(which(names(nda17)=="nihtbx_picvocab_uncorrected"),which(names(nda17)=="nihtbx_flanker_uncorrected"),which(names(nda17)=="nihtbx_cardsort_uncorrected"), which(names(nda17)=="beh_sst_ssrt_mean_total"))
+for(j in 1:length(predictor_iv)) nda17[,predictor_iv[j]] = scale(as.numeric(nda17[,predictor_iv[j]])) # standardize the DV to get standardized betas
+
+#outcome variables, demographic variables with non-zero missing cases
+
+outcome_dv = c(which(names(nda17)=="race.ethnicity"), which(names(nda17)=="high.educ"), which(names(nda17)=="household.income"))
+
+pause = function()
+{
+    if (interactive())
+    {
+        invisible(readline(prompt = "Press <Enter> to continue..."))
+    }
+    else
+    {
+        cat("Press <Enter> to continue...")
+        invisible(readLines(file("stdin"), 1))
+    }
+}
+for(i in 1:length(outcome_dv)){
+print(i)
+	for(j in 1:length(predictor_iv)){
+	print(j)
+		pred<-paste(names(nda17)[predictor_iv[j]],"-->",names(nda17)[outcome_dv[i]])
+		print(pred)
+		form = paste("as.numeric(","is.na(",names(nda17)[outcome_dv[i]],")",")","~" ,names(nda17)[predictor_iv[j]], sep = "")
+		print(form)
+		log.mod<-glm(form, family = binomial, data = nda17)
+    print(summary(log.mod))
+	}
+}
+
+###End Missing Data Summary############
+#######################################
+
 #create a smaller data frame that is easier to manage
 
 abcd_subset<-data.frame(completedData$src_subject_id, as.numeric(nda17$nihtbx_reading_uncorrected), as.numeric(nda17$nihtbx_picture_uncorrected), as.numeric(nda17$nihtbx_list_uncorrected),
@@ -103,9 +160,10 @@ colnames(abcd_subset)<-c("src_subject_id", "nihtbx_reading_uncorrected", "nihtbx
 "nihtbx_picvocab_uncorrected", "nihtbx_flanker_uncorrected", "nihtbx_cardsort_uncorrected", "beh_sst_ssrt_mean_total", "SSRTr", "abcd_site","rel_family_id","accult_q1_y","accult_q2_y",
 "accult_q3_dropdwn_y","accult_q3_other_y","accult_q4_y","accult_q5_y")
 
-#saveRDS(abcd_subset, "<filename>.Rds") # save this dataframe for future use
-#abcd_subset<-readRDS("<filename>.Rds") # load the dataframe if you come back to the analysis at a later point so you can start here
+#saveRDS(abcd_subset, "abcd_subset_7_8_18.Rds") # save this dataframe for future use
+#abcd_subset<-readRDS("abcd_subset_7_8_18.Rds") # load the dataframe
 
+#summarize data
 
 names(abcd_subset)
 dim(abcd_subset)
@@ -120,12 +178,13 @@ count(accult_q4_y)
 count(accult_q5_y)
 count(nda17$demo_origin_v2)
 count(nda17$demo_prnt_16)
-count(nda17$demo_biofather_v2) #
-count(nda17$demo_biomother_v2) #
+count(nda17$demo_biofather_v2) # 814  0.1799293
+count(nda17$demo_biomother_v2) # 540  0.1193634
 
 #recode the accult_q2_y variable into a binary "Bilingual Status", 0 = not bilingual; 1 = bilingual
 
 bilingual_status <- accult_q2_y
+sum(is.na(bilingual_status))
 
 #dimension a 'bilingual degree' variable, where 1 = participant said they were bilingual, and they speak the other language with friends all the time, most of the time,
 #or equally, OR they speak the other language with family all the time, most of the time, or equally.
@@ -133,11 +192,13 @@ bilingual_status <- accult_q2_y
 bilingual_degree <- ifelse(bilingual_status == 0, 0, ifelse(bilingual_status == 1 & (as.numeric(accult_q4_y) <= 3 | as.numeric(accult_q5_y) <= 3), 1, NA))
 
 count(bilingual_degree) #check the data
+sum(is.na(bilingual_degree))
 
 #dimension a continuous 'bilingual use' variable, and reverse-score so that if participants speak the other language with friends all the time, most of the time...,
 #they will receive high scores on this measure (range 0-8, with 8 indicating a high-degree of other language use)
 
 bilingual_use<-10-(as.numeric(abcd_subset$accult_q4_y)+as.numeric(abcd_subset$accult_q5_y))
+sum(is.na(bilingual_use))
 
 count(bilingual_use) # check the data
 for.hist.use <- melt(bilingual_use)
@@ -155,6 +216,7 @@ dev.off()
 
 fluidIQ<-scale((nihtbx_picture_uncorrected+nihtbx_list_uncorrected+nihtbx_pattern_uncorrected)/3)
 hist(fluidIQ)
+sum(is.na(fluidIQ))
 
 ###################################################
 ###################################################
@@ -162,24 +224,59 @@ hist(fluidIQ)
 ###################################################
 ###################################################
 
-#abcd_subset<-cbind.data.frame(abcd_subset, fluidIQ, bilingual_status, bilingual_degree, bilingual_use) # create a new data frame
 abcd_subset<-cbind.data.frame(abcd_subset, fluidIQ, bilingual_status, bilingual_degree, bilingual_use) # create a new data frame
+#abcd_subset.noimp<-cbind.data.frame(abcd_subset, fluidIQ, bilingual_status, bilingual_degree, bilingual_use) # create a new data frame
+
+#Descriptive Statistics by bilingual grouping
+stat_vocab<-ddply(abcd_subset,~bilingual_status,summarise,mean=mean(nihtbx_picvocab_uncorrected, na.rm = TRUE),sd=sd(nihtbx_picvocab_uncorrected, na.rm = TRUE))
+stat_flanker<-ddply(abcd_subset,~bilingual_status,summarise,mean=mean(nihtbx_flanker_uncorrected, na.rm = TRUE),sd=sd(nihtbx_flanker_uncorrected, na.rm = TRUE))
+stat_card<-ddply(abcd_subset,~bilingual_status,summarise,mean=mean(nihtbx_cardsort_uncorrected, na.rm = TRUE),sd=sd(nihtbx_cardsort_uncorrected, na.rm = TRUE))
+stat_SST<-ddply(abcd_subset,~bilingual_status,summarise,mean=mean(SSRTr, na.rm = TRUE),sd=sd(SSRTr, na.rm = TRUE))
+
+degree_vocab<-ddply(abcd_subset,~bilingual_degree,summarise,mean=mean(nihtbx_picvocab_uncorrected, na.rm = TRUE),sd=sd(nihtbx_picvocab_uncorrected, na.rm = TRUE))
+degree_flanker<-ddply(abcd_subset,~bilingual_degree,summarise,mean=mean(nihtbx_flanker_uncorrected, na.rm = TRUE),sd=sd(nihtbx_flanker_uncorrected, na.rm = TRUE))
+degree_card<-ddply(abcd_subset,~bilingual_degree,summarise,mean=mean(nihtbx_cardsort_uncorrected, na.rm = TRUE),sd=sd(nihtbx_cardsort_uncorrected, na.rm = TRUE))
+degree_SST<-ddply(abcd_subset,~bilingual_degree,summarise,mean=mean(SSRTr, na.rm = TRUE),sd=sd(SSRTr, na.rm = TRUE))
+
+use_vocab<-ddply(abcd_subset,~bilingual_use,summarise,mean=mean(nihtbx_picvocab_uncorrected, na.rm = TRUE),sd=sd(nihtbx_picvocab_uncorrected, na.rm = TRUE))
+use_flanker<-ddply(abcd_subset,~bilingual_use,summarise,mean=mean(nihtbx_flanker_uncorrected, na.rm = TRUE),sd=sd(nihtbx_flanker_uncorrected, na.rm = TRUE))
+use_card<-ddply(abcd_subset,~bilingual_use,summarise,mean=mean(nihtbx_cardsort_uncorrected, na.rm = TRUE),sd=sd(nihtbx_cardsort_uncorrected, na.rm = TRUE))
+use_SST<-ddply(abcd_subset,~bilingual_use,summarise,mean=mean(SSRTr, na.rm = TRUE),sd=sd(SSRTr, na.rm = TRUE))
+
+col1<-c(round(stat_vocab$mean[1], 1), round(stat_flanker$mean[1], 1), round(stat_card$mean[1], 1), round(stat_SST$mean[1], 1),round(degree_vocab$mean[1], 1), round(degree_flanker$mean[1], 1), round(degree_card$mean[1], 1), round(degree_SST$mean[1], 1))
+col2<-c(round(stat_vocab$sd[1], 1), round(stat_flanker$sd[1], 1), round(stat_card$sd[1], 1), round(stat_SST$sd[1], 1),round(degree_vocab$sd[1], 1), round(degree_flanker$sd[1], 1), round(degree_card$sd[1], 1), round(degree_SST$sd[1], 1))
+col3<-c(round(stat_vocab$mean[2], 1), round(stat_flanker$mean[2], 1), round(stat_card$mean[2], 1), round(stat_SST$mean[2], 1),round(degree_vocab$mean[2], 1), round(degree_flanker$mean[2], 1), round(degree_card$mean[2], 1), round(degree_SST$mean[2], 1))
+col4<-c(round(stat_vocab$sd[2], 1), round(stat_flanker$sd[2], 1), round(stat_card$sd[2], 1), round(stat_SST$sd[2], 1),round(degree_vocab$sd[2], 1), round(degree_flanker$sd[2], 1), round(degree_card$sd[2], 1), round(degree_SST$sd[2], 1))
+descr.m<-data.frame(col1,col2,col3,col4)
+descr.use.m<-data.frame(use_vocab, use_flanker, use_card, use_SST)
+
+colnames(descr.m)<-c("Mean Monolingual","Standard Deviation Monolingual", "Mean Bilingual", "Standard Deviation Bilingual")
+write.table(descr.m, file = "descriptives.csv", sep = ",")
+
 
 ###### Select covariates for model  ######
 
-ind_cov = c(which(names(abcd_subset)=="age"),which(names(abcd_subset)=="female"),which(names(abcd_subset)=="race.ethnicity"), which(names(abcd_subset)=="high.educ"),which(names(abcd_subset)=="married"),which(names(abcd_subset)=="household.income"))
+##Vocabulary Only##
+
+ind_cov = c(which(names(abcd_subset)=="nihtbx_picvocab_uncorrected"))
+names(abcd_subset)[ind_cov]
+summary(abcd_subset[,ind_cov])
+
+###### Select demographic covariates for model  ######
+
+ind_cov = c(which(names(abcd_subset)=="age"), which(names(abcd_subset)=="female"), which(names(abcd_subset)=="race.ethnicity"), which(names(abcd_subset)=="high.educ"), which(names(abcd_subset)=="married"), which(names(abcd_subset)=="household.income"))
 names(abcd_subset)[ind_cov]
 summary(abcd_subset[,ind_cov])
 
 ###### Select covariates for model, including IQ  ######
 
-#ind_cov = c(which(names(abcd_subset)=="nihtbx_reading_uncorrected"), which(names(abcd_subset)=="fluidIQ"), which(names(abcd_subset)=="age"),which(names(abcd_subset)=="female"),which(names(abcd_subset)=="race.thnicity"), which(names(abcd_subset)=="high.educ"),which(names(abcd_subset)=="married"),which(names(abcd_subset)=="household.income"))
+ind_cov = c(which(names(abcd_subset)=="nihtbx_reading_uncorrected"), which(names(abcd_subset)=="fluidIQ"), which(names(abcd_subset)=="age"),which(names(abcd_subset)=="female"),which(names(abcd_subset)=="race.thnicity"), which(names(abcd_subset)=="high.educ"),which(names(abcd_subset)=="married"),which(names(abcd_subset)=="household.income"))
 
 ## Select covariates for model, including Picture Vocabulary and IQ
 
-#ind_cov = c(which(names(abcd_subset)=="nihtbx_reading_uncorrected"), which(names(abcd_subset)=="fluidIQ"), which(names(abcd_subset)=="nihtbx_picvocab_uncorrected"), which(names(abcd_subset)=="age"),which(names(abcd_subset)=="female"),which(names(abcd_subset)=="race.thnicity"), which(names(abcd_subset)=="high.educ"),which(names(abcd_subset)=="married"),which(names(abcd_subset)=="household.income"))
+ind_cov = c(which(names(abcd_subset)=="nihtbx_reading_uncorrected"), which(names(abcd_subset)=="fluidIQ"), which(names(abcd_subset)=="nihtbx_picvocab_uncorrected"), which(names(abcd_subset)=="age"),which(names(abcd_subset)=="female"),which(names(abcd_subset)=="race.thnicity"), which(names(abcd_subset)=="high.educ"),which(names(abcd_subset)=="married"),which(names(abcd_subset)=="household.income"))
 
-## Select nesting variables
+## Select nesting variables (doesn't change)
 
 ind_nest = c(which(names(abcd_subset)=="abcd_site"),which(names(abcd_subset)=="rel_family_id"));names(abcd_subset)[ind_nest]
 
@@ -189,7 +286,8 @@ ind_dv = c(which(names(abcd_subset)=="nihtbx_picvocab_uncorrected"),which(names(
 
 #ind_dv = c(which(names(abcd_subset)=="nihtbx_flanker_uncorrected"),which(names(abcd_subset)=="nihtbx_cardsort_uncorrected"), which(names(abcd_subset)=="SSRTr"))
 
-#for(j in 1:length(ind_dv)) abcd_subset[,ind_dv[j]] = scale(as.numeric(abcd_subset[,ind_dv[j]])) # standardize the DV to get standardized betas
+abcd_subset_scale<-abcd_subset
+for(j in 1:length(ind_dv)) abcd_subset_scale[,ind_dv[j]] = scale(as.numeric(abcd_subset_scale[,ind_dv[j]])) # standardize the DV to get standardized betas
 
 names(abcd_subset)[ind_dv]
 boxplot(abcd_subset[ind_dv])
@@ -197,9 +295,9 @@ summary(abcd_subset[ind_dv])
 
 ###### Select IVs  ######
 
-ind_iv = c(which(names(abcd_subset)=="bilingual_status"),which(names(abcd_subset)=="bilingual_degree"), which(names(abcd_subset)=="bilingual_use"))
+ind_iv = c(which(names(abcd_subset)=="bilingual_status"), which(names(abcd_subset)=="bilingual_degree"), which(names(abcd_subset)=="bilingual_use"))
 
-#for(k in 1:length(ind_iv)) abcd_subset[,ind_iv[k]] = scale(as.numeric(abcd_subset[,ind_iv[k]])) # standardize the IV to get standardized betas
+for(k in 1:length(ind_iv)) abcd_subset_scale[,ind_iv[k]] = scale(as.numeric(abcd_subset_scale[,ind_iv[k]])) # standardize the IV to get standardized betas
 
 names(abcd_subset)[ind_iv]
 summary(abcd_subset[,ind_iv])
@@ -225,8 +323,8 @@ pause = function()
         invisible(readLines(file("stdin"), 1))
     }
 }
-m <- matrix(nrow=33, ncol=9) # set the output matrix size to hold the results
-mat_row = 25 # to fill out the m matrix, start at row 13, row 25 when running with different covariates
+m <- matrix(nrow=42, ncol=15) # set the output matrix size to hold the results
+mat_row = 34 # to fill out the m matrix, start at row 1, row 10, row 22, row 34  when running with different covariates
 for(i in 1:length(ind_dv)){
 print(i)
 	for(j in 1:length(ind_iv)){
@@ -248,31 +346,32 @@ print(i)
 		ran = formula(ran)
 			mod0 = gamm4(formula = form.simple, random = ran, data = abcd_subset)
 			print(summary(mod0$gam))
+			mod1 = gamm4(formula = form.simple, random = ran, data = abcd_subset_scale)
+			#print(summary(mod1$gam))
 			#pause()
 			gamm0 = gamm4(formula = form, random = ran, data = abcd_subset)
 			#pause()
 			print(summary(gamm0$gam))
-			output<-c(pred,summary(mod0$gam)$p.coeff[2],summary(mod0$gam)$se[2],summary(mod0$gam)$p.t[2],summary(mod0$gam)$p.pv[2],summary(gamm0$gam)$p.coeff[2],summary(gamm0$gam)$se[2],summary(gamm0$gam)$p.t[2],summary(gamm0$gam)$p.pv[2])
+			gamm1 = gamm4(formula = form, random = ran, data = abcd_subset_scale)
+			#pause()
+			print(summary(gamm1$gam))
+			output<-c(pred,mod0$gam$df.residual, summary(mod0$gam)$p.coeff[2],summary(mod0$gam)$se[2],summary(mod1$gam)$p.coeff[2],summary(mod1$gam)$se[2], summary(mod0$gam)$p.t[2],summary(mod0$gam)$p.pv[2], gamm0$gam$df.residual, summary(gamm0$gam)$p.coeff[2],summary(gamm0$gam)$se[2], summary(gamm1$gam)$p.coeff[2], summary(gamm1$gam)$se[2], summary(gamm0$gam)$p.t[2],summary(gamm0$gam)$p.pv[2])
 			m[mat_row,]<-output # place the output into the matrix that was dimensioned above
 			mat_row = mat_row + 1 # step to the next matrix row
-		#print("########################################")
-		#print("########################################")
 	}
 }
-colnames(m)<-c("prediction","nocov_b","nocov_se","nocov_tval","nocov_pval","gamm_b","gamm_se","gamm_tval","gamm_pval")
+colnames(m)<-c("prediction","mod0_df","nocov_b","nocov_se","nocov_beta","nocov_se_beta","nocov_tval","nocov_pval","gamm_df", "gamm_b","gamm_se","gamm_beta","gamm_se_beta","gamm_tval","gamm_pval")
 print(m)
-write.table(m,file="filename.txt") # save for filling out tables, figures
+write.table(m,file="m_bvals_11_27_18.txt")
+#m<-read.csv(file="m_bvals_11_15_18.csv")
 
 ################################################
 ################################################
-## Tests of Equivalence                       ##
+## Two-One-Sided Tests of Equivalence (TOSTs) ##
 ################################################
 ################################################
 
-## The function below uses regression slopes not raw data
-## Based on Counsell, A., & Cribbie, R. A. (2014). Equivalence tests for comparing correlation and regression coefficients.
-## British Journal of Mathematical and Statistical Psychology, 68, 292-309.
-## Special thanks to Alyssa Counsell for sharing the code
+## function below uses regression slopes not raw data
 
 # b1 is the regression coefficient for group 1
 # b2 is the regression coefficient for group 2
@@ -298,203 +397,249 @@ equivbs<-function(b1,b2,se1,se2,equiv_int,alpha=.05){
 	names(out)<-c("P value", "CI", "Decision")
 print(out)
 }
-########
-#use the m matrix from the regressions to fill out the figure
 
-y<-rev(c(NA, NA, as.numeric(m[1,2]),
-as.numeric(m[2,2]),
-as.numeric(m[3,2]),
+
+y<-rev(c(NA, NA, as.numeric(m[10,4]),
+as.numeric(m[11,4]),
+as.numeric(m[12,4]),
 NA,
-as.numeric(m[4,2]),
-as.numeric(m[5,2]),
-as.numeric(m[6,2]),
+as.numeric(m[13,4]),
+as.numeric(m[14,4]),
+as.numeric(m[15,4]),
 NA,
-as.numeric(m[7,2]),
-as.numeric(m[8,2]),
-as.numeric(m[9,2]),
+as.numeric(m[16,4]),
+as.numeric(m[17,4]),
+as.numeric(m[18,4]),
 NA,
-as.numeric(m[10,2]),
-as.numeric(m[11,2]),
-as.numeric(m[12,2]),
+as.numeric(m[19,4]),
+as.numeric(m[20,4]),
+as.numeric(m[21,4]),
 NA, NA, NA,
-as.numeric(m[1,6]),
-as.numeric(m[2,6]),
-as.numeric(m[3,6]),
+as.numeric(m[1,10]),
+as.numeric(m[2,10]),
+as.numeric(m[3,10]),
 NA,
-as.numeric(m[4,6]),
-as.numeric(m[5,6]),
-as.numeric(m[6,6]),
+as.numeric(m[4,10]),
+as.numeric(m[5,10]),
+as.numeric(m[6,10]),
 NA,
-as.numeric(m[7,6]),
-as.numeric(m[8,6]),
-as.numeric(m[9,6]),
+as.numeric(m[7,10]),
+as.numeric(m[8,10]),
+as.numeric(m[9,10]),
+NA, NA, NA,
+as.numeric(m[10,10]),
+as.numeric(m[11,10]),
+as.numeric(m[12,10]),
+NA,
+as.numeric(m[13,10]),
+as.numeric(m[14,10]),
+as.numeric(m[15,10]),
 NA,
-as.numeric(m[10,6]),
-as.numeric(m[11,6]),
-as.numeric(m[12,6]),
+as.numeric(m[16,10]),
+as.numeric(m[17,10]),
+as.numeric(m[18,10]),
+NA,
+as.numeric(m[19,10]),
+as.numeric(m[20,10]),
+as.numeric(m[21,10]),
 NA, NA, NA,
-as.numeric(m[13,6]),
-as.numeric(m[14,6]),
-as.numeric(m[15,6]),
+as.numeric(m[22,10]),
+as.numeric(m[23,10]),
+as.numeric(m[24,10]),
 NA,
-as.numeric(m[16,6]),
-as.numeric(m[17,6]),
-as.numeric(m[18,6]),
+as.numeric(m[25,10]),
+as.numeric(m[26,10]),
+as.numeric(m[27,10]),
 NA,
-as.numeric(m[19,6]),
-as.numeric(m[20,6]),
-as.numeric(m[21,6]),
+as.numeric(m[28,10]),
+as.numeric(m[29,10]),
+as.numeric(m[30,10]),
 NA,
-as.numeric(m[22,6]),
-as.numeric(m[23,6]),
-as.numeric(m[24,6]),
+as.numeric(m[31,10]),
+as.numeric(m[32,10]),
+as.numeric(m[33,10]),
 NA, NA, NA,
-as.numeric(m[25,6]),
-as.numeric(m[26,6]),
-as.numeric(m[27,6]),
+as.numeric(m[34,10]),
+as.numeric(m[35,10]),
+as.numeric(m[36,10]),
 NA,
-as.numeric(m[28,6]),
-as.numeric(m[29,6]),
-as.numeric(m[30,6]),
+as.numeric(m[37,10]),
+as.numeric(m[38,10]),
+as.numeric(m[39,10]),
 NA,
-as.numeric(m[31,6]),
-as.numeric(m[32,6]),
-as.numeric(m[33,6])))
+as.numeric(m[40,10]),
+as.numeric(m[41,10]),
+as.numeric(m[42,10])))
+
 
-ylo<-rev(c(NA, NA, equivbs(as.numeric(m[1,2]),0,as.numeric(m[1,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[2,2]),0,as.numeric(m[2,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[3,2]),0,as.numeric(m[3,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+ylo<-rev(c(NA, NA, equivbs(as.numeric(m[10,4]),0,as.numeric(m[10,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[11,4]),0,as.numeric(m[11,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[12,4]),0,as.numeric(m[12,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[4,2]),0,as.numeric(m[4,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[5,2]),0,as.numeric(m[5,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[6,2]),0,as.numeric(m[6,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[13,4]),0,as.numeric(m[13,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[14,4]),0,as.numeric(m[14,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[15,4]),0,as.numeric(m[15,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+NA,
+equivbs(as.numeric(m[16,4]),0,as.numeric(m[16,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[17,4]),0,as.numeric(m[17,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[18,4]),0,as.numeric(m[18,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+NA,
+equivbs(as.numeric(m[19,4]),0,as.numeric(m[19,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[20,4]),0,as.numeric(m[20,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[21,4]),0,as.numeric(m[21,5]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+NA, NA, NA,
+equivbs(as.numeric(m[1,10]),0,as.numeric(m[1,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[2,10]),0,as.numeric(m[2,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[3,10]),0,as.numeric(m[3,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[7,2]),0,as.numeric(m[7,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[8,2]),0,as.numeric(m[8,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[9,2]),0,as.numeric(m[9,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[4,10]),0,as.numeric(m[4,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[5,10]),0,as.numeric(m[5,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[6,10]),0,as.numeric(m[6,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[10,2]),0,as.numeric(m[10,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[11,2]),0,as.numeric(m[11,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[12,2]),0,as.numeric(m[12,3]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[7,10]),0,as.numeric(m[7,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[8,10]),0,as.numeric(m[8,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[9,10]),0,as.numeric(m[9,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA, NA, NA,
-equivbs(as.numeric(m[1,6]),0,as.numeric(m[1,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[2,6]),0,as.numeric(m[2,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[3,6]),0,as.numeric(m[3,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[10,10]),0,as.numeric(m[10,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[11,10]),0,as.numeric(m[11,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[12,10]),0,as.numeric(m[12,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[4,6]),0,as.numeric(m[4,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[5,6]),0,as.numeric(m[5,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[6,6]),0,as.numeric(m[6,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[13,10]),0,as.numeric(m[13,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[14,10]),0,as.numeric(m[14,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[15,10]),0,as.numeric(m[15,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[7,6]),0,as.numeric(m[7,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[8,6]),0,as.numeric(m[8,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[9,6]),0,as.numeric(m[9,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[16,10]),0,as.numeric(m[16,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[17,10]),0,as.numeric(m[17,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[18,10]),0,as.numeric(m[18,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[10,6]),0,as.numeric(m[10,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[11,6]),0,as.numeric(m[11,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[12,6]),0,as.numeric(m[12,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[19,10]),0,as.numeric(m[19,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[20,10]),0,as.numeric(m[20,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[21,10]),0,as.numeric(m[21,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA, NA, NA,
-equivbs(as.numeric(m[13,6]),0,as.numeric(m[13,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[14,6]),0,as.numeric(m[14,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[15,6]),0,as.numeric(m[15,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[22,10]),0,as.numeric(m[22,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[23,10]),0,as.numeric(m[23,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[24,10]),0,as.numeric(m[24,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[16,6]),0,as.numeric(m[16,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[17,6]),0,as.numeric(m[17,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[18,6]),0,as.numeric(m[18,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[25,10]),0,as.numeric(m[25,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[26,10]),0,as.numeric(m[26,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[27,10]),0,as.numeric(m[27,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[19,6]),0,as.numeric(m[19,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[20,6]),0,as.numeric(m[20,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[21,6]),0,as.numeric(m[21,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[28,10]),0,as.numeric(m[28,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[29,10]),0,as.numeric(m[29,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[30,10]),0,as.numeric(m[30,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[22,6]),0,as.numeric(m[22,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[23,6]),0,as.numeric(m[23,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[24,6]),0,as.numeric(m[24,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[31,10]),0,as.numeric(m[31,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[32,10]),0,as.numeric(m[32,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[33,10]),0,as.numeric(m[33,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA, NA, NA,
-equivbs(as.numeric(m[25,6]),0,as.numeric(m[25,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[26,6]),0,as.numeric(m[26,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[27,6]),0,as.numeric(m[27,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[34,10]),0,as.numeric(m[34,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[35,10]),0,as.numeric(m[35,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[36,10]),0,as.numeric(m[36,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[28,6]),0,as.numeric(m[28,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[29,6]),0,as.numeric(m[29,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[30,6]),0,as.numeric(m[30,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[37,10]),0,as.numeric(m[37,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[38,10]),0,as.numeric(m[38,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[39,10]),0,as.numeric(m[39,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
 NA,
-equivbs(as.numeric(m[31,6]),0,as.numeric(m[31,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[32,6]),0,as.numeric(m[32,7]),0,c(-.1, .1),alpha = .05)$`CI`[1],
-equivbs(as.numeric(m[33,6]),0,as.numeric(m[33,7]),0,c(-.1, .1),alpha = .05)$`CI`[1]))
+equivbs(as.numeric(m[40,10]),0,as.numeric(m[40,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[41,10]),0,as.numeric(m[41,11]),0,c(-.1, .1),alpha = .05)$`CI`[1],
+equivbs(as.numeric(m[42,10]),0,as.numeric(m[42,11]),0,c(-.1, .1),alpha = .05)$`CI`[1]))
 
-yhi<-rev(c(NA, NA, equivbs(as.numeric(m[1,2]),0,as.numeric(m[1,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[2,2]),0,as.numeric(m[2,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[3,2]),0,as.numeric(m[3,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+
+yhi<-rev(c(NA, NA, equivbs(as.numeric(m[10,4]),0,as.numeric(m[10,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[11,4]),0,as.numeric(m[11,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[12,4]),0,as.numeric(m[12,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+NA,
+equivbs(as.numeric(m[13,4]),0,as.numeric(m[13,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[14,4]),0,as.numeric(m[14,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[15,4]),0,as.numeric(m[15,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+NA,
+equivbs(as.numeric(m[16,4]),0,as.numeric(m[16,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[17,4]),0,as.numeric(m[17,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[18,4]),0,as.numeric(m[18,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[4,2]),0,as.numeric(m[4,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[5,2]),0,as.numeric(m[5,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[6,2]),0,as.numeric(m[6,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[19,4]),0,as.numeric(m[19,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[20,4]),0,as.numeric(m[20,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[21,4]),0,as.numeric(m[21,5]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+NA, NA, NA,
+equivbs(as.numeric(m[1,10]),0,as.numeric(m[1,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[2,10]),0,as.numeric(m[2,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[3,10]),0,as.numeric(m[3,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[7,2]),0,as.numeric(m[7,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[8,2]),0,as.numeric(m[8,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[9,2]),0,as.numeric(m[9,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[4,10]),0,as.numeric(m[4,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[5,10]),0,as.numeric(m[5,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[6,10]),0,as.numeric(m[6,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[10,2]),0,as.numeric(m[10,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[11,2]),0,as.numeric(m[11,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[12,2]),0,as.numeric(m[12,3]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[7,10]),0,as.numeric(m[7,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[8,10]),0,as.numeric(m[8,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[9,10]),0,as.numeric(m[9,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA, NA, NA,
-equivbs(as.numeric(m[1,6]),0,as.numeric(m[1,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[2,6]),0,as.numeric(m[2,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[3,6]),0,as.numeric(m[3,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[10,10]),0,as.numeric(m[10,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[11,10]),0,as.numeric(m[11,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[12,10]),0,as.numeric(m[12,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[4,6]),0,as.numeric(m[4,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[5,6]),0,as.numeric(m[5,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[6,6]),0,as.numeric(m[6,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[13,10]),0,as.numeric(m[13,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[14,10]),0,as.numeric(m[14,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[15,10]),0,as.numeric(m[15,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[7,6]),0,as.numeric(m[7,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[8,6]),0,as.numeric(m[8,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[9,6]),0,as.numeric(m[9,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[16,10]),0,as.numeric(m[16,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[17,10]),0,as.numeric(m[17,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[18,10]),0,as.numeric(m[18,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[10,6]),0,as.numeric(m[10,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[11,6]),0,as.numeric(m[11,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[12,6]),0,as.numeric(m[12,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[19,10]),0,as.numeric(m[19,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[20,10]),0,as.numeric(m[20,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[21,10]),0,as.numeric(m[21,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA, NA, NA,
-equivbs(as.numeric(m[13,6]),0,as.numeric(m[13,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[14,6]),0,as.numeric(m[14,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[15,6]),0,as.numeric(m[15,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[22,10]),0,as.numeric(m[22,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[23,10]),0,as.numeric(m[23,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[24,10]),0,as.numeric(m[24,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[16,6]),0,as.numeric(m[16,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[17,6]),0,as.numeric(m[17,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[18,6]),0,as.numeric(m[18,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[25,10]),0,as.numeric(m[25,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[26,10]),0,as.numeric(m[26,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[27,10]),0,as.numeric(m[27,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[19,6]),0,as.numeric(m[19,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[20,6]),0,as.numeric(m[20,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[21,6]),0,as.numeric(m[21,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[28,10]),0,as.numeric(m[28,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[29,10]),0,as.numeric(m[29,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[30,10]),0,as.numeric(m[30,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[22,6]),0,as.numeric(m[22,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[23,6]),0,as.numeric(m[23,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[24,6]),0,as.numeric(m[24,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[31,10]),0,as.numeric(m[31,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[32,10]),0,as.numeric(m[32,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[33,10]),0,as.numeric(m[33,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA, NA, NA,
-equivbs(as.numeric(m[25,6]),0,as.numeric(m[25,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[26,6]),0,as.numeric(m[26,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[27,6]),0,as.numeric(m[27,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[34,10]),0,as.numeric(m[34,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[35,10]),0,as.numeric(m[35,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[36,10]),0,as.numeric(m[36,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[28,6]),0,as.numeric(m[28,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[29,6]),0,as.numeric(m[29,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[30,6]),0,as.numeric(m[30,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[37,10]),0,as.numeric(m[37,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[38,10]),0,as.numeric(m[38,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[39,10]),0,as.numeric(m[39,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
 NA,
-equivbs(as.numeric(m[31,6]),0,as.numeric(m[31,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[32,6]),0,as.numeric(m[32,7]),0,c(-.1, .1),alpha = .05)$`CI`[2],
-equivbs(as.numeric(m[33,6]),0,as.numeric(m[33,7]),0,c(-.1, .1),alpha = .05)$`CI`[2]))
+equivbs(as.numeric(m[40,10]),0,as.numeric(m[40,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[41,10]),0,as.numeric(m[41,11]),0,c(-.1, .1),alpha = .05)$`CI`[2],
+equivbs(as.numeric(m[42,10]),0,as.numeric(m[42,11]),0,c(-.1, .1),alpha = .05)$`CI`[2]))
 
 x<-as_factor(c("Bilingual Use -> SSRTa", "Bilingual Degree -> SSRTa", "Bilingual Status -> SSRTa", "",
-"Bilingual Use -> Card Sorta", "Bilingual Degree -> Card Sorta", "Bilingual Status -> Card Sorta", " ", "Bilingual Use -> Flankera", "Bilingual Degree -> Flankera",
-"Bilingual Status -> Flankera", "   ", "GAMM Model (w covariates 1-9)", "    ",
-"Bilingual Use -> SSRTb", "Bilingual Degree -> SSRTb", "Bilingual Status -> SSRTb", "     ",
-"Bilingual Use -> Card Sortb", "Bilingual Degree -> Card Sortb", "Bilingual Status -> Card Sortb", "      ",
-"Bilingual Use -> Flankerb", "Bilingual Degree -> Flankerb", "Bilingual Status -> Flankerb", "       ", "Bilingual Use -> Vocabularyb",
-"Bilingual Degree -> Vocabularyb", "Bilingual Status -> Vocabularyb", "        ","GAMM Model (w covariates 1-8)",
-"         ", "Bilingual Use -> SSRTc", "Bilingual Degree -> SSRTc", "Bilingual Status -> SSRTc", "          ", "Bilingual Use -> Card Sortc", "Bilingual Degree -> Card Sortc",
-"Bilingual Status -> Card Sortc", "           ",
-"Bilingual Use -> Flankerc", "Bilingual Degree -> Flankerc", "Bilingual Status -> Flankerc", "            ", "Bilingual Use -> Vocabularyc", "Bilingual Degree -> Vocabularyc",
-"Bilingual Status -> Vocabularyc","             ", "GAMM Model (w covariates 1-6)", "              ", "Bilingual Use -> SSRTd", "Bilingual Degree -> SSRTd",
-"Bilingual Status -> SSRTd", "               ", "Bilingual Use -> Card Sortd", "Bilingual Degree -> Card Sortd",
-"Bilingual Status -> Card Sortd", "                ", "Bilingual Use -> Flankerd", "Bilingual Degree -> Flankerd", "Bilingual Status -> Flankerd",
-"                  ", "Bilingual Use -> Vocabularyd", "Bilingual Degree -> Vocabularyd", "Bilingual Status -> Vocabularyd", "                   ", "OLS Model"))
+"Bilingual Use -> Card Sorta", "Bilingual Degree -> Card Sorta", "Bilingual Status -> Card Sorta", " ",
+"Bilingual Use -> Flankera", "Bilingual Degree -> Flankera","Bilingual Status -> Flankera","  ",
+"GAMM Model (w covariates 1-9)", "   ",
+"Bilingual Use -> SSRTb", "Bilingual Degree -> SSRTb", "Bilingual Status -> SSRTb", "    ",
+"Bilingual Use -> Card Sortb", "Bilingual Degree -> Card Sortb", "Bilingual Status -> Card Sortb", "     ",
+"Bilingual Use -> Flankerb", "Bilingual Degree -> Flankerb", "Bilingual Status -> Flankerb", "      ",
+"Bilingual Use -> Vocabularyb","Bilingual Degree -> Vocabularyb", "Bilingual Status -> Vocabularyb", "       ",
+"GAMM Model (w covariates 1-8)","        ",
+"Bilingual Use -> SSRTc", "Bilingual Degree -> SSRTc", "Bilingual Status -> SSRTc", "         ",
+"Bilingual Use -> Card Sortc", "Bilingual Degree -> Card Sortc","Bilingual Status -> Card Sortc", "          ",
+"Bilingual Use -> Flankerc", "Bilingual Degree -> Flankerc", "Bilingual Status -> Flankerc", "            ",
+"Bilingual Use -> Vocabularyc", "Bilingual Degree -> Vocabularyc","Bilingual Status -> Vocabularyc","              ",
+"GAMM Model (w covariates 1-6)","               ",
+"Bilingual Use -> SSRTd", "Bilingual Degree -> SSRTd","Bilingual Status -> SSRTd", "                ",
+"Bilingual Use -> Card Sortd", "Bilingual Degree -> Card Sortd","Bilingual Status -> Card Sortd", "                 ",
+"Bilingual Use -> Flankerd", "Bilingual Degree -> Flankerd", "Bilingual Status -> Flankerd","                  ",
+"GAMM Model (w Covariate 9 only)", "                   ",
+"Bilingual Use -> SSRTe", "Bilingual Degree -> SSRTe","Bilingual Status -> SSRTe", "                    ",
+"Bilingual Use -> Card Sorte", "Bilingual Degree -> Card Sorte","Bilingual Status -> Card Sorte", "                     ",
+"Bilingual Use -> Flankere", "Bilingual Degree -> Flankere", "Bilingual Status -> Flankere","                      ",
+"Bilingual Use -> Vocabularye", "Bilingual Degree -> Vocabularye", "Bilingual Status -> Vocabularye", "                       ",
+"OLS Model"))
 
 d<-data.frame(x, y, ylo, yhi)
 
@@ -514,19 +659,28 @@ credplot.gg <- function(d){
  coord_flip()+scale_y_continuous(breaks=seq(-4,4,.1)) +
  ylab('Effect Size') +
  scale_x_discrete('Predictor -> Outcome', waiver(), labels = c("Bilingual Use -> SSRT", "Bilingual Degree -> SSRT", "Bilingual Status -> SSRT", "",
-"Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort", "Bilingual Status -> Card Sort", " ", "Bilingual Use -> Flanker", "Bilingual Degree -> Flanker",
-"Bilingual Status -> Flanker", "   ", "GAMM Model (w covariates 1-9)", "    ",
-"Bilingual Use -> SSRT", "Bilingual Degree -> SSRT", "Bilingual Status -> SSRT", "     ",
-"Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort", "Bilingual Status -> Card Sort", "      ",
-"Bilingual Use -> Flanker", "Bilingual Degree -> Flanker", "Bilingual Status -> Flanker", "       ", "Bilingual Use -> Vocabulary",
-"Bilingual Degree -> Vocabulary", "Bilingual Status -> Vocabulary", "        ","GAMM Model (w covariates 1-8)",
-"         ", "Bilingual Use -> SSRT", "Bilingual Degree -> SSRT", "Bilingual Status -> SSRT", "          ", "Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort",
-"Bilingual Status -> Card Sort", "           ",
-"Bilingual Use -> Flanker", "Bilingual Degree -> Flanker", "Bilingual Status -> Flanker", "            ", "Bilingual Use -> Vocabulary", "Bilingual Degree -> Vocabulary",
-"Bilingual Status -> Vocabulary","             ", "GAMM Model (w covariates 1-6)", "              ", "Bilingual Use -> SSRT", "Bilingual Degree -> SSRT",
-"Bilingual Status -> SSRT", "               ", "Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort",
-"Bilingual Status -> Card Sort", "                ", "Bilingual Use -> Flanker", "Bilingual Degree -> Flanker", "Bilingual Status -> Flanker",
-"                  ", "Bilingual Use -> Vocabulary", "Bilingual Degree -> Vocabulary", "Bilingual Status -> Vocabulary", "                   ", "OLS Model")) +
+ "Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort", "Bilingual Status -> Card Sort", " ",
+ "Bilingual Use -> Flanker", "Bilingual Degree -> Flanker","Bilingual Status -> Flanker","  ",
+ "GAMM Model (w covariates 1-9)", "   ",
+ "Bilingual Use -> SSRT", "Bilingual Degree -> SSRT", "Bilingual Status -> SSRT", "    ",
+ "Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort", "Bilingual Status -> Card Sort", "     ",
+ "Bilingual Use -> Flanker", "Bilingual Degree -> Flanker", "Bilingual Status -> Flanker", "      ",
+ "Bilingual Use -> Vocabulary","Bilingual Degree -> Vocabulary", "Bilingual Status -> Vocabulary", "       ",
+ "GAMM Model (w covariates 1-8)","        ",
+ "Bilingual Use -> SSRT", "Bilingual Degree -> SSRT", "Bilingual Status -> SSRT", "         ",
+ "Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort","Bilingual Status -> Card Sort", "          ",
+ "Bilingual Use -> Flanker", "Bilingual Degree -> Flanker", "Bilingual Status -> Flanker", "            ",
+ "Bilingual Use -> Vocabulary", "Bilingual Degree -> Vocabulary","Bilingual Status -> Vocabulary","              ",
+ "GAMM Model (w covariates 1-6)","               ",
+ "Bilingual Use -> SSRT", "Bilingual Degree -> SSRT","Bilingual Status -> SSRT", "                ",
+ "Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort","Bilingual Status -> Card Sort", "                 ",
+ "Bilingual Use -> Flanker", "Bilingual Degree -> Flanker", "Bilingual Status -> Flanker","                  ",
+ "GAMM Model (w Covariate 9 only)", "                   ",
+ "Bilingual Use -> SSRT", "Bilingual Degree -> SSRT","Bilingual Status -> SSRT", "                    ",
+ "Bilingual Use -> Card Sort", "Bilingual Degree -> Card Sort","Bilingual Status -> Card Sort", "                     ",
+ "Bilingual Use -> Flanker", "Bilingual Degree -> Flanker", "Bilingual Status -> Flanker","                      ",
+ "Bilingual Use -> Vocabulary", "Bilingual Degree -> Vocabulary", "Bilingual Status -> Vocabulary", "                       ",
+ "OLS Model")) +
  ggtitle(expression(paste("Effect Sizes (", italic(beta), ") and Confidence Intervals Plotted Against the Interval of Equivalence"))) +
  theme(panel.border = element_blank(),
           panel.grid.major = element_blank(),
@@ -536,7 +690,7 @@ credplot.gg <- function(d){
  return(p)
 }
 
-tiff("figure.tiff", units = 'in', width = 20, height = 14, res = 200, compression = "lzw")
+tiff("figure1_newestestest.tiff", units = 'in', width = 20, height = 20, res = 200, compression = "lzw")
 credplot.gg(d)
 dev.off()
 
@@ -545,7 +699,9 @@ dev.off()
 ## Bootstrapping at Small Samples	  ##
 ########################################
 ########################################
-
+abcd_subset<-cbind.data.frame(abcd_subset, random_data=rnorm(4524, mean=0, sd=1))
+ind_iv = c(which(names(abcd_subset)=="bilingual_status"), which(names(abcd_subset)=="bilingual_degree"), which(names(abcd_subset)=="bilingual_use"))
+ind_dv = c(which(names(abcd_subset)=="nihtbx_picvocab_uncorrected"),which(names(abcd_subset)=="nihtbx_flanker_uncorrected"),which(names(abcd_subset)=="nihtbx_cardsort_uncorrected"), which(names(abcd_subset)=="SSRTr"), which(names(abcd_subset)=="random_data"))
 
 samp.n = 30
 mat_row = 1
@@ -576,19 +732,21 @@ for(b in 1:nboot){
 		} # go to the next step in the loop
 		mat_col = 1
 		mat_row = mat_row + 1 # step to the next matrix column
-	colnames(mp)<-c("Bilingual Status --> Vocabulary", "Bilingual Degree --> Vocabulary", "Bilingual Use --> Vocabulary", "Bilingual Status --> Flanker", "Bilingual Degree --> Flanker", "Bilingual Use --> Flanker", "Bilingual Status --> Card Sort", "Bilingual Degree --> Card Sort", "Bilingual Use --> Card Sort", "Bilingual Status --> SSRT", "Bilingual Degree --> SSRT", "Bilingual Use --> SSRT")
+	colnames(mp)<-c("Bilingual Status --> Vocabulary", "Bilingual Degree --> Vocabulary", "Bilingual Use --> Vocabulary", "Bilingual Status --> Flanker",
+  "Bilingual Degree --> Flanker", "Bilingual Use --> Flanker", "Bilingual Status --> Card Sort", "Bilingual Degree --> Card Sort",
+  "Bilingual Use --> Card Sort", "Bilingual Status --> SSRT", "Bilingual Degree --> SSRT", "Bilingual Use --> SSRT", "Bilingual Status --> Random Data",
+  "Bilingual Degree --> Random Data", "Bilingual Use --> Random Data")
 	#print(mp)
 	#print(mdf)
 }
 
-
 for.hist <- melt(mp)
 for.hist<-for.hist[which(for.hist$value < .05),]
 
-tiff("figure.tiff", units = 'in', width = 12, height = 10, res = 300, compression = "lzw")
+tiff("figure_hist.tiff", units = 'in', width = 12, height = 10, res = 300, compression = "lzw")
 #par(mar = c(5, 5, 8, 8), xpd=FALSE) #set figure boundaries
 ggplot(for.hist,aes(x = for.hist$value)) +
-    facet_wrap(~Var2,scales = "free_x", nrow = 4, ncol = 3) + theme_bw() + geom_histogram(
+    facet_wrap(~Var2,scales = "free_x", nrow = 5, ncol = 3) + theme_bw() + geom_histogram(
     col = "grey", bins = 20,
     fill="light blue") +
     labs(x=expression(paste(italic("p "), "Value")), y="Frequency") +