# R0 and the exponential growth of a pandemic, an update

A few days ago, I wrote a blog post – R0 and the exponential growth of a pandemic – where I was trying to generate some visualization of some exponential growth, in the context of a pandemic. After giving some thoughts, the previous graph might not be the best one to see an exponential based contagion.

Having graphs evolving, from the left to the right, gives us the (false) idea of some temporal evolution. Which is no necessarily correct. It simply means that contaminated people will contaminate other people, and we look at the number of iterations here. So maybe some concentric dots would look better.

And from a technical perspective, what I did was fun, but probably too complicated. In my previous post, I wanted to pack optimally $k$ identical disks intro a unit circle. On http://hydra.nat.uni-magdeburg.de/packing, it was possible to get the “best known packings of equal circles in a circle”, with the coordinates. But as we will see, we can use something much more simple here.

My idea is now to create some picture like one below, with concentric colored dot. In the center, we have the first people that were contaminated, and then, we can see the transmission, somehow

From a technical perspective, here, I use a different strategy. I decided to draw random points, uniformly. The problem with randomness is the natural high discrepancy, with monte carlo methods: it is very likely that some disks will overlap. It is not a major issue, but it might distort the message. So I decided to use some low-discrepancy sequences, such as Halton‘s sequence.

library(randtoolbox) S = halton(n=5000, dim = 2)*2-1

Here, I have disk coordinates in $[-1,+1]^2$. Then, to get disks in a circle, I simply compute the distance to the origin $(0,0)$,

D0 = S[,1]^2+S[,2]^2

and take the ranks. If I want to visualize $k=200$ people, I consider the 200 smaller ranks. To get concentric circles, each part having $k_i$ individuals, I use as thresholds $R_0^{\bar k_{i-1}},R_0^{\bar k_{i}},R_0^{\bar k_{i+1}}$, etc, where $\bar k_i=\bar k_{i-1}+k_i$,

R0 = rank(D0,ties.method = "random") C0 = as.numeric(cut(R0,c(0,cumsum(k)+.5)),100000)

where

R0=1.8 k=round(R0^(seq(1,9,by=2)))

Then we can plot the dots, with appropriate colors,

points(S,pch=19,col=colrpal[C0],cex=.75)

And of course, we can try that with different values, for $R_0$

R0=2.2 k=round(R0^(seq(1,9,by=2))) kmax=max(k) S = halton(n=5000, dim = 2)*2-1 plot(S,col="light yellow",axes=FALSE,xlab="",ylab="",xlim=c(-1.3,1),ylim=c(-1,1),cex=.75,pch=19) D0 = S[,1]^2+S[,2]^2 R0 = rank(D0,ties.method = "random") C0 = as.numeric(cut(R0,c(0,cumsum(k)+.5)),100000) points(S,pch=19,col=colrpal[C0],cex=.75)

# Modeling Joint Lives within Families

With Olivier Cabrignac and Ewen Gallic, we recently uploaded a research paper, entitled “Modeling Joint Lives within Families

Family history is usually seen as a significant factor insurance companies look at when applying for a life insurance policy. Where it is used, family history of cardiovascular diseases, death by cancer, or family history of high blood pressure and diabetes could result in higher premiums or no coverage at all. In this article, we use massive (historical) data to study dependencies between life length within families. If joint life contracts (between a husband and a wife) have been long studied in actuarial literature, little is known about child and parents dependencies. We illustrate those dependencies using 19th century family trees in France, and quantify implications in annuities computations. For parents and children, we observe a modest but significant positive association between life lengths. It yields different estimates for remaining life expectancy, present values of annuities, or whole life insurance guarantee, given information about the parents (such as the number of parents alive). A similar but weaker pattern is observed when using information on grandparents.

The paper is online on https://arxiv.org/abs/2006.08446.

# Slides 7 – modèle de Poisson

On va maintenant parler des modèles de comptage, en revenant sur le modèle de Poisson. Les slides sont en ligne (slides 7) et la vidéo aussi (slides 7)

Again, this post is related to my MAT7381 course, where we will see that it is actually possible to write our own code to compute Lasso regression, $$\min\left\lbrace\frac{1}{2}\|\mathbf{y}-\mathbf{X}\mathbf{\beta}\|_{\ell_2}^2+\lambda\|\mathbf{\beta}\|_{\ell_1}\right\rbrace$$We have to define the soft-thresholding function$$S(z,\gamma)=\text{sign}(z)\cdot(|z|-\gamma)_+=\begin{cases}z-\gamma&\text{ if }\gamma>|z|\text{ and }z<0\\z+\gamma&\text{ if }\gamma<|z|\text{ and }z<0 \\0&\text{ if }\gamma\geq|z|\end{cases}$$The R function would be

soft_thresholding = function(x,a){ sign(x) * pmax(abs(x)-a,0) }

