Need bayesfactor means from all the runs for all, west, east: bfmeansWithScafPosition_all, east and west
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
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
combine scaffold, positions, all,east and west bfmeans and all, east and west afmeans to get a final datatable.
FILE LOG:
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
effectsSurvProbitTotal.txt = file with survival data from the experiments
effectsWgtTotal.txt = file with weight data from the experiments
bf_survivaldata_combine.py #used to combine survival data with bayesfactors from baypass
bf_wghtdata_combine.py #used to combine weight data with bayesfactors from baypass
all_popallelefreqs.txt, west_popallelefreqs.txt, east_popallelefreqs.txt = population allele frequencies for each of these sets
all_afmeans,east_afmeans,west_afmeans : allele frequency means from the previous files for each of these runs
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
final_wghtdata_combine.py #main script for combining weight data and bayesfactors in final_datatable for all the test sets
survfile_tab.txt - Tabdelimited files for running the python script
wghtfile_tab.txt - same as above
Outcombine_surv_final.csv #FINAL FILE USED FOR ANALYSIS
Outcombine_wgt_final.csv #FINAL FILE USED FOR ANALYSIS
quantilesSource_glawght.R,quantilesSource_slasurv.R, quantilesSource_slawght.R = source code files to do quantiles for Experimental comparisons.
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