Méthode de Monte Carlo appliquée au biais des variables omises

Version imprimable

Bonjour, je suis nouveau sur le forum (et néophyte de R) et j'ai besoin d'aide concernant un DM que je suis censé rendre dans 1 mois (on ne peut pas me reprocher de ne pas m'y être pris à l'avance en tout cas :mrgreen:)

Le principe est simple en apparence : on nous pose un premier modèle du type Yi = β0 + β1 · Xi + β2 · Zi + ui ,
où β0 = 0.75, β1 = 0.50 et β2 = 0.25. On imagine qu'on pense pouvoir estimer β1 sans biais par le sous-modèle suivant : Yi = b0 + b1 · Xi + vi

Je dois écrire un code R qui montre que :
1. Quand Xi et Zi sont dé-corrélés, l'estimateur des moindres carrés ordinaires pour b1 est non biaisé et convergent.
2. Quand Xi et Zi sont positivement corrélés, l'estimateur des moindres carrés ordinaires pour b1 est biaisé vers le haut et divergent
3. Quand Xi et Zi sont négativement corrélés, l'estimateur des moindres carrés pour b1 est biaisé vers le bas et divergent

Par "montrer" on entend obtenir une distribution de l'estimateur des moindre carrés pour b1, calculer les statistiques pertinentes et tracer un graphique dans l'espace approprié

Voici où mon travail m'a mené jusqu'ici :
- je sais que, pour être dans une situation de variable aléatoire omise, X doit être corrélé avec la variable omise et la variable omise doit pouvoir expliquer Y
- j'ai trouvé une façon de générer des simulations d'expériences de Monte Carlo par le programme suivant :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 set.seed(1) n<-X; r<-X b0<-X;b1<-X;su<-X b0hat<-numeric(r) b1hat<-numeric(r) x1<-rnorm(n,0,X) for (j in 1:r) { print(j) u<-rnorm(n,0,su) y<-b0+b1*x1+u bhat<-coefficients(lm(y~x1)) b0hat[j]<-bhat["(Intercept)"] b1hat[j]<-bhat["x1"]
Voici le problème que je rencontre pour compléter le programme : je ne comprends pas la logique de l'exercice, parce que je suis supposé prendre en compte des cas généraux et des conditions sur des variables (donc des commandes du type if/then que j'ai déjà eu l'occasion d'utiliser sur Scilab me semblent appropriées, j'ai trouvé qu'on utilisait if/else sur R, mais impossible d'aller plus loin et de savoir quoi y mettre dedans). Mais avec la commande que je viens de vous montrer, on ne prend pas en compte de cas généraux mais un cas spécifique (certes généré aléatoirement, c'est le principe de la méthode de Monte Carlo) qui donne des résultats spécifiques... Je ne vois pas en quoi Xi et Zi peuvent varier dans une simulation de Monte Carlo qui donne des résultats figés...

Merci d'avance pour votre aide. Je demande surtout à comprendre, davantage qu'à avoir les solutions au DM sans que je ne comprenne rien à ce qui se passe. Je n'ai jamais utilisé R, mais il faut que j'apprenne alors je vous en appelle à l'aide ! :D

28/11/2019, 10h21
tototode

Bonjour,

A mon avis pour répondre à ces questions tu vas devoir simuler des données pour x1 et x2 tout d'abord sans que ces variables soient corrélées puis avec des corrélations positives et négatives que tu auras choisie arbitrairement. Tu peux aussi le faire pour différentes valeurs de corrélations. Une fois que tu as les valeurs pour chaque variable, tu les multiplies par tes coefficients (sans oublier b0) et tu calcules ta régression linéaire entre y et x1 (a voir s'il faut ajouter ou non de la variabilité sur y pour avoir une erreur résiduelle). Tu répètes ça au moins 999 fois pour avoir une distribution de b1 et tu peux faire des histogrammes (ou des courbes de densité) et voir ou se situe la "vraie" valeur par rapport à cette distribution.
Tu verras que dans le 1er cas, b1 est proche du centre de la distribution et que dans les deux autres cas, b1 est < à la moyenne des b1 estimés et que dans le dernier cas b1 est > à la moyenne des b1 estimés. Tu peux faire des tests, comme pour les tests de permutation ou un test de student.

cdlt

1 pièce(s) jointe(s)

Merci tototode pour ton aide :D

Je ne suis pas sûr d'avoir compris, mais de ce que j'ai compris dans ce que tu m'as dit, je vais devoir utiliser le code que j'ai donné mais pour le modèle d'origine avec Xi et Zi (alors que dans ce que j'ai mis dans mon premier message, je l'avais utilisé pour le sous-modèle) ? Comment le faire, et surtout comment faire pour contrôler la corrélation que les deux variables ont entre elles ? C'est le cœur de l'exercice et surtout le centre de mon incompréhension, puisque par définition une simulation Monte Carlo génère une situation aléatoire non? Comment je peux donc avoir un impact sur les variables?

