#install have package if not installed 
#package(haven)
library(haven)
dir()

if (sum(dim(adva)==c(114365,126)) != 2) {
adva<-read_xpt("c4591001-A-D-adva.zip")
adva<-as.data.frame(adva)
 #check size
dim(adva)
 #check number of subjects
length(unique(adva$SUBJID))
}


cohorts<-c("Stage 1, Age 65 to 85, Medium dose level (30mcg), 21 Day (BNT162b2 or PBO)" , "Stage 1, Age 18 to 55, Medium dose level (30mcg), 21 Day (BNT162b2 or PBO)" )
titles<-c("Stage 1, Age 65 - 85, BNT162b2 or placebo", "Stage 1, Age 18 - 55, BNT162b2 or placebo")

 #set two graph areas
par(mar=c(7,4,4,2), mfrow = c(2,2))

for (cohort in cohorts) {

 #select 65+ age cohort
data<-adva[adva$COHORT==cohort &  adva$PARAM == "COVID-19 S1 IgG (U/mL) - Luminex Immunoassay",]

 #remove Placebos
data<-data[data$TRTP!="Placebo",]
data<-data[data$AVISIT!="",]


 #redorder the factors so the plot is in the correct order
reordered<-factor(data$AVISIT,levels(as.factor(data$AVISIT))[c(6,4,2,3,1,5)])

 #create means and standard errors for plots
stderrs<-tapply(data$AVALC,reordered,sd)/sqrt(tapply(data$AVALC, reordered, length))
means<-(tapply(as.numeric(data$AVALC),reordered,mean))
print(cohort)
print(tapply(as.numeric(data$AVALC),reordered,length))
max<-max(tapply(as.numeric(data$AVALC),reordered,length))
title<-titles[match(cohort, cohorts)]

 #plot linear with rotated text
barplot(means,xaxt="n",col='steelblue', ylim=c(0,max(means+stderrs)*1.1))->plt
title(main=paste(title,"N=", max, "(linear scale)"), cex.main=1)
arrows(plt, y0=means + stderrs, y1 = means - stderrs, angle = 90, code = 3, length = 0.1)
text(plt, par("usr")[3], labels = levels(as.factor(reordered)), srt=45, adj=c(1.1,1.1),xpd=TRUE, cex=.7)
 
 #plot log with rotated text (add 10 to values because pretreatment group log level is <0)
barplot(means+10,xaxt="n",col='steelblue', ylim=c(10,max(means+stderrs)*1.1), log="y")->plt
title(main=paste(title,"N=", max, "(log scale)"), cex.main=1)
arrows(plt, y0=means + stderrs +10, y1 = means - stderrs+10, angle = 90, code = 3, length = 0.1)
text(plt, 10^par("usr")[3], labels = levels(as.factor(reordered)), srt=45, adj=c(1.1,1.1),xpd=TRUE, cex=.7)

}