# Testing for Covid-19 in the U.S.

For almost a month, on a daily basis, we are working with colleagues (Romuald, Chi and Mathieu) on modeling the dynamics of the recent pandemic. I learn of lot of things discussing with them, but we keep struggling with the tests. Paul, in Montréal, helped me a little bit, but I think we will still have to more to get a better understand. To but honest, we stuggle with two very simple questions

• how many people are tested on a daily basis ?

Recently, I discovered Modelling COVID-19 exit strategies for policy makers in the United Kingdom, which is very close to what we try to do… and in the document two interesting scenarios are discussed, with, for the first one, “1 million ‘reliable’ daily tests are deployed” (in the U.K.) and “5 million ‘useless’ daily tests are deployed”. There are about 65 millions unhabitants in the U.K. so we talk here about 1.5% people tested, on a daily basis, or 7.69% people ! It could make sense, but our question was, at some point, is that realistic ? where are we today with testing ? In the U.S. https://covidtracking.com/ collects interesting data, on a daily basis, per state.

url = "https://raw.githubusercontent.com/COVID19Tracking/covid-tracking-data/master/data/states_daily_4pm_et.csv" download.file(url,destfile="covid.csv") base = read.csv("covid.csv")

Unfortunately, there is no information about the population. That we can find on wikipedia. But in that table, the state is given by its full name (and the symbol in the previous dataset). So we new also to match the two datasets properly,

url="https://en.wikipedia.org/wiki/List_of_states_and_territories_of_the_United_States_by_population" download.file(url,destfile = "popUS.html") #pas contaminé 2/3 R=3 library(XML) tables=readHTMLTable("popUS.html") T=tables[][3:54,c("V3","V4")] names(T)=c("state","pop") url="https://en.wikipedia.org/wiki/List_of_U.S._state_abbreviations" download.file(url,destfile = "nameUS.html") tables=readHTMLTable("nameUS.html") T2=tables[][13:63,c(1,4)] names(T2)=c("state","symbol") T=merge(T,T2) T$population = as.numeric(gsub(",", "", T$pop, fixed = TRUE)) names(base)="symbol" base = merge(base,T[,c("symbol","population")])

Now our dataset is fine… and we can get a function to plot the number of people tested in the U.S. (cumulated). Here, we distinguish between the positive and the negative,

