library(stringr)
noncompleters <-read.csv("enrolled_but_not_completed.csv")
noncompleters$id <-str_trim(noncompleters$id, side ="both")
##  [1] "record_id"                        "subject"                         
##  [3] "dob"                              "medical_history"                 
##  [5] "daily_meds"                       "daily_dose_of_prednisone"        
##  [7] "prednisone_dose"                  "recent_surgery"                  
##  [9] "other_vaccines"                   "date_of_vaccines"                
## [11] "covid_history"                    "covid_diagnosis_date"            
## [13] "covid_vaccine"                    "covid_19_vaccination_first_dose" 
## [15] "covid_19_vaccination_second_dose" "first_blood_dose"                
## [17] "blood_draw_second_dose"           "days_after_second_dose"          
## [19] "legal_gender"                     "daily_dose_of_prednisone.factor" 
## [21] "recent_surgery.factor"            "other_vaccines.factor"           
## [23] "covid_history.factor"             "covid_vaccine.factor"            
## [25] "legal_gender.factor"              "collection_date"                 
## [27] "collection_time"                  "cancer_type"                     
## [29] "cancer_diagnosis_date"            "cancer_surgery"                  
## [31] "surgery_date"                     "radiation"                       
## [33] "radiation_date"                   "chemotherapy_regimen"            
## [35] "treatment_date"                   "g_csf_use"                       
## [37] "g_csf_dose"                       "po_cancer_treatments"            
## [39] "targeted_treatments"              "hormonal_rx"                     
## [41] "hormonal_rx_injections"           "hormonal_rx_po"                  
## [43] "covid19_history"                  "covid19_diagnosis_date"          
## [45] "cancer_cohort_complete"           "cancer_surgery.factor"           
## [47] "radiation.factor"                 "g_csf_use.factor"                
## [49] "hormonal_rx.factor"               "cancer_cohort_complete.factor"   
## [51] "group"                            "charnum"                         
## [53] "num"                              "id"                              
## [55] "DOB"                              "datefirstshot"                   
## [57] "Age"                              "time.from.last.treatment"
library(readxl)
## Warning: package 'readxl' was built under R version 3.6.2
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.6.3
## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
## v tibble  3.0.4     v purrr   0.3.4
## v readr   1.3.1     v forcats 0.5.1
## Warning: package 'tibble' was built under R version 3.6.3
## Warning: package 'purrr' was built under R version 3.6.3
## Warning: package 'forcats' was built under R version 3.6.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date()        masks base::date()
## x dplyr::filter()          masks stats::filter()
## x lubridate::intersect()   masks base::intersect()
## x dplyr::lag()             masks stats::lag()
## x lubridate::setdiff()     masks base::setdiff()
## x Hmisc::src()             masks dplyr::src()
## x Hmisc::summarize()       masks dplyr::summarize()
## x lubridate::union()       masks base::union()
library(zoo)
## Warning: package 'zoo' was built under R version 3.6.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
####  need to remove /ml from excel file for this code to run

raw_data <- read_excel("Polished table_050921.xlsx",
                        sheet = 1)  %>%
  rename("id" = 1) %>%
  as.data.frame()
## New names:
## * `` -> ...1
## * `` -> ...3
## * `` -> ...4
## * `` -> ...5
## * `` -> ...6
## * ...
# Get draw names from first row
draw_num <- colnames(raw_data)[-1]
draw_num <- zoo::na.locf(replace(draw_num, grep("\\.{3}", draw_num), NA)) %>%
  str_replace("[:space:]", "\\_")

# Rename columns
new_names <- as.character(raw_data[1,-1]) %>%
  str_replace("\\( x10\\^6\\)", "") %>% # remove scientific notation
  str_replace("freq", "") %>% # "S1 NTD IgM MBCs" columns were inconsistently labeled
  str_trim(side = "both") %>%
  str_replace_all("[:blank:]", "_") %>% # remove blanks from names
  tolower() # all lowercase colnames so faster to type

# Create long dataset
colnames(raw_data) <- c("id", paste0(new_names, "_", draw_num))

labdata3 <- raw_data %>%
  filter(!is.na(id)) %>%
  gather(column, value, -id) %>%
  mutate(draw = case_when(str_detect(column, "Draw_1") ~ "Draw 1",
                          str_detect(column, "Draw_2") ~ "Draw 2",
                          str_detect(column, "Draw_3") ~ "Draw 3"),
         column = str_sub(column, end = -8)) %>%
  spread(column, value) %>%
  mutate(pbmcs = str_replace(pbmcs, " \\(half  v. of blood \\)|\\*", ""),
         across(unique(new_names), as.numeric)) %>%
  dplyr::select(id, draw, unique(new_names))




labdata.draw1 <- labdata3 %>% filter(
  draw=="Draw 1")