To solve our optimization problem, set$$\mathbf{r}_j=\mathbf{y} - \left(\beta_0\mathbf{1}+\sum_{k\neq j}\beta_k\mathbf{x}_k\right)=\mathbf{y}-\widehat{\mathbf{y}}^{(j)}$$
so that the optimization problem can be written, equivalently
$$\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p [\mathbf{r}_j-\beta_j\mathbf{x}_j]^2+\lambda |\beta_j|\right\rbrace$$
hence$$\min\left\lbrace\frac{1}{2n}\sum_{j=1}^p \beta_j^2\|\mathbf{x}_j\|-2\beta_j\mathbf{r}_j^T\mathbf{x}_j+\lambda |\beta_j|\right\rbrace$$
and one gets
$$\beta_{j,\lambda} = \frac{1}{\|\mathbf{x}_j\|^2}S(\mathbf{r}_j^T\mathbf{x}_j,n\lambda)$$
or, if we develop
$$\beta_{j,\lambda} = \frac{1}{\sum_i x_{ij}^2}S\left(\sum_ix_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)$$
Again, if there are weights $\mathbf{\omega}=(\omega_i)$, the coordinate-wise update becomes
$$\beta_{j,\lambda,{\color{red}{\omega}}} = \frac{1}{\sum_i {\color{red}{\omega_i}}x_{ij}^2}S\left(\sum_i{\color{red}{\omega_i}}x_{i,j}[y_i-\widehat{y}_i^{(j)}],n\lambda\right)$$
The code to compute this componentwise descent is

lasso_coord_desc = function(X,y,beta,lambda,tol=1e-6,maxiter=1000){ beta = as.matrix(beta) X = as.matrix(X) omega = rep(1/length(y),length(y)) obj = numeric(length=(maxiter+1)) betalist = list(length(maxiter+1)) betalist[[1]] = beta beta0list = numeric(length(maxiter+1)) beta0 = sum(y-X%*%beta)/(length(y)) beta0list[1] = beta0 for (j in 1:maxiter){ for (k in 1:length(beta)){ r = y - X[,-k]%*%beta[-k] - beta0*rep(1,length(y)) beta[k] = (1/sum(omega*X[,k]^2))* soft_thresholding(t(omega*r)%*%X[,k],length(y)*lambda) } beta0 = sum(y-X%*%beta)/(length(y)) beta0list[j+1] = beta0 betalist[[j+1]] = beta obj[j] = (1/2)*(1/length(y))*norm(omega*(y - X%*%beta - beta0*rep(1,length(y))),'F')^2 + lambda*sum(abs(beta)) if (norm(rbind(beta0list[j],betalist[[j]]) - rbind(beta0,beta),'F') &lt; tol) { break } } return(list(obj=obj[1:j],beta=beta,intercept=beta0)) }

For instance, consider the following (simple) dataset, with three covariates

chicago = read.table("http://freakonometrics.free.fr/chicago.txt",header=TRUE,sep=";")

that we can “normalize” (or “standardize“)

X = model.matrix(lm(Fire~.,data=chicago))[,2:4] for(j in 1:3) X[,j] = (X[,j]-mean(X[,j]))/sd(X[,j]) y = chicago$Fire y = (y-mean(y))/sd(y) To initialize the algorithm, use the OLS estimate beta_init = lm(Fire~0+.,data=chicago)$coef

For instance

lasso_coord_desc(X,y,beta_init,lambda=.001) $obj [1] 0.001014426 0.001008009 0.001009558 0.001011094 0.001011119 0.001011119$beta [,1] X_1 0.0000000 X_2 0.3836087 X_3 -0.5026137   $intercept [1] 2.060999e-16 and we can get the standard Lasso plot by looping, # October, grant proposal season In 2012, Danielle Herbert, Adrian Barnett, Philip Clarke and Nicholas Graves published an article entitled “on the time spent preparing grant proposals: an observational study of Australian researchers“, whose conclusions had been included in Nature under a more explicit title, “Australia’s grant system wastes time” ! In this study, they included 3700 grant applications sent to the National Health and Medical Research Council, and showed that each application represented 37 working days: “Extrapolating this to all 3,727 submitted proposals gives an estimated 550 working years of researchers’ time (95% confidence interval, 513-589)“. But in these times when I have to write my funding application, I find that losing 37 days of work is huge. Because it’s become the norm! And somehow, it’s sad. Forget about the crazy idea that I would rather, in fact, spend more time doing my research. In fact, the thought I had this morning was that it is rather sad that in the Faculty of Science, mathematicians are asked to spend a considerable amount of time, comparable to that required of physicists or chemists, for often smaller amounts of funding… And I thought it could be easily verified. We start by retrieving the discipline codes url="http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSC-ResultatsCSS_eng.asp" download.file(url,destfile = "GSC.html") library(XML) tables=readHTMLTable("GSC.html") GSC=tables[[1]]$V1 GSC=as.character(GSC[-(1:2)]) namesGSC=tables[[1]]$V2 namesGSC=as.character(namesGSC[-(1:2)]) We’re going to need a small function, to remove the$ and other symbols that pollute the data (and prevent them from being treated as numbers)