drawing = function(st ="NY"){ sbase=base[base$symbol==st,c("date","positive","negative","population")] sbase$DATE = as.Date(as.character(sbase$date),"%Y%m%d") sbase=sbase[order(sbase$DATE),] par(mfrow=c(1,2)) plot(sbase$DATE,(sbase$positive+sbase$negative)/sbase$population,ylab="Proportion Test (/population of state)",type="l",xlab="",col="blue",lwd=3) lines(sbase$DATE,sbase$positive/sbase$population,col="red",lwd=2) legend("topleft",c("negative","positive"),lwd=2,col=c("blue","red"),bty="n") title(st) plot(sbase$DATE,sbase$positive/(sbase$positive+sbase$negative),ylab="Ratio of positive tests",ylim=c(0,1),type="l",xlab="",col="black",lwd=3) title(st)} Let us start with New York drawing("NY") As at now, 4% of the entiere population got tested… over 6 weeks…. The graph on the right is the proportion of people who tested positive… I won’t get back on that one here today, I keep it for our work. In New Jersey, we got about 2.5% of the entiere population tested, overall, drawing("NJ") Let us try a last one, Florida drawing("FL") As at today, it is 1.5% of the population, over 6 weeks. Overall, in the U.S. less than 0.1% people are tested, on a daily basis. Which is far from the 1.5% in the U.K. scenarios. Now, here come the second question, • what are we actually testing for ? On that one, my experience in biology is… very limited, and Paul helped me. He mentioned this morning a nice report, from a lab in UC Berkeley One of my question was for instance, if you get tested positive, and you do it again, can you test negative ? Or, in the context of our data, do we test different people ? are some people tested on a regular basis (perhaps every week) ? For instance, with antigen tests (Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR) – also called molecular or PCR – Polymerase Chain Reaction – test) we test if someone is infectious, while with antibody test (using serological immunoassays that detect viral-specific antibodies — Immunoglobin M (IgM) and G (IgG) — also called serology test), we test for immunity. Which is rather different… I have no idea what we have in our database, to be honest… and for the past six weeks, I have seen a lot of databases, and most of the time, I don’t know how to interpret, I don’t know what is measured… and it is scary. So, so far, we try to do some maths, to test dynamics by tuning parameters “the best we can” (and not estimate them). But if anyone has good references on testing, in the context of Covid-19 (for instance on specificity, sensitivity of all those tests) I would love to hear about it ! # “Family History” and Life Insurance Today and tomorrow, I will attend the Online International Conference in Actuarial science, data science and finance, organised by colleagues in Lyon. But I won’t be in Lyon, I will be at home, in Montréal… I will give a talk on wednesday afternoon, on a joint paper with Ewen Gallic and Olivier Cabrignac. Slides are available here, and if I can get a copy of the video, I will share it… # Telematics: a revolution ? Tomorrow morning, Laurence Barry will present some joint work at the Online International Conference in Actuarial science, data science and finance, organised by colleagues in Lyon. The paper, “Personalization as a Promise: Can Big Data Change the Practice of Insurance?” is online, as a Working Paper of the PARI chair, The purpose of this paper is to measure the impact of technologies from the Big Bang. data on thehe pricing ofhe products car insurance. The first part describes how the aggregated view buildsuit by statistics enables highlighting invisible regularities at the individual level. Despite a very granular segmentations in automobile insurance, the approach remained classificatory, hypothesizing the risk identity of individuals from the same class. The second part highlights the reversal of big data-induced perspective in the’analysis ofgiven ; awith theur volume and the new algorithms, the aggregate viewpoint is questioned.. The hypothesis of class homogeneity is becoming increasingly difficult to test. maintain, especially since predictive analysis boasts the ability to predict the rs results at the individual level. The third part is studying the’influence of telematics boxes able to import the new pinsurance aradigm automobile. However, a reading of the most recent research articles on a pricing automobile including this new monter that the epistemological leap, at least for now, has not taken place. # INF7100, à distance Suite au contexte de la covid-19, la session d’été commencera une semaine plus tard. Le cours INF7100 Initiation à la science des données et à l’intelligence artificielle – qui devait commencer mardi 28 avril commencera le 5 mai prochain. Les étudiants inscrits recevront un lien pour se connecter, pour une présentation globale. On le donnera conjointement avec Marie-Jean Meurs (du département d’informatique) et Jean-Hugues Roy (de l’école des médias). Je posterais des informations dans les semaines à venir sur les points que je vais aborder (essentiellement sur la science des données). # Estimer la surmortalité Ce matin, Baptiste Coulmont m’envoyait un tweet avec un joli graphique, présentant le nombre quotidien de décès en France. Comme les données sont en ligne, je me suis dit que je pourrais jouer un peu avec. Pour les plus curieux, on a la liste de tous les décès depuis… longtemps ! (on a plus de 25 millions de décès). La seule chose qui nous intéresse, c’est la date, alors on la récupère. Puis on compte, par jour, combien de décès on a eu, D = read.csv("insee_deces.csv",header=TRUE) Vecteur_Dates = as.character(D[,8]) TV = table(Vecteur_Dates) B = data.frame(dateC=names(TV),nb=as.numeric(TV)) B$date=as.Date(B$dateC,"%Y-%m-%d") B$year = format(B$date,"%Y") B$month = format(B$date,"%m") B$day = format(B$date,"%d") B$deb=as.Date(paste(B$year,"01-01",sep="-"),"%Y-%m-%d") B$fin=as.Date(paste(B$year,"12-31",sep="-"),"%Y-%m-%d") B$dif=as.numeric(B$date-B$deb)/as.numeric(B$fin-B$deb)

Je crée ici une une variable qui me dit où je suis dans l’année en cours (0 pour le 1er janvier, et 1 pour pour 31 décembre). Ça me permet de contourner un petit problème technique : les années bissextiles. On peut ensuite représenter toutes les années,

plot(B2$dif[B2$year == "2020"],B2$NormCpte[B2$year == "2020"],type="l",xlim=0:1,ylim=c(1200,4000),col="white") for(i in as.character(2000:2019)) lines(B2$dif[B2$year == i],B2$NormCpte[B2$year == i],col="light blue")

On veut ensuite avoir la tendance moyenne. Baptiste prenait la moyenne par jour (la courbe noire). Ici, on va lisser avec des splines

