Daily Archives: 29/05/2014

Demography, Pension, Statistics

Your Life in Weeks

1 Comment

This week, I discovered a picture on http://waitbutwhy.com/, which represent a (so-called) typical human life, in weeks,

I found that interesting. But the first problem is that I don’t understand the limit, below: 90 years, that’s not the average life length. That’s not what you should expect to live when you get born. The second problem is that it cannot be as static as it might seem, when you look at the picture. I mean, life expectancy at age 0 is not the same as life expectancy at age 30, or 50. So I did try to make an animated graph, using prospective life tables. Here a code to generate life tables, at different period, for a French population (I distinguish, here male and female)

library(demography)
france.LC1 <- lca(fr.mort,adjust="e0",series="female",years=c(1900,2100))
france.fcast <- forecast(france.LC1,h=100)
L2 <- lifetable(france.fcast)
ex2=L2$ex
L1=lifetable(fr.mort,series="female")
ex1=L1$ex
exF=cbind(ex1,ex2)
france.LC1 <- lca(fr.mort,adjust="e0",series="male",years=c(1900,2100))
france.fcast <- forecast(france.LC1,h=100)
L2 <- lifetable(france.fcast)
ex2=L2$ex
L1=lifetable(fr.mort,series="male")
ex1=L1$ex
exM=cbind(ex1,ex2)
Y=colnames(exF)

Based on those lifetables, we can extract remaining life expectancy, at various ages (say, for instance 50, 51, 52, etc), for someone born on some given year (say 1950). Based on those expected remaining lifetimes, we can plot

picture=function(yearborn=1950,age=50){
k=which(Y==yearborn)
M=diag(exM[,k+0:100])
F=diag(exF[,k+0:100])
par(mfrow=c(1,2))
va=0:(52*100-1)
plot(va%%52,va%/%52,cex=.6,pch=15,col=c("light yellow","light blue","white")[1+
(va>=age*52)*1+(va>(age+M[age+1])*52)*1],ylim=c(100,0),axes=FALSE,xlab="Week",
ylab="Age",main=paste("Man, born on ",yearborn,
", age ",age,sep=""))
axis(1)
axis(2)
plot(va%%52,va%/%52,cex=.6,pch=15,col=c("light yellow","pink","white")[1+
(va>=age*52)*1+(va>(age+F[age+1])*52)*1],ylim=c(100,0),axes=FALSE,xlab="Week",
ylab="Age",main=paste("Woman, born on ",yearborn,
", age ",age,sep=""))
axis(1)
axis(2)}

For instance, if we want the graph above, for someone age 30, born in 1980, we use

picture(1980,30)

Now, if we run a code to get an animated gif, we can get, for someone born in 1950,

and for someone born in 2000

Now, if I could get historical datasets, with the average time spent in schools, ages of retirement, etc, I guess I could add it on the graph. But that’s another story…

2014-E

Régression linéaire, quelques codes

Leave a comment

Un rapide billet pour mettre en ligne les codes utilisés la semaine passée, complétant les codes des transparents. On travaille toujours sur la même base, ou on cherche à prévoir une distance de freinage d’un véhicule, tenant compte de la vitesse du véhicule.

> plot(cars)
> reg=lm(dist~speed,data=cars)
> summary(reg)

