# From a random generator to a sample function

This week-end, I wrote a post since I had some trouble to generate a sample random sample with R, to reproduce one obtained by a co-author, with SAS (generated using Fishman and Moore (1982) used in function RANUNI). I was lucky since another contributor for that book, Christrophe Dutang, got the anwer to the last question I asked: is it possible to reproduce the random generator ? Yes, we can. And it is quite simple, if you use the appropriate library and the appropriate function,

```> library(randtoolbox)
> a <- 397204094
> b <- 0
> m <- 2^(31)-1
> set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123)
> runif(10)
[1] 0.7503961 0.3209120 0.1783896 0.9060334 0.3571171
[6] 0.2211140 0.7864383 0.3980819 0.1246652 0.1876858```

If you check in the previous post this is exactly what SAS gave us (and that I could not reproduce by myself). But that was only one part of my problems, since the goal was actually to reproduce indices for a training subsample for credit scoring issues.

I have to admit that I had never though about it before: how should we write a sample function? If values can be replaced, that is fine, we just have to split the unit interval correctly. Like

```> set.seed(95)
> (U=runif(10))
[1] 0.15171155 0.57584087 0.05309844 0.07044485 0.48887914 0.15276707
[7] 0.37405684 0.30006096 0.96997126 0.30373498
> set.seed(95)
> (R=sample(0:99,size=10,replace=TRUE))
[1] 15 57  5  7 48 15 37 30 96 30```

Here, we just truncate from the values obtained from the random generator. And that is just fine. But how do we write a code to sample without replacements? I mean, how do you get that :

```> (S=sample(0:99,size=10,replace=FALSE))
[1] 15 57  5  6 46 14 35 27 89 92```

My initial idea was to use the following technique. The first value is easy to get: we just split the unit interval into 100 subdivision (as before since for the first value, replacement or not, we don’t care) and see in which interval the random value is. And we remove that value from our sample. Then, we split the unit intervall into 99 subdivision, and see in which interval the random value is. It is the 10th? Fine, then our second value is the 10th from our sample (the first value has been removed). Then we split the unit interval in 98 subdivision, etc. The code I wrote to produce that algorithm was the following,

```> mysample1=function(N,unif){
+  n=length(N)
+  size=length(unif)
+  V0=N[trunc(unif[1]*n)+1]
+  N=N[-which(N==V0)]
+  V=V0
+  for(i in 2:length(unif)){
+    V0=N[trunc(unif[i]*(n-i+1))+1]
+    N=N[-which(N==V0)]
+    V=c(V,V0)}
+ return(V)}```

Unfortunetely, I could not reproduce the sample obtained with the R function,

```> mysample1(0:99,unif=U)
[1] 15 58  5  7 49 17 39 31 97 32
> S
[1] 15 57  5  6 46 14 35 27 89 92```

Since Christrophe is an expert on random generators, I did ask him, one more time. And he came up with the following code,

```> mysample2=function(N,unif){
+   integerset=1:length(N)
+   result=rep(NA,length(unif))
+   for(i in 1:length(unif)){
+     intchosen=integerset[ceiling(U[i]*(length(N)-i+1))]
+     integerset[intchosen]=length(integerset)
+     integerset=integerset[-length(integerset)]
+     result[i]=intchosen}
+   return(N[result])}```

which works just fine !

```> mysample2(0:99,unif=U)
[1] 15 57  5  6 46 14 35 27 89 92
> S
[1] 15 57  5  6 46 14 35 27 89 92```

So now, not only can we reproduce random numbers obtained with other software, we can also obtain the same samples indices, with or without replacement ! Thanks Christophe !

[May, 15th] Note that this is note the generator used in SAS, actually. In order to reproduce the sample function of SAS, the algorithm is much more simple, by clearly not that efficicient since we generate a random sample of size 100 (if we have 100 observations), and then, we keep the values associated to the indices of the 10 smallest (if we want a sample of size 10). The code could be

```> mysample3=function(N,unif,size){
+ q=sort(unif)[size]
+ return(N[U<=q])}

> library(randtoolbox)
> a <- 397204094
> b <- 0
> m <- 2^(31)-1
> set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123) #OK

> U=runif(100)
> mysample3(1:100,U,size=10)
[1] 27 37 47 59 60 71 75 82 84 87```

Thanks Jean-Philippe for the idea (which works).

# Reproducibility and randomness