library(gam) reg = gam(NormCpte~bs(dif,40),data=B2[B2$year != "2020",]) vx = seq(0,1,length=501) vy = predict(reg,newdata=data.frame(dif=vx)) lines(vx,vy,col="blue",lwd=3) ce qui donne la courbe bleue. C’est l’écart à cette courbe qui donne la surmortalité. Par exemple (prenons un exemple bien connu, et largement visible sur le graphique précédant) l’année 2003, et les quinze premiers jours d’août, correspondant (plus ou moins) à la période de canicule. i = 2003 x = B2$dif[(B2$year == i)&amp;((B2$month == "08"))&amp;(B2$day %in% c(paste("0",1:9,sep=""),10:15))] y = B2$NormCpte[(B2$year == i)&amp;((B2$month == "08"))&amp;(B2$day %in% c(paste("0",1:9,sep=""),10:15))] yp = predict(reg,newdata=data.frame(dif=x)) e = y-yp sum(e)  14294.39 Plusieurs chiffres sont mentionnés sur wikipédia, dont le rapport de l’INSERM qui annonçait 14,802 morts. On n’est pas trop loin… On peut visualiser cette surmortalité sur le graphique ci-dessous plot(B2$dif[B2$year == "2020"],B2$NormCpte[B2$year == "2020"],type="l",col="white",xlim=0:1,ylim=c(1200,4000)) for(i in as.character(2000:2019)) lines(B2$dif[B2$year == i],B2$NormCpte[B2$year == i],col="light blue") lines(vx,vy,col="blue") lines(B2$dif[B2$year == 2003],B2$NormCpte[B2$year == 2003],col="red") for(u in 1:length(x)) segments(x[u],y[u],x[u],yp[u],lwd=3,col="red") Prenons un autre exemple qu’on voit bien sur le graphique de Baptiste, la grippe de 2017 (dite de 2016-2017). Si on prend les 45 premiers jours de 2017, on obtient les chiffres suivants i = 2017 x=B2$dif[(B2$year == i)&amp;((B2$month == "01")| ((B2$month == "02")&amp;(B2$day %in% c(paste("0",1:9,sep=""),10:15))))] y=B2$NormCpte[(B2$year == i)&amp;((B2$month == "01")| ((B2$month == "02")&amp;(B2$day %in% c(paste("0",1:9,sep=""),10:15))))] yp=predict(reg,newdata=data.frame(dif=x)) yp = predict(reg,newdata=data.frame(dif=x)) e = y-yp sum(e)  21177.33 Le rapport de Santé Publique France parle de 21,200 décès (pour la “surmortalité toute cause”) ce qui est là aussi comparable… Maintenant, avant de conclure, et avant de me faire troller pendant quelques jours, je tiens à préciser que c’est juste un exercice de modélisation, et oui, j’ai exclu 2020. Je ne tiens pas à faire de comparaison avec la Covid-19, à qui 21,000 morts sont attribués, en France, ce soir. Si l’ordre de grandeur est le même, qu’on en me fasse pas dire qu’on est dans une situation comparable, et que la Covid-19 est finalement une simple grippe. Je ne me prononcerais pas sur ce point car je ne suis pas virologue. Mais surtout, on a 21,000 morts malgré des mesures incroyables de confinement. Je reviendrais probablement dans les jours qui viennent sur l’impact du contrôle par confinement, mais l’objet du court billet de ce soir était de voir comment on pouvait rapidement quantifier une surmortalité. • petite note technique : pour les 3 graphiques (mais pas pour les codes et les calculs) j’avais initialement normalisé les données. En effet, la population française a relativement pas mal changé en 20 ans, augmentant de 10%… Aussi, 1000 morts en 2003, ça n’est pas autant (relativement à la taille de la population) qu’en 2020. # Une courte histoire des expériences randomisées On a beaucoup parlé ces derniers mois des expériences randomisées. Cet automne pour commencer, lorsque le prix Nobel d’économie a été attribué à Esther Duflo, Michael Kremer et Abhijit Banerjee. Plus récemment, lors de la polémique, en France, où, dans une conférence donnée le 13 février 2020, intitulée contre la méthode, un médecin marseillais affirmait « moi je n’ai jamais fait d’essais randomisés […] faire ça sur des maladies infectieuses, ça n’a pas de sens ». # Les expériences randomisées Quand on se demande l’effet d’un traitement (par rapport à un autre), la bonne méthode consiste à se poser la question « que ce serait-il passé si la personne avait choisi l’autre option ? » – ce que John Stuart Mill parlait de la « méthode des différences ». On parle ici d’inférence causale car on ne peut pas savoir ce qu’aurait donné l’autre option, qui n’a pas été retenue. L’effet causal d’une option, ou d’un traitement, n’est jamais observé car il s’agit de la différence potentielle entre deux alternatives. La méthode la plus simple consiste à prendre deux patients proches (à défaut d’être identiques), et de voir comment les deux répondent aux deux options, la proximité se mesurant en fonction de covariables (comme l’âge, le genre, etc), et les patients n’ayant pas choisi les options proposées. Il y a trois points importants dans cette approche : (1) la notion de groupe de contrôle, autrement dit il faut qu’il y ait des personnes « non traitées » (2) la notion d’aveuglement, ou d’absence d’information, autrement dit, les patients ne doivent pas savoir quel traitement ils ont eu (3) la notion d’affectation aléatoire, autrement dit, le traitement doit être choisi de manière totalement indépendante des patients. Ce sont des points qui opposent fondamentalement les expériences randomisées aux données observationnelles, comme les oppose par exemple Paul Rosenbaum dans Observation and Experiment. Si la méthode d’expérience randomisée a été formalisée par Ronald Fisher en 1934, nous allons revenir ici sur deux moments importants dans l’histoire de cette méthodologie. La première est l’utilisation de la randomisation pour faire des tests autour de 1890, en particulier avec les travaux de Peirce puis Richet, sur le spiritualisme et la transmission de pensée (qui ont très probablement influencé Fisher). Et la seconde est l’expérience qui a été menée en 1954 aux États-Unis, pour tester l’efficacité d’un vaccin contre la polio (et qui montre au passage que faire des expériences randomisées dans un contexte de maladies infectieuses a tout son sens). # STT5100, quiz (régression de Poisson #1) Pour la fin de la semaine, j’avais donné un petit quiz, basé sur la base suivante, qui donne le nombre de cyclistes à une intersection, à New York City, download.file("http://freakonometrics.free.fr/NYCVelo.RData","velo.RData") load("velo.RData") str(base) 'data.frame': 214 obs. of 7 variables:$ Date : chr "1-Apr-17" "2-Apr-17" "3-Apr-17" "4-Apr-17" ... $HIGH_T : num 46 62.1 63 51.1 63 48.9 48 55.9 66 73.9 ...$ LOW_T : num 37 41 50 46 46 41 43 39.9 45 55 ... $PRECIP : num 0 0 0.03 1.18 0 0.73 0.01 0 0 0 ...$ BB_COUNT: int 606 2021 2470 723 2807 461 1222 1674 2375 3324 ... $DAY : chr "Sam" "Dim" "Lun" "Mar" ...$ DIFF_T : num 9 21.1 13 5.1 17 7.9 5 16 21 18.9 ......