Call:
lm(formula = dist ~ speed, data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.069  -9.525  -2.272   9.215  43.201 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -17.5791     6.7584  -2.601   0.0123 *  
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared:  0.6511,	Adjusted R-squared:  0.6438 
F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Pour faire plusieurs prévisions, à la main, on peut utiliser le code suivant (la boucle permet de faire des prévisions pour plusieurs valeurs)

> for(x in seq(3,30,by=.25)){
+ b0=coef(reg)[1]
+ b1=coef(reg)[2]
+ Yx=b0+b1*x
+ V=vcov(reg)
+ Vx=V[1,1]+2*V[1,2]*x+V[2,2]*x^2
+ IC1=Yx+c(-1,+1)*1.96*sqrt(Vx)
+ s=summary(reg)$sigma
+ IC2=Yx+c(-1,+1)*1.96*s
+ points(x,Yx,pch=19,col="red")
+ points(c(x,x),IC1,pch=3,col="blue")
+ points(c(x,x),IC2,pch=3,col="purple")}

On avait ensuite fait une régression linéaire sur une sous-base, avec 20 observations tirées au hasard

> I=sample(1:50,size=20)
> reg=lm(dist~speed,data=cars[I,])

Le but était de visualiser l’impact du nombre d’observation sur la qualité de la régression

> summary(reg)

Call:
lm(formula = dist ~ speed, data = cars[I, ])

Residuals:
    Min      1Q  Median      3Q     Max 
-23.529  -7.998  -5.394  11.634  39.348 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -20.7408     9.4639  -2.192   0.0418 *  
speed         4.2247     0.6129   6.893 1.91e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.62 on 18 degrees of freedom
Multiple R-squared:  0.7252,	Adjusted R-squared:   0.71 
F-statistic: 47.51 on 1 and 18 DF,  p-value: 1.91e-06

> for(x in seq(3,30,by=.25)){
+   b0=coef(reg)[1]
+   b1=coef(reg)[2]
+   Yx=b0+b1*x
+   V=vcov(reg)
+   Vx=V[1,1]+2*V[1,2]*x+V[2,2]*x^2
+   IC=Yx+c(-1,+1)*1.96*sqrt(Vx)
+   points(x,Yx,pch=19,col="purple")
+   points(c(x,x),IC,pch=3,col="green")}

Notons qu’il est possible d’utiliser des fonctions de R pour faire des prévisions, avec des intervalles de confiance

> predict(reg,
+ newdata=data.frame(speed=c(15,25)),interval= "confidence")
       fit      lwr       upr
1 42.62976 34.75450  50.50502
2 84.87677 68.92746 100.82607
> predict(reg,
+ newdata=data.frame(speed=15),interval= "prediction")
       fit      lwr      upr
1 42.62976 6.836077 78.42344

Quand on a plus d’une variable explicative, c’est plus compliqué de “visualiser” la régression

>  chicago=read.table("http://freakonometrics.free.fr/chicago.txt",
+  header=TRUE,sep=";")
>  Y=chicago$Fire
>  X1=chicago$X_1
>  X2=chicago$X_2
>  X3=chicago$X_3
>  base=data.frame(Y,X1,X2,X3)
> plot(X2,X3)
> reg=lm(Y~X2+X3,data=base)
> y=function(x2,x3) predict(reg,newdata=data.frame(X2=x2,X3=x3))
> VX2=seq(0,80)
> VX3=seq(5,25)
> VY=outer(VX2,VX3,y)
> image(VX2,VX3,VY)
> contour(VX2,VX3,VY,add=TRUE)

qui correspond à un plan de régression

> persp(VX2,VX3,VY,theta=30,ticktype=detailed)

On reviendra plus en détails sur ce point, mais il est possible de faire des régressions non linéaires assez facilement, à partir de ce modèle linéaire. On avait commencé par un modèle linéaire sur le logarithme de la distance

> plot(cars$speed,log(cars$dist))
> reg1=lm(log(dist)~speed,data=cars)
> abline(reg1,col="red")

(on le verra, ce n’est pas fini, car on n’a pas ici de prévision sur la distance, juste sur son logarithme… mais promis, on en reparlera) ou sur la racine carrée

> plot(cars$speed,sqrt(cars$dist))
> reg1=lm(sqrt(dist)~speed,data=cars)
> abline(reg1,col="red")

Au lieu de transformer la variable d’intérêt, on peut aussi transformer la variable explicative. On peut pendre des puissances, ou des fonctions simples, mais aussi mettre des ruptures. On avait commencé par une variable indicatrice,

> plot(cars$speed,cars$dist)
> s=10
> abline(v=s,col="green")
> regs=lm(dist~speed+I(speed>s),data=cars)
> summary(regs)

Call:
lm(formula = dist ~ speed + I(speed > s), data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.472  -9.559  -2.088   7.456  44.412 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      -17.2964     6.7709  -2.555   0.0139 *  
speed              4.3140     0.5762   7.487  1.5e-09 ***
I(speed > s)TRUE  -7.5116     7.8511  -0.957   0.3436    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.39 on 47 degrees of freedom
Multiple R-squared:  0.6577,	Adjusted R-squared:  0.6432 
F-statistic: 45.16 on 2 and 47 DF,  p-value: 1.141e-11

Mais on peut aussi mettre des fonctions afin d’avoir un modèle linéaire par morceaux, tout en étant continu

> plot(cars)
> s=15
> abline(v=s,col="green")
> positive=function(x) ifelse(x>0,x,0)
> regs=lm(dist~speed+positive(speed-s),data=cars)
> summary(regs)

Call:
lm(formula = dist ~ speed + positive(speed - s), data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.502  -9.513  -2.413   5.195  45.391 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)   
(Intercept)          -7.6519    10.6254  -0.720  0.47500   
speed                 3.0186     0.8627   3.499  0.00103 **
positive(speed - s)   1.7562     1.4551   1.207  0.23350   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.31 on 47 degrees of freedom
Multiple R-squared:  0.6616,	Adjusted R-squared:  0.6472 
F-statistic: 45.94 on 2 and 47 DF,  p-value: 8.761e-12

On a ici une rupture, mais on pourrait imaginer en avoir plusieurs

> nouvellebase=data.frame(speed=5:25)
> y=predict(regs,newdata=nouvellebase)
> lines(5:25,y,col="red")
> 
> plot(cars$speed,cars$dist)
> s1=10
> s2=20
> abline(v=c(s1,s2),col="green")
> positive=function(x) ifelse(x>0,x,0)
> regs=lm(dist~speed+positive(speed-s1)+positive(speed-s2),data=cars)
> summary(regs)

Call:
lm(formula = dist ~ speed + positive(speed - s1) + positive(speed - s2), data = cars)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.374  -9.475  -2.625   6.639  43.914 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)  
(Intercept)           -7.6305    16.2941  -0.468   0.6418  
speed                  3.0630     1.8238   1.679   0.0998 .
positive(speed - s1)   0.2087     2.2453   0.093   0.9263  
positive(speed - s2)   4.2812     2.2843   1.874   0.0673 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15 on 46 degrees of freedom
Multiple R-squared:  0.6821,	Adjusted R-squared:  0.6613 
F-statistic: 32.89 on 3 and 46 DF,  p-value: 1.643e-11

Comme vu en cours, le test de significativité des deux derniers coefficients ne veut pas dire que la pente est nulle, mais qu’elle est significativement différente de cette obtenue sur la zone de gauche (avant les deux seuils).