With Stéphane Tufféry, we were working this week on a chapter of a book, entitled Statistical Learning in Actuarial Science. The chapter should be based on R functions, and we wanted to reproduce some outputs he previously obtained with SAS. The good thing is that even complex functions (logistic regression, regression trees, etc) produce the same kind of outputs. But we found a problem that we could not fix: generating identical training subsets of observations… Out of 1,000 lines, we subsample about 600 lines. The problem is that we could not generate the same sets of indexes with R, and SAS. Even using similar random generators… (execpt if we want to extract 1 or 2 lines, no more).

Let us try to explain what’s going on (based on code produced by Stéphane). According to Eubank (2010), there are (at least) two generators of random numbers in SAS,

For instance, for the RAND function, if we generate a Gaussian sample – with Mersenne-Twister algorithm – the code should be

```%LET SEED =6;
%LET NREP=10;
DATA TESTRANDOM ;
CALL STREAMINIT(&SEED);
DO REP = 1 TO &NREP;
X = RAND ('NORMAL');
OUTPUT;
END;
RUN;
PROC PRINT DATA = TESTRANDOM ;
RUN ;```

And we get here

 Obs. REP X 1 1 2.10680 2 2 -0.25604 3 3 0.28692 4 4 -0.22806 5 5 1.34569 6 6 0.16341 7 7 -0.27788 8 8 0.02133 9 9 1.24050 10 10 1.01054

If we want a Uniform sample, it should be

```%LET SEED =6;
%LET NREP=10;
DATA TESTRANDM ;
CALL STREAMINIT(&SEED);
DO REP = 1 TO &NREP;
X = RAND ('UNIFORM');
OUTPUT;
END;
RUN;
PROC PRINT DATA = TESTRANDOM ;
RUN ;```
 Obs. REP X 1 1 0.66097 2 2 0.48044 3 3 0.87849 4 4 0.19916 5 5 0.04838 6 6 0.19966 7 7 0.81353 8 8 0.53807 9 9 0.01105 10 10 0.53753

On good thing (so far) about the latest, is that Mersenne-Twister has been coded in R, in the RNGkind function

```> RNGkind("Mersenne")
> set.seed(6)
> runif(10)
[1] 0.64357277 0.91590888 0.09523258 0.29537280
[5] 0.76993170 0.25589353 0.51789573 0.67784993
[9] 0.14722782 0.70052604```

But the output is different, even if we’re supposed to start, here, with the same seed. Now, if we want to make sure about what is done here, let us write our own codes of the Fishman and Moore (1982) algorithm (in order to reproduce the SAS output). The R version of

```> a = 397204094      # RANUNI multiplier
> seed = 123         # seed
> n = 10             # sample size
> m = (2^31) - 1     # period
> for (i in (1:n-1)) {
+ seed = (a*seed)%%m
+ u = seed / m
+ print(u)
+ }
[1] 0.7503961
[1] 0.3209121
[1] 0.3453204
[1] 0.2683455
[1] 0.241798
[1] 0.9888181
[1] 0.3943279
[1] 0.9710172
[1] 0.001632214
[1] 0.921537```

Let us now run a similar code with SAS,

```%LET SEED =123;
%LET NREP=10;
DATA FRANUNI (KEEP = x) ;
seed = &SEED ;
DO REP = 1 TO &NREP;
CALL RANUNI(seed, x);
OUTPUT;
END;
RUN;
PROC PRINT DATA = FRANUNI ;
RUN ;```

and we get the following output

 Obs. x 1 0.75040 2 0.32091 3 0.17839 4 0.90603 5 0.35712 6 0.22111 7 0.78644 8 0.39808 9 0.12467 10 0.18769

It looks like here, indeed, we start with the same seed, since the first two numbers generated are similar. But then, it looks like we really have random numbers… If we change the seed, the first two numbers are similar, but that’s all.

We might be missing something trivial here, but we did not see it. So if anyone has a clue about reproducibility issues when generating random samples, with R and SAS, we are interested !

# La tarification avec SAS

En tarification, il est possible d’utiliser d’autres logiciels que R, en particulier, il semble que l’on puisse faire deux ou trois choses avec SAS…. J’en parle un peu car il semble  que, paradoxalement, les asssureurs préfèrent encore SAS à R (par exemple). Et comme plusieurs étudiants m’avaient demandé “et comment on fait avec SAS ?“. Bon, par contre je ne mets que les choses de base, parce que SAS est assez limité sur ce qu’il peut faire….

Pour suivre un peu le plan du cours, la première étape est de définir une variable d’exposition dans la table,

```DATA contrats;
SET lib.contrats;
lnexpo = log(expo);
run;```

Pour faire une régression de Poisson, ce n’est pas forcément compliqué,

```PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1_Param
Type3=Genmod1_Var
Modelfit=Genmod1_InfoModele;
MODEL nbsin = ageconducteur /
dist = poisson
offset = lnexpo
type3;
RUN; QUIT;```

La sortie SAS a alors l’allure suivante

```                                  The GENMOD Procedure
Critère pour évaluer la qualité de l'ajustement
Critère                   DF          Valeur       Valeur/DF
Deviance                63E3      26872.5334          0.4237
Scaled Deviance         63E3      26872.5334          0.4237
Pearson Chi-Square      63E3      73275.5362          1.1553
Scaled Pearson X2       63E3      73275.5362          1.1553
Log Likelihood                   -18474.2667

Algorithm converged.
Analyse des résultats estimés de paramètres

Erreur      Wald 95Limites
Paramètre    DF   Estimation   standard      de confiance %    Khi 2   Pr > Khi 2
Intercept     1      -3.5164     0.0851    -3.6832  -3.3496   1708.02       <.0001
ageconducteur 1       0.0168     0.0014     0.0141   0.0195    146.73       <.0001
Scale         0       1.0000     0.0000     1.0000   1.0000
NOTE: The scale parameter was held fixed.

Statistiques LR pour Analyse de Type 3
Source           DF      Khi 2    Pr > Khi 2
ageconducteur     1     148.72        <.0001```

Il est aussi possible de faire des GAM (i.e. du lissage de la variable explicative – continue – avec des fonctions splines)

```PROC GAM DATA = base;
MODEL nbsin = spline(ageconducteur) / dist = Poisson;
OUTPUT OUT=gam PREDICTED;
RUN; QUIT;
PROC SORT DATA = gam NODUPKEY; BY age_cond; RUN; QUIT;```

et on peut faire des prédictions avec ce modèle (la sortie n’apporte pas grand chose, en pratique),

```DATA gam;
SET gam;
pred_nbsin_gam = exp(P_nbsin);
KEEP ageconducteur pred_nbsin_gam;
RUN;```

Enfin, on peut tenter de faire un joli graphique. Pour cela, on calcule les prédictions de trois modèles, le premier étant des nombres moyens de sinistres par âge

```PROC SORT DATA = base; BY ageconducteur; RUN; QUIT;
PROC MEANS DATA = base NOPRINT;
BY ageconducteur;
VAR nbsin;
WEIGHT expo;
OUTPUT OUT = nbsin_age (DROP = _TYPE_ _FREQ_) MEAN=mo
y_uni_nbsin;
RUN; QUIT;```

ensuite, on fait un modèle GLM, et  un modèle GAM, et on récupère les sorties

```PROC SORT DATA = nbsin_age; BY age_cond; RUN; QUIT;
PROC SORT DATA = gam; BY age_cond; RUN; QUIT;
DATA nbsin_age;
MERGE nbsin_age
gam;
BY age_cond;
RUN;```

On essaye de faire le dessin (je passe les lignes de commande, il y en a une vingtaine)

Pour faire une régression quasiPoisson, le code a l’allure suivante,s

```PROC GENMOD DATA = base;
ODS OUTPUT ParameterEstimates=Genmod1bis_Param
Type3=Genmod1bis_Var
Modelfit=Genmod1bis_InfoModele;
MODEL nbsin = ageconducteur /
dist = poisson
offset = lnexpo
type3
scale = deviance;
RUN; QUIT;```

La sortie donne alors l’estimation du paramètre de surdispersion (ou sur cet exemple de sousdispersion)

```                      Analyse des résultats estimés de paramètres

Erreur    Wald 95Limites
Paramètre      DF   Estimation   standard    de confiance %     Khi 2   Pr > Khi 2

Intercept        1     -3.5164     0.0554  -3.6249  -3.4078   4031.29       <.0001
ageconducteur    1      0.0168     0.0009   0.0150   0.0186    346.32       <.0001
Scale            0      0.6509     0.0000   0.6509   0.6509```

On notera que pour calculer le critère d’Akaike, ça n’est pas forcément trivial,

