# Analysis of the dataset from Elgazzar et al. (2020), 10.21203/rs.3.rs-100956/v3.
# By Nick Brown (nicholasjlbrown@gmail.com), July 2021.
# Licence: CC-0.

# The preprint sends the reader here to get the data file: https://filetransfer.io/data-package/qGiU0mw6#link
# If you put that file in the R working directory and run this script, that should be all you need.

library(adagio)
library(ez)
library(readxl)
library(stringr)

filename <- "Copy of covid_19 final master sheet12-12-2020 (1).xlsx"
rawdata <- as.data.frame(readxl::read_excel(filename))

# Remove the trailing notes, plus row 203 (Excel row 204) which has the name "MM" but no other data.
data <- rawdata[-c(203, 602:608),]

data$GroupNumber <- c(rep(1, 100), rep(2, 100), rep(3, 100), rep(4, 100), rep(5, 100), rep(6, 100))
data$pID <- 1:600

# Function to convert a row number from our data frame to the corresponding line number in the Excel sheet.
xlrow <- function (datarow)
{
  return(datarow + (if (datarow > 202) 2 else 1))    # 1 for the header, 1 for row 203 which we removed
}

# Excel file line 5, data frame row 4: The Excel cell K4 contains "06-Jan", because it is formatted
#  as a date. But read_excel() manages to see the underlying number (6.3), so that doesn't need work here.

# Excel file line 300, data frame row 298: The TLC value appears to be missing a decimal point.
glitch <- as.numeric(data[298, 11])
if (glitch > 50) {                    # this if() provides a sanity check
  data[298, 11] <- glitch / 10
}

# Excel file line 302, data frame row 300: It looks like two records have got smashed together here.
# Perhaps that's why the authors reduced the height of the line to almost nothing.
# I've chosen to ignore the values in Excel columns O through Q (77.40%, 1,	11.80%),
#  move the 214 from column R to column P ("serum ferritin one week after treament"),
#  and impute the value in column 14 ("serum ferritin before treament") from the mean of the other 99
#  patients in the same group. Other principled choices would have been possible here.
# I hope the authors will be able to tell us what they did when converting the data to SPSS format.
data[300, 16] <- data[300, 18]
data[300, c(15, 17, 18)] <- NA
data[300, 14] <- mean(data[301:399, 14])

# Excel file line 403, data frame row 401: There are two extra columns, apparently caused by
#  whoever entered the data adding the spurious keywords "comorbidities" and "prognosis" as cells.
data[401, 5] <- data[401, 6]
data[401, 7] <- data[401, 8]
data[401, c(6, 8)] <- NA

# Excel file line 105, data frame row 104: HGB seems to be missing its decimal point.
glitch <- as.numeric(data[104, 19])
if (glitch > 1) {                    # this if() provides a sanity check
  data[104, 19] <- glitch / 100
}

# Excel file line 170, data frame row 169: The discharge date 15/5/2020 is after the admission date.
# Looking at the other data it seems that this should probably say 15/6/2020.
data[169, 24] <- sub("15/5/2020", "15/6/2020", data[169, 24])

# Excel file line 449, data frame row 447: Sex is recorded as "A". Nothing in the paper suggests
#  that patients identified as anything other than M or F. We will set it to "M" as this gives the
#  authors the benefit of the doubt (the resulting M/F numbers for this group match the paper).
data[447, 3] <- sub("A", "M", data[447, 3])

# In a few cases (e.g., Excel row 8), columns X and Y are merged. We split them here.
for (i in 1:nrow(data)) {
  colX <- data[i, 24]
  if (grepl("/2020 +\\d+ +[A-Za-z]", colX)) {
    data[i, 24] <- sub("/2020.*", "/2020", colX)
    data[i, 25] <- sub(".*/2020", "", colX)
  }
}

# Identify comorbidities from the various strings in the file.
data[401:600, 10] <- data[401:600, 5]   # outpatient comorbidities are in a different column
data[401:600, 5] <- NA

data.comorb <- unlist(data["comorbidities"])
data$Como_BA <- ( !is.na(str_extract(data.comorb, "B[\\.,]?A"))
                  | !is.na(str_extract(data.comorb, "Bronchial asthma"))
                )
