JH Aug 28, 2020

https://www.inspq.qc.ca/covid-19/donnees

https://experience.arcgis.com/experience/a6f23959a8b14bfa989e3cda29297ded

https://covid19.who.int/explorer

https://gisanddata.maps.arcgis.com/apps/opsdashboard/index.html#/bda7594740fd40299423467b48e9ecf6

1a. Immunogenicity of the ChAdOx1 nCoV-19 vaccine against SARS-CoV-2

1 What LEVELS of immunity are found in patients who have recovered from COVID-19? (PANEL B)
.
2 Relative to these, at 28 days say, what LEVELS of immunity are found in persons who have received the ChAdOx1 nCoV-19 vaccine? CONTRAST Panel A (prime, 28 days) vs PANEL B.
.
3 What is the TIME COURSE of immunity LEVELS in persons who have received the ChAdOx1 nCoV-19 vaccine?: Panel A (prime) 0, 7, 14, 28, 56 days
.
http://www.biostat.mcgill.ca/hanley/statbook/CovisVaccineOxford.pdf

1b. Aerosol and Surface Stability of SARS-CoV-2 as Compared with SARS-CoV-1

http://www.biostat.mcgill.ca/hanley/statbook/AerosolStabilitySARS2vsSARS1.pdf

Viability of SARS-CoV-1 and SARS-CoV-2 in Aerosols and on Various Surfaces. As shown in Panel A, the titer of aerosolized viable virus is expressed in 50% tissue-culture infectious dose (TCID50) per liter of air. Viruses were applied to copper, cardboard, stainless steel, and plastic main- tained at 21 to 23°C and 40% relative humidity over 7 days. The titer of viable virus is expressed as TCID50 per milliliter of collection medium. All samples were quantified by end-point titration on Vero E6 cells. Plots show the means and standard errors (I bars) across three replicates. As shown in Panel B, regression plots indicate the predicted decay of virus titer over time; the titer is plotted on a logarithmic scale. Points show measured titers and are slightly jittered (i.e., they show small rapid variations in the amplitude or timing of a waveform arising from fluctuations) along the time axis to avoid overplotting. Lines are random draws from the joint posterior distribu- tion of the exponential decay rate (negative of the slope) and intercept (initial virus titer) to show the range of possible decay patterns for each experimental condition. There were 150 lines per panel, including 50 lines from each plotted replicate. As shown in Panel C, violin plots indicate posterior distribution for the half-life of viable virus based on the estimated ex- ponential decay rates of the virus titer. The dots indicate the posterior median estimates, and the black lines indicate a 95% credible interval. Experimental conditions are ordered according to the posterior median half-life of SARS-CoV-2. The dashed lines indicate the limit of detection, which was 3.33×100.5 TCID50 per liter of air for aerosols, 100.5 TCID50 per milliliter of medium for plastic, steel, and cardboard, and 101.5 TCID50 per milliliter of medium for copper.

2a. Case Fatality Rates

Full report: https://www.publichealthontario.ca/-/media/documents/ncov/epi/2020/06/covid19-epi-case-identification-age-only-template.pdf

1 Case Fatality Rate, Ontario overall [grey]
.
3 Gradient: Case Fatality Rate, Age specific Case fatality [shown]

2b. Cases within 1 mo. of Stay-at-Home Order [Illinois] or not [Iowa]

.
http://www.biostat.mcgill.ca/hanley/statbook/IowaIllinoisBorder.pdf

3. (Excess) deaths

http://www.biostat.mcgill.ca/hanley/statbook/AnalysisExcessAllCauseMortalityIrelandDuringCOVID19epidemic.pdf

See also https://www.rip.ie and http://www.medicine.mcgill.ca/epidemiology/hanley/Pandemics/ [bottom of page]

https://www.medrxiv.org/content/10.1101/2020.07.22.20159913v2

Statistical Models

1a. LEVELS of immunity: patients who recovered from COVID-19