En utiisant une régression de Poisson, il fallait prédire combien de cyclistes passeront un dimanche, s’il fait une température maximale de 85F, minimale de 70F, et s’il ne pleut pas. Et il fallait voir ce que donnerait la prévision pour un lundi.

newbase = data.frame(DAY=as.factor(c("Lun","Dim")), HIGH_T=c(85,85),LOW_T=c(70,70), PRECIP=c(0,0))

Faisons un modèle avec toutes les variables explicatives.

reg0 = glm(BB_COUNT~HIGH_T+LOW_T+PRECIP+DAY,data=base,family=poisson)

Rajoutons aussi une indicatrice pour indiquer les jours où il ne pleut pas,

reg = glm(BB_COUNT~HIGH_T+LOW_T+PRECIP+I(PRECIP==0)+DAY,data=base,family=poisson) summary(reg)   Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 6.8844970 0.0110463 623.241 &lt;2e-16 *** HIGH_T 0.0210950 0.0003133 67.328 &lt;2e-16 *** LOW_T -0.0114006 0.0003351 -34.024 &lt;2e-16 *** PRECIP -0.6570450 0.0071899 -91.384 &lt;2e-16 *** I(PRECIP == 0)TRUE 0.1303908 0.0033283 39.176 &lt;2e-16 *** DAYJeu 0.1683475 0.0049690 33.880 &lt;2e-16 *** DAYLun 0.1144129 0.0050480 22.665 &lt;2e-16 *** DAYMar 0.1655886 0.0049936 33.160 &lt;2e-16 *** DAYMer 0.1688035 0.0049190 34.317 &lt;2e-16 *** DAYSam 0.0466447 0.0051838 8.998 &lt;2e-16 *** DAYVen 0.1003536 0.0050919 19.709 &lt;2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1   (Dispersion parameter for poisson family taken to be 1)   Null deviance: 70021 on 213 degrees of freedom Residual deviance: 26493 on 203 degrees of freedom AIC: 28580   Number of Fisher Scoring iterations: 4