data$Como_CHO <- !is.na(str_extract(data.comorb, "cholecystitis"))
data$Como_CKD <- !is.na(str_extract(data.comorb, "CKD"))
data$Como_DM <- !is.na(str_extract(data.comorb, "D[\\.,]?M"))
data$Como_HBV <- !is.na(str_extract(data.comorb, "HBV"))
data$Como_HCV <- !is.na(str_extract(data.comorb, "HCV"))
data$Como_HTN <- !is.na(str_extract(data.comorb, "HTN"))
data$Como_IHD <- !is.na(str_extract(data.comorb, "IHD"))
data$Como_OHS <- !is.na(str_extract(data.comorb, "Open HT surgery"))
data$Como_PRG <- !is.na(str_extract(data.comorb, "pregnant"))
data$Como_PE <- !is.na(str_extract(data.comorb, "pulmonary embolism"))

# Identify symptoms from the various strings in the file.
data.symp <- unlist(data["other symptoms"])
data$Symp_ABDO <- !is.na(str_extract(data.symp, "abd[a-z]* pain"))
data$Symp_ANOS <- !is.na(str_extract(data.symp, "anosmia"))
data$Symp_BODY <- !is.na(str_extract(data.symp, "body *ache"))
data$Symp_CHEST <- !is.na(str_extract(data.symp, "chest pain"))
# Interestingly, "cough" spelt "coguh" when and only when it accompanies "sputum".
data$Symp_COUGH <- ( !is.na(str_extract(data.symp, regex("cough", ignore_case=TRUE)))
                   | !is.na(str_extract(data.symp, "coguh"))
                   )
data$Symp_DIARR <- !is.na(str_extract(data.symp, regex("diarrhea", ignore_case=TRUE)))
data$Symp_DROW <- !is.na(str_extract(data.symp, "drowsy"))
data$Symp_DYSP <- !is.na(str_extract(data.symp, "dyspnea"))
data$Symp_EPIG <- !is.na(str_extract(data.symp, "epigastric pain"))
data$Symp_HEAD <- !is.na(str_extract(data.symp, "headache"))
data$Symp_JOINT <- !is.na(str_extract(data.symp, "joint pain"))
data$Symp_LOIN <- !is.na(str_extract(data.symp, "loin pain"))
data$Symp_MYALGIA <- !is.na(str_extract(data.symp, "myal?gia"))
data$Symp_SPUT <- !is.na(str_extract(data.symp, "sput"))
data$Symp_TASTE <- !is.na(str_extract(data.symp, "loss of taste"))
data$Symp_VOM <- !is.na(str_extract(data.symp, regex("vomiti[mn]g", ignore_case=TRUE)))

# Convert "yes" and "no" for symptoms that have their own column to numbers
data$Num_Fever <- !is.na(str_extract(unlist(data[5]), regex("yes", ignore_case=TRUE)))
data$Num_Fatigue <- !is.na(str_extract(unlist(data[6]), regex("yes", ignore_case=TRUE)))
data$Num_Dyspnea <- !is.na(str_extract(unlist(data[7]), regex("yes", ignore_case=TRUE)))
data$Num_SoreThr <- !is.na(str_extract(unlist(data[8]), regex("yes", ignore_case=TRUE)))

# Identify non-numeric items that should be numbers.
for (col in c(
    "TLC (X 103"
  , "lympocytes absolute"
  , "lymph. %"
  , "serum ferritin before"
  , "serum ferritin one week after treament"
  , "HGB (gm/dl)"
  , "CRP before"
  , "CO- RADS"
  , "CRP at discharge"
)) {
  need.header <- TRUE

  vec <- unlist(data[col])[1:400]       # only include hospital cases
  notnum <- which(is.na(suppressWarnings(as.numeric(vec))))
  for (row in notnum) {
    value <- vec[row]
    if (is.na(value)) {
      next
    }

    if (need.header) {
      cat("Non-numeric items in column \"", col, "\"\n", sep="")
      cat("Data frame row\tExcel row\tValue\n")
      need.header <- FALSE
    }

    cat(row, "\t\t", xlrow(row), "\t\t", "\"", value, "\"", "\n", sep="")
  }

  if (! need.header) {
    cat("\n")
  }
}