Data (imperfectly!) scraped from Postscript file ‘behind’ the pdf file.

ds=read.table("http://www.biostat.mcgill.ca/hanley/statbook/immunogenicityChAdOx1.nCoV-19vaccine.txt")
str(ds)
## 'data.frame':    307 obs. of  2 variables:
##  $ RefIndexCategory            : chr  "Convalescent" "Convalescent" "Convalescent" "Convalescent" ...
##  $ IgGResponse.log10.ElisaUnits: num  2.56 2.74 2.79 3.32 3.15 2.35 2.72 2.95 2.42 2.64 ...
tail(ds)
##             RefIndexCategory IgGResponse.log10.ElisaUnits
## 302 Day28PostChAdOx1 nCoV-19                         1.99
## 303 Day28PostChAdOx1 nCoV-19                         1.99
## 304 Day28PostChAdOx1 nCoV-19                         2.42
## 305 Day28PostChAdOx1 nCoV-19                         2.46
## 306 Day28PostChAdOx1 nCoV-19                         2.42
## 307 Day28PostChAdOx1 nCoV-19                         1.17
# just the naturally induced responses

natural = ds[ds$RefIndexCategory=="Convalescent",]
str(natural)
## 'data.frame':    180 obs. of  2 variables:
##  $ RefIndexCategory            : chr  "Convalescent" "Convalescent" "Convalescent" "Convalescent" ...
##  $ IgGResponse.log10.ElisaUnits: num  2.56 2.74 2.79 3.32 3.15 2.35 2.72 2.95 2.42 2.64 ...
hist(natural$IgGResponse.log10.ElisaUnits,breaks=20)

summary(natural$IgGResponse.log10.ElisaUnits,breaks=20)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.417   2.570   2.577   2.780   3.860
summary( glm(IgGResponse.log10.ElisaUnits ~ 1,
             data = natural) )
## 
## Call:
## glm(formula = IgGResponse.log10.ElisaUnits ~ 1, data = natural)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.57733  -0.15983  -0.00733   0.20267   1.28267  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.57733    0.03432   75.09   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.2120565)
## 
##     Null deviance: 37.958  on 179  degrees of freedom
## Residual deviance: 37.958  on 179  degrees of freedom
## AIC: 234.65
## 
## Number of Fisher Scoring iterations: 2

COMPARISON of naturally vs. vaccine-induced response LEVELS

# in pictures
boxplot(ds$IgGResponse.log10.ElisaUnits ~ ds$RefIndexCategory)

# in numbers
by(ds$IgGResponse.log10.ElisaUnits,ds$RefIndexCategory,summary)
## ds$RefIndexCategory: Convalescent
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.417   2.570   2.577   2.780   3.860 
## ------------------------------------------------------------ 
## ds$RefIndexCategory: Day28PostChAdOx1 nCoV-19
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.170   1.985   2.050   2.047   2.120   2.850
# via statistical model
ds$vaccine.induced = ds$RefIndexCategory =="Day28PostChAdOx1 nCoV-19"

fitted.model = glm(IgGResponse.log10.ElisaUnits ~
               ds$vaccine.induced, data = ds) 
summary(fitted.model)
## 
## Call:
## glm(formula = IgGResponse.log10.ElisaUnits ~ ds$vaccine.induced, 
##     data = ds)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.57733  -0.14654   0.00346   0.14807   1.28267  
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             2.57733    0.02874   89.67   <2e-16 ***
## ds$vaccine.inducedTRUE -0.53080    0.04469  -11.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1487187)
## 
##     Null deviance: 66.339  on 306  degrees of freedom
## Residual deviance: 45.359  on 305  degrees of freedom
## AIC: 290.17
## 
## Number of Fisher Scoring iterations: 2
plot(ds$vaccine.induced,ds$IgGResponse.log10.ElisaUnits,pch=19,cex=0.5)
abline(h=seq(0,4,0.5),col="lightblue")
lines(ds$vaccine.induced,fitted.model$fitted.values,col="red")

