library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union
library(ggplot2)
library(ggpubr)
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
library(tibble)
library(tidyr)
pops<-read.table("../data/french_region_pop.csv", sep=",")
d<- read.table("../results/mixture_Rt_foldchanges.csv", sep=",", h=T)
pops
head(d)
summary(d)
                     region    region_id  Rt_foldchange          R0       
 auvergne.rhone.alpes   :1   Min.   : 1   Min.   :0.1742   Min.   :2.030  
 bourgogne.franche.comte:1   1st Qu.: 4   1st Qu.:0.1925   1st Qu.:2.989  
 bretagne               :1   Median : 7   Median :0.1977   Median :3.338  
 centre.val.de.loire    :1   Mean   : 7   Mean   :0.1957   Mean   :3.284  
 corse                  :1   3rd Qu.:10   3rd Qu.:0.2008   3rd Qu.:3.555  
 grand.est              :1   Max.   :13   Max.   :0.2068   Max.   :3.971  
 (Other)                :7                                                

Adding the area to the populations

Source: wikipedia : https://fr.wikipedia.org/wiki/R%C3%A9gion_fran%C3%A7aise

pops$area <- c(69711, 47784, 27208, 39151, 8680, 57441, 31806, 12011, 29907, 84036, 72724, 32082, 31400)
colnames(pops)<-c("region", "census", "area")
pops$density <- pops$census / pops$area
pops

Pre-post lockdown population

prepost <- read.table("../data/CP_Donnees_Population-confinement/CP_Tableaux_Population-confinement.csv", h=T, sep=";", skip=0)[,1:5]
prepost$Code.département <- as.numeric(as.character(prepost$Code.département))
NAs introduced by coercion
depToReg <- read.table("../data/departements-region.csv", h=T, sep=",")
depToReg$num_dep <- as.numeric(as.character(depToReg$num_dep))
NAs introduced by coercion
prePostReg <-merge(prepost, depToReg, by.x="Code.département", by.y="num_dep")
ppReg <- group_by(prePostReg, region_name) %>% 
  summarise(diff = sum(Ecart.moyen.en.nuitées.après...avant.le.confinement))

plots and correlations

plotAndCor <- function(x, y, xlab, ylab) {
  res<-cor.test(x,y)
  plot(x, y, pch=20, ylab=ylab, xlab=xlab, main=paste("cor=",res$estimate, "; pvalue=", res$p.value, sep=""))
}
plotAndCor(pops$census, d$Rt_foldchange, "Population size", "Median R_t fold change")

plotAndCor(pops$density, d$Rt_foldchange, "Population density", "Median R_t fold change")

plotAndCor(pops$area, d$Rt_foldchange, "Area", "Median R_t fold change")

plotAndCor(d$R0, d$Rt_foldchange, "Median R0", "Median R_t fold change")

plotAndCor(pops$census, d$R0, "Population size", "Median R_0")

plotAndCor(pops$density, d$R0, "Population density", "Median R_0")

plotAndCor(ppReg$diff, d$Rt_foldchange, "Pre-post lockdown difference", "Median R_t fold change")

Linear model

full_data <- cbind(d, ppReg, pops)
mod <- lm(Rt_foldchange~R0+diff+census+density, data=full_data)
summary(mod)