### apply Ab test cut-offs 
labdata.draw1$RBD.upperbound <- ifelse(labdata.draw1 $rbd_od >= 0.4160, "Positive", "Negative")

labdata.draw1$RBD.cutpoint <- ifelse(labdata.draw1$rbd_od < 0.1489, "Negative", "Positive")
labdata.draw1$S2.cutpoint <-ifelse(labdata.draw1$s2_od < 0.195, "Negative", "Positive")

labdata.draw1$TestResult <- ifelse(labdata.draw1$RBD.upperbound =="Negative" & labdata.draw1$RBD.cutpoint == "Positive" & labdata.draw1$S2.cutpoint == "Positive", "Positive", labdata.draw1$RBD.upperbound)
labdata.draw1$TestResult <-toupper(labdata.draw1$TestResult)

seropos <- subset(labdata.draw1, TestResult=="POSITIVE")
sero.pos <- seropos[, 1]
seroneg <- subset(labdata.draw1, TestResult == "NEGATIVE")
sero.neg <-seroneg[, 1]

#### only one person who reported being positive for COVID-19 in the past tested positive using ####  UA/AZ State Ab assay cut-points
#######  LAB TABLE OF DATA

labdata2 <- labdata3 %>% filter(id %notin% sero.pos) %>% filter(id %notin% moderna.id) %>% filter(id %notin% noncompleters)

##labdata2 <- labdata3 %>% filter(id %in% sero.neg) %>% filter(id %in% pfizer.id) why doesn't 
labdata <-as.data.frame(labdata2)
labdata$group <- substr(labdata$id, 1, 2)

colnames(labdata)[c(4,10:13, 18:21)] <- c("per_b_cells", "spike_t_cells", "total_t_cells", "dn2_rbd", "dn2_s1", "sw_neg_rbd", "sw_neg_s1","sw_pos_rbd", "sw_pos_s1" )


lab.v1 <-subset(labdata, draw == "Draw 1")
lab.v2 <-subset(labdata, draw == "Draw 2")
lab.v3 <-subset(labdata, draw == "Draw 3")


#no.miss <- labdata %>% drop_na(rbd_od)

no.miss.v1 <- (lab.v1%>% drop_na(rbd_od))[, c(1:2, 22)]
no.miss.v2 <- (lab.v2 %>% drop_na(rbd_od))[, c(1:2, 22)]
no.miss.v3 <- (lab.v3 %>% drop_na(rbd_od))[, c(1:2, 22)]

### 113, 107, 101

total.nomiss1 <- inner_join(no.miss.v1, no.miss.v2, by="id")
total.nomiss <- inner_join(total.nomiss1, no.miss.v3, by="id")



with(total.nomiss,table(draw, group))
##         group
## draw     CC HC
##   Draw 3 51 47
missing.data <- labdata %>% filter(id %notin% total.nomiss$id)
unique(missing.data$id)
##  [1] "CC22" "CC29" "CC45" "CC48" "CC55" "CC56" "CC57" "CC59" "HC2"  "HC22"
## [11] "HC29" "HC32" "HC44" "HC46" "HC54" "HC57" "HC59" "HC60" "HC64" "HC68"
## [21] "HC69" "HC74" "HC77"
complete.labdata <-labdata %>% filter(id %in% total.nomiss$id)
complete.labdata.noid <- complete.labdata[, -1]

