experimental comparison_aug14

1. Need bayesfactor means from all the runs for all, west, east: bfmeansWithScafPosition_all, east and west

2. Combined these bayesfactors with the survival and weight data using bf_sruv_combine.py scripts and got the Outcombine_surv_bf files for all,east and west. After combining these files I have 95217 SNPs for survival data and 97082 SNPs for weight data

3. Need allele frequency averages from all, west, east : these were the main population allele freq files from which I generated allele counts for the baypass analysis. For this analysis I just took averages of the allele frequencies from all, west and east populations in R and created the all_afmeans. west_afmeans and east_afmeans files

4. combine scaffold, positions, all,east and west bfmeans and all, east and west afmeans to get a final datatable.

FILE LOG:

1. bfmeansWithScafpositions_all, east, west (206,029 SNPs)- bayPass file for all, east and western populations bayesfactor means from all baypass runs with scaffold and positions

2. effectsSurvProbitTotal.txt = file with survival data from the experiments

3. effectsWgtTotal.txt = file with weight data from the experiments

4. bf_survivaldata_combine.py #used to combine survival data with bayesfactors from baypass

5. bf_wghtdata_combine.py #used to combine weight data with bayesfactors from baypass

6. all_popallelefreqs.txt, west_popallelefreqs.txt, east_popallelefreqs.txt = population allele frequencies for each of these sets

7. all_afmeans,east_afmeans,west_afmeans : allele frequency means from the previous files for each of these runs

8. final_datatable = (206029 SNPs) this is the final table used in the analysis. contains the following:scaffold,position,all_bf,west_bf,east_bf,all_af,west_af,east_af,gla_af,sla_affinal_survdata_combine.py #main script for combining survival data and bayesfactors in final_datatable for all the test sets

9. final_wghtdata_combine.py #main script for combining weight data and bayesfactors in final_datatable for all the test sets

10. survfile_tab.txt - Tabdelimited files for running the python script

11. wghtfile_tab.txt - same as above

12. Outcombine_surv_final.csv #FINAL FILE USED FOR ANALYSIS

13. Outcombine_wgt_final.csv #FINAL FILE USED FOR ANALYSIS

14. quantilesSource_glawght.R,quantilesSource_slasurv.R, quantilesSource_slawght.R = source code files to do quantiles for Experimental comparisons.

15. effSigns_top58Snps.R = code for getting effect signs for top 58 shared snps

Here is the link to the table with results: https://docs.google.com/spreadsheets/d/11MNP_wrZeUs0ynDdeCMNGVrM-yUHvsFSMZ8_DrXAIzI/edit#gid=0

survcounts = c(9,15,8,4,17,17,13,9,18,15,5,8)

pdf("survival_shared.pdf", width = 10, height = 10)

par(mar=c(6,8,8,6),mgp=c(4,1,0))

barplot(survcounts, col = c("darkgreen","darkolivegreen1","slateblue1","skyblue","darkgreen","darkolivegreen1","slateblue1","skyblue","darkgreen","darkolivegreen1","slateblue1","skyblue"), xlab = "Population and host plants", ylab = "Number of observed common SNVs", cex.lab = 2, ylim = c(0,20))

dev.off()

wgtcounts = c(12,16,7,11,12,13,9,18,13,14,12,8)

pdf("wght_shared.pdf", width = 10, height = 10)

par(mar=c(6,8,8,6),mgp=c(4,1,0))

barplot(survcounts, col = c("darkgreen","darkolivegreen1","slateblue1","skyblue","darkgreen","darkolivegreen1","slateblue1","skyblue","darkgreen","darkolivegreen1","slateblue1","skyblue"), xlab = "Population and host plants", ylab = "Number of observed common SNVs", cex.lab = 2, ylim = c(0,20))

dev.off()

On August 14th

## I have switched east and west in the main files. I redid this again because before I had not taken the absolute values for the effect sizes.

################ SURVIVAL #####################################

#glaMs filtering#############

#get absolute values for effect sizes for glams and remove rows with NA

surv_glams = survival[!is.na(survival$glaMs) & !is.na(survival$gla_af),]

glams_fil = abs(surv_glams$glaMs)

#glm for effect sizes

fit_glams = glm(glams_fil ~ surv_glams$gla_af, family = gaussian)

##all populations

#glm for bayes factors

fit_all = glm(surv_glams$all_bf ~ surv_glams$all_af, family=gaussian)

ms = fit_glams$residuals

bf = fit_all$residuals

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(surv_glams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glams)[1],0)

msnull<-sample(1:dim(surv_glams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

#raw values

ms_raw = glams_fil

bf_raw = surv_glams$all_bf

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glams)[1],0)

