For these results I have correctly switched east and west data. So these are correct data. However, in the code its still in the old way--east is west.
*********** survival ********************************
#all slaMs
#all glaAc
#west slaAc
#west slaMs
#east slaAc( this is coded as west in the code. I called east as west and then ZG corrected me at evolution)
#east slaMs
#results
#all glaMs
#all glaAc
#east_glaMs
#east_glaAc
#west_glaMs
#west_glaAc
*******************************************************
##plots
#enrichments for all ms
enrich_ms = c(1.121,1.213,1.16,1.31,1.34,1.33,1.29,1.25,1.22,1.19)
setEPS()
postscript("allms_Exp_enrichments.eps", width=12, height=9)
par(mar=c(7,8,8,7), mgp=c(4,1,0))
barplot(enrich_ms, xlab = "top SNPs %", ylab = "xfold enrichment", ylim = c(0,1.5), col = "mistyrose4", cex.lab=2, cex.axis=1.5, names.arg=c("1%","2%","3%","4%","5%","6%","7%","8%","9%" ,"10%"), cex.names=0.8)
dev.off()
#enrichments for all ac
enrich_ac = c(1.820,1.428,1.502,1.467,1.464,1.431,1.337,1.301,1.259,1.266)
setEPS()
postscript("allac_Exp_enrichments.eps", width=12, height=9)
par(mar=c(7,8,8,7), mgp=c(4,1,0))
barplot(enrich_ac, xlab = "top SNPs %", ylab = "xfold enrichment", ylim = c(0,2), col = "mistyrose4", cex.lab=2, cex.axis=1.5, names.arg=c("1%","2%","3%","4%","5%","6%","7%","8%","9%" ,"10%"), cex.names=0.8)
dev.off()
#all ac and ms
#glm and filtering
fil_glams = survival[!is.na(survival$glaMs) & !is.na(survival$gla_af),]
fit_glams = glm(fil_glams$glaMs ~ fil_glams$gla_af, family = gaussian)
fit_all = glm(fil_glams$all_bf ~ fil_glams$all_af , family=gaussian)
#doing correlations and randomization ####
ms = fit_glams$residuals
bf = fit_all$residuals
bndsb = quantile(bf,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndsm = quantile(ms,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
mstop1 = which(ms > quantile(ms,probs=0.99,na.rm=T))
dim(fil_glams)[1]
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.01*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop1 %in% bftop1)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop2 = which(bf > quantile(bf,probs=0.98,na.rm=T))
mstop2 = which(ms > quantile(ms,probs=0.98,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.02*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop2 %in% bftop2)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop3 = which(bf > quantile(bf,probs=0.97,na.rm=T))
mstop3 = which(ms > quantile(ms,probs=0.97,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.03*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop3 %in% bftop3)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop4 = which(bf > quantile(bf,probs=0.96,na.rm=T))
mstop4 = which(ms > quantile(ms,probs=0.96,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.04*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop4 %in% bftop4)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop5 = which(bf > quantile(bf,probs=0.95,na.rm=T))
mstop5 = which(ms > quantile(ms,probs=0.95,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.05*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop5 %in% bftop5)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop6 = which(bf > quantile(bf,probs=0.94,na.rm=T))
mstop6 = which(ms > quantile(ms,probs=0.94,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.06*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop6 %in% bftop6)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop7 = which(bf > quantile(bf,probs=0.93,na.rm=T))
mstop7 = which(ms > quantile(ms,probs=0.93,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.07*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop7 %in% bftop7)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop8 = which(bf > quantile(bf,probs=0.92,na.rm=T))
mstop8 = which(ms > quantile(ms,probs=0.92,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.08*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop8 %in% bftop8)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop9 = which(bf > quantile(bf,probs=0.91,na.rm=T))
mstop9 = which(ms > quantile(ms,probs=0.91,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.09*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop9 %in% bftop9)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop10 = which(bf > quantile(bf,probs=0.9,na.rm=T))
mstop10 = which(ms > quantile(ms,probs=0.9,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.1*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop10 %in% bftop10)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
##all glaAc
#glm and filtering
fil_glaac = survival[!is.na(survival$glaAc) & !is.na(survival$gla_af),]
fit_glaac = glm(fil_glaac$glaAc ~ fil_glaac$gla_af, family = gaussian)
fit_all = glm(fil_glaac$all_bf ~ fil_glaac$all_af , family=gaussian)
#doing correlations and randomization ####
ac = fit_glaac$residuals
bf = fit_all$residuals
bndsb = quantile(bf,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndsa = quantile(ac,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
actop1 = which(ac > quantile(ac,probs=0.99,na.rm=T))
dim(fil_glaac)[1]
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.01*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop1 %in% bftop1)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop2 = which(bf > quantile(bf,probs=0.98,na.rm=T))
actop2 = which(ac > quantile(ac,probs=0.98,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.02*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop2 %in% bftop2)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop3 = which(bf > quantile(bf,probs=0.97,na.rm=T))
actop3 = which(ac > quantile(ac,probs=0.97,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.03*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop3 %in% bftop3)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop4 = which(bf > quantile(bf,probs=0.96,na.rm=T))
actop4 = which(ac > quantile(ac,probs=0.96,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.04*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop4 %in% bftop4)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop5 = which(bf > quantile(bf,probs=0.95,na.rm=T))
actop5 = which(ac > quantile(ac,probs=0.95,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.05*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop5 %in% bftop5)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop6 = which(bf > quantile(bf,probs=0.94,na.rm=T))
actop6 = which(ac > quantile(ac,probs=0.94,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.06*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop6 %in% bftop6)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop7 = which(bf > quantile(bf,probs=0.93,na.rm=T))
actop7 = which(ac > quantile(ac,probs=0.93,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.07*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop7 %in% bftop7)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop8 = which(bf > quantile(bf,probs=0.92,na.rm=T))
actop8 = which(ac > quantile(ac,probs=0.92,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.08*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop8 %in% bftop8)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop9 = which(bf > quantile(bf,probs=0.91,na.rm=T))
actop9 = which(ac > quantile(ac,probs=0.91,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.09*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop9 %in% bftop9)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop10 = which(bf > quantile(bf,probs=0.9,na.rm=T))
actop10 = which(ac > quantile(ac,probs=0.9,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.1*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop10 %in% bftop10)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
#glm and filtering west ac
fil_glaac = survival[!is.na(survival$glaAc) & !is.na(survival$gla_af),]
fit_glaac = glm(fil_glaac$glaAc ~ fil_glaac$gla_af, family = gaussian)
fit_west = glm(fil_glaac$west_bf ~ fil_glaac$west_af , family = gaussian)
#doing correlations and randomization ####
ac = fit_glaac$residuals
bf = fit_west$residuals
bndsb = quantile(bf,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndsa = quantile(ac,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
actop1 = which(ac > quantile(ac,probs=0.99,na.rm=T))
dim(fil_glaac)[1]
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.01*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop1 %in% bftop1)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop2 = which(bf > quantile(bf,probs=0.98,na.rm=T))
actop2 = which(ac > quantile(ac,probs=0.98,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.02*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop2 %in% bftop2)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop3 = which(bf > quantile(bf,probs=0.97,na.rm=T))
actop3 = which(ac > quantile(ac,probs=0.97,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.03*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop3 %in% bftop3)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop4 = which(bf > quantile(bf,probs=0.96,na.rm=T))
actop4 = which(ac > quantile(ac,probs=0.96,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.04*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop4 %in% bftop4)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop5 = which(bf > quantile(bf,probs=0.95,na.rm=T))
actop5 = which(ac > quantile(ac,probs=0.95,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.05*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop5 %in% bftop5)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop6 = which(bf > quantile(bf,probs=0.94,na.rm=T))
actop6 = which(ac > quantile(ac,probs=0.94,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.06*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop6 %in% bftop6)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop7 = which(bf > quantile(bf,probs=0.93,na.rm=T))
actop7 = which(ac > quantile(ac,probs=0.93,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.07*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop7 %in% bftop7)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop8 = which(bf > quantile(bf,probs=0.92,na.rm=T))
actop8 = which(ac > quantile(ac,probs=0.92,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.08*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop8 %in% bftop8)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop9 = which(bf > quantile(bf,probs=0.91,na.rm=T))
actop9 = which(ac > quantile(ac,probs=0.91,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.09*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop9 %in% bftop9)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop10 = which(bf > quantile(bf,probs=0.9,na.rm=T))
actop10 = which(ac > quantile(ac,probs=0.9,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.1*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop10 %in% bftop10)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
##west glaMs
#glm and filtering
fil_glams = survival[!is.na(survival$glaMs) & !is.na(survival$gla_af),]
fit_glams = glm(fil_glams$glaMs ~ fil_glams$gla_af, family = gaussian)
fit_west = glm(fil_glams$west_bf ~ fil_glams$west_af , family = gaussian)
#doing correlations and randomization ####
ms = fit_glams$residuals
bf = fit_west$residuals
bndsb = quantile(bf,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndsm = quantile(ms,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
mstop1 = which(ms > quantile(ms,probs=0.99,na.rm=T))
dim(fil_glams)[1]
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.01*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop1 %in% bftop1)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop2 = which(bf > quantile(bf,probs=0.98,na.rm=T))
mstop2 = which(ms > quantile(ms,probs=0.98,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.02*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop2 %in% bftop2)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop3 = which(bf > quantile(bf,probs=0.97,na.rm=T))
mstop3 = which(ms > quantile(ms,probs=0.97,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.03*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop3 %in% bftop3)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop4 = which(bf > quantile(bf,probs=0.96,na.rm=T))
mstop4 = which(ms > quantile(ms,probs=0.96,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.04*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop4 %in% bftop4)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop5 = which(bf > quantile(bf,probs=0.95,na.rm=T))
mstop5 = which(ms > quantile(ms,probs=0.95,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.05*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop5 %in% bftop5)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop6 = which(bf > quantile(bf,probs=0.94,na.rm=T))
mstop6 = which(ms > quantile(ms,probs=0.94,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.06*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop6 %in% bftop6)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop7 = which(bf > quantile(bf,probs=0.93,na.rm=T))
mstop7 = which(ms > quantile(ms,probs=0.93,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.07*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop7 %in% bftop7)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop8 = which(bf > quantile(bf,probs=0.92,na.rm=T))
mstop8 = which(ms > quantile(ms,probs=0.92,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.08*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop8 %in% bftop8)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop9 = which(bf > quantile(bf,probs=0.91,na.rm=T))
mstop9 = which(ms > quantile(ms,probs=0.91,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.09*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop9 %in% bftop9)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop10 = which(bf > quantile(bf,probs=0.9,na.rm=T))
mstop10 = which(ms > quantile(ms,probs=0.9,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.1*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop10 %in% bftop10)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
#glm and filtering east ac
fil_glaac = survival[!is.na(survival$glaAc) & !is.na(survival$gla_af),]
fit_glaac = glm(fil_glaac$glaAc ~ fil_glaac$gla_af, family = gaussian)
fit_east = glm(fil_glaac$east_bf ~ fil_glaac$east_af , family = gaussian)
#doing correlations and randomization ####
ac = fit_glaac$residuals
bf = fit_east$residuals
bndsb = quantile(bf,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndsa = quantile(ac,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
actop1 = which(ac > quantile(ac,probs=0.99,na.rm=T))
dim(fil_glaac)[1]
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.01*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop1 %in% bftop1)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop2 = which(bf > quantile(bf,probs=0.98,na.rm=T))
actop2 = which(ac > quantile(ac,probs=0.98,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.02*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop2 %in% bftop2)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop3 = which(bf > quantile(bf,probs=0.97,na.rm=T))
actop3 = which(ac > quantile(ac,probs=0.97,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.03*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop3 %in% bftop3)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop4 = which(bf > quantile(bf,probs=0.96,na.rm=T))
actop4 = which(ac > quantile(ac,probs=0.96,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.04*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop4 %in% bftop4)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop5 = which(bf > quantile(bf,probs=0.95,na.rm=T))
actop5 = which(ac > quantile(ac,probs=0.95,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.05*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop5 %in% bftop5)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop6 = which(bf > quantile(bf,probs=0.94,na.rm=T))
actop6 = which(ac > quantile(ac,probs=0.94,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.06*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop6 %in% bftop6)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop7 = which(bf > quantile(bf,probs=0.93,na.rm=T))
actop7 = which(ac > quantile(ac,probs=0.93,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.07*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop7 %in% bftop7)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop8 = which(bf > quantile(bf,probs=0.92,na.rm=T))
actop8 = which(ac > quantile(ac,probs=0.92,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.08*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop8 %in% bftop8)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop9 = which(bf > quantile(bf,probs=0.91,na.rm=T))
actop9 = which(ac > quantile(ac,probs=0.91,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.09*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop9 %in% bftop9)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop10 = which(bf > quantile(bf,probs=0.9,na.rm=T))
actop10 = which(ac > quantile(ac,probs=0.9,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.1*dim(fil_glaac)[1],0)
acnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glaac)[1],n,replace=F)
null[i]<-sum(acnull %in% bfnull)
}
obs = sum(actop10 %in% bftop10)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
##east glaMs
#glm and filtering
fil_glams = survival[!is.na(survival$glaMs) & !is.na(survival$gla_af),]
fit_glams = glm(fil_glams$glaMs ~ fil_glams$gla_af, family = gaussian)
fit_east = glm(fil_glams$east_bf ~ fil_glams$east_af , family = gaussian)
#doing correlations and randomization ####
ms = fit_glams$residuals
bf = fit_east$residuals
bndsb = quantile(bf,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bndsm = quantile(ms,probs=c(0.99,0.98, 0.97, 0.96, 0.95,0.94,0.93,0.92,0.91,0.9))
bftop1 = which(bf > quantile(bf,probs=0.99,na.rm=T))
mstop1 = which(ms > quantile(ms,probs=0.99,na.rm=T))
dim(fil_glams)[1]
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.01*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop1 %in% bftop1)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop2 = which(bf > quantile(bf,probs=0.98,na.rm=T))
mstop2 = which(ms > quantile(ms,probs=0.98,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.02*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop2 %in% bftop2)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop3 = which(bf > quantile(bf,probs=0.97,na.rm=T))
mstop3 = which(ms > quantile(ms,probs=0.97,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.03*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop3 %in% bftop3)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop4 = which(bf > quantile(bf,probs=0.96,na.rm=T))
mstop4 = which(ms > quantile(ms,probs=0.96,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.04*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop4 %in% bftop4)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop5 = which(bf > quantile(bf,probs=0.95,na.rm=T))
mstop5 = which(ms > quantile(ms,probs=0.95,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.05*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop5 %in% bftop5)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop6 = which(bf > quantile(bf,probs=0.94,na.rm=T))
mstop6 = which(ms > quantile(ms,probs=0.94,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.06*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop6 %in% bftop6)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop7 = which(bf > quantile(bf,probs=0.93,na.rm=T))
mstop7 = which(ms > quantile(ms,probs=0.93,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.07*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop7 %in% bftop7)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop8 = which(bf > quantile(bf,probs=0.92,na.rm=T))
mstop8 = which(ms > quantile(ms,probs=0.92,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.08*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop8 %in% bftop8)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop9 = which(bf > quantile(bf,probs=0.91,na.rm=T))
mstop9 = which(ms > quantile(ms,probs=0.91,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.09*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop9 %in% bftop9)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold
bftop10 = which(bf > quantile(bf,probs=0.9,na.rm=T))
mstop10 = which(ms > quantile(ms,probs=0.9,na.rm=T))
null<-rep(NA,10000)
for(i in 1:10000){
n<-round(0.1*dim(fil_glams)[1],0)
msnull<-sample(1:dim(fil_glams)[1],n,replace=F)
bfnull<-sample(1:dim(fil_glams)[1],n,replace=F)
null[i]<-sum(msnull %in% bfnull)
}
obs = sum(mstop10 %in% bftop10)
mean(null)
p = mean(null >= obs)
xfold = obs / mean(null)
obs
p
xfold