Quand je cherche sur Internet un moyen d'imposer une corrélation a priori à un code de modèle de régression linéaire sur Internet, je ne trouve rien. La seule chose que je peux faire, c'est avoir la corrélation après le calcul du modèle, donc a posteriori. Sachant que le modèle est aléatoire, la corrélation l'est donc aussi. Je n'ai donc aucun contrôle dessus...

J'ai considéré les options que je connais : boucle while, boucle if/else, mais ce n'est pas exactement ce que l'on veut non ? Dans une boucle while, on répète une action tant qu'elle est satisfaite. Mais si elle n'est pas satisfaite, le programme s'arrête... ce n'est pas ce que l'on veut, on veut partir du principe que corr(x1,x2)=0 (pour la question 1) et pas continuer le code si et seulement si la condition est vérifiée. Pareil pour if/else.

J'ai considéré une autre façon de penser : modifier le modèle en lui-même (par exemple, l'échantillon). Par exemple, si je prends un échantillon de 250 personnes, j'ai une corrélation positive. Si j'ai un échantillon de 25 personnes, une corrélation négative.

Déjà, si j'essaie d'adapter le code de mon message introductif à un modèle à plusieurs variables aléatoires pour 250personnes, voici ce que j'obtiens :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 ## Expérience de Monte Carlo} set.seed(1) ## Taille de l'échantillon et nombre de répétitions n<-25000; r<-20 ## Les vraies valeurs du modèle b0<-0.75;b1<-0.5;b2<-0.25;su<-100 ## Initialiser les 3 vecteurs pour sauvegarder les résultats b0hat<-numeric(r) b1hat<-numeric(r) b2hat<-numeric(r) ## Générer X} x1<-rnorm(n,0,10) x2<-rnorm(n,0,10) ## Faire la boucle pour répéter l'opération r fois for (j in 1:r) { print(j) # Dessiner l'échantillon u<-rnorm(n,0,su) # Générer le modèle y<-b0+b1*x1+b2*x2+u # Estimer le modèle par MCO bhat<-coefficients(lm(y~x1*x2)) ## Enregistrer les résultats b0hat[j]<-bhat["(Intercept)"] b1hat[j]<-bhat["x1"] b2hat[j]<-bhat["x2"] }
J'ai ensuite fait ceci :
Code:

1 2 3 4 cor(x1,x2) plot(b1hat,b2hat) summary(b1hat)) summary(b2hat)
J'ai obtenu une corrélation de 0.04, ce qui est très faible. Le graphique donne ce que j'ai montré en bas de page.
On constate effectivement que les b1hat sont proches de 0.5 (sa vraie valeur) et b2hat de 0.25 (sa vraie valeur)

La commande summary donne ceci, respectivement pour b1hat et b2hat
Code:

1 2 3 4 5 Min. 1st Qu. Median Mean 3rd Qu. Max. -1.85639 0.05913 0.47651 0.48726 0.94473 3.15970 Min. 1st Qu. Median Mean 3rd Qu. Max. -1.5532 -0.1799 0.2540 0.2419 0.6708 1.9633
On a b1hat < b1 et b2hat < b2, ce qui est l'inverse de ce que l'on recherche puisque si Xi et Zi sont corrélés positivement, il y a biais vers le haut, autrement dit b1hat > b1 et b2hat > b2

Si je teste le même modèle avec 25 personnes, je tombe sur une corrélation négative avec un biais vers le haut cette fois-ci pour b1 (b1hat > b1) mais bien négatif comme on doit le trouver pour b2hat (b2hat < b2)

