projects/04_aim4/r_old/v_18.05.30_chmi_aim_04.R
In mvazquezs/chmitools: Package to manage & obtain results of the CHMI Study
##############################################################
### CODE TO ANALYZE CHMI AB DATA FOLLOWING THE PAPER: AIM4 ###
##############################################################
# Start: 2017-06-16
rm(list=ls())
library(Hmisc)
library(plyr)
library(compareGroups)
manual.theme <- theme(panel.background = element_rect(fill = "white",
colour = NA), panel.border = element_rect(fill = NA,
colour = "grey20"), panel.grid.major = element_line(colour = "grey92"),
panel.grid.minor = element_line(colour = "grey92",
size = 0.25), strip.background = element_rect(fill = "grey85",
colour = "grey20"), legend.key = element_rect(fill = "white",
colour = NA), complete = TRUE, plot.title = element_text(hjust = 0.5))
fcn.pvalue.display<-function(p.vec, show.p = FALSE){
p.vec.plot<-as.character(rep(NA,length(p.vec)))
if(any((p.vec >= .01 & p.vec <= .04) | p.vec >= .06)==TRUE) {
p.vec.plot<-as.character(round(p.vec,2))
p.vec.plot[nchar(as.character(p.vec.plot))==3 & p.vec.plot > 0.1]<-paste(as.character(p.vec.plot)[nchar(as.character(p.vec.plot))==3 & p.vec.plot > 0.1],"0",sep="")
}
p.vec.plot<-ifelse ((p.vec >= 0.001 & p.vec < 0.01) | (p.vec > 0.04 & p.vec < 0.06), as.character(round(p.vec,3)),p.vec.plot)
if(show.p){
p.vec.plot[p.vec < 0.001]<-c("P<0.001")
p.vec.plot[p.vec >= 0.001]<-c(paste0("P=", p.vec.plot))
} else {
p.vec.plot[p.vec < 0.001]<-c("<0.001")
}
return(p.vec.plot)
}
## localpath <- "/Users/hsanz/Dropbox (M067)/m067_team/working_areas/MALIMMUNO Team Support/3. CHMI study Antibodies"
############################################################################################################################################
## load("/Users/hsanz/Dropbox (M067)/m067_team/working_areas/MALIMMUNO Team Support/3. CHMI study Antibodies/Data/1. All isotypes data 20171107/All_data_20171218.Rdata")
load('~/datasets_ori/v_17.12.18_chmi/v_17.12.18_chmi_dat.RData')
### SOME RECODES
#################
dat <- chmi_ab_l
names(dat) <- tolower(names(dat))
dim(dat)
head(dat)
xxx <- unique(dat[,c("original_id", "dataset", "immune_status", "status")])
with(xxx, table(status, dataset))
with(xxx, table(immune_status, dataset))
xxx <- unique(dat[,c("original_id", "dataset", "hb_status","status")])
with(xxx[xxx$dataset=="L1",], table(hb_status, status))
## Some recodes
dat$timepoint <- as.character(dat$timepoint)
dat$timepoint <- with(dat, ifelse(dataset=="T2" & timepoint=="DM", "D11", timepoint))
dat$hb <- as.character(dat$hb)
dat$hb <- with(dat, ifelse(dataset=="L1", "", hb))
dat$cvac_dose <- with(dat, ifelse(dataset=="T2", cvac_dose, NA))
dat$hb_status_benjamin <- NULL
dat$pre_patent_period <- as.character(dat$pre_patent_period)
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-005" & dataset=="T1", "9.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-018" & dataset=="T1", "7",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-039" & dataset=="T1", "11",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-040" & dataset=="T1", "8",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-041" & dataset=="T1", "9",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-043" & dataset=="T1", "7.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-044" & dataset=="T1", "10",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-045" & dataset=="T1", "11.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-048" & dataset=="T1", "12.5",pre_patent_period))
dat$gender <- ""
dat$gender <- with(dat, ifelse(original_id=="T1-005" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-018" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-039" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-040" & dataset=="T1", "F",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-041" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-043" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-044" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-045" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-048" & dataset=="T1", "F",gender))
datx <- dat
datx$log10_mfi <- NULL
datx$batch <- NULL
datx$hb <- NULL
datx$cvac_dose <- NULL
datx$pre_patent_period <- with(datx, ifelse(is.na(pre_patent_period), "NA", pre_patent_period))
datx$pcr.pos <- with(datx, ifelse(is.na(pcr.pos), 999, pcr.pos))
datx <- datx[, c("original_id", "timepoint","ig","study_number",
"pre_patent_period", "dataset", "antigen","mfi","pcr.pos","hb_status",
"immune_status", "status", "gender", "malaria")]
datx <- aggregate(mfi ~ original_id + timepoint + ig +
+ study_number + malaria
+ pre_patent_period + dataset + antigen + pcr.pos + immune_status
+ status + gender + hb_status, FUN = mean, data=datx)
xxx <- as.data.frame(with(datx[datx$antigen=="CSP",], table(original_id,ig,study_number, timepoint)))
xxx[xxx$Freq>1,]
dat <- datx
dat$log10_mfi <- log10(dat$mfi)
#### PREPARE DATASET FOR ANALYSIS
##################################
tmp.dat <- subset(dat, status!="Naive" | dataset!="T2")
tmp.dat <- subset(tmp.dat, antigen%nin%c("BSA", "agal"))
tmp.dat$status <- with(tmp.dat, ifelse(status=="Naive", "Naive",
ifelse(status=="PfSPZ_CVac_vaccinated", "Vaccinated",
ifelse(status=="Semi_immune", "Semi-immune","XXX"))))
tmp.dat$status <- factor(tmp.dat$status, levels=c("Naive", "Vaccinated", "Semi-immune"))
lantigen <- sort(unique(tmp.dat$antigen))
with(dat, table(status, dataset))
with(tmp.dat, table(status, dataset))
tmp.dat$timepoint2c <- with(tmp.dat, ifelse(timepoint=="D13" & dataset=="L1", "D11-D13",
ifelse(timepoint=="DM" & dataset=="T1", "D11-D13",
ifelse(timepoint=="D11" & dataset=="T2", "D11-D13",timepoint))))
## Read Extra variables.
extra.vars <- read.xls(file.path(localpath,"/Data/20180530/chmi_pid_info_PG.xlsx"))
dim(tmp.dat)
tmp.dat$pcr.pos <- NULL
tmp.dat$gender <- NULL
tmp.dat <- merge(tmp.dat, extra.vars[,c("original_id", "dataset", "pcr.pos", "gender", "height", "weight")],
by=c("original_id", "dataset"), all.x = TRUE, sort=FALSE)
dim(tmp.dat)
## Add mal.pos variable only for LaCHMI-001
extra.vars <- read.xls(file.path(localpath,"/Data/20180530/qPCRdata_CHMI_Patricia.xlsx"), sheet = "LACHMI-001")
extra.vars$original_id <- extra.vars$vid
extra.vars <- extra.vars[,c("original_id", "mal.pos")]
dim(tmp.dat)
tmp.dat <- merge(tmp.dat, extra.vars, by=c("original_id"), all.x = TRUE, sort=FALSE)
dim(tmp.dat)
# ok
with(extra.vars, table(mal.pos, useNA = "always"))
with(unique(tmp.dat[,c("original_id", "dataset","mal.pos")]), table(mal.pos, useNA = "always"))
### CHECK 1
table(tmp.dat[tmp.dat$original_id=="L1-023", "hb_status"])
table(tmp.dat$hb_status)
x <- read.xls(file.path(localpath,"/Data/20180530/qPCRdata_CHMI_Patricia.xlsx"), sheet = "LACHMI-001")
x$original_id <- x$vid
x <- x[,c("original_id", "ppp.tbs")]
dim(tmp.dat)
y <- merge(tmp.dat, x, by=c("original_id"), all.x = TRUE, sort=FALSE)
dim(y)
x <- y[!is.na(y$ppp.tbs),c("original_id", "dataset","ppp.tbs","pre_patent_period")]
with(x, table(pre_patent_period, ppp.tbs))
### CHECK 2
x <- read.xls(file.path(localpath,"/Data/20180530/qPCRdata_CHMI_Patricia.xlsx"), sheet = "LACHMI-001")
x$original_id <- x$vid
x <- x[,c("original_id", "tbs.pos")]
dim(tmp.dat)
y <- merge(tmp.dat, x, by=c("original_id"), all.x = TRUE, sort=FALSE)
dim(y)
x <- y[!is.na(y$tbs.pos),c("original_id", "dataset","tbs.pos","malaria")]
with(x, table(malaria, tbs.pos))
## Export dataset for Patricia
# export.example <- unique(tmp.dat[,c("original_id","study_number","dataset", "malaria","pcr.pos",
# "immune_status","status")])
# wh <- with(export.example, order(original_id))
# export.example <- export.example[wh,]
# export.example
# library(xlsx)
# write.xlsx(export.example, file="/Users/hsanz/Downloads/chmi_data_pid_20180510.xlsx", row.names = FALSE)
########################################
########################################
########################################
### AIM 4 Analysis: MALARIA OUTCOME (malaria variable)
########################################
########################################
########################################
### Confidence Interval
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x){
model <- glm(y ~ x,family=binomial(link='logit'))
ans <- confint(model)
ans <- c(coef(model)[2], ans[2,])
ans <- round(exp(ans),2)
ans <- paste0(ans[1], " (", ans[2],";", ans[3],")")
#ans <- data.frame(ans=ans)
return(ans)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise, OR=ff(malaria, log10_mfi))
ans.summ.wide <- reshape2::dcast(ans.summ, ig + antigen ~ timepoint2c, value.var="OR")
ans.summ.wide <- ans.summ.wide[, c(1,2,3,5,4)]
### P value
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x, rr = "pval"){
model <- glm(y ~ x,family=binomial(link='logit'))
or <- round(exp(coef(model)[2]),2)
raw.pval <- summary(model)$coefficients[2,4]
if(rr=="pval") return(raw.pval)
if(rr=="or") return(or)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise,
OR=ff(malaria, log10_mfi, "or"), RAW.P=ff(malaria, log10_mfi, "pval"))
xiso <- unique(ans.summ$ig)
adj.pvals <- sapply(xiso, function(x) p.adjust(ans.summ[ans.summ$ig==x,"RAW.P"], "BH"))
adj.pvals <- unlist(adj.pvals, use.names=FALSE)
ans.summ$ADJ.PVAL <- fcn.pvalue.display(adj.pvals)
ans.summ$RAW.P <- fcn.pvalue.display(ans.summ$RAW.P)
ans.summ.wide <- data.table::dcast(data.table::setDT(ans.summ),
ig + antigen ~ timepoint2c,
value.var=c("OR","RAW.P", "ADJ.PVAL"))
ans.summ.wide <- ans.summ.wide[, c(1,2,3,6,9,5,8,11,4,7,10)]
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_malaria_with_pvals.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Logistic regression between malaria outcome (defined by clinical malaria - 'malaria' variable in the dataset) and antigen (log10[MFI] scale) by isotype and timepoint independently. Odd ratio, raw p value and isotype specific adjusted p value are shown for ", xvector[i] , " isotype."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.summ.wide[ans.summ.wide$ig==xvector[i],] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 8))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Timepoint"),
colspan = c(1,1,9), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","C-1","D7","D11-D13"),
colspan = c(1,1,3,3,3), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen",rep(c("OR","RAW.P","ADJ.P"),3)),
colspan = rep(1,11), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_malaria_with_pvals.docx"))
########################################
########################################
########################################
### AIM 4 Analysis: MALARIA OUTCOME (pcr.pos variable)
########################################
########################################
########################################
### Confidence Interval
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x){
model <- glm(y ~ x,family=binomial(link='logit'))
ans <- confint(model)
ans <- c(coef(model)[2], ans[2,])
ans <- round(exp(ans),2)
ans <- paste0(ans[1], " (", ans[2],";", ans[3],")")
#ans <- data.frame(ans=ans)
return(ans)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise, OR=ff(pcr.pos, log10_mfi))
ans.summ.wide <- reshape2::dcast(ans.summ, ig + antigen ~ timepoint2c, value.var="OR")
ans.summ.wide <- ans.summ.wide[, c(1,2,3,5,4)]
### P value
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x, rr = "pval"){
model <- glm(y ~ x,family=binomial(link='logit'))
or <- round(exp(coef(model)[2]),2)
raw.pval <- summary(model)$coefficients[2,4]
if(rr=="pval") return(raw.pval)
if(rr=="or") return(or)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise,
OR=ff(pcr.pos, log10_mfi, "or"), RAW.P=ff(pcr.pos, log10_mfi, "pval"))
xiso <- unique(ans.summ$ig)
adj.pvals <- sapply(xiso, function(x) p.adjust(ans.summ[ans.summ$ig==x,"RAW.P"], "BH"))
adj.pvals <- unlist(adj.pvals, use.names=FALSE)
ans.summ$ADJ.PVAL <- fcn.pvalue.display(adj.pvals)
ans.summ$RAW.P <- fcn.pvalue.display(ans.summ$RAW.P)
ans.summ.wide <- data.table::dcast(data.table::setDT(ans.summ),
ig + antigen ~ timepoint2c,
value.var=c("OR","RAW.P", "ADJ.PVAL"))
ans.summ.wide <- ans.summ.wide[, c(1,2,3,6,9,5,8,11,4,7,10)]
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_pcr_with_pvals.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Logistic regression between malaria outcome (defined by PCR - 'pcr.pos' variable in the dataset) and antigen (log10[MFI] scale) by isotype and timepoint independently. Odd ratio, raw p value and isotype specific adjusted p value are shown for ", xvector[i] , " isotype."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.summ.wide[ans.summ.wide$ig==xvector[i],] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 8))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Timepoint"),
colspan = c(1,1,9), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","C-1","D7","D11-D13"),
colspan = c(1,1,3,3,3), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen",rep(c("OR","RAW.P","ADJ.P"),3)),
colspan = rep(1,11), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_pcr_with_pvals.docx"))
########################################
########################################
########################################
### AIM 4 Analysis: MALARIA OUTCOME (pcr.pos variable)
########################################
########################################
########################################
### Confidence Interval
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt & dataset=="L1")
ff <- function(y, x){
model <- glm(y ~ x,family=binomial(link='logit'))
ans <- confint(model)
ans <- c(coef(model)[2], ans[2,])
ans <- round(exp(ans),2)
ans <- paste0(ans[1], " (", ans[2],";", ans[3],")")
#ans <- data.frame(ans=ans)
return(ans)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise, OR=ff(mal.pos, log10_mfi))
ans.summ.wide <- reshape2::dcast(ans.summ, ig + antigen ~ timepoint2c, value.var="OR")
ans.summ.wide <- ans.summ.wide[, c(1,2,3,5,4)]
### P value
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x, rr = "pval"){
model <- glm(y ~ x,family=binomial(link='logit'))
or <- round(exp(coef(model)[2]),2)
raw.pval <- summary(model)$coefficients[2,4]
if(rr=="pval") return(raw.pval)
if(rr=="or") return(or)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise,
OR=ff(mal.pos, log10_mfi, "or"), RAW.P=ff(mal.pos, log10_mfi, "pval"))
xiso <- unique(ans.summ$ig)
adj.pvals <- sapply(xiso, function(x) p.adjust(ans.summ[ans.summ$ig==x,"RAW.P"], "BH"))
adj.pvals <- unlist(adj.pvals, use.names=FALSE)
ans.summ$ADJ.PVAL <- fcn.pvalue.display(adj.pvals)
ans.summ$RAW.P <- fcn.pvalue.display(ans.summ$RAW.P)
ans.summ.wide <- data.table::dcast(data.table::setDT(ans.summ),
ig + antigen ~ timepoint2c,
value.var=c("OR","RAW.P", "ADJ.PVAL"))
ans.summ.wide <- ans.summ.wide[, c(1,2,3,6,9,5,8,11,4,7,10)]
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_lachmi_with_pvals.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Logistic regression between malaria outcome (defined by symptomatic malaria only on LACHMI obs. - 'mal.pos' variable in the dataset) and antigen (log10[MFI] scale) by isotype and timepoint independently. Odd ratio, raw p value and isotype specific adjusted p value are shown for ", xvector[i] , " isotype."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.summ.wide[ans.summ.wide$ig==xvector[i],] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 8))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Timepoint"),
colspan = c(1,1,9), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","C-1","D7","D11-D13"),
colspan = c(1,1,3,3,3), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen",rep(c("OR","RAW.P","ADJ.P"),3)),
colspan = rep(1,11), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_lachmi_with_pvals.docx"))
################################################
################################################
### AIM 4 Analysis: Pre-Patent period outcome
################################################
################################################
## CInterval
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
aux.data <- subset(aux.data, malaria==1)
aux.data$pre_patent_period <- with(aux.data, ifelse(pre_patent_period%in%c("missing", "NA"), NA, pre_patent_period))
aux.data$pre_patent_period <- as.numeric(aux.data$pre_patent_period)
summary(aux.data$pre_patent_period)
sum(!is.na(aux.data$pre_patent_period))
ff <- function(y, x){
model <- lm(y ~ x)
ans <- confint(model)
ans <- c(coef(model)[2], ans[2,])
ans <- round(ans,2)
ans <- paste0(ans[1], " (", ans[2],";", ans[3],")")
#ans <- data.frame(ans=ans)
return(ans)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise, Beta=ff(pre_patent_period, log10_mfi))
ans.summ.wide <- reshape2::dcast(ans.summ, ig + antigen ~ timepoint2c, value.var="Beta")
ans.summ.wide <- ans.summ.wide[, c(1,2,3,5,4)]
### P value
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
aux.data <- subset(aux.data, malaria==1)
aux.data$pre_patent_period <- with(aux.data, ifelse(pre_patent_period%in%c("missing", "NA"), NA, pre_patent_period))
aux.data$pre_patent_period <- as.numeric(aux.data$pre_patent_period)
ff <- function(y, x, rr = "pval"){
model <- lm(y ~ x)
beta <- round((coef(model)[2]),2)
raw.pval <- summary(model)$coefficients[2,4]
if(rr=="pval") return(raw.pval)
if(rr=="beta") return(beta)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise,
BETA=ff(pre_patent_period, log10_mfi, "beta"), RAW.P=ff(pre_patent_period, log10_mfi, "pval"))
xiso <- unique(ans.summ$ig)
adj.pvals <- sapply(xiso, function(x) p.adjust(ans.summ[ans.summ$ig==x,"RAW.P"], "BH"))
adj.pvals <- unlist(adj.pvals, use.names=FALSE)
ans.summ$ADJ.PVAL <- fcn.pvalue.display(adj.pvals)
ans.summ$RAW.P <- fcn.pvalue.display(ans.summ$RAW.P)
ans.summ.wide <- data.table::dcast(data.table::setDT(ans.summ),
ig + antigen ~ timepoint2c,
value.var=c("BETA","RAW.P", "ADJ.PVAL"))
ans.summ.wide <- ans.summ.wide[, c(1,2,3,6,9,5,8,11,4,7,10)]
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_prepatentperiod_with_pvals.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Linear regression between pre-patent period and antigen (log10[MFI] scale) by isotype and timepoint independently. Only malaria cases are selected. Beta coefficient, raw p value and isotype specific adjusted p value are shown for ", xvector[i] , " isotype."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.summ.wide[ans.summ.wide$ig==xvector[i],] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 8))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Timepoint"),
colspan = c(1,1,9), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","C-1","D7","D11-D13"),
colspan = c(1,1,3,3,3), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen",rep(c("Beta","RAW.P","ADJ.P"),3)),
colspan = rep(1,11), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_prepatentperiod_with_pvals.docx"))
mvazquezs/chmitools documentation built on May 1, 2020, 2:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
##############################################################
### CODE TO ANALYZE CHMI AB DATA FOLLOWING THE PAPER: AIM4 ###
##############################################################
# Start: 2017-06-16
rm(list=ls())
library(Hmisc)
library(plyr)
library(compareGroups)
manual.theme <- theme(panel.background = element_rect(fill = "white",
colour = NA), panel.border = element_rect(fill = NA,
colour = "grey20"), panel.grid.major = element_line(colour = "grey92"),
panel.grid.minor = element_line(colour = "grey92",
size = 0.25), strip.background = element_rect(fill = "grey85",
colour = "grey20"), legend.key = element_rect(fill = "white",
colour = NA), complete = TRUE, plot.title = element_text(hjust = 0.5))
fcn.pvalue.display<-function(p.vec, show.p = FALSE){
p.vec.plot<-as.character(rep(NA,length(p.vec)))
if(any((p.vec >= .01 & p.vec <= .04) | p.vec >= .06)==TRUE) {
p.vec.plot<-as.character(round(p.vec,2))
p.vec.plot[nchar(as.character(p.vec.plot))==3 & p.vec.plot > 0.1]<-paste(as.character(p.vec.plot)[nchar(as.character(p.vec.plot))==3 & p.vec.plot > 0.1],"0",sep="")
}
p.vec.plot<-ifelse ((p.vec >= 0.001 & p.vec < 0.01) | (p.vec > 0.04 & p.vec < 0.06), as.character(round(p.vec,3)),p.vec.plot)
if(show.p){
p.vec.plot[p.vec < 0.001]<-c("P<0.001")
p.vec.plot[p.vec >= 0.001]<-c(paste0("P=", p.vec.plot))
} else {
p.vec.plot[p.vec < 0.001]<-c("<0.001")
}
return(p.vec.plot)
}
## localpath <- "/Users/hsanz/Dropbox (M067)/m067_team/working_areas/MALIMMUNO Team Support/3. CHMI study Antibodies"
############################################################################################################################################
## load("/Users/hsanz/Dropbox (M067)/m067_team/working_areas/MALIMMUNO Team Support/3. CHMI study Antibodies/Data/1. All isotypes data 20171107/All_data_20171218.Rdata")
load('~/datasets_ori/v_17.12.18_chmi/v_17.12.18_chmi_dat.RData')
### SOME RECODES
#################
dat <- chmi_ab_l
names(dat) <- tolower(names(dat))
dim(dat)
head(dat)
xxx <- unique(dat[,c("original_id", "dataset", "immune_status", "status")])
with(xxx, table(status, dataset))
with(xxx, table(immune_status, dataset))
xxx <- unique(dat[,c("original_id", "dataset", "hb_status","status")])
with(xxx[xxx$dataset=="L1",], table(hb_status, status))
## Some recodes
dat$timepoint <- as.character(dat$timepoint)
dat$timepoint <- with(dat, ifelse(dataset=="T2" & timepoint=="DM", "D11", timepoint))
dat$hb <- as.character(dat$hb)
dat$hb <- with(dat, ifelse(dataset=="L1", "", hb))
dat$cvac_dose <- with(dat, ifelse(dataset=="T2", cvac_dose, NA))
dat$hb_status_benjamin <- NULL
dat$pre_patent_period <- as.character(dat$pre_patent_period)
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-005" & dataset=="T1", "9.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-018" & dataset=="T1", "7",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-039" & dataset=="T1", "11",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-040" & dataset=="T1", "8",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-041" & dataset=="T1", "9",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-043" & dataset=="T1", "7.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-044" & dataset=="T1", "10",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-045" & dataset=="T1", "11.5",pre_patent_period))
dat$pre_patent_period <- with(dat, ifelse(original_id=="T1-048" & dataset=="T1", "12.5",pre_patent_period))
dat$gender <- ""
dat$gender <- with(dat, ifelse(original_id=="T1-005" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-018" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-039" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-040" & dataset=="T1", "F",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-041" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-043" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-044" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-045" & dataset=="T1", "M",gender))
dat$gender <- with(dat, ifelse(original_id=="T1-048" & dataset=="T1", "F",gender))
datx <- dat
datx$log10_mfi <- NULL
datx$batch <- NULL
datx$hb <- NULL
datx$cvac_dose <- NULL
datx$pre_patent_period <- with(datx, ifelse(is.na(pre_patent_period), "NA", pre_patent_period))
datx$pcr.pos <- with(datx, ifelse(is.na(pcr.pos), 999, pcr.pos))
datx <- datx[, c("original_id", "timepoint","ig","study_number",
"pre_patent_period", "dataset", "antigen","mfi","pcr.pos","hb_status",
"immune_status", "status", "gender", "malaria")]
datx <- aggregate(mfi ~ original_id + timepoint + ig +
+ study_number + malaria
+ pre_patent_period + dataset + antigen + pcr.pos + immune_status
+ status + gender + hb_status, FUN = mean, data=datx)
xxx <- as.data.frame(with(datx[datx$antigen=="CSP",], table(original_id,ig,study_number, timepoint)))
xxx[xxx$Freq>1,]
dat <- datx
dat$log10_mfi <- log10(dat$mfi)
#### PREPARE DATASET FOR ANALYSIS
##################################
tmp.dat <- subset(dat, status!="Naive" | dataset!="T2")
tmp.dat <- subset(tmp.dat, antigen%nin%c("BSA", "agal"))
tmp.dat$status <- with(tmp.dat, ifelse(status=="Naive", "Naive",
ifelse(status=="PfSPZ_CVac_vaccinated", "Vaccinated",
ifelse(status=="Semi_immune", "Semi-immune","XXX"))))
tmp.dat$status <- factor(tmp.dat$status, levels=c("Naive", "Vaccinated", "Semi-immune"))
lantigen <- sort(unique(tmp.dat$antigen))
with(dat, table(status, dataset))
with(tmp.dat, table(status, dataset))
tmp.dat$timepoint2c <- with(tmp.dat, ifelse(timepoint=="D13" & dataset=="L1", "D11-D13",
ifelse(timepoint=="DM" & dataset=="T1", "D11-D13",
ifelse(timepoint=="D11" & dataset=="T2", "D11-D13",timepoint))))
## Read Extra variables.
extra.vars <- read.xls(file.path(localpath,"/Data/20180530/chmi_pid_info_PG.xlsx"))
dim(tmp.dat)
tmp.dat$pcr.pos <- NULL
tmp.dat$gender <- NULL
tmp.dat <- merge(tmp.dat, extra.vars[,c("original_id", "dataset", "pcr.pos", "gender", "height", "weight")],
by=c("original_id", "dataset"), all.x = TRUE, sort=FALSE)
dim(tmp.dat)
## Add mal.pos variable only for LaCHMI-001
extra.vars <- read.xls(file.path(localpath,"/Data/20180530/qPCRdata_CHMI_Patricia.xlsx"), sheet = "LACHMI-001")
extra.vars$original_id <- extra.vars$vid
extra.vars <- extra.vars[,c("original_id", "mal.pos")]
dim(tmp.dat)
tmp.dat <- merge(tmp.dat, extra.vars, by=c("original_id"), all.x = TRUE, sort=FALSE)
dim(tmp.dat)
# ok
with(extra.vars, table(mal.pos, useNA = "always"))
with(unique(tmp.dat[,c("original_id", "dataset","mal.pos")]), table(mal.pos, useNA = "always"))
### CHECK 1
table(tmp.dat[tmp.dat$original_id=="L1-023", "hb_status"])
table(tmp.dat$hb_status)
x <- read.xls(file.path(localpath,"/Data/20180530/qPCRdata_CHMI_Patricia.xlsx"), sheet = "LACHMI-001")
x$original_id <- x$vid
x <- x[,c("original_id", "ppp.tbs")]
dim(tmp.dat)
y <- merge(tmp.dat, x, by=c("original_id"), all.x = TRUE, sort=FALSE)
dim(y)
x <- y[!is.na(y$ppp.tbs),c("original_id", "dataset","ppp.tbs","pre_patent_period")]
with(x, table(pre_patent_period, ppp.tbs))
### CHECK 2
x <- read.xls(file.path(localpath,"/Data/20180530/qPCRdata_CHMI_Patricia.xlsx"), sheet = "LACHMI-001")
x$original_id <- x$vid
x <- x[,c("original_id", "tbs.pos")]
dim(tmp.dat)
y <- merge(tmp.dat, x, by=c("original_id"), all.x = TRUE, sort=FALSE)
dim(y)
x <- y[!is.na(y$tbs.pos),c("original_id", "dataset","tbs.pos","malaria")]
with(x, table(malaria, tbs.pos))
## Export dataset for Patricia
# export.example <- unique(tmp.dat[,c("original_id","study_number","dataset", "malaria","pcr.pos",
# "immune_status","status")])
# wh <- with(export.example, order(original_id))
# export.example <- export.example[wh,]
# export.example
# library(xlsx)
# write.xlsx(export.example, file="/Users/hsanz/Downloads/chmi_data_pid_20180510.xlsx", row.names = FALSE)
########################################
########################################
########################################
### AIM 4 Analysis: MALARIA OUTCOME (malaria variable)
########################################
########################################
########################################
### Confidence Interval
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x){
model <- glm(y ~ x,family=binomial(link='logit'))
ans <- confint(model)
ans <- c(coef(model)[2], ans[2,])
ans <- round(exp(ans),2)
ans <- paste0(ans[1], " (", ans[2],";", ans[3],")")
#ans <- data.frame(ans=ans)
return(ans)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise, OR=ff(malaria, log10_mfi))
ans.summ.wide <- reshape2::dcast(ans.summ, ig + antigen ~ timepoint2c, value.var="OR")
ans.summ.wide <- ans.summ.wide[, c(1,2,3,5,4)]
### P value
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x, rr = "pval"){
model <- glm(y ~ x,family=binomial(link='logit'))
or <- round(exp(coef(model)[2]),2)
raw.pval <- summary(model)$coefficients[2,4]
if(rr=="pval") return(raw.pval)
if(rr=="or") return(or)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise,
OR=ff(malaria, log10_mfi, "or"), RAW.P=ff(malaria, log10_mfi, "pval"))
xiso <- unique(ans.summ$ig)
adj.pvals <- sapply(xiso, function(x) p.adjust(ans.summ[ans.summ$ig==x,"RAW.P"], "BH"))
adj.pvals <- unlist(adj.pvals, use.names=FALSE)
ans.summ$ADJ.PVAL <- fcn.pvalue.display(adj.pvals)
ans.summ$RAW.P <- fcn.pvalue.display(ans.summ$RAW.P)
ans.summ.wide <- data.table::dcast(data.table::setDT(ans.summ),
ig + antigen ~ timepoint2c,
value.var=c("OR","RAW.P", "ADJ.PVAL"))
ans.summ.wide <- ans.summ.wide[, c(1,2,3,6,9,5,8,11,4,7,10)]
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_malaria_with_pvals.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Logistic regression between malaria outcome (defined by clinical malaria - 'malaria' variable in the dataset) and antigen (log10[MFI] scale) by isotype and timepoint independently. Odd ratio, raw p value and isotype specific adjusted p value are shown for ", xvector[i] , " isotype."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.summ.wide[ans.summ.wide$ig==xvector[i],] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 8))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Timepoint"),
colspan = c(1,1,9), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","C-1","D7","D11-D13"),
colspan = c(1,1,3,3,3), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen",rep(c("OR","RAW.P","ADJ.P"),3)),
colspan = rep(1,11), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_malaria_with_pvals.docx"))
########################################
########################################
########################################
### AIM 4 Analysis: MALARIA OUTCOME (pcr.pos variable)
########################################
########################################
########################################
### Confidence Interval
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x){
model <- glm(y ~ x,family=binomial(link='logit'))
ans <- confint(model)
ans <- c(coef(model)[2], ans[2,])
ans <- round(exp(ans),2)
ans <- paste0(ans[1], " (", ans[2],";", ans[3],")")
#ans <- data.frame(ans=ans)
return(ans)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise, OR=ff(pcr.pos, log10_mfi))
ans.summ.wide <- reshape2::dcast(ans.summ, ig + antigen ~ timepoint2c, value.var="OR")
ans.summ.wide <- ans.summ.wide[, c(1,2,3,5,4)]
### P value
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x, rr = "pval"){
model <- glm(y ~ x,family=binomial(link='logit'))
or <- round(exp(coef(model)[2]),2)
raw.pval <- summary(model)$coefficients[2,4]
if(rr=="pval") return(raw.pval)
if(rr=="or") return(or)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise,
OR=ff(pcr.pos, log10_mfi, "or"), RAW.P=ff(pcr.pos, log10_mfi, "pval"))
xiso <- unique(ans.summ$ig)
adj.pvals <- sapply(xiso, function(x) p.adjust(ans.summ[ans.summ$ig==x,"RAW.P"], "BH"))
adj.pvals <- unlist(adj.pvals, use.names=FALSE)
ans.summ$ADJ.PVAL <- fcn.pvalue.display(adj.pvals)
ans.summ$RAW.P <- fcn.pvalue.display(ans.summ$RAW.P)
ans.summ.wide <- data.table::dcast(data.table::setDT(ans.summ),
ig + antigen ~ timepoint2c,
value.var=c("OR","RAW.P", "ADJ.PVAL"))
ans.summ.wide <- ans.summ.wide[, c(1,2,3,6,9,5,8,11,4,7,10)]
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_pcr_with_pvals.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Logistic regression between malaria outcome (defined by PCR - 'pcr.pos' variable in the dataset) and antigen (log10[MFI] scale) by isotype and timepoint independently. Odd ratio, raw p value and isotype specific adjusted p value are shown for ", xvector[i] , " isotype."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.summ.wide[ans.summ.wide$ig==xvector[i],] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 8))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Timepoint"),
colspan = c(1,1,9), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","C-1","D7","D11-D13"),
colspan = c(1,1,3,3,3), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen",rep(c("OR","RAW.P","ADJ.P"),3)),
colspan = rep(1,11), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_pcr_with_pvals.docx"))
########################################
########################################
########################################
### AIM 4 Analysis: MALARIA OUTCOME (pcr.pos variable)
########################################
########################################
########################################
### Confidence Interval
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt & dataset=="L1")
ff <- function(y, x){
model <- glm(y ~ x,family=binomial(link='logit'))
ans <- confint(model)
ans <- c(coef(model)[2], ans[2,])
ans <- round(exp(ans),2)
ans <- paste0(ans[1], " (", ans[2],";", ans[3],")")
#ans <- data.frame(ans=ans)
return(ans)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise, OR=ff(mal.pos, log10_mfi))
ans.summ.wide <- reshape2::dcast(ans.summ, ig + antigen ~ timepoint2c, value.var="OR")
ans.summ.wide <- ans.summ.wide[, c(1,2,3,5,4)]
### P value
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
ff <- function(y, x, rr = "pval"){
model <- glm(y ~ x,family=binomial(link='logit'))
or <- round(exp(coef(model)[2]),2)
raw.pval <- summary(model)$coefficients[2,4]
if(rr=="pval") return(raw.pval)
if(rr=="or") return(or)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise,
OR=ff(mal.pos, log10_mfi, "or"), RAW.P=ff(mal.pos, log10_mfi, "pval"))
xiso <- unique(ans.summ$ig)
adj.pvals <- sapply(xiso, function(x) p.adjust(ans.summ[ans.summ$ig==x,"RAW.P"], "BH"))
adj.pvals <- unlist(adj.pvals, use.names=FALSE)
ans.summ$ADJ.PVAL <- fcn.pvalue.display(adj.pvals)
ans.summ$RAW.P <- fcn.pvalue.display(ans.summ$RAW.P)
ans.summ.wide <- data.table::dcast(data.table::setDT(ans.summ),
ig + antigen ~ timepoint2c,
value.var=c("OR","RAW.P", "ADJ.PVAL"))
ans.summ.wide <- ans.summ.wide[, c(1,2,3,6,9,5,8,11,4,7,10)]
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_lachmi_with_pvals.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Logistic regression between malaria outcome (defined by symptomatic malaria only on LACHMI obs. - 'mal.pos' variable in the dataset) and antigen (log10[MFI] scale) by isotype and timepoint independently. Odd ratio, raw p value and isotype specific adjusted p value are shown for ", xvector[i] , " isotype."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.summ.wide[ans.summ.wide$ig==xvector[i],] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 8))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Timepoint"),
colspan = c(1,1,9), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","C-1","D7","D11-D13"),
colspan = c(1,1,3,3,3), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen",rep(c("OR","RAW.P","ADJ.P"),3)),
colspan = rep(1,11), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_logistic_lachmi_with_pvals.docx"))
################################################
################################################
### AIM 4 Analysis: Pre-Patent period outcome
################################################
################################################
## CInterval
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
aux.data <- subset(aux.data, malaria==1)
aux.data$pre_patent_period <- with(aux.data, ifelse(pre_patent_period%in%c("missing", "NA"), NA, pre_patent_period))
aux.data$pre_patent_period <- as.numeric(aux.data$pre_patent_period)
summary(aux.data$pre_patent_period)
sum(!is.na(aux.data$pre_patent_period))
ff <- function(y, x){
model <- lm(y ~ x)
ans <- confint(model)
ans <- c(coef(model)[2], ans[2,])
ans <- round(ans,2)
ans <- paste0(ans[1], " (", ans[2],";", ans[3],")")
#ans <- data.frame(ans=ans)
return(ans)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise, Beta=ff(pre_patent_period, log10_mfi))
ans.summ.wide <- reshape2::dcast(ans.summ, ig + antigen ~ timepoint2c, value.var="Beta")
ans.summ.wide <- ans.summ.wide[, c(1,2,3,5,4)]
### P value
tt <- c("C-1", "D7", "D11-D13")
aux.data <- subset(tmp.dat, timepoint2c%in%tt)
aux.data <- subset(aux.data, malaria==1)
aux.data$pre_patent_period <- with(aux.data, ifelse(pre_patent_period%in%c("missing", "NA"), NA, pre_patent_period))
aux.data$pre_patent_period <- as.numeric(aux.data$pre_patent_period)
ff <- function(y, x, rr = "pval"){
model <- lm(y ~ x)
beta <- round((coef(model)[2]),2)
raw.pval <- summary(model)$coefficients[2,4]
if(rr=="pval") return(raw.pval)
if(rr=="beta") return(beta)
}
ans.summ <- ddply(aux.data, .(ig, antigen, timepoint2c), summarise,
BETA=ff(pre_patent_period, log10_mfi, "beta"), RAW.P=ff(pre_patent_period, log10_mfi, "pval"))
xiso <- unique(ans.summ$ig)
adj.pvals <- sapply(xiso, function(x) p.adjust(ans.summ[ans.summ$ig==x,"RAW.P"], "BH"))
adj.pvals <- unlist(adj.pvals, use.names=FALSE)
ans.summ$ADJ.PVAL <- fcn.pvalue.display(adj.pvals)
ans.summ$RAW.P <- fcn.pvalue.display(ans.summ$RAW.P)
ans.summ.wide <- data.table::dcast(data.table::setDT(ans.summ),
ig + antigen ~ timepoint2c,
value.var=c("BETA","RAW.P", "ADJ.PVAL"))
ans.summ.wide <- ans.summ.wide[, c(1,2,3,6,9,5,8,11,4,7,10)]
#################################
### Export models to word
#################################
library(ReporteRs)
#########################################################
## Create Tables in Word describing
#########################################################
mydoc = docx(file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_prepatentperiod_with_pvals.docx"))
xvector <- c("IgG", "IgG1", "IgG2", "IgG3", "IgG4","IgM")
for(i in 1:6){
my_text = pot(paste0("Table", i, ". Linear regression between pre-patent period and antigen (log10[MFI] scale) by isotype and timepoint independently. Only malaria cases are selected. Beta coefficient, raw p value and isotype specific adjusted p value are shown for ", xvector[i] , " isotype."))
mydoc = addParagraph( mydoc, value = my_text, stylename = "Normal")
mydoc = addParagraph( mydoc, value = "", stylename = "Normal")
ftab = FlexTable(ans.summ.wide[ans.summ.wide$ig==xvector[i],] ,add.rownames = FALSE, header.columns=FALSE,
body.par.props = parCenter(),
header.par.props = parCenter(),
body.text.props = textProperties( font.size = 8))
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","Timepoint"),
colspan = c(1,1,9), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("", "","C-1","D7","D11-D13"),
colspan = c(1,1,3,3,3), par.properties = parCenter())
ftab = addHeaderRow( ftab, text.properties = textProperties( font.size = 8, font.weight='bold'),
value = c("Isotype", "Antigen",rep(c("Beta","RAW.P","ADJ.P"),3)),
colspan = rep(1,11), par.properties = parCenter())
mydoc = addFlexTable( mydoc, ftab, par.properties = parCenter())
mydoc <- addPageBreak(mydoc)
}
writeDoc( mydoc, file = file.path(localpath, "Results", "Paper", "aim4", "mfi", "independent_timepoint_prepatentperiod_with_pvals.docx"))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.