# Clean up numeric items.
tlc <- unlist(data["TLC (X 103"])[1:400]
tlc <- gsub("o", "0", tlc)
tlc <- as.numeric(tlc)
data$Clean_TLC <- c(tlc, rep(NA, 200))

lym.a <- unlist(data["lympocytes absolute"])[1:400]
lym.a <- gsub("o", "0", lym.a)
lym.a <- as.numeric(lym.a)
lym.a[which(lym.a > 5)] <- lym.a[which(lym.a > 5)] / 10   # fix the biologically implausible numbers in cells L65, L132, L277, and L290
data$Clean_LYM.a <- c(lym.a, rep(NA, 200))

lym.p <- unlist(data["lymph. %"])[1:400]
lym.p <- gsub(",", ".", lym.p)
lym.p <- gsub("%", "", lym.p)
lym.p <- gsub("o", "0", lym.p)
lym.p <- as.numeric(lym.p)
lym.p[which(lym.p > 1)] <- lym.p[which(lym.p > 1)] / 100    # fix the numbers we converted from strings, plus cell S105 (1420%)
data$Clean_LYM.p <- c(lym.p, rep(NA, 200))

crp.t0 <- unlist(data["CRP before"])[1:400]
data$Clean_CRP.t0 <- c(crp.t0, rep(NA, 200))   # it seems no cleaning is needed, if we allow for the NAs.
crp.t1 <- unlist(data["CRP at discharge"])[1:400]
data$Clean_CRP.t1 <- c(crp.t1, rep(NA, 200))   # it seems no cleaning is needed, if we allow for the NAs.

sf.t0 <- unlist(data["serum ferritin before"])[1:400]
data$Clean_SF.t0 <- c(sf.t0, rep(NA, 200))  # I don't know what to do about cells N302-R302
sf.t1 <- unlist(data["serum ferritin one week after treament"])[1:400]
data$Clean_SF.t1 <- c(sf.t1, rep(NA, 200))

hgb <- unlist(data["HGB (gm/dl)"])[1:400]
hgb <- gsub(",", ".", hgb)
hgb <- gsub("%", "", hgb)
hgb <- gsub("o", "0", hgb)
hgb <- as.numeric(hgb)
# The HGB column contains numbers expressed as percentages.
# However, it seems that these should all be numbers in the general range 10-20.
# At this point we have some numbers in that range which we converted from strings, plus cell S105 (1420%).
# But the ones formatted as Excel percentages are numbers below 1.
# We convert those to more reasonable haemoglobin g/dl values in the next line.
hgb[which(hgb < 1)] <- hgb[which(hgb < 1)] * 100
data$Clean_HGB <- c(hgb, rep(NA, 200))

# Let's try to make sense of the date fields.
#
# It appears that the authors pasted strings in the format dd/mm/yyyy into Excel, but that
#  unfortunately Excel was in US date mode. So only dates that also matched mm/dd/yyyy
#  (i.e., those where dd and mm were both 12 or less) became "proper" Excel dates,
#  and these are now incorrect and need to have the date and month switched.
#
# The rest are still dd/mm/yyyy strings, some of which also need to be cleaned up.
#
# Possible exception: The "symptoms date&+ ve PCR" (entry to study) field in Excel cell W170
#  has the numeric date 2020-06-05, which is 21 days after the (string) value in cell X170,
#  corresponding to "recovery date & -ve PCR" (exiting study). However, based on the chronology
#  of column X, I think the most likely explanation is that the recovery data (the string "15/5/2020")
#  is an error and should read "15/06/2020", which would correspond to the reported hospital stay
#  of 11 days in cell Y170. That is, I do not think that W170 is an exception to my interpretation
#  of how the mixed date formats in columns W and X came about.