1b. Aerosol Stability (SARS-CoV-2)

Item=rep("Aerosols",15)
Hour=rep(c(0,0.5,1,2,3),each=3)
CoVType=rep(2.15)
Titer.log10 = c(3,  3.45,3.55, 
                3.5,3,  2.8,
                3.45,3.55,2.8,
                3,  2.8,2.3,
                3,  2.5,2.35)

plot(Hour,Titer.log10,pch=19,cex=1,
     ylim=c(0,4),cex.lab=1.5)
abline(h=0:4,col="lightblue")

fitted.line = glm(Titer.log10 ~ Hour)

summary(fitted.line)
## 
## Call:
## glm(formula = Titer.log10 ~ Hour)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -0.52773  -0.25273  -0.02773   0.25826   0.47141  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.32945    0.13504  24.656 2.67e-12 ***
## Hour        -0.25086    0.07999  -3.136  0.00788 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.1113323)
## 
##     Null deviance: 2.5423  on 14  degrees of freedom
## Residual deviance: 1.4473  on 13  degrees of freedom
## AIC: 13.493
## 
## Number of Fisher Scoring iterations: 2
lines(Hour,fitted.line$fitted.values,col="red")

fitted.slope = fitted.line$coefficients[2]
fitted.half.life = log10(0.5)/fitted.slope

points(fitted.half.life,
    fitted.line$coefficients[1]+log10(0.5),pch=18,
    col="red",cex=4)
abline(v=fitted.half.life,col="red",lty="dotted")

2b. Comparison of Estimated Rates of Coronavirus Disease 2019 (COVID-19) in Border Counties in Iowa Without a Stay-at-Home Order and Border Counties in Illinois With a Stay-at-Home Order.

State=c("Iowa","Illinois")
Lockdown = c(0,1)
Population  =c(462445,272385)

Rate.per10000.to.April.20 = c(15.5,  10)
Ave.Daily.Rate.per10000.to.March.21 = c( 0.024, 0.026)

( Cases.to.April.20 = round(Rate.per10000.to.April.20 * 
                              Population/10000) )
## [1] 717 272
( Cases.to.March.21 = round(
  Ave.Daily.Rate.per10000.to.March.21 * 6 * Population/10000) )
## [1] 7 4
( Cases.March.22.to.April.20 = Cases.to.April.20  - Cases.to.March.21 )
## [1] 710 268
( ds = data.frame(State,Lockdown,Population,Cases.March.22.to.April.20) )
##      State Lockdown Population Cases.March.22.to.April.20
## 1     Iowa        0     462445                        710
## 2 Illinois        1     272385                        268
ds
##      State Lockdown Population Cases.March.22.to.April.20
## 1     Iowa        0     462445                        710
## 2 Illinois        1     272385                        268
sum(ds$Cases.March.22.to.April.20)/sum(ds$Pop)
## [1] 0.00133092
summary(glm(Cases.March.22.to.April.20 ~ Lockdown)) 
## 
## Call:
## glm(formula = Cases.March.22.to.April.20 ~ Lockdown)
## 
## Deviance Residuals: 
## [1]  0  0
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)      710         NA      NA       NA
## Lockdown        -442         NA      NA       NA
## 
## (Dispersion parameter for gaussian family taken to be NaN)
## 
##     Null deviance: 9.7682e+04  on 1  degrees of freedom
## Residual deviance: 5.4930e-26  on 0  degrees of freedom
## AIC: -106.04
## 
## Number of Fisher Scoring iterations: 1
summary(glm(Cases.March.22.to.April.20 ~ -1 + Population),
        data = ds)