Call:
lm(formula = Rt_foldchange ~ R0 + diff + census + density, data = full_data)

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0106742 -0.0015679 -0.0009788  0.0025482  0.0121430 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.883e-01  1.210e-02  15.568 2.89e-07 ***
R0           2.854e-03  4.187e-03   0.682   0.5147    
diff        -1.117e-08  2.349e-08  -0.476   0.6471    
census       9.765e-10  1.206e-09   0.810   0.4413    
density     -3.944e-05  1.368e-05  -2.884   0.0204 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.00624 on 8 degrees of freedom
Multiple R-squared:  0.6302,    Adjusted R-squared:  0.4453 
F-statistic: 3.409 on 4 and 8 DF,  p-value: 0.06583
LS0tCnRpdGxlOiAiTWl4dHVyZSBhbW9uZyByZWdpb25zIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7cn0KbGlicmFyeShkcGx5cikKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGdncHVicikKbGlicmFyeSh0aWJibGUpCmxpYnJhcnkodGlkeXIpCmBgYAoKCmBgYHtyfQpwb3BzPC1yZWFkLnRhYmxlKCIuLi9kYXRhL2ZyZW5jaF9yZWdpb25fcG9wLmNzdiIsIHNlcD0iLCIpCmQ8LSByZWFkLnRhYmxlKCIuLi9yZXN1bHRzL21peHR1cmVfUnRfZm9sZGNoYW5nZXMuY3N2Iiwgc2VwPSIsIiwgaD1UKQpgYGAKCgpgYGB7cn0KcG9wcwpgYGAKYGBge3J9CmhlYWQoZCkKYGBgCgpgYGB7cn0Kc3VtbWFyeShkKQpgYGAKCiMgQWRkaW5nIHRoZSBhcmVhIHRvIHRoZSBwb3B1bGF0aW9ucwpTb3VyY2U6IHdpa2lwZWRpYSA6IGh0dHBzOi8vZnIud2lraXBlZGlhLm9yZy93aWtpL1IlQzMlQTlnaW9uX2ZyYW4lQzMlQTdhaXNlCmBgYHtyfQpwb3BzJGFyZWEgPC0gYyg2OTcxMSwgNDc3ODQsIDI3MjA4LCAzOTE1MSwgODY4MCwgNTc0NDEsIDMxODA2LCAxMjAxMSwgMjk5MDcsIDg0MDM2LCA3MjcyNCwgMzIwODIsIDMxNDAwKQpjb2xuYW1lcyhwb3BzKTwtYygicmVnaW9uIiwgImNlbnN1cyIsICJhcmVhIikKcG9wcyRkZW5zaXR5IDwtIHBvcHMkY2Vuc3VzIC8gcG9wcyRhcmVhCmBgYAoKYGBge3J9CnBvcHMKYGBgCgoKIyBQcmUtcG9zdCBsb2NrZG93biBwb3B1bGF0aW9uCmBgYHtyfQpwcmVwb3N0IDwtIHJlYWQudGFibGUoIi4uL2RhdGEvQ1BfRG9ubmVlc19Qb3B1bGF0aW9uLWNvbmZpbmVtZW50L0NQX1RhYmxlYXV4X1BvcHVsYXRpb24tY29uZmluZW1lbnQuY3N2IiwgaD1ULCBzZXA9IjsiLCBza2lwPTApWywxOjVdCnByZXBvc3QkQ29kZS5kw6lwYXJ0ZW1lbnQgPC0gYXMubnVtZXJpYyhhcy5jaGFyYWN0ZXIocHJlcG9zdCRDb2RlLmTDqXBhcnRlbWVudCkpCmBgYAoKYGBge3J9CmRlcFRvUmVnIDwtIHJlYWQudGFibGUoIi4uL2RhdGEvZGVwYXJ0ZW1lbnRzLXJlZ2lvbi5jc3YiLCBoPVQsIHNlcD0iLCIpCmRlcFRvUmVnJG51bV9kZXAgPC0gYXMubnVtZXJpYyhhcy5jaGFyYWN0ZXIoZGVwVG9SZWckbnVtX2RlcCkpCmBgYAoKYGBge3J9CnByZVBvc3RSZWcgPC1tZXJnZShwcmVwb3N0LCBkZXBUb1JlZywgYnkueD0iQ29kZS5kw6lwYXJ0ZW1lbnQiLCBieS55PSJudW1fZGVwIikKCmBgYAoKYGBge3J9CnBwUmVnIDwtIGdyb3VwX2J5KHByZVBvc3RSZWcsIHJlZ2lvbl9uYW1lKSAlPiUgCiAgc3VtbWFyaXNlKGRpZmYgPSBzdW0oRWNhcnQubW95ZW4uZW4ubnVpdMOpZXMuYXByw6hzLi4uYXZhbnQubGUuY29uZmluZW1lbnQpKQoKYGBgCgoKIyBwbG90cyBhbmQgY29ycmVsYXRpb25zCmBgYHtyfQpwbG90QW5kQ29yIDwtIGZ1bmN0aW9uKHgsIHksIHhsYWIsIHlsYWIpIHsKICByZXM8LWNvci50ZXN0KHgseSkKICBwbG90KHgsIHksIHBjaD0yMCwgeWxhYj15bGFiLCB4bGFiPXhsYWIsIG1haW49cGFzdGUoImNvcj0iLHJlcyRlc3RpbWF0ZSwgIjsgcHZhbHVlPSIsIHJlcyRwLnZhbHVlLCBzZXA9IiIpKQp9CmBgYAoKYGBge3J9CnBsb3RBbmRDb3IocG9wcyRjZW5zdXMsIGQkUnRfZm9sZGNoYW5nZSwgIlBvcHVsYXRpb24gc2l6ZSIsICJNZWRpYW4gUl90IGZvbGQgY2hhbmdlIikKYGBgCgoKYGBge3J9CnBsb3RBbmRDb3IocG9wcyRkZW5zaXR5LCBkJFJ0X2ZvbGRjaGFuZ2UsICJQb3B1bGF0aW9uIGRlbnNpdHkiLCAiTWVkaWFuIFJfdCBmb2xkIGNoYW5nZSIpCgpgYGAKCmBgYHtyfQpwbG90QW5kQ29yKHBvcHMkYXJlYSwgZCRSdF9mb2xkY2hhbmdlLCAiQXJlYSIsICJNZWRpYW4gUl90IGZvbGQgY2hhbmdlIikKCmBgYAoKCgpgYGB7cn0KcGxvdEFuZENvcihkJFIwLCBkJFJ0X2ZvbGRjaGFuZ2UsICJNZWRpYW4gUjAiLCAiTWVkaWFuIFJfdCBmb2xkIGNoYW5nZSIpCmBgYAoKYGBge3J9CnBsb3RBbmRDb3IocG9wcyRjZW5zdXMsIGQkUjAsICJQb3B1bGF0aW9uIHNpemUiLCAiTWVkaWFuIFJfMCIpCmBgYAoKYGBge3J9CnBsb3RBbmRDb3IocG9wcyRkZW5zaXR5LCBkJFIwLCAiUG9wdWxhdGlvbiBkZW5zaXR5IiwgIk1lZGlhbiBSXzAiKQpgYGAKCmBgYHtyfQpwbG90QW5kQ29yKHBwUmVnJGRpZmYsIGQkUnRfZm9sZGNoYW5nZSwgIlByZS1wb3N0IGxvY2tkb3duIGRpZmZlcmVuY2UiLCAiTWVkaWFuIFJfdCBmb2xkIGNoYW5nZSIpCmBgYAoKIyBMaW5lYXIgbW9kZWwKYGBge3J9CmZ1bGxfZGF0YSA8LSBjYmluZChkLCBwcFJlZywgcG9wcykKYGBgCgpgYGB7cn0KbW9kIDwtIGxtKFJ0X2ZvbGRjaGFuZ2V+UjArZGlmZitjZW5zdXMrZGVuc2l0eSwgZGF0YT1mdWxsX2RhdGEpCmBgYAoKYGBge3J9CnN1bW1hcnkobW9kKQpgYGAKCg==