for (col in c("symptoms date&+ ve PCR", "recovery date & -ve PCR")) {
  vec <- unlist(data[col])

# Reduce the various strings to either an Excel day number (since 1/1/1900), or dd/mm.
  dv <- str_extract(unlist(vec), "[0-9][^ ]*[0-9]")   # get the first 1 or 2 numbers
  dv <- gsub("l", "/", dv)            # there is a case of "1l6l2020" in column 23, Excel row 110
  dv <- gsub("//", "/", dv)           # there are case of "dd/mm//2020" in column 23 and 24
  dv <- gsub("\\.", "/", dv)          # there is a case of "19.8/2020" in column 24, Excel row 290
  dv <- gsub("28/2020", "28/05", dv)  # the month is missing in "28/2020" in column 24, Excel row 167
  dv <- gsub("23/2020", "23/05", dv)  # the month is missing in "23/2020" in column 24, Excel row 168
  dv <- gsub("2020$", "", dv)         # remove year, whether it is preceded by / or not (e.g., "06/82020")
  dv <- gsub("/$", "", dv)            # remove trailing slashes after the previous operation
  dv <- gsub("/20", "", dv)           # deal with dates such as "7/8/20/20"
  dv <- gsub("31/6", "1/7", dv)       # there is a case of 31 June (!) in column 24, Excel row 155

  ddmm <- str_split_fixed(dv, "/", 2)
  datevec <- rep(NA, length(vec))

  for (i in 1:length(vec)) {
    pair <- ddmm[i, 1:2]
    if (pair[1] == "") {              # date is blank
      date <- NA
    }
    else if (pair[2] == "") {         # date is an Excel number
      usdate <- as.character(as.Date(as.integer(pair[1]), origin="1899-12-30"))
      date <- paste0("2020-", substr(usdate, 9, 10), "-", substr(usdate, 6, 7))
    }
    else {
      date <- as.character(as.Date(paste0("2020-", pair[2], "-", pair[1])))
    }

    datevec[i] <- date
  }

  new.col <- paste0("Clean_", col)
  data[new.col] <- datevec
}

# Calculate the difference between the dates of admission and discharge.
# We add 1, to include both endpoints; it seems the authors did this for groups 1 and 3 but not 2 and 4.
admitted <- as.Date(unlist(data["Clean_symptoms date&+ ve PCR"]))
discharged <- as.Date(unlist(data["Clean_recovery date & -ve PCR"]))
data$CalculatedStay <- discharged - admitted + 1

# Convert the authors' recorded text version of the hospital stay length to a number.
vec <- unlist(data["Hospital stay"])
vec <- gsub(" ", "", vec)
vec <- gsub("[A-Za-z]+", "", vec)
data$RecordedStay <- as.numeric(vec)

# Calculate the difference between the two versions of the stay length
data$StayDiscrepancy <- data$RecordedStay - data$CalculatedStay

cat("Comparing comorbidites across all 4 hosptialised groups\n")
comonames <- c("BA", "CHO", "CKD", "DM", "HBV", "HCV", "HTN", "IHD", "OHS", "PRG", "PE")
for (cn in comonames) {
  cat(cn)
  for (g in 1:4) {
    gr1 <- ((g - 1) * 100) + 1
    gr100 <- gr1 + 99
    grdata <- data[gr1:gr100,]
    col <- paste0("Como_", cn)
    cat("\t", sum(unlist(grdata[col])), sep="")
  }
  cat("\n")
}
cat("\n")

cat("Comparing symptoms across all 4 hosptialised groups\n")
sympnames <- c("ABDO", "ANOS", "BODY", "CHEST", "COUGH", "DIARR", "DROW", "DYSP", "EPIG", "HEAD", "JOINT", "LOIN", "MYALGIA", "SPUT", "TASTE", "VOM")
for (sn in sympnames) {
  cat(sn)
  for (g in 1:4) {
    grdata <- data[data$GroupNumber == g,]
    col <- paste0("Symp_", sn)
    cat("\t", sum(unlist(grdata[col])), sep="")
  }
  cat("\n")
}
cat("\n")

cat("Recalculating numbers that appeared in the text of the Results section\n")
res_como <- c("DM", "HTN", "IHD", "BA")
res_symp <- c("Fatigue", "Dyspnea")

cat("Group\tmAge\tsdAge\tMale\tFemale")
for (title in c(res_como, res_symp)) {
  cat("\t", title, sep="")
}
cat("\n")