table1( ~ . |group*draw, data=complete.labdata.noid, overall=FALSE)
CC
HC
Draw 1
(n=51)		 Draw 2
(n=51)		 Draw 3
(n=51)		 Draw 1
(n=47)		 Draw 2
(n=47)		 Draw 3
(n=47)
draw
Draw 1 51 (100%) 0 (0%) 0 (0%) 47 (100%) 0 (0%) 0 (0%)
Draw 2 0 (0%) 51 (100%) 0 (0%) 0 (0%) 47 (100%) 0 (0%)
Draw 3 0 (0%) 0 (0%) 51 (100%) 0 (0%) 0 (0%) 47 (100%)
pbmcs
Mean (SD) 2.74 (2.35) 3.70 (3.94) 2.87 (1.83) 2.81 (0.838) 2.74 (0.830) 2.71 (0.995)
Median [Min, Max] 2.20 [0.625, 16.8] 2.40 [0.950, 24.6] 2.53 [0.750, 11.5] 2.60 [1.58, 5.01] 2.62 [1.66, 5.55] 2.54 [1.43, 6.10]
Missing 0 (0%) 0 (0%) 1 (2.0%) 0 (0%) 0 (0%) 0 (0%)
per_b_cells
Mean (SD) 10.3 (7.00) 9.71 (7.06) 10.0 (7.10) 18.8 (5.73) 16.8 (5.10) 15.4 (5.35)
Median [Min, Max] 9.28 [0.280, 26.3] 7.66 [0.250, 31.8] 8.48 [0.150, 28.0] 18.9 [7.29, 33.0] 16.1 [8.04, 28.4] 14.2 [8.51, 30.6]
Missing 3 (5.9%) 2 (3.9%) 2 (3.9%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
%_myeloid_cells
Mean (SD) 8.42 (5.88) 10.5 (10.1) 9.84 (8.32) 4.91 (3.02) 3.99 (2.11) 4.04 (2.94)
Median [Min, Max] 6.60 [1.86, 25.0] 8.10 [0.710, 46.2] 6.11 [1.01, 31.5] 4.33 [1.01, 17.6] 4.12 [0.840, 8.96] 3.32 [0.500, 14.2]
Missing 3 (5.9%) 2 (3.9%) 2 (3.9%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
rbd_od
Mean (SD) 0.0581 (0.00939) 0.171 (0.164) 0.408 (0.317) 0.0654 (0.0426) 0.457 (0.188) 0.819 (0.152)
Median [Min, Max] 0.0550 [0.0440, 0.0920] 0.0930 [0.0480, 0.702] 0.326 [0.0500, 1.08] 0.0540 [0.0450, 0.320] 0.491 [0.0500, 0.842] 0.854 [0.379, 1.08]
s2_od
Mean (SD) 0.0864 (0.0364) 0.156 (0.0983) 0.236 (0.164) 0.109 (0.0638) 0.287 (0.117) 0.400 (0.115)
Median [Min, Max] 0.0770 [0.0510, 0.238] 0.131 [0.0520, 0.545] 0.176 [0.0560, 0.800] 0.0920 [0.0530, 0.342] 0.273 [0.0790, 0.595] 0.404 [0.163, 0.700]
rbd_auc
Mean (SD) NA (NA) 0.175 (0.100) 0.338 (0.291) NA (NA) 0.388 (0.209) 0.895 (0.273)
Median [Min, Max] NA [NA, NA] 0.133 [0.106, 0.620] 0.194 [0.104, 1.31] NA [NA, NA] 0.312 [0.137, 0.958] 0.941 [0.226, 1.31]
Missing 51 (100%) 0 (0%) 1 (2.0%) 47 (100%) 0 (0%) 0 (0%)
neut_titer
Mean (SD) NA (NA) 40.8 (81.7) 159 (267) NA (NA) 153 (262) 1490 (2410)
Median [Min, Max] NA [NA, NA] 20.0 [0.00, 540] 60.0 [0.00, 1620] NA [NA, NA] 60.0 [0.00, 1620] 540 [60.0, 14600]
Missing 51 (100%) 1 (2.0%) 3 (5.9%) 47 (100%) 0 (0%) 0 (0%)
spike_t_cells
Mean (SD) 3.91 (7.04) 7.24 (17.9) 20.4 (29.2) 5.94 (11.7) 14.4 (17.7) 31.2 (31.4)
Median [Min, Max] 0.400 [0.00, 30.0] 1.60 [0.00, 91.2] 6.40 [0.00, 145] 2.40 [0.00, 64.4] 7.00 [0.00, 66.4] 18.4 [0.00, 116]
Missing 0 (0%) 0 (0%) 0 (0%) 0 (0%) 1 (2.1%) 0 (0%)
total_t_cells
Mean (SD) 179 (105) 164 (72.2) 166 (103) 211 (122) 175 (67.6) 205 (111)
Median [Min, Max] 172 [0.400, 437] 161 [24.0, 327] 146 [0.00, 407] 190 [0.00, 604] 179 [31.2, 327] 182 [0.00, 574]
dn2_rbd
Mean (SD) 0.00 (0.00) 6.34e-06 (3.33e-05) 0.000128 (0.000253) 1.02e-05 (4.08e-05) 5.15e-05 (0.000110) 0.000902 (0.00110)
Median [Min, Max] 0.00 [0.00, 0.00] 0.00 [0.00, 0.000213] 0.00 [0.00, 0.00102] 0.00 [0.00, 0.000214] 0.00 [0.00, 0.000544] 0.000406 [0.00, 0.00491]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
dn2_s1
Mean (SD) 2.55e-05 (0.000146) 3.57e-05 (8.35e-05) 9.8e-05 (0.000227) 1.26e-05 (4.32e-05) 8.52e-05 (0.000123) 0.000638 (0.000899)
Median [Min, Max] 0.00 [0.00, 0.00101] 0.00 [0.00, 0.000344] 0.00 [0.00, 0.00118] 0.00 [0.00, 0.000214] 0.00 [0.00, 0.000544] 0.000280 [0.00, 0.00460]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
%_unswitched_memory__cd21-_rbd+
Mean (SD) 5.87e-05 (0.000130) 7.33e-05 (0.000151) 0.000109 (0.000289) 4.12e-05 (9e-05) 9.33e-05 (0.000117) 0.000263 (0.000543)
Median [Min, Max] 0.00 [0.00, 0.000563] 0.00 [0.00, 0.000665] 0.00 [0.00, 0.00189] 0.00 [0.00, 0.000445] 5.93e-05 [0.00, 0.000514] 0.000139 [0.00, 0.00321]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
%_unswitched_memory__cd21-_s1+
Mean (SD) 0.000308 (0.000481) 0.000431 (0.000952) 0.000378 (0.000613) 0.000431 (0.000510) 0.000548 (0.000550) 0.000865 (0.00101)
Median [Min, Max] 0.000140 [0.00, 0.00204] 0.000151 [0.00, 0.00616] 0.000126 [0.00, 0.00248] 0.000244 [0.00, 0.00205] 0.000424 [0.00, 0.00216] 0.000524 [0.00, 0.00371]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
%_unswitched_memory_cd21+_rbd+
Mean (SD) 5.93e-05 (0.000146) 9.52e-05 (0.000213) 0.000118 (0.000252) 0.000261 (0.000226) 0.000445 (0.000474) 0.000488 (0.000547)
Median [Min, Max] 0.00 [0.00, 0.000640] 0.00 [0.00, 0.00123] 0.00 [0.00, 0.00126] 0.000213 [0.00, 0.000769] 0.000267 [0.00, 0.00193] 0.000352 [0.00, 0.00275]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
%__unswitched_memory__cd21+_s1+
Mean (SD) 0.00268 (0.0135) 0.00218 (0.0103) 0.00525 (0.0324) 0.00188 (0.00104) 0.00184 (0.00121) 0.00183 (0.00149)
Median [Min, Max] 0.000455 [0.00, 0.0950] 0.000430 [0.00, 0.0730] 0.000393 [0.00, 0.230] 0.00186 [0.000174, 0.00428] 0.00187 [0.000118, 0.00731] 0.00155 [0.000169, 0.00836]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
sw_neg_rbd
Mean (SD) 4.54e-05 (0.000221) 3.32e-06 (2.35e-05) 0.000309 (0.000485) 1.75e-05 (6.94e-05) 0.000156 (0.000256) 0.00310 (0.00411)
Median [Min, Max] 0.00 [0.00, 0.00149] 0.00 [0.00, 0.000166] 0.000107 [0.00, 0.00201] 0.00 [0.00, 0.000428] 0.00 [0.00, 0.00120] 0.00157 [0.00, 0.0200]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
sw_neg_s1
Mean (SD) 3.83e-05 (0.000117) 4.45e-05 (0.000106) 0.000236 (0.000474) 4.02e-05 (6.8e-05) 0.000225 (0.000482) 0.00251 (0.00407)
Median [Min, Max] 0.00 [0.00, 0.000691] 0.00 [0.00, 0.000532] 0.00 [0.00, 0.00259] 0.00 [0.00, 0.000276] 7.93e-05 [0.00, 0.00286] 0.000986 [0.00, 0.0200]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
sw_pos_rbd
Mean (SD) 2.11e-05 (5.59e-05) 5.81e-05 (0.000168) 0.000360 (0.00132) 8.49e-05 (0.000102) 0.000569 (0.000798) 0.00205 (0.00398)
Median [Min, Max] 0.00 [0.00, 0.000249] 0.00 [0.00, 0.00109] 0.00 [0.00, 0.00918] 7.86e-05 [0.00, 0.000454] 0.000361 [0.00, 0.00453] 6e-04 [0.00, 0.0190]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
sw_pos_s1
Mean (SD) 0.000932 (0.00426) 0.000677 (0.00280) 0.00162 (0.00846) 0.000977 (0.000901) 0.00116 (0.000784) 0.00120 (0.00186)
Median [Min, Max] 0.000174 [0.00, 0.0300] 0.000243 [0.00, 0.0200] 0.000224 [0.00, 0.0600] 0.000638 [0.00, 0.00400] 0.00110 [0.00, 0.00288] 0.000633 [0.00, 0.0120]
Missing 2 (3.9%) 1 (2.0%) 1 (2.0%) 1 (2.1%) 1 (2.1%) 1 (2.1%)
group
CC 51 (100%) 51 (100%) 51 (100%) 0 (0%) 0 (0%) 0 (0%)
HC 0 (0%) 0 (0%) 0 (0%) 47 (100%) 47 (100%) 47 (100%) 
######  CLINICAL TABLE 1 



All.Clinical.completers <- All.Clinical %>% filter(id %notin% sero.pos) %>% filter(id %notin% moderna.id) %>% filter(id %notin% noncompleters$id)

table(All.Clinical.completers$group)
## 
## CC HC 
## 57 50
table(All.Clinical$group)
## 
## CC HC 
## 60 73
#### TABLE OF ALL COMPLETERS CLINICAL

#table1(~ Age + legal_gender.factor + daily_dose_of_prednisone.factor + recent_surgery.factor + other_vaccines.factor + covid_history.factor + time.from.last.treatment + radiation.factor + hormonal_rx.factor  | group, data=All.Clinical.completers)
####### Repeated Measures Analysis
complete.labdata$draw.cont <- as.numeric(substr(complete.labdata$draw, 6, 7))
draw2HC.draw3CC <- subset(complete.labdata, 
                          (group=="CC" & draw=="Draw 3") | 
                            (group=="HC" & draw=="Draw 2"))



library(nlme)
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(lme4)
## Warning: package 'lme4' was built under R version 3.6.3
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
library(emmeans)
## Warning: package 'emmeans' was built under R version 3.6.3
## 
## Attaching package: 'emmeans'
## The following object is masked from 'package:rms':
## 
##     contrast
library(effects)
## Warning: package 'effects' was built under R version 3.6.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 3.6.3
## Use the command
##     lattice::trellis.par.set(effectsTheme())
##   to customize lattice options for effects plots.
## See ?effectsTheme for details.
##### Analysis for RBD.OD ######

rbd_od <-lme(rbd_od ~ group*draw, random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(rbd_od)
##             numDF denDF  F-value p-value
## (Intercept)     1   192 657.9413  <.0001
## group           1    96  85.7435  <.0001
## draw            2   192 299.5311  <.0001
## group:draw      2   192  43.1376  <.0001
rbd.od.lme <- emmeans(rbd_od, ~ group|draw)
contrasts.rbd.od <- contrast(rbd.od.lme, interaction = "pairwise")
pairs(contrasts.rbd.od, by = NULL)
##  contrast                            estimate     SE  df t.ratio p.value
##  (CC - HC Draw 1) - (CC - HC Draw 2)    0.279 0.0445 192 6.269   <.0001 
##  (CC - HC Draw 1) - (CC - HC Draw 3)    0.404 0.0445 192 9.070   <.0001 
##  (CC - HC Draw 2) - (CC - HC Draw 3)    0.125 0.0445 192 2.801   0.0154 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot.rbd.od_lme <- effect(term="group*draw",rbd_od, multiline=FALSE)
plot(plot.rbd.od_lme)

t.test(rbd_od ~ factor(group), data=draw2HC.draw3CC)
## 
##  Welch Two Sample t-test
## 
## data:  rbd_od by factor(group)
## t = -0.94283, df = 82.446, p-value = 0.3485
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.15325337  0.05469101
## sample estimates:
## mean in group CC mean in group HC 
##        0.4078039        0.4570851
######  RBD AUC

rbd_auc <-lme(rbd_auc ~ group*draw,random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(rbd_auc)
##             numDF denDF  F-value p-value
## (Intercept)     1    96 479.9436  <.0001
## group           1    96  91.3413  <.0001
## draw            1    95 202.1294  <.0001
## group:draw      1    95  54.4792  <.0001
rbd.auc.lme <- emmeans(rbd_auc, ~ group|draw)
contrasts.rbd.auc <- contrast(rbd.auc.lme, interaction = "pairwise")
pairs(contrasts.rbd.auc, by = NULL)
##  contrast                            estimate     SE df t.ratio p.value
##  (CC - HC Draw 2) - (CC - HC Draw 3)    0.343 0.0465 95 7.381   <.0001 
## 
## Degrees-of-freedom method: containment
plot.rbd.auc_lme <- effect(term="group*draw", rbd_auc, multiline=TRUE)
plot(plot.rbd.auc_lme)

#####  S2_OD

s2_od <-lme(s2_od ~ group*draw, random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(s2_od)
##             numDF denDF  F-value p-value
## (Intercept)     1   192 701.7810  <.0001
## group           1    96  44.7025  <.0001
## draw            2   192 145.9409  <.0001
## group:draw      2   192  16.9892  <.0001
s2_od.lme <- emmeans(s2_od, ~ group|draw)
contrasts.s2_od <- contrast(s2_od.lme, interaction = "pairwise")
pairs(contrasts.s2_od, by = NULL)
##  contrast                            estimate     SE  df t.ratio p.value
##  (CC - HC Draw 1) - (CC - HC Draw 2)   0.1092 0.0256 192 4.272   0.0001 
##  (CC - HC Draw 1) - (CC - HC Draw 3)   0.1424 0.0256 192 5.571   <.0001 
##  (CC - HC Draw 2) - (CC - HC Draw 3)   0.0332 0.0256 192 1.298   0.3978 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot.s2_od_lme <- effect(term="group*draw", s2_od, multiline=FALSE)
plot(plot.s2_od_lme)

t.test(s2_od ~ factor(group), data=draw2HC.draw3CC)
## 
##  Welch Two Sample t-test
## 
## data:  s2_od by factor(group)
## t = -1.8039, df = 90.649, p-value = 0.07457
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.108319990  0.005217779
## sample estimates:
## mean in group CC mean in group HC 
##        0.2357255        0.2872766
#####



#y <- as.data.frame(ef2)
#ggplot(y , aes(draw.cont, fit, color=group)) + geom_point() + geom_errorbar(aes(ymin=fit-se, #ymax=fit+se), width=0.4) + 

#####  Spike t cells



spike.t.cells.model <-lme(log(spike_t_cells + 1) ~ group*draw, random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(spike.t.cells.model)
##             numDF denDF  F-value p-value
## (Intercept)     1   191 337.7918  <.0001
## group           1    96  14.4453  0.0003
## draw            2   191  54.8096  <.0001
## group:draw      2   191   3.5996  0.0292
spike.t.cells.lme <- emmeans(spike.t.cells.model, ~ group|draw)
contrasts.spike.t.cells <- contrast(spike.t.cells.lme, interaction = "pairwise")
pairs(contrasts.spike.t.cells, by = NULL)
##  contrast                            estimate    SE  df t.ratio p.value
##  (CC - HC Draw 1) - (CC - HC Draw 2)    0.711 0.276 191  2.573  0.0290 
##  (CC - HC Draw 1) - (CC - HC Draw 3)    0.534 0.275 191  1.940  0.1302 
##  (CC - HC Draw 2) - (CC - HC Draw 3)   -0.177 0.276 191 -0.640  0.7983 
## 
## Degrees-of-freedom method: containment 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot.spike.t.cells_lme <- effect(term="group*draw", spike.t.cells.model, multiline=FALSE)
plot(plot.spike.t.cells_lme)

spike.t.cells.model
## Linear mixed-effects model fit by REML
##   Data: complete.labdata 
##   Log-restricted-likelihood: -457.1185
##   Fixed: log(spike_t_cells + 1) ~ group * draw 
##        (Intercept)            groupHC         drawDraw 2         drawDraw 3 
##          0.9121113          0.2989558          0.2352406          1.1749299 
## groupHC:drawDraw 2 groupHC:drawDraw 3 
##          0.7112180          0.5343735 
## 
## Random effects:
##  Formula: ~1 | id
##         (Intercept)  Residual
## StdDev:    0.742721 0.9630254
## 
## Number of Observations: 293
## Number of Groups: 98
draw2HC.draw3CC$log.spike.t.cells <- log(draw2HC.draw3CC$spike_t_cells + 1)
t.test(log(spike_t_cells + 1)~ factor(group), data=draw2HC.draw3CC)
## 
##  Welch Two Sample t-test
## 
## data:  log(spike_t_cells + 1) by factor(group)
## t = -0.22273, df = 91.428, p-value = 0.8242
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.595416  0.475347
## sample estimates:
## mean in group CC mean in group HC 
##         2.087041         2.147076
###### FIGURE 4 Memory B cell responses

time.3.data <-subset(complete.labdata, draw=="Draw 3")

plot(asin(sqrt(sw_neg_rbd)) ~ factor(group), data=time.3.data)

plot(asin(sqrt(sw_neg_s1)) ~ factor(group), data=time.3.data)

plot(asin(sqrt(sw_pos_rbd)) ~ factor(group), data=time.3.data)

plot(asin(sqrt(sw_pos_s1)) ~ factor(group), data=time.3.data)

### SW neg rbd
sw.neg.rbd <-lme(asin(sqrt(sw_neg_rbd)) ~ group*draw, random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(sw.neg.rbd)
##             numDF denDF   F-value p-value
## (Intercept)     1   187 144.54433  <.0001
## group           1    94  53.54893  <.0001
## draw            2   187 102.47160  <.0001
## group:draw      2   187  37.73968  <.0001
sw.neg.rbd.lme <- emmeans(sw.neg.rbd, ~ group|draw)
## Warning in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : There are unevaluated constants in the response formula
## Auto-detection of the response transformation may be incorrect
contrasts.sw.neg.rbd<- contrast(sw.neg.rbd.lme, interaction = "pairwise")
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
pairs(contrasts.sw.neg.rbd, by = NULL)
##  contrast                            estimate      SE  df t.ratio p.value
##  (CC - HC Draw 1) - (CC - HC Draw 2)  0.00825 0.00414 187 1.991   0.1171 
##  (CC - HC Draw 1) - (CC - HC Draw 3)  0.03445 0.00414 187 8.312   <.0001 
##  (CC - HC Draw 2) - (CC - HC Draw 3)  0.02620 0.00413 187 6.337   <.0001 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot.sw.neg.rbd_lme <- effect(term="group*draw", sw.neg.rbd, multiline=FALSE)
plot(plot.sw.neg.rbd_lme)

### SW neg s1
sw.neg.s1 <-lme(asin(sqrt(sw_neg_s1)) ~ group*draw, random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(sw.neg.s1)
##             numDF denDF   F-value p-value
## (Intercept)     1   187 131.62761  <.0001
## group           1    94  45.35522  <.0001
## draw            2   187  69.73331  <.0001
## group:draw      2   187  34.20721  <.0001
sw.neg.s1.lme <- emmeans(sw.neg.s1, ~ group|draw)
## Warning in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : There are unevaluated constants in the response formula
## Auto-detection of the response transformation may be incorrect
contrasts.sw.neg.s1<- contrast(sw.neg.s1.lme, interaction = "pairwise")
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
pairs(contrasts.sw.neg.s1, by = NULL)
##  contrast                            estimate      SE  df t.ratio p.value
##  (CC - HC Draw 1) - (CC - HC Draw 2)  0.00533 0.00388 187 1.374   0.3567 
##  (CC - HC Draw 1) - (CC - HC Draw 3)  0.03006 0.00388 187 7.741   <.0001 
##  (CC - HC Draw 2) - (CC - HC Draw 3)  0.02472 0.00387 187 6.385   <.0001 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot.sw.neg.s1_lme <- effect(term="group*draw", sw.neg.s1, multiline=FALSE)
plot(plot.sw.neg.s1_lme)

### SW pos rbd
sw.pos.rbd <-lme(asin(sqrt(sw_pos_rbd)) ~ group*draw, random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(sw.pos.rbd)
##             numDF denDF   F-value p-value
## (Intercept)     1   187 117.06297  <.0001
## group           1    94  44.64199  <.0001
## draw            2   187  38.27652  <.0001
## group:draw      2   187  11.84580  <.0001
sw.pos.rbd.lme <- emmeans(sw.pos.rbd, ~ group|draw)
## Warning in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : There are unevaluated constants in the response formula
## Auto-detection of the response transformation may be incorrect
contrasts.sw.pos.rbd<- contrast(sw.pos.rbd.lme, interaction = "pairwise")
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
pairs(contrasts.sw.pos.rbd, by = NULL)
##  contrast                            estimate      SE  df t.ratio p.value
##  (CC - HC Draw 1) - (CC - HC Draw 2)  0.01092 0.00399 187 2.737   0.0186 
##  (CC - HC Draw 1) - (CC - HC Draw 3)  0.01937 0.00399 187 4.855   <.0001 
##  (CC - HC Draw 2) - (CC - HC Draw 3)  0.00845 0.00398 187 2.124   0.0877 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot.sw.pos.rbd_lme <- effect(term="group*draw", sw.pos.rbd, multiline=FALSE)
plot(plot.sw.pos.rbd_lme)

### SW pos s1
sw.pos.s1 <-lme(asin(sqrt(sw_pos_s1)) ~ group*draw, random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(sw.pos.s1)
##             numDF denDF   F-value p-value
## (Intercept)     1   187 120.06110  <.0001
## group           1    94   9.06439  0.0033
## draw            2   187   0.53450  0.5869
## group:draw      2   187   1.71273  0.1832
sw.pos.s1.lme <- emmeans(sw.pos.s1, ~ group|draw)
## Warning in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : There are unevaluated constants in the response formula
## Auto-detection of the response transformation may be incorrect
contrasts.sw.pos.s1<- contrast(sw.pos.s1.lme, interaction = "pairwise")
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
pairs(contrasts.sw.pos.s1, by = NULL)
##  contrast                            estimate      SE  df t.ratio p.value
##  (CC - HC Draw 1) - (CC - HC Draw 2)  0.00471 0.00355 187  1.326  0.3826 
##  (CC - HC Draw 1) - (CC - HC Draw 3) -0.00159 0.00355 187 -0.448  0.8953 
##  (CC - HC Draw 2) - (CC - HC Draw 3) -0.00630 0.00354 187 -1.780  0.1790 
## 
## Degrees-of-freedom method: containment 
## P value adjustment: tukey method for comparing a family of 3 estimates
plot.sw.pos.s1_lme <- effect(term="group*draw", sw.pos.s1, multiline=FALSE)
plot(plot.sw.pos.s1_lme)

#####  neut titer
#plot(log10(neut_titer + 1) ~ factor(group), data=time.3.data)


neut.titer.model <-lme(log10(neut_titer + 1) ~ group*draw, random= ~1|id, data=complete.labdata, na.action = na.omit)
anova(neut.titer.model)
##             numDF denDF  F-value p-value
## (Intercept)     1    95 702.0269  <.0001
## group           1    95  48.4325  <.0001
## draw            1    93 225.7382  <.0001
## group:draw      1    93  14.5866   2e-04
neut.titer.lme <- emmeans(neut.titer.model, ~ group|draw)
contrasts.neut.titer <- contrast(neut.titer.lme, interaction = "pairwise")
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates
pairs(contrasts.neut.titer, by = NULL)
##  contrast                            estimate    SE df t.ratio p.value
##  (CC - HC Draw 2) - (CC - HC Draw 3)    0.398 0.104 93 3.819   0.0002 
## 
## Degrees-of-freedom method: containment
plot.contrasts.neut.titer_lme <- effect(term="group*draw", neut.titer.model, multiline=FALSE)
plot(plot.contrasts.neut.titer_lme)

#stargazer(s2_od, rbd_od, rbd_auc, type="html", align=TRUE)
unique.lab <-unique(complete.labdata$id)

nomatch <- All.Clinical.completers %>% filter(id %notin%unique.lab)


with(complete.labdata, table(draw,group))
##         group
## draw     CC HC
##   Draw 1 51 47
##   Draw 2 51 47
##   Draw 3 51 47
data.all <-left_join(complete.labdata, All.Clinical.completers, by="id")
 
dim(data.all)
## [1] 297  80
data.all <- data.all[-which(data.all$id == "CC19")[c(1,3,5)],]

dim(data.all)
## [1] 294  80
with(data.all, table(draw,group.y))
##         group.y
## draw     CC HC
##   Draw 1 51 47
##   Draw 2 51 47
##   Draw 3 51 47
#### TABLE WITH ALL DATA COMPLETE PLUS LAB

for.table <-subset(data.all, draw=="Draw 3")

table1(~ Age + legal_gender.factor + daily_dose_of_prednisone.factor + recent_surgery.factor + other_vaccines.factor + covid_history.factor + time.from.last.treatment + radiation.factor + hormonal_rx.factor  | group.y, data=for.table)
CC
(n=51)		 HC
(n=47)		 Overall
(n=98)
Age
Mean (SD) 63.2 (9.63) 40.8 (16.8) 52.4 (17.6)
Median [Min, Max] 62.9 [38.5, 85.0] 36.6 [19.5, 81.5] 55.2 [19.5, 85.0]
legal_gender.factor
Female 41 (80.4%) 31 (66.0%) 72 (73.5%)
Male 10 (19.6%) 16 (34.0%) 26 (26.5%)
daily_dose_of_prednisone.factor
Yes 1 (2.0%) 0 (0%) 1 (1.0%)
No 50 (98.0%) 47 (100%) 97 (99.0%)
recent_surgery.factor
Yes 2 (3.9%) 0 (0%) 2 (2.0%)
No 49 (96.1%) 47 (100%) 96 (98.0%)
other_vaccines.factor
Yes 0 (0%) 0 (0%) 0 (0%)
No 51 (100%) 47 (100%) 98 (100%)
covid_history.factor
Yes 3 (5.9%) 1 (2.1%) 4 (4.1%)
No 48 (94.1%) 46 (97.9%) 94 (95.9%)
time.from.last.treatment
Mean (SD) 4.76 (7.01) NA (NA) 4.76 (7.01)
Median [Min, Max] 2.00 [-8.00, 23.0] NA [NA, NA] 2.00 [-8.00, 23.0]
Missing 18 (35.3%) 47 (100%) 65 (66.3%)
radiation.factor
Yes 18 (35.3%) 0 (0%) 18 (18.4%)
No 32 (62.7%) 0 (0%) 32 (32.7%)
Missing 1 (2.0%) 47 (100%) 48 (49.0%)
hormonal_rx.factor
Yes 8 (15.7%) 0 (0%) 8 (8.2%)
No 43 (84.3%) 0 (0%) 43 (43.9%)
Missing 0 (0%) 47 (100%) 47 (48.0%) 
data.all <- subset(data.all, covid_vaccine.factor != "NA")

#pdf("pfizer versus moderna.pdf")

p <- ggplot(data=data.all, aes(x = draw, y = rbd_od, group=id ))
p +  stat_smooth(aes(group = 1)) + stat_summary(aes(group = 1),
                                                                           geom = "point", fun.y = mean, shape = 17, size = 3) + facet_grid(. ~ covid_vaccine.factor)+ xlab("Blood Draws (shot 1, shot 2, 7 days post shot 2)") +
  labs(caption = "RBD.OD by Vaccine (HCW)") + theme(plot.caption = element_text(hjust = 0, face= "italic"), panel.background = element_rect(fill = 'white', colour = 'black'))
## Warning: `fun.y` is deprecated. Use `fun` instead.
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at 0.99
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 2.01
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.4271e-015
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 4.0401
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : pseudoinverse used at
## 0.99
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : neighborhood radius 2.01
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : reciprocal condition
## number 1.4271e-015
## Warning in predLoess(object$y, object$x, newx = if
## (is.null(newdata)) object$x else if (is.data.frame(newdata))
## as.matrix(model.frame(delete.response(terms(object)), : There are other near
## singularities as well. 4.0401

#dev.off()