Tout semble vraiment significatif. Et même avoir du sens (si on regarde les signes, par exemple). Et si on a peur de rater un effet non-linéaire, on peut mettre des splines sur toutes les variables continues

library(gam) reggam = gam(BB_COUNT~bs(HIGH_T)+bs(LOW_T)+bs(PRECIP)+I(PRECIP==0)+DAY,data=base,family=poisson) plot(reggam, se=TRUE) pour la température maximale, ou pour la température minimale et la courbe suivante pour la précipitation, avec un lissage linéaire entre l’observation maximale (3) et celle juste avant (vers 1.8) On peut aussi régresser sur la température minimale, et l’écart entre le maximum et le minimum (dans un modèle linéaire les modèles sont équivalents, mais avec de transformation nonlinéaire, la différence pourrait donner un modèle plus simple)

library(gam) reggam2 = gam(BB_COUNT~bs(DIFF_T)+bs(LOW_T)+bs(PRECIP)+I(PRECIP==0)+DAY,data=base,family=poisson) plot(reggam2, se=TRUE) On peut maintenant comparer les quatre modèles, et leurs prévision. Par exemple, pour le modèle linéaire (avec la variable indicatrice indiquant qu’il ne pleut pas),

P = predict(reg,newdata=newbase,type="response",se.fit=TRUE)

pour le lundi, on obtient l’intervalle de confiance à 95% pour $\widehat{\lambda}$

P$fit+c(-2,2)*P$se.fit  3349.842 3401.395

et pour le dimanche l’intervalle de confiance à 95% est

P$fit+c(-2,2)*P$se.fit  2987.497 3033.861

On peut visualiser les quatre intervalles de confiance (avec, de haut en bas, le second modèle gam, le premier, le modèle linéaire avec l’indicatrice d’absence de pluie, puis le premier modèle linéaire) alors que pour le dimanche, on a Autrement dit, en changeant de modèle, on change beaucoup les intervalles de confiance sur la prévision (et les intervalles sont parfois complètement disjoints). Ce qui n’est pas forcément une bonne nouvelle.

# De la démarche scientifique en période de crise

Dans une conférence donnée le 13 février 2020[i], intitulée contre la méthode, Didier Raoult affirmait « moi je n’ai jamais fait d’essais randomisés […] faire ça sur des maladies infectieuses, ça n’a pas de sens ». Cette vision était reprise dans une tribune plus détaillée, où face à « la méthode » (et « aux mathématiques »), Didier Raoult défendait (ce qu’il appelait) « la morale [et] l’humanisme » du serment d’Hippocrate. Comme il le rappelle, faire des groupes de contrôle, c’est « dire au malade qu’on va lui donner au hasard soit le médicament dont on sait qu’il marche, soit le médicament dont on ne sait pas s’il marche » (Raoult (2020a, 2020b)). Alors que cette méthode d’expériences randomisées est aujourd’hui saluée dans toutes les disciplines – comme le rappelle le prix Nobel d’économie attribué en 2019 à Esther Duflo, Michael Kremer et Abhijit Banerjee – comment un chercheur peut-il prendre une telle position, aujourd’hui ? Continue reading De la démarche scientifique en période de crise

# On the “correlation” between a continuous and a categorical variable

Let us get back on the Titanic dataset,

