Reanalysis of Laird et al. (2020) https://pubmed.ncbi.nlm.nih.gov/32603576/

Laird-VitD.R

# Reanalysis of Laird et al. (2020) https://pubmed.ncbi.nlm.nih.gov/32603576/
# By Nick Brown, January 2021.
# Licence: CC-0.

# Laird et al. did not state how they went from four (Italy) or three (Spain) measurements
#  of Vitamin D levels to a single figure.
# When I calculate the unweighted mean of their samples, this seems to correspond to the
#  Y-axis values in their figure 1, but I have included the relevant numbers from their
#  Table 1 here to allow a weighted average to be calculated.

library(ggplot2)

rawVitD <- data.frame(
  "Country"=c("Italy", "Italy", "Italy", "Italy", "Spain", "Spain", "Spain")
  , "N"=c(104, 1005, 570, 700, 171, 237, 256)
  , "VitaminD"=c(15.0, 48.5, 45.7, 27.2, 33.5, 43.0, 45.7)
)

vitD <- data.frame(
  "Country"=c(
    "Italy", "Spain", "Sweden", "Norway", "Finland"
  , "United Kingdom", "Ireland", "Scotland"
  , "Germany", "France", "Portugal", "Netherlands"
  )
, "VitaminD"=c(
    NA, NA, 73.5, 67.7, 65.7
  , 48.7, 51.3, 28.0
  , 43.2, 35.9, 42.3, 51.2
  )
, "DeathsApr20"=c(
    281.0, 300.0, 59.1, 13.2, 6.2
  , 93.0, 43.7, 65.5
  , 22.0, 156.0, 32.8, 124.0
  )
# January 2021 death rates calculated from Johns Hopkins data, 29 January 2021, except for
#  UK and Scotland, which are taken from the UK ONS web site and gov.uk estimates of population.
, "DeathsJan21"=c(
    1445.23, 1236.37, 1137.31, 102.74, 119.84
  , 1519.86, 641.38, 1092.75
  , 666.99, 1142.90, 1138.41, 812.67
  )
)

# Merge the multiple studies from Italy and Spain into a single record.
# You can change the formula here if you want to use, say, a weighted average.
for (cc in unique(rawVitD$Country)) {
  vitD[(vitD$Country == cc),]$VitaminD <- mean(rawVitD[(rawVitD$Country == cc),]$VitaminD)
}

ymax <- (floor(max(vitD$VitaminD) / 20) + 1) * 20

xmax.apr20 <- (floor(max(vitD$DeathsApr20) / 50) + 1) * 50

p.apr20 <- cor.test(vitD$VitaminD, vitD$DeathsApr20, method="spearman")$p.value
p.apr20.label <- paste0("p=", sprintf("%.3f", p.apr20))
p.apr20.x <- xmax.apr20 * 0.9
p.y <- ymax * 0.9

ggplot(vitD, aes(x=DeathsApr20, y=VitaminD)) +
  geom_point(colour="red") +
  geom_text(label=vitD$Country, size=3, position=position_nudge(x=5, y=-1.5)) +
  geom_smooth(method="lm", formula=y~x, se=FALSE) +
  annotate(geom="label", label=p.apr20.label, x=p.apr20.x, y=p.y, fill="white") +
  xlim(c(0, xmax.apr20)) +
  ylim(c(0, ymax)) +
  labs(x="Deaths/100k April 2020")

xmax.jan21 <- (floor(max(vitD$DeathsJan21) / 100) + 1) * 100
p.jan21 <- cor.test(vitD$VitaminD, vitD$DeathsJan21, method="spearman")$p.value
p.jan21.label <- paste0("p=", sprintf("%.3f", p.jan21))
p.jan21.x <- xmax.jan21 * 0.9

ggplot(vitD, aes(x=DeathsJan21, y=VitaminD)) +
  geom_point(colour="red") +
  geom_text(label=vitD$Country, size=3, position=position_nudge(x=5, y=-1.5)) +
  geom_smooth(method="lm", formula=y~x, se=FALSE) +
  annotate(geom="label", label=p.jan21.label, x=p.jan21.x, y=p.y, fill="white") +
  xlim(c(0, xmax.jan21)) +
  ylim(c(0, ymax)) +
  labs(x="Deaths/100k January 2021")