msnull<-sample(1:dim(surv_glams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

############### west ####################

fit_glams = glm(glams_fil ~ surv_glams$gla_af, family = gaussian)

fit_west = glm(surv_glams$west_bf ~ surv_glams$west_af, family = gaussian)

ms = fit_glams$residuals

bf = fit_west$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(surv_glams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glams)[1],0)

msnull<-sample(1:dim(surv_glams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

#raw values

ms_raw = glams_fil

bf_raw = surv_glams$west_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glams)[1],0)

msnull<-sample(1:dim(surv_glams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

########## east glaMs ##################

#raw values

ms_raw = glams_fil

bf_raw = surv_glams$east_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glams)[1],0)

msnull<-sample(1:dim(surv_glams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

#glm residuals

fit_glams = glm(glams_fil ~ surv_glams$gla_af, family = gaussian)

fit_east = glm(surv_glams$east_bf ~ surv_glams$east_af, family = gaussian)

ms = fit_glams$residuals

bf = fit_east$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(surv_glams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glams)[1],0)

msnull<-sample(1:dim(surv_glams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

################glaAc######################################

#get absolute values for effect sizes for glaac and remove rows with NA

surv_glaac = survival[!is.na(survival$glaAc) & !is.na(survival$gla_af),]

glaac_fil = abs(surv_glaac$glaAc)

#raw values

ac_raw = glaac_fil

bf_raw = surv_glaac$all_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glaac)[1],0)

acnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#glm for residuals

fit_glaac = glm(glaac_fil ~ surv_glaac$gla_af, family = gaussian)

##all populations

#glm for bayes factors

fit_all = glm(surv_glaac$all_bf ~ surv_glaac$all_af, family=gaussian)

ac = fit_glaac$residuals

bf = fit_all$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(surv_glaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glaac)[1],0)

acnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

############### west ####################

#raw values

ac_raw = glaac_fil

bf_raw = surv_glaac$west_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glaac)[1],0)

acnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#residuals

fit_glaac = glm(glaac_fil ~ surv_glaac$gla_af, family = gaussian)

fit_west = glm(surv_glaac$west_bf ~ surv_glaac$west_af, family = gaussian)

ac = fit_glaac$residuals

bf = fit_west$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(surv_glaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glaac)[1],0)

acnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

########## east glaac ##################

#raw values

ac_raw = glaac_fil

bf_raw = surv_glaac$east_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glaac)[1],0)

acnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#glm residuals

fit_glaac = glm(glaac_fil ~ surv_glaac$gla_af, family = gaussian)

fit_east = glm(surv_glaac$east_bf ~ surv_glaac$east_af, family = gaussian)

ac = fit_glaac$residuals

bf = fit_east$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(surv_glaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_glaac)[1],0)

acnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

##########slaAc

################slaAc######################################

#get absolute values for effect sizes for slaac and remove rows with NA

surv_slaac = survival[!is.na(survival$slaAc) & !is.na(survival$sla_af),]

slaac_fil = abs(surv_slaac$slaAc)

#raw values

ac_raw = slaac_fil

bf_raw = surv_slaac$all_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slaac)[1],0)

acnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#glm for residuals

fit_slaac = glm(slaac_fil ~ surv_slaac$sla_af, family = gaussian)

##all populations

#glm for bayes factors

fit_all = glm(surv_slaac$all_bf ~ surv_slaac$all_af, family=gaussian)

ac = fit_slaac$residuals

bf = fit_all$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(surv_slaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slaac)[1],0)

acnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

############### west ####################

#raw values

ac_raw = slaac_fil

bf_raw = surv_slaac$west_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slaac)[1],0)

acnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#residuals

fit_slaac = glm(slaac_fil ~ surv_slaac$sla_af, family = gaussian)

fit_west = glm(surv_slaac$west_bf ~ surv_slaac$west_af, family = gaussian)

ac = fit_slaac$residuals

bf = fit_west$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(surv_slaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slaac)[1],0)

acnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

########## east slaac ##################

#raw values

ac_raw = slaac_fil

bf_raw = surv_slaac$east_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slaac)[1],0)

acnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#glm residuals

fit_slaac = glm(slaac_fil ~ surv_slaac$sla_af, family = gaussian)

fit_east = glm(surv_slaac$east_bf ~ surv_slaac$east_af, family = gaussian)

ac = fit_slaac$residuals

bf = fit_east$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(surv_slaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slaac)[1],0)

acnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

################slaMs######################################

#get absolute values for effect sizes for slams and remove rows with NA

surv_slams = survival[!is.na(survival$slaMs) & !is.na(survival$sla_af),]

slams_fil = abs(surv_slams$slaMs)

#raw values

ms_raw = slams_fil

bf_raw = surv_slams$all_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slams)[1],0)