for (g in 1:6) {
  grdata <- data[data$GroupNumber == g,]
  m.age <- mean(grdata$Age)
  sd.age <- sd(grdata$Age)

  cat(g
      , "\t", sprintf("%.1f", m.age)
      , "\t", sprintf("%.1f", sd.age)
      , "\t", sum(grdata$Sex == "M")
      , "\t", sum(grdata$Sex == "F")
      , sep="")

  for (como in res_como) {
    var <- paste0("Como_", como)
    n <- sum(unlist(grdata[var]), na.rm=TRUE)
    cat("\t", n, sep="")
  }

  for (symp in res_symp) {
    var <- paste0("Num_", symp)
    n <- if (g <= 4) sum(unlist(grdata[var]), na.rm=TRUE) else "n/a"
    cat("\t", n, sep="")
  }

  cat("\n")
}
cat("\n")

cat("Recalculating results from Table 1\n")
cat("\t\t\tGroup I\t\tGroup II\tGroup III\tGroup IV\tANOVA\n")
for (col in c(19, 11, 13, 20, 14)) {
  colname <- names(data)[col]
  cat(colname)
  if (nchar(colname) < 16) {
    cat("\t")
  }

  t1coldata <- suppressWarnings(as.numeric(unlist(data[col])))
  if ((col == 19) || (col == 13)) {
    t1coldata <- t1coldata * 100    # convert percentage (HGB probably shouldn't even be a percentage...)
  }

  for (g in 1:4) {
    gr1 <- ((g - 1) * 100) + 1
    gr100 <- gr1 + 99
    grdata <- t1coldata[gr1:gr100]
    gm <- mean(grdata, na.rm=TRUE)
    gs <- sd(grdata, na.rm=TRUE)
    ms <- sprintf("%.1f ± %.1f", gm, gs)
    cat("\t", ms, sep="")
  }

  cat("\t")
  dano <- data[1:400, c(colname, "pID", "GroupNumber")]    # data for ANOVA
  names(dano)[1] <- "DV"
  dano$DV <- suppressWarnings(as.numeric(unlist(dano$DV)))
  dano <- dano[!is.na(dano$DV),]
  dano$pID <- as.factor(dano$pID)
  dano$GroupNumber <- as.factor(dano$GroupNumber)
  eza <- suppressMessages(suppressWarnings(ezANOVA(data=dano, dv=DV, wid=pID, between=GroupNumber)))
  anova.result <- paste0("F(", eza$ANOVA$DFn, ",", eza$ANOVA$DFd, ")=", sprintf("%.1f", eza$ANOVA$F))
  cat(anova.result)
  cat("\n")

  cat("\t\tRange\t")

  for (g in 1:4) {
    if (g > 1) {
      cat("\t\t")
    }
    gr1 <- ((g - 1) * 100) + 1
    gr100 <- gr1 + 99
    grdata <- t1coldata[gr1:gr100]
    gmin <- floor(min(grdata, na.rm=TRUE))
    gmax <- ceiling(max(grdata, na.rm=TRUE))
    cat(gmin, "–", gmax, sep="")
  }

  cat("\n")
}
cat("\n")

for (tnum in 2:3) {
  cat("Recalculating results from Table ", tnum, "\n", sep="")
  cat("\t\t\t")
  cat(if (tnum==2) "Group I\t\tGroup II" else "Group III\tGroup IV")
  cat("\tt stat\tp value\n")
  for (colname in c("CRP at discharge", "serum ferritin one week after treament", "CalculatedStay", "RecordedStay")) {
    cat(sub(" one week after treament", "\t", colname))
    if (nchar(colname) < 16) {
      cat("\t")
    }

    t1coldata <- suppressWarnings(as.numeric(unlist(data[colname])))
    t1halves <- list()
    do.t.test <- TRUE

    g1 <- ((tnum - 2) * 2) + 1
    for (g in g1:(g1 + 1)) {
      gr1 <- ((g - 1) * 100) + 1
      gr100 <- gr1 + 99
      grdata <- t1coldata[gr1:gr100]

      if (any(!is.na(grdata))) {
        gm <- mean(grdata, na.rm=TRUE)
        gs <- sd(grdata, na.rm=TRUE)
        ms <- sprintf("%.1f ± %.1f", gm, gs)
      }
      else {
        ms <- "n/a ± n/a"
        do.t.test <- FALSE
      }
      cat("\t", ms, sep="")
      
      t1halves[[g]] <- grdata
    }

# If the authors used SPSS then they probably used Student's t test (SPSS default) rather than Welch's,
#  but I will leave the default for t.test() (i.e., Welch's test, equal variance not assumed) here.
    if (do.t.test) {
      ttest <- t.test(t1halves[[g1]], t1halves[[(g1 + 1)]])
      cat("\t", sprintf("%.2f", ttest$statistic), "\t", sprintf("%.2f", ttest$p.value), sep="")
    }

    cat("\n")
  }

  cat("\n")
}