library(stringr) Correction = function(x) as.numeric(gsub('[$,]', '', x)) We will now read the 12 pages, and harvest (we will just take the 2017 data, but we could go back a few years before) grants= function(gsc){ url=paste("http://www.nserc-crsng.gc.ca/NSERC-CRSNG/FundingDecisions-DecisionsFinancement/ResearchGrants-SubventionsDeRecherche/ResultsGSCDetail-ResultatsCSSDetails_eng.asp?Year=2017&amp;GSC=",gsc,sep="") download.file(url,destfile = "GSC.html") library(XML) tables=readHTMLTable("GSC.html") X=as.character(tables[[1]]$"Awarded Amount") A=as.numeric(Vectorize(Correction)(X)) return(c(median(A),mean(A),as.numeric(quantile(A,(1:99)/100)))) } M=Vectorize(grants)(GSC[1:12])

The average amounts of individual grants can be compared,

barplot(M[2,])

In mathematics, the average grant amount is $24400. If we normalize by this quantity, we obtain barplot(M[2,]/M[2,8]) In other words, the average amount of a (individual) grant in chemistry (to pay for students, conferences, etc.) is twice that in mathematics, 60% higher in physics than in maths… We can also look at the median values (rather than the averages) barplot(M[1,]) Here again, it is in mathematics that it is the weakest…. barplot(M[1,]/M[1,8]) in comparable proportions. If we think that the time spent writing should be proportional to the amount allocated, we should spend half as much time in math as in chemistry. Cumulative functions can also be ploted, plot(M[3:101,8],(1:99)/100,type="s",xlim=range(M)) lines(M[3:101,5],(1:99)/100,type="s",col="red") lines(M[3:101,4],(1:99)/100,type="s",col="blue") with math in black, physics in red, and chemistry in blue. What is surprising is the bottom part: a “bad” researcher in chemistry or physics will earn more than the median researcher in mathematics… Now that my intuition is confirmed, I have to go back, writing my proposal… and explain to my coauthors that I have to postpone some research projects because, well, you know… # Tests non-paramétriques et simulations Lors du dernier cours de statistique, nous avons présenter les tests d’ajustment de lois. Nous avions illustré ces tests à partir de la taille d’individus (déjà utilisé pour présenter l’ajustement de lois, et l’estimation de densité), correspondant à $\boldsymbol{x}=\{x_1,\cdots,x_n\}$. > Davis=read.table( + "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt") > Davis[12,c(2,3)]=Davis[12,c(3,2)] > X=Davis$height > n=length(X)

On notera $(x_{i:n})$ la statistique d’ordre, au sens où

$x_{1:n}\leq x_{2:n}\leq\cdots\leq x_{n-1:n}\leq x_{n:n}$

Parmi les outils graphiques, nous avons vu le PP plot (graphique probabilité-probabilité) et le QQ plot (graphique quantile). Le code pour créer un PP plot peut être le suivant

> PP=function(Y,F=pnorm){ +   n=length(Y) +   x=F(sort(Y)) +   y=seq(1/n/2,1-1/n/2,by=1/n) +   return(data.frame(x=x,y=y)) + }

qui représente (à un détail près) le nuage de points

$\left\{F_0(x_{i:n}),\frac{i}{n}\right\}$

et celui pour le QQ plot serait

> QQ=function(Y,Q=qnorm){ +   n=length(Y) +   x=Q(seq(1/n/2,1-1/n/2,by=1/n)) +   y=sort(Y) +   return(data.frame(x=x,y=y)) + }

qui représente (toujours à un détail près) le nuage de points

$\left\{F_0^{-1}\left(\frac{i}{n}\right),x_{i:n}\right\}$

où $F_0$ est la loi que l’on cherche à tester, au sens où on a $H_0:F=F_0$ avec comme hypothèse alternative $H_1:F\neq F_0$.

# Soutenance de thèse sur les inégalités

Ce jeudi, je participais au jury de thèse de Fattouma Souissi en tant que rapporteur, à l’université de Montpellier, sur les régressions PLS-Gini. Le manuscript de la thèse est désormais en ligne.

# La tarification en assurance, et la sélection à l’université

