Simulation de maladie (Petit exemple avec r)

Dans une population la probabilité de loccurrence d'une maladie (Y) depend du sexe  (S) et de la presence d'une mutation (M). en plus on sait qe sexe et mutation sont des fcteurs independants,

P(S=M) = 0.5
P(M=1) = 0.25,

on note      M = 1 quand la maladie est presente sinon M=0
                  Y = 1 si l'individu est malade Y=0 sinon
                   S = 1 si ll'individu est un male S=0 sinon

On suppose que la facon dont  probabilite de la maladie est relie aux autres variable peut s'ecrire

Logit(P(Y=1/ S=m,M=1) ) = log(1.1) + log(1.5) * s + log(2) * m

logit(p) = $\quad \log { \frac { p }{ 1-p }  } $

Nous allons coder sur r et obtenir un echantillon aleatoir de 1000 000  (Si. Mi. Yi )
 Nous allons obtenir des graphiques et verifier les resultats theoriques avec ca


CODE

set.seed(100)


n = 1000000

# le vecteur (1,0) est cree pour nous permettre de tirer ces valeurs
# On cree la matrice M avec 3 colonnes et n lignes

x = c(1,0)
M =  matrix (nrow = n, ncol = 3)






# on genere les valeurs de S  et m a l'aide la loi binomiale


s = rbinom(n,1,0.5)
m = rbinom(n,1,0.25)
#sachant ces 2 valeur on genere p et y


logit = log(1.1) + log(1.5) * s + log(2) * m

# p est la probabilite d etre malade sachant le sexe et m

p = exp(logit) / (1 + exp(logit))

y =  rbinom(n,1,p)


# bij  on cherce les lignes de M dont s = i et m = j
M= as.matrix(cbind(s,m,y))
b11 =which( M[,1] ==1 & M[,2] == 1 )
b10 =which( M[,1] ==1 & M[,2] == 0 )
b00 =which( M[,1] ==0 & M[,2] == 0 )
b01 =which( M[,1] ==0 & M[,2] == 1 )


# Mij sous matrice de M dont s = i et m = j

M11 = as.matrix(M [b11,])
M10 =as.matrix(M [b10,])
M00 =as.matrix(M [b00,])
M01 =as.matrix(M [b01,])

# on fait la somme des colonnes de nos sous matrices

nb11 = as.matrix(colSums(M11))
nb10 = as.matrix(colSums(M10))
nb00 = as.matrix(colSums(M00))
nb01 = as.matrix(colSums(M01))

# nbobsij donne le nombre de non malade dans chaque Mij

nbObs11 = length(M11[ ,3]) - nb11
nbObs10 = length(M10[ ,3]) - nb10
nbObs00 = length(M00[ ,3]) - nb00
nbObs01 = length(M01[ ,3]) - nb01

# proba est un vecteur qui donne la probabilite de tombe mlade ou pas sachant le couple s et m


proba = c(nb11[3,] /125710 ,nb10[3,] /374533,nb00[3,] /374804,nb01[3,] /124953,
          nbObs11[3,1]  /125710, nbObs10[3,1]  /374533,nbObs00[3,1]  /374804,nbObs01[3,1]  /124953)

#   SiMj est le label y = 1 / s=i,m=j, no est y = 0

classe = c("S1M1","S1M0","S2M0","S2M1","noS1M1","noS1M0","noS2M0","noS2M1")

barplot(proba,names.arg = classe)

proba
        y         y         y         y         y         y         y         y
0.7658818 0.6223964 0.5245355 0.6868823 0.2341182 0.3776036 0.4754645 0.3131177