cat("Recalculating results from Table 4\n")
cat("\t\t\tGroup I\t\tGroup II\tGroup III\tGroup IV\n")
for (colname in c("CalculatedStay", "RecordedStay")) {
  cat(colname, "\n\tRange\t\t", sep="")
  t1coldata <- unlist(data[colname])

  for (g in 1:4) {
    if (g > 1) {
      cat("\t\t")
    }
    gr1 <- ((g - 1) * 100) + 1
    gr100 <- gr1 + 99
    grdata <- t1coldata[gr1:gr100]
    gmin <- floor(min(grdata, na.rm=TRUE))
    gmax <- ceiling(max(grdata, na.rm=TRUE))
    cat(gmin, "–", gmax, sep="")
  }
  
  cat("\n\tMean/SD\t")
  
  for (g in 1:4) {
    gr1 <- ((g - 1) * 100) + 1
    gr100 <- gr1 + 99
    grdata <- t1coldata[gr1:gr100]
    gm <- mean(grdata, na.rm=TRUE)
    gs <- sd(grdata, na.rm=TRUE)
    ms <- sprintf("%.1f ± %.1f", gm, gs)
    cat("\t", ms, sep="")
  }
  
  cat("\n")
}
cat("\n")

start.date <- "2020-06-08"
cat("Identifying study entry/exit dates before reported start date (", start.date, ")\n", sep="")
cat("Group\tEntry\tExit\n")
for (g in 1:4) {
  grdata <- data[data$GroupNumber == g,]
  entry <- sum(grdata["Clean_symptoms date&+ ve PCR"] < start.date)
  exit <- sum(grdata["Clean_recovery date & -ve PCR"] < start.date)

  cat(g, "\t", entry, "\t", exit, "\n", sep="")
}
cat("\n")

cat("See plots for histogram of patient ages\n");
hist(data$Age, breaks=100)
cat("\n")

# Table of odd/even patient ages.
cat("Even (0) / odd (1) patient ages");     # omit trailing newline, as print() will supply it
print(table(data$Age %% 2))
cat("\n")

cat("Trailing digits of numeric values\n")
for (col in list(
    c("TLC (X 103", "Clean_TLC")
  , c("lympocytes absolute", "Clean_LYM.a")
  , c("lymph. %", "Clean_LYM.p")
  , c("serum ferritin before", "Clean_SF.t0")
  , c("serum ferritin one week after treament", "Clean_SF.t1")
  , c("HGB (gm/dl)", "Clean_HGB")
  , c("CRP before", "Clean_CRP.t0")
  , c("CRP at discharge", "Clean_CRP.t1")
  )) {
  colname <- col[1]
  cleancol <- col[2]
  vec <- as.vector(unlist(data[cleancol])[1:400])
  if (    (colname == "CRP at discharge")
       || (colname == "lympocytes absolute")
     ){
    vec <- vec * 10                       # convert decimal number to integer
  }
  else if (colname == "lymph. %") {
    vec <- vec * 100                      # convert decimal number to integer
  }
  tab <- table(round(vec) %% 10)
  csq <- chisq.test(tab)
  cat(colname)
  print(tab)                  # print() is a bit clunky but at least it works when some numbers are absent
  cat("Chi-sq(9)=", round(csq$statistic, 2), " p=", csq$p.value, "\n", sep="")
  cat("\n")
}