En assurance, même si les primes sont calculées ex ante, on distingue la tarification a priori et a posteriori. Dans la tarification a priori on va utiliser des variables dites explicatives pour affecter dans une classe tarifaire. Par exemple en assurance automobile, on va utiliser l’âge, le genre, le lieu d’habitation, la puissance du véhicule pour prédire, a priori, la charge sinistre pour l’année à venir. Le principe de la tarification repose sur la loi des grands nombres qui consiste à faire demande un prime égale à ce que l’on pense payer – en moyenne – par classe de risque (la compagnie est alors – en moyenne – à l’équilibre financier). Mais on peut utiliser une tarification dite a posteriori où on utilise de l’information passé (par exemple le nombre de sinistres survenus les 3 dernières années, information que l’on peut retrouver, en France, dans le mécanisme bonus-malus). Si j’oppose ici les deux, on imagine que ces deux informations sont intéressantes, et peut être complémentaires. En revanche, il n’y a pas de tarification ex post avec ce mécanisme. Il y aura éventuellement un bonus qui sera offert l’année prochaine (et non pas de manière rétroactive), mais l’assureur ne peut pas rétropédaler en disant “finalement, vu la manière dont vous conduisez, la prime que vous nous devez pour l’année en cours n’est pas de 600 euros, mais de 900 euros“. La seule possibilité est de rediscuter, au moment du renouvellement, le montant de la prime, au vu de l’expérience.

# 2015, October 21

I am almost sure I have already lived that day before,

# “A 99% TVaR is generally a 99.6% VaR”

Almost 6 years ago, I posted a brief comment on a sentence I found surprising, by that time, discovered in a report claiming that

the expected shortfall […] at the 99 % level corresponds quite closely to the […] value-at-risk at a 99.6% level

which was inspired by a remark in Swiss Experience report,

expected shortfall […] on a 99% confidence level […} corresponds to approximately 99.6% to 99.8% Value at Risk

# Another Interactive Map for the Cholera Dataset

Following my previous post, François (aka @FrancoisKeck) posted a comment mentionning another package I could use to get an interactive map, the rleafmap package. And the heatmap was here easy to include.

The first part is still the same, to get the data,

> require(rleafmap)
> library(sp)
> library(rgdal)
> library(maptools)
> library(KernSmooth)
> setwd("/home/arthur/Documents/")
> df_deaths <- data.frame(deaths@coords)
> coordinates(df_deaths)=~coords.x1+coords.x2
> proj4string(df_deaths)=CRS("+init=epsg:27700")
> df_deaths = spTransform(df_deaths,CRS("+proj=longlat +datum=WGS84"))
> df=data.frame(df_deaths@coords)

To get a first visualisation, use

> stamen_bm <- basemap("stamen.toner")
> j_snow <- spLayer(df_deaths, stroke = FALSE)
> writeMap(stamen_bm, j_snow, width = 1000, height = 750, setView = c( mean(df[,1]),mean(df[,2])), setZoom = 14)

and again, using the + and the – in the top left area, we can zoom in, or out. Or we can do it manually,

> writeMap(stamen_bm, j_snow, width = 1000, height = 750, setView = c( mean(df[,1]),mean(df[,2])), setZoom = 16)

To get the heatmap, use

> library(spatstat)
> library(maptools)

> win <- owin(xrange = bbox(df_deaths)[1,] + c(-0.01,0.01), yrange = bbox(df_deaths)[2,] + c(-0.01,0.01))
> df_deaths_ppp <- ppp(coordinates(df_deaths)[,1],  coordinates(df_deaths)[,2], window = win)
>
> df_deaths_ppp_d <- density.ppp(df_deaths_ppp,
sigma = min(bw.ucv(df[,1]),bw.ucv(df[,2])))

> df_deaths_d <- as.SpatialGridDataFrame.im(df_deaths_ppp_d)
> df_deaths_d$v[df_deaths_d$v < 10^3] <- NA

> stamen_bm <- basemap("stamen.toner")
> mapquest_bm <- basemap("mapquest.map")

> j_snow <- spLayer(df_deaths, stroke = FALSE)
> df_deaths_den <- spLayer(df_deaths_d, layer = "v", cells.alpha = seq(0.1, 0.8, length.out = 12))
> my_ui <- ui(layers = "topright")

> writeMap(stamen_bm, mapquest_bm, j_snow, df_deaths_den, width = 1000, height = 750, interface = my_ui, setView = c( mean(df[,1]),mean(df[,2])), setZoom = 16)

The amazing thing here are the options in the top right corner. For instance, we can remove some layers, e.g. to remove the points

or to change the background

To get an html file, instead of a standard visualisation in RStudio, use

> writeMap(stamen_bm, mapquest_bm, j_snow, df_deaths_den, width = 450, height = 350, interface = my_ui, setView = c( mean(df[,1]),mean(df[,2])), setZoom = 16, directView ="browser")

which will generate the html table (as well as some additional files actually) above. Awesome, isn’t it?

# Somewhere else, part 176

yes, extreme winds  cause a waterfall (in England) to blow upward