```%MACRO CALCUL_AIC_BIC(infomodel=, param=);
DATA _null_;
SET &infomodel.;
IF Criterion = "Log Likelihood" THEN CALL SYMPUT("Loglike", Value);
IF Criterion = "Deviance" THEN CALL SYMPUT("n_etoile", Df);
RUN;
DATA _null_;
SET &param.  end=fin;
RETAIN nb_df 0;
nb_df = nb_df + df;
IF fin THEN CALL SYMPUT("k", nb_df);
RUN;
DATA genmod_aic_bic;
SET &param.;
FORMAT Loglike 12.2 K 10. N 10. AIC_CALC 12.2 BIC_CALC 12.2;
Loglike = 0; K = 0; N = 0; AIC_CALC = 0; BIC_CALC = 0;
IF Parameter = "Intercept";
KEEP Loglike K N AIC_CALC BIC_CALC;
RUN;
DATA genmod_aic_bic;
SET genmod_aic_bic;
Loglike = &loglike.;
K = &k.;
N = %eval(&n_etoile. + &k.);
AIC_CALC = 2 * Loglike + 2 * K;
BIC_CALC = 2 * Loglike + K * log(N);
RUN;
PROC PRINT DATA = genmod_aic_bic;
RUN; QUIT;
%MEND CALCUL_AIC_BIC;
%CALCUL_AIC_BIC(infomodel=Genmod2_InfoModele, param=Genmod2_Param);```

# Qui peut m’aider à comprendre les sorties de SAS ?

Je m’étais promis que j’évoquerais une bizarrerie rencontrée avec SAS lors d’une formation…. Écrire ce billet permettra à ceux qui auraient des éléments d’explication de poster un commentaire.
Pour cela, comparons une régression logistique faite avec deux outils différents, sous SAS,

• avec la procédure logistique

Le code pour faire une régression logistique ressemble à ça

```PROC LOGISTIC DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / selection=stepwise;
RUN; QUIT;```

ce qui donne la sortie suivante (je passe l’introduction pour insister sur les coefficients)

```                                 The LOGISTIC Procedure

Analyse des estimations de la vraisemblance maximum
Erreur         Khi 2
Paramètre                    DF    Estimation         std       de Wald    Pr > Khi 2

Intercept                     1        1.7833      0.0676      696.9022        <.0001
sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001
age_soc      1_AGESOC_-60     1        0.4578      0.0667       47.1020        <.0001
fract_paiemt Annuel           1        0.6021      0.0997       36.4862        <.0001
fract_paiemt Mensuel          1       -0.5410      0.0842       41.2342        <.0001```
• avec la procédure genmod (car la régression logistique est un glm)

On peut faire exactement la même chose (théoriquement) en ajustement un modèle GLM,

```PROC GENMOD DATA=base_logistq;
FORMAT age_soc f2_ageso.;
CLASS sexe_soc age_soc fract_paiemt;
MODEL SPOCAM = sexe_soc age_soc fract_paiemt / dist = binomial;
RUN;```

et la sortie ressemble à ça

```                                  The GENMOD Procedure
Analyse des résultats estimés de paramètres

Erreur      Wald 95Limites
Paramètre                     DF   Estimation   standard      de confiance %       Khi 2
Intercept                      1       1.5073     0.1501     1.2131     1.8014    100.85
sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000       .
age_soc        1_AGESOC_-60    1       0.9156     0.1334     0.6542     1.1771     47.10
age_soc        Z_AGESOC_+60    0       0.0000     0.0000     0.0000     0.0000       .
fract_paiemt   Annuel          1       0.6634     0.1770     0.3165     1.0104     14.05
fract_paiemt   Mensuel         1      -0.4798     0.1510    -0.7759    -0.1838     10.09
fract_paiemt   Semestriel      0       0.0000     0.0000     0.0000     0.0000       .
Scale                          0       1.0000     0.0000     1.0000     1.0000```
• comparaison des deux sorties

Si on regarde l’impact du sexe par exemple, dans la première sortie on peut lire

`sexe_soc     Femme            1       -0.2429      0.0619       15.4237        <.0001`
alors que dans la seconde sortie, on a
```sexe_soc       Femme           1      -0.4859     0.1237    -0.7284    -0.2434     15.42
sexe_soc       Homme           0       0.0000     0.0000     0.0000     0.0000```

On dira ce qu’on veut, mais moi je trouve cette différence troublante…. Dans la seconde sortie, le coefficient vaut le double de l’autre….
Alors SAS semble s’y retrouver car si on lui demande d’afficher le score prédit pour un individu au hasard (le premier de la base par exemple), les prédictions sont très proches,

```                           fract_                                 proba1_       proba1_
Obs  sexe_soc   age_soc  paiemt      SPOCAM  proba1_logit
1      Homme          71  Annuel         0      0.10242637    0.10241302```

Si quelqu’un sait interpréter ce qui est fait avec cette procédure logistique (car R donne la même chose que la sortie GLM), je suis preneur…..