Edit: en choisissant un échantillon beaucoup plus grand (j'ai testé n=2500 et n=25000) les résultats redeviennent logiques. Quand n=25000, la corrélation est de 0.005 et b1hat>b1 et b2hat>b2. Quand n=2500, la corrélation est de -0.005 et b1hat<b1 et b2hat<b2. Sauf que c'est un peu facile de choisir l'échantillon qui nous arrange, comment faire ça de façon plus générale ? Et surtout, comment prendre en compte le cas où corr(x1,x2)=0 ?

En attendant, merci beaucoup pour ton aide.

16/12/2019, 21h07
tototode

Bonjour,

si tu n'as pas trouvé le moyen de fixer une corrélation a priori c'est que tu n'as pas très bien cherché :-)
Utilises la fonction mvrnorm du package MASS et tu pourras via l'argument Sigma spécifier une matrice de variance-covariance et donc des corrélations entre Zi et Xi et ensuite fixé des paramètres et faire tes simulations.

cdlt

Merci pour ce conseil tototode, malheureusement je n'ai pas réussi à utiliser cette commande malgré mes recherches (j'ai réussi à l'installer mais pas à l'appliquer à l'exercice).

J'ai trouvé un autre moyen pour imposer une corrélation aux variables que je génère :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 set.seed(1) n<-25000; r<-1000; b0<-0.75;b1<-0.5;b2<-0.25;su<-100;p<--0.9 b0hat<-numeric(r) b1hat<-numeric(r) b2hat<-numeric(r) x1<-rnorm(n,0,10) x2<-10*p*(x1-0)/10 + 0 + 10*rnorm(n,0,sqrt(1-p^2)) for (j in 1:r) { print(j) u<-rnorm(n,0,su) y<-b0+b1*x1+b2*x2+u bhat<-coefficients(lm(y~x1 + x2)) b0hat[j]<-bhat["(Intercept)"] b1hat[j]<-bhat["x1"] b2hat[j]<-bhat["x2"] } plot(b1hat,b2hat) summary(b1hat) cor(x1,x2)
Je trouve bel et bien, lorsque cor(x1,x2) s'affiche, une corrélation identique à celle que j'ai imposé dans p.

Le souci, c'est que je trouve le biais inverse à celui que je devrais trouver (quand cor(x1,x2)>0, je trouve un biais vers le bas et inversement, alors que ça devrait être le contraire).
Entre p = 0 et p = 0.5, le biais est vers le haut mais il diminue jusqu'à passer sous les b1hat < 0.5

En revanche, si je change le code et que je fais ceci :
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 ## Monte Carlo experiment r=20} set.seed(1) n<-25000; r<-1000; b0<-0.75;b1<-0.5;b2<-0.25;su<-100;p<-0.5 b0hat<-numeric(r) b1hat<-numeric(r) b2hat<-numeric(r) x1<-rnorm(n,0,10) x2<-10*p*(x1-0)/10 + 0 + 10*rnorm(n,0,sqrt(1-p^2)) for (j in 1:r) { print(j) u<-rnorm(n,0,su) y<-b0+b1*x1+b2*x2+u bhat<-lm(y~x1 + x2) b0hat[j]<-bhat["(Intercept)"] b1hat[j]<-bhat["x1"] b2hat[j]<-bhat["x2"] } plot(x1,x2) summary(bhat) cor(x1,x2)
J'obtiens, si p=0.5 par exemple (donc supérieur à 0)
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 Call: lm(formula = y ~ x1 + x2) Residuals: Min 1Q Median 3Q Max -407.05 -67.63 0.21 67.34 379.28 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.31854 0.63048 2.091 0.0365 * x1 0.55669 0.07278 7.649 2.09e-14 *** x2 0.16982 0.07200 2.359 0.0184 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 99.69 on 24997 degrees of freedom Multiple R-squared: 0.004358, Adjusted R-squared: 0.004279 F-statistic: 54.71 on 2 and 24997 DF, p-value: < 2.2e-16
Un biais positif pour b1hat (0.55 > 0.5) ce que l'on cherche !

Je ne sais pas d'où vient le souci. Est-ce que la technique que j'ai employée est mauvaise ?
Merci d'avance pour votre aide.

Voilà ce que j'aurai fait avec du monte-carlo. Après je n'ai pas mis du bruit sur y mais tu peux le faire, ça ne change rien au
Une autre possibilité serait de générer des données et de faire des tirages aléatoires dans le jeu de données crées (considéré comme la pop) et refaire des modèles sur ces jeux de données.
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 library(MASS) b0 <- 0.75 b1 <- 0.5 b2 <- 0.25 bs <- c(b0, b1, b2) bsm <- matrix(bs, ncol = 1) # sans corrélation sim1 <- replicate(1000, { tab <- cbind(1, mvrnorm(1000, rep(0,2), diag(1, 2))) y <- tab %*% bsm # y <- y + rnorm(1000, 0, 0.5) cfs <- lm.fit(tab[,1:2], y)$coef cfs } ) sim1 <- t(sim1) # avec corrélation positive matp <- matrix(c(1,0.5, 0.5, 1), 2) sim2 <- replicate(1000, { tab <- cbind(1, mvrnorm(1000, rep(0,2), matp)) y <- tab %*% bsm # y <- y + rnorm(1000, 0, 0.5) cfs <- lm.fit(tab[,1:2], y)$coef cfs } ) sim2 <- t(sim2) # avec corrélation négative matn <- matrix(c(1,-0.5, -0.5, 1), 2) sim3 <- replicate(1000, { tab <- cbind(1, mvrnorm(1000, rep(0,2), matn)) y <- tab %*% bsm # y <- y + rnorm(1000, 0, 0.5) cfs <- lm.fit(tab[,1:2], y)$coef cfs } ) sim3 <- t(sim3)
Pour chaque cas de figure tu peux représenter la densité des coefficients estimés et calculer les fréquences des coefficients > ou < à b1 selon les hypothèses.
cdlt