loc_fichier = "http://freakonometrics.free.fr/titanic.RData" download.file(loc_fichier, "titanic.RData") load("titanic.RData") base = base[!is.na(base$Age),] On consider two variables, the age $x$ (the continuous one) and the survivor indicator $y$ (the qualitative one) X = base$Age Y = base$Survived It looks like the age might be a valid explanatory variable in the logistic regression, summary(glm(Survived~Age,data=base,family=binomial)) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) -0.05672 0.17358 -0.327 0.7438 Age -0.01096 0.00533 -2.057 0.0397 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 964.52 on 713 degrees of freedom Residual deviance: 960.23 on 712 degrees of freedom AIC: 964.23 The significance test here has a $p$-value just below $4%$. Actually, one can relate it with the value of the deviance (the null deviance and the residual deviance). Recall that$$D=2\big(\log\mathcal{L}(\boldsymbol{y})-\log\mathcal{L}(\widehat{\boldsymbol{\mu}})\big)$$while$$D_0=2\big(\log\mathcal{L}(\boldsymbol{y})-\log\mathcal{L}(\overline{y})\big)$$Under the assumption that $x$ is worthless, $D_0-D$ tends to a $\chi^2$ distribution with 1 degree of freedom. And we can compute the $p$-value dof that likelihood ratio test, 1-pchisq(964.52-960.23,1)  0.03833717 (which is consistent with a Gaussian test). But if we consider a nonlinear transformation summary(glm(Survived~bs(Age),data=base,family=binomial)) Coefficients: Estimate Std. Error z value Pr(&gt;|z|) (Intercept) 0.8648 0.3460 2.500 0.012433 * bs(Age)1 -3.6772 1.0458 -3.516 0.000438 *** bs(Age)2 1.7430 1.1068 1.575 0.115299 bs(Age)3 -3.9251 1.4544 -2.699 0.006961 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 964.52 on 713 degrees of freedom Residual deviance: 948.69 on 710 degrees of freedom which seems to be “more significant” 1-pchisq(964.52-948.69,3)  0.001228712 So it looks like the variable $x$ is interesting here. To visualize the non-null correlation, one can consider the condition distribution of $x$ given $y=1$, and compare it with the condition distribution of $x$ given $y=0$, ks.test(X[Y==0],X[Y==1]) Two-sample Kolmogorov-Smirnov test data: X[Y == 0] and X[Y == 1] D = 0.088777, p-value = 0.1324 alternative hypothesis: two-sided i.e. with a $p$-value above $10%$, the two distributions are not significatly different. F0 = function(x) mean(X[Y==0]&lt;=x) F1 = function(x) mean(X[Y==1]&lt;=x) vx = seq(0,80,by=.1) vy0 = Vectorize(F0)(vx) vy1 = Vectorize(F1)(vx) plot(vx,vy0,col="red",type="s") lines(vx,vy1,col="blue",type="s") (we can also look at the density, but it looks like that there is not much to see) An alternative is discretize variable $x$ and to use Pearson’s independence test, k=5 LV = quantile(X,(0:k)/k) LV = 0 Xc = cut(X,LV) table(Xc,Y) Y Xc 0 1 (0,19] 85 79 (19,25] 92 45 (25,31.8] 77 50 (31.8,41] 81 63 (41,80] 89 53 chisq.test(table(Xc,Y)) Pearson's Chi-squared test data: table(Xc, Y) X-squared = 8.6155, df = 4, p-value = 0.07146 The $p$-value is here $7%$, with five categories for the age. And actually, we can compare the $p$-value pvalue = function(k=5){ LV = quantile(X,(0:k)/k) LV = 0 Xc = cut(X,LV) chisq.test(table(Xc,Y))$p.value} vk = 2:20 vp = Vectorize(pvalue)(vk) plot(vk,vp,type="l") abline(h=.05,col="red",lty=2) which gives a $p$-value close to $5$%, as soon as we have enough categories. In the slides of the course (STT5100), I claim that actually, the age is an important variable when trying to predict if a passenger survived. Test mentioned here are not as conclusive, nevertheless…

# Slides 21 – Poisson vs. Binomiale

Poisson ou Binomiale ? Application sur un problème de démographie. Les slides sont en ligne (slides 21) et la vidéo aussi (slides 21)

# Slides 20 – GLM et sélection de variables (stepwise)

Rapide discussion sur l’utilisation des méthodes de simulations en petite dimension dans les GLM. Les slides sont en ligne (slides 20) et la vidéo aussi (slides 20)

# Slides 19 – GLM et résultats non-asymptotiques

Rapide discussion sur l’utilisation des méthodes de simulations en petite dimension dans les GLM. Les slides sont en ligne (slides 19) et la vidéo aussi (slides 19)

# Slides 18 – Tests et GLM

Un peu de tests et d’intervalles de confiance pour les GLM. Les slides sont en ligne (slides 18) et la vidéo aussi (slides 18)

# Slides 17 – Sur-dispersion

Une série de quelques slides sur la dispersion. Les slides sont en ligne (slides 17) et la vidéo aussi (slides 17)

# Slides 16 – Poids et modèle Tweedie

Petite digression sur l’utilisation des poids en régression, et quelques slides sur la régression Tweedie. Les slides sont en ligne (slides 16) et la vidéo aussi (slides 16)