tanu/LSHTM_analysis
1
0
Fork 0
Code Issues Pull requests Releases Wiki Activity

extracting other params from logistic

Browse source
This commit is contained in:
Tanushree Tunstall 2020-06-22 14:11:16 +01:00
parent 28e52d4194
commit 6b5ced65e5
3 changed files with 94 additions and 31 deletions
Download patch file Download diff file Expand all files Collapse all files

1
scripts/af_or_calcs.R
View file

@ -388,7 +388,6 @@ pvals_logistic = sapply(gene_snps_unique,function(m){
mut = grepl(m,raw_data$all_muts_gene)
#print(table(dst, mut))
model<-glm(dst ~ mut , family = binomial)
#or_logistic = exp(summary(model)$coefficients[2,1])
pval_logistic = summary(model)$coefficients[2,4]
})

108
scripts/af_or_calcs_scratch.R
View file

@ -163,7 +163,7 @@ cat(paste0('Total no. of distinct comp snps to perform OR calcs: ', length(gene_
# Define OR function
#x = as.numeric(mut)
#y = dst
my_chisq_or = function(x,y){
logistic_chisq_or = function(x,y){
tab = as.matrix(table(x,y))
a = tab[2,2]
if (a==0){ a<-0.5}
@ -214,7 +214,7 @@ chisq.test(table(mut,dst))
fisher.test(table(mut, dst))
fisher.test(table(mut, dst))$p.value
fisher.test(table(mut, dst))$estimate
my_chisq_or(mut,dst)
logistic_chisq_or(mut,dst)
# logistic or
summary(model<-glm(dst ~ mut, family = binomial))
@ -230,7 +230,7 @@ pval_logistic2 = summary(model2)$coefficients[2,4]; print(pval_logistic2)
ors = sapply(gene_snps_unique,function(m){
mut = grepl(m,raw_data$all_muts_gene)
my_chisq_or(mut,dst)
logistic_chisq_or(mut,dst)
})
ors
@ -249,15 +249,69 @@ afs = sapply(gene_snps_unique,function(m){
afs
# logistic
logistic_ors = sapply(gene_snps_unique,function(m){
# logistic or
ors_logistic = sapply(gene_snps_unique,function(m){
mut = grepl(m,raw_data$all_muts_gene)
model<-glm(dst ~ mut, family = binomial)
or_logistic = exp(summary(model)$coefficients[2,1])
#pval_logistic = summary(model)$coefficients[2,4]
#logistic_se = summary(model)$coefficients[2,2]
#logistic_zval = summary(model)$coefficients[2,3]
#ci_mod = exp(confint(model))[2,]
#logistic_ci_lower = ci_mod[["2.5 %"]]
#logistic_ci_upper = ci_mod[["97.5 %"]]
})
logistic_ors
ors_logistic
head(ors_logistic); head(names(ors_logistic))
## logistic p-value
pvals_logistic = sapply(gene_snps_unique,function(m){
mut = grepl(m,raw_data$all_muts_gene)
model<-glm(dst ~ mut , family = binomial)
pval_logistic = summary(model)$coefficients[2,4]
})
head(pvals_logistic); head(names(pvals_logistic))
## logistic se
se_logistic = sapply(gene_snps_unique,function(m){
mut = grepl(m,raw_data$all_muts_gene)
model<-glm(dst ~ mut , family = binomial)
logistic_se = summary(model)$coefficients[2,2]
})
head(se_logistic); head(names(se_logistic))
## logistic z-value
zval_logistic = sapply(gene_snps_unique,function(m){
mut = grepl(m,raw_data$all_muts_gene)
model<-glm(dst ~ mut , family = binomial)
logistic_zval = summary(model)$coefficients[2,3]
})
head(zval_logistic); head(names(zval_logistic))
## logistic ci - lower bound
ci_lb_logistic = sapply(gene_snps_unique,function(m){
mut = grepl(m,raw_data$all_muts_gene)
model<-glm(dst ~ mut , family = binomial)
ci_mod = exp(confint(model))[2,]
logistic_ci_lower = ci_mod[["2.5 %"]]
})
head(ci_lb_logistic); head(names(ci_lb_logistic))
## logistic ci - upper bound
ci_ub_logistic = sapply(gene_snps_unique,function(m){
mut = grepl(m,raw_data$all_muts_gene)
model<-glm(dst ~ mut , family = binomial)
ci_mod = exp(confint(model))[2,]
logistic_ci_upper = ci_mod[["97.5 %"]]
})
head(ci_ub_logistic); head(names(ci_ub_logistic))
# logistic adj # Doesn't seem to make a difference
logistic_ors2 = sapply(gene_snps_unique,function(m){
@ -276,16 +330,35 @@ or_logistic2; pval_logistic2
head(logistic_ors)
#====================================
# logistic
summary(model<-glm(dst ~ mut
, family = binomial
#, control = glm.control(maxit = 1)
#, options(warn = 1)
))
or_logistic_maxit = exp(summary(model)$coefficients[2,1]); print(or_logistic_maxit)
or_logistic = exp(summary(model)$coefficients[2,1]); print(or_logistic)
pval_logistic_maxit = summary(model)$coefficients[2,4]; print(pval_logistic_maxit)
# extract SE of the logistic model for a given snp
logistic_se = summary(model)$coefficients[2,2]
print(paste0('SE:', logistic_se))
# extract Z of the logistic model for a given snp
logistic_zval = summary(model)$coefficients[2,3]
print(paste0('Z-value:', logistic_zval))
# extract confint interval of snp (2 steps, since the output is a named number)
ci_mod = exp(confint(model))[2,]
print(paste0('CI:', ci_mod))
#logistic_ci = paste(ci_mod[["2.5 %"]], ",", ci_mod[["97.5 %"]])
logistic_ci_lower = ci_mod[["2.5 %"]]
logistic_ci_upper = ci_mod[["97.5 %"]]
print(paste0('CI_lower:', logistic_ci_lower))
print(paste0('CI_upper:', logistic_ci_upper))
#####################
# iterate: subset
@ -296,10 +369,6 @@ snps_test = c("pnca_p.trp68gly", "pnca_p.leu4ser", "pnca_p.leu159arg","pnca_p.hi
data = snps_test[1:2]
data
################# start loop
for (i in data){
@ -328,9 +397,14 @@ for (i in data){
#, control = glm.control(maxit = 1)
#, options(warn = 1)
))
or_logistic_maxit = exp(summary(model)$coefficients[2,1]); print(or_logistic_maxit)
#, warning = perfectSeparation))
or_logistic = exp(summary(model)$coefficients[2,1]); print(or_logistic)
pval_logistic_maxit = summary(model)$coefficients[2,4]; print(pval_logistic_maxit)
logistic_se = summary(model)$coefficients[2,2]
logistic_zval = summary(model)$coefficients[2,3]
ci_mod = exp(confint(model))[2,]
logistic_ci_lower = ci_mod[["2.5 %"]]
logistic_ci_upper = ci_mod[["97.5 %"]]
#=====================
# fishers test
#=====================
@ -353,11 +427,9 @@ for (i in data){
# all output
writeLines(c(paste0("mutation:", i)
, paste0("=========================")
, paste0("OR_logistic_maxit:", or_logistic_maxit,"--->", "P-val_logistic_maxit:", pval_logistic_maxit )
, paste0("or_logistic:", or_logistic,"--->", "P-val_logistic_maxit:", pval_logistic_maxit )
, paste0("OR_fisher:", or_fisher, "--->","P-val_fisher:", pval_fisher )
, paste0("Chi_sq_estimate:", est_chisq, "--->","P-val_chisq:", pval_chisq)))
}
@ -382,4 +454,4 @@ table(dst, mut)
# https://stats.stackexchange.com/questions/259635/what-is-the-difference-using-a-fishers-exact-test-vs-a-logistic-regression-for
exact2x2(table(dst, mut),tsmethod="central")

16
scripts/combine_afs_ors.py
View file

@ -308,33 +308,25 @@ print('Finished writing file:', outfile
, '\nNo. of cols:', len(combined_or_df.columns)
, '\n=========================================================')
#%%
#%% practice
df = pd.DataFrame()
column_names = ['x','y','z','mean']
for col in column_names:
df[col] = np.random.randint(0,100, size=10000)
df.head()
# drop duplicate col with dup values not necessarily colnames
df['xdup'] = df['x']
df
df.T.drop_duplicates().T
df = df.T.drop_duplicates().T
import math
#import math
math.exp(0)
df['expX'] = np.exp(df['x']) # math doesn't understand series dtype
df
#%%