## 
## Call:
## glm(formula = Cases.March.22.to.April.20 ~ -1 + Population)
## 
## Deviance Residuals: 
##       1        2  
##   65.68  -111.51  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## Population 0.0013933  0.0002411   5.778    0.109
## 
## (Dispersion parameter for gaussian family taken to be 16748.63)
## 
##     Null deviance: 575924  on 2  degrees of freedom
## Residual deviance:  16749  on 1  degrees of freedom
## AIC: 27.742
## 
## Number of Fisher Scoring iterations: 2
ds$Product = ds$Population * ds$Lockdown

fitted.model = glm(Cases.March.22.to.April.20 ~ -1 + Population + Product,
            data = ds)

summary(fitted.model)
## 
## Call:
## glm(formula = Cases.March.22.to.April.20 ~ -1 + Population + 
##     Product, data = ds)
## 
## Deviance Residuals: 
## [1]  0  0
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## Population  0.0015353         NA      NA       NA
## Product    -0.0005514         NA      NA       NA
## 
## (Dispersion parameter for gaussian family taken to be NaN)
## 
##     Null deviance: 575924  on 2  degrees of freedom
## Residual deviance:      0  on 0  degrees of freedom
## AIC: -Inf
## 
## Number of Fisher Scoring iterations: 1
plot(ds$Population,ds$Cases.March.22.to.April.20,
     ylim = c(0,max(ds$Cases.March.22.to.April.20)),
     xlim = c(0,1.1*max(ds$Population)))
segments( 0,0, ds$Population,ds$Cases.March.22.to.April.20)
text(ds$Population,
     ds$Cases.March.22.to.April.20,
     ds$State,adj=c(-0.25,0.5))

slope = toString(
  round(fitted.model$coefficients[1],5) )
text(ds$Population[1]/1.95,
     ds$Cases.March.22.to.April.20[1]/2,
     paste("slope (cases/capita):", slope ),
     adj=c(0,0.5))

slope = toString(
  round(fitted.model$coefficients[2],5) )
text(ds$Population[2]/1.95,
     ds$Cases.March.22.to.April.20[2]/2,
     paste("difference in.slope (in cases/capita):", slope ),
     adj=c(0,0.5))

fitted.poisson.model = 
  glm(Cases.March.22.to.April.20 ~
             -1 + Population + Product, 
            family =  poisson(link="identity"),
            data = ds)

summary(fitted.poisson.model)
## 
## Call:
## glm(formula = Cases.March.22.to.April.20 ~ -1 + Population + 
##     Product, family = poisson(link = "identity"), data = ds)
## 
## Deviance Residuals: 
## [1]  0  0
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## Population  1.535e-03  5.762e-05  26.646  < 2e-16 ***
## Product    -5.514e-04  8.326e-05  -6.623 3.52e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:         Inf  on 2  degrees of freedom
## Residual deviance: -4.0856e-14  on 0  degrees of freedom
## AIC: 19.833
## 
## Number of Fisher Scoring iterations: 2
Stay.at.Home = c(0,1)
Fitted.rates = c(fitted.poisson.model$coefficients[1],
       fitted.poisson.model$coefficients[1]+
       fitted.poisson.model$coefficients[2])
plot(Stay.at.Home,Fitted.rates,
     ylab="Fitted Rates (per capita)",
     cex.lab=1.5,
     ylim=c(0,fitted.poisson.model$coefficients[1]),
     cex=3,pch=19,col="red")
abline(h=(0:3)*5/10000,col="lightblue")

lines(Stay.at.Home,Fitted.rates,col="red")

When (as here) proportions (risks) are low, and sample sizes large, and everyone is observed for the same (here 30 day) period, the binomial and Poisson models are just about interchangeable.

BUT, since the cases probably cluster by workplace/household, the binomial/Poisson model is probably ‘too tight’

In the Ontario study, the case-fatality rates (risks) in older persons are quite high; if all cases in an age band were followed to death or recovery, the binomial might be a reasonable first choice of model. But since some of the more recently diagnosed cases were not followed until the ultimate 0/1 outcome, ie. their observations were ‘censored,’ am (actuarial) correction (effectively a reduction in the ‘effective sample size’, was applied.