msnull<-sample(1:dim(surv_slams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

#glm for residuals

fit_slams = glm(slams_fil ~ surv_slams$sla_af, family = gaussian)

##all populations

#glm for bayes fmstors

fit_all = glm(surv_slams$all_bf ~ surv_slams$all_af, family=gaussian)

ms = fit_slams$residuals

bf = fit_all$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(surv_slams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slams)[1],0)

msnull<-sample(1:dim(surv_slams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

############### west ####################

#raw values

ms_raw = slams_fil

bf_raw = surv_slams$west_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slams)[1],0)

msnull<-sample(1:dim(surv_slams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

#residuals

fit_slams = glm(slams_fil ~ surv_slams$sla_af, family = gaussian)

fit_west = glm(surv_slams$west_bf ~ surv_slams$west_af, family = gaussian)

ms = fit_slams$residuals

bf = fit_west$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(surv_slams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slams)[1],0)

msnull<-sample(1:dim(surv_slams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

########## east slams ##################

#raw values

ms_raw = slams_fil

bf_raw = surv_slams$east_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slams)[1],0)

msnull<-sample(1:dim(surv_slams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

#glm residuals

fit_slams = glm(slams_fil ~ surv_slams$sla_af, family = gaussian)

fit_east = glm(surv_slams$east_bf ~ surv_slams$east_af, family = gaussian)

ms = fit_slams$residuals

bf = fit_east$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(surv_slams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(surv_slams)[1],0)

msnull<-sample(1:dim(surv_slams)[1],n,replace=F)

bfnull<-sample(1:dim(surv_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

################ WEIGHT #####################################

weight = read.table("Outcombine_wgt_final.csv", header=TRUE)

################## glaMs #############

#get absolute values for effect sizes for glams and remove rows with NA

wght_glams = weight[!is.na(weight$glaMs) & !is.na(weight$gla_af),]

glams_fil = abs(wght_glams$glaMs)

#glm for effect sizes

fit_glams = glm(glams_fil ~ wght_glams$gla_af, family = gaussian)

############## all glaMs ###################

#raw values

ms_raw = glams_fil

bf_raw = wght_glams$all_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glams)[1],0)

msnull<-sample(1:dim(wght_glams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

#glm residuals

fit_all = glm(wght_glams$all_bf ~ wght_glams$all_af, family=gaussian)

ms = fit_glams$residuals

bf = fit_all$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(wght_glams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glams)[1],0)

msnull<-sample(1:dim(wght_glams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

############### west glaMs ####################

#raw values

ms_raw = glams_fil

bf_raw = wght_glams$west_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glams)[1],0)

msnull<-sample(1:dim(wght_glams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

#residuals

fit_glams = glm(glams_fil ~ wght_glams$gla_af, family = gaussian)

fit_west = glm(wght_glams$west_bf ~ wght_glams$west_af, family = gaussian)

ms = fit_glams$residuals

bf = fit_west$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(wght_glams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glams)[1],0)

msnull<-sample(1:dim(wght_glams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

########## east glaMs ##################

#raw values

ms_raw = glams_fil

bf_raw = wght_glams$east_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glams)[1],0)

msnull<-sample(1:dim(wght_glams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

#glm residuals

fit_glams = glm(glams_fil ~ wght_glams$gla_af, family = gaussian)

fit_east = glm(wght_glams$east_bf ~ wght_glams$east_af, family = gaussian)

ms = fit_glams$residuals

bf = fit_east$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(wght_glams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glams)[1],0)

msnull<-sample(1:dim(wght_glams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

p

xfold

obs

################glaAc######################################

#get absolute values for effect sizes for glaac and remove rows with NA

wght_glaac = weight[!is.na(weight$glaAc) & !is.na(weight$gla_af),]

glaac_fil = abs(wght_glaac$glaAc)

#raw values

ac_raw = glaac_fil

bf_raw = wght_glaac$all_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glaac)[1],0)

acnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#glm for residuals

fit_glaac = glm(glaac_fil ~ wght_glaac$gla_af, family = gaussian)

##all populations

#glm for bayes factors

fit_all = glm(wght_glaac$all_bf ~ wght_glaac$all_af, family=gaussian)

ac = fit_glaac$residuals

bf = fit_all$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(wght_glaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glaac)[1],0)

acnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

############### west ####################

#raw values

ac_raw = glaac_fil

bf_raw = wght_glaac$west_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glaac)[1],0)

acnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#residuals

fit_glaac = glm(glaac_fil ~ wght_glaac$gla_af, family = gaussian)

fit_west = glm(wght_glaac$west_bf ~ wght_glaac$west_af, family = gaussian)

ac = fit_glaac$residuals

bf = fit_west$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(wght_glaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glaac)[1],0)

acnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

########## east glaac ##################

#raw values

ac_raw = glaac_fil

bf_raw = wght_glaac$east_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glaac)[1],0)

acnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#glm residuals

fit_glaac = glm(glaac_fil ~ wght_glaac$gla_af, family = gaussian)

fit_east = glm(wght_glaac$east_bf ~ wght_glaac$east_af, family = gaussian)

ac = fit_glaac$residuals

bf = fit_east$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(wght_glaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_glaac)[1],0)

acnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_glaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

##########slaAc

################slaAc######################################

#get absolute values for effect sizes for slaac and remove rows with NA

wght_slaac = weight[!is.na(weight$slaAc) & !is.na(weight$sla_af),]

slaac_fil = abs(wght_slaac$slaAc)

#raw values

ac_raw = slaac_fil

bf_raw = wght_slaac$all_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slaac)[1],0)

acnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#glm for residuals

fit_slaac = glm(slaac_fil ~ wght_slaac$sla_af, family = gaussian)

##all populations

#glm for bayes factors

fit_all = glm(wght_slaac$all_bf ~ wght_slaac$all_af, family=gaussian)

ac = fit_slaac$residuals

bf = fit_all$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(wght_slaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slaac)[1],0)

acnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

############### west ####################

#raw values

ac_raw = slaac_fil

bf_raw = wght_slaac$west_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slaac)[1],0)

acnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#residuals

fit_slaac = glm(slaac_fil ~ wght_slaac$sla_af, family = gaussian)

fit_west = glm(wght_slaac$west_bf ~ wght_slaac$west_af, family = gaussian)

ac = fit_slaac$residuals

bf = fit_west$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(wght_slaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slaac)[1],0)

acnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

########## east slaac ##################

#raw values

ac_raw = slaac_fil

bf_raw = wght_slaac$east_bf

cor.test(ac_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac_raw > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slaac)[1],0)

acnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

#glm residuals

fit_slaac = glm(slaac_fil ~ wght_slaac$sla_af, family = gaussian)

fit_east = glm(wght_slaac$east_bf ~ wght_slaac$east_af, family = gaussian)

ac = fit_slaac$residuals

bf = fit_east$residuals

cor.test(ac,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

actop = which(ac > quantile(ac,probs=0.99,na.rm=T))

obs = sum(actop %in% bftop)

dim(wght_slaac)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slaac)[1],0)

acnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slaac)[1],n,replace=F)

null[i]<-sum(acnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(actop %in% bftop)

p

xfold

obs

################slaMs######################################

#get absolute values for effect sizes for slams and remove rows with NA

wght_slams = weight[!is.na(weight$slaMs) & !is.na(weight$sla_af),]

slams_fil = abs(wght_slams$slaMs)

#raw values

ms_raw = slams_fil

bf_raw = wght_slams$all_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slams)[1],0)

msnull<-sample(1:dim(wght_slams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

#glm for residuals

fit_slams = glm(slams_fil ~ wght_slams$sla_af, family = gaussian)

##all populations

#glm for bayes fmstors

fit_all = glm(wght_slams$all_bf ~ wght_slams$all_af, family=gaussian)

ms = fit_slams$residuals

bf = fit_all$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(wght_slams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slams)[1],0)

msnull<-sample(1:dim(wght_slams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

############### west ####################

#raw values

ms_raw = slams_fil

bf_raw = wght_slams$west_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slams)[1],0)

msnull<-sample(1:dim(wght_slams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

#residuals

fit_slams = glm(slams_fil ~ wght_slams$sla_af, family = gaussian)

fit_west = glm(wght_slams$west_bf ~ wght_slams$west_af, family = gaussian)

ms = fit_slams$residuals

bf = fit_west$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(wght_slams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slams)[1],0)

msnull<-sample(1:dim(wght_slams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

########## east slams ##################

#raw values

ms_raw = slams_fil

bf_raw = wght_slams$east_bf

cor.test(ms_raw,bf_raw)

bftop = which(bf_raw > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms_raw > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slams)[1],0)

msnull<-sample(1:dim(wght_slams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold

#glm residuals

fit_slams = glm(slams_fil ~ wght_slams$sla_af, family = gaussian)

fit_east = glm(wght_slams$east_bf ~ wght_slams$east_af, family = gaussian)

ms = fit_slams$residuals

bf = fit_east$residuals

cor.test(ms,bf)

#quantiles

bftop = which(bf > quantile(bf,probs=0.99,na.rm=T))

mstop = which(ms > quantile(ms,probs=0.99,na.rm=T))

obs = sum(mstop %in% bftop)

dim(wght_slams)[1]

#randomization

null<-rep(NA,10000)

for(i in 1:10000){

n<-round(0.01*dim(wght_slams)[1],0)

msnull<-sample(1:dim(wght_slams)[1],n,replace=F)

bfnull<-sample(1:dim(wght_slams)[1],n,replace=F)

null[i]<-sum(msnull %in% bfnull)

}

p = mean(null >= obs)

xfold = obs / mean(null)

obs = sum(mstop %in% bftop)

obs

p

xfold