*********************************************************** * STT-4100/STT-7230 * Planification des expériences * * Exemple du chapitre 4: Analyse de la covariance * (inspiré de l'exemple 5.3, SAS System for Mixed Models, * Seulement les traitements alimentaires 1 et 3) ***********************************************************; **************************** * Lecture des données * ****************************; data bouffe; input obs trt cov rep; label trt="Diète"; label cov="Poids initial"; label rep="Gain de poids"; drop obs; cards; 1 1 338 1.03 2 1 403 1.31 3 1 394 1.59 4 1 499 2.09 5 1 371 1.66 6 1 395 1.42 7 1 414 1.41 8 1 315 0.18 9 3 444 1.82 10 3 450 2.13 11 3 482 2.33 12 3 391 2.21 13 3 486 2.65 14 3 316 1.58 15 3 309 1.08 16 3 308 0.76 ; proc print data=bouffe; run; proc means data=bouffe mean std printall; class trt; var rep cov ; run; proc sgplot data=bouffe; vbox rep /category=trt; run; *********************************************** * i) Modèle 1 à 2 droites différentes * * Les pentes sont-elles nulles simultanément? * ***********************************************; * Regardons le test sur cov*trt ; * On rejette l'hypothèse que les 2 pentes sont nulles; proc glm data=bouffe; class trt ; model rep = trt cov*trt/noint solution; run;quit; ********************************************* * ii) Modèle 1 à 2 droites différentes * * Les pentes sont-elles égales entre elles? * *********************************************; * On a ajouté cov dans le modèle, une pente commune. Le terme cov*trt représente la déviation de chaque pente par rapport à la pente commune. On veut tester si ces déviations sont toutes nulles. On regarde le test sur cov*trt; * On ne rejette pas l'hypothèse que les 2 pentes sont égales; proc glm data=bouffe; class trt ; model rep = trt cov cov*trt/ clparm solution; lsmeans trt / cl stderr at cov = 330; estimate "moy aj. trt1" intercept 1 trt 1 0 cov 330 cov*trt 330 0; estimate "moy aj. trt3" intercept 1 trt 0 1 cov 330 cov*trt 0 330; estimate "moy aj. trt1-trt3" intercept 0 trt 1 -1 cov 0 cov*trt 330 -330; run;quit; *********************************************************************************** * iii) Modèle 2 à 2 droites de même pente, mais d'ordonnées à l'origine différentes * * La pente commune est-elle nulle? * * Les ordonnées à l'origine sont-elles égales entre elles? * ***********************************************************************************; * On rejette l'hypothèse que la pente commune est nulle (test sur cov). Donc la covariable réduit sensiblement la variabilité de la réponse; * On rejette l'hypothèse que les 2 ordonnées à l'origine sont égales (test sur trt). Donc les trt ont des moyennes différentes. Cette différence est considérée constante pour toutes les valeurs de la covariable, car les droites ajustées sont parallèles. Par contre la valeur estimée de la moyenne (et son erreur-type) varie selon la valeur de la covariable (voir les 3 énoncés lsmeans ci-dessous); proc glm data=bouffe; class trt ; model rep= trt cov/solution; lsmeans trt/stderr pdiff; *(qui correspond à la valeur moyenne de la covariable); lsmeans trt/stderr pdiff at cov=330; lsmeans trt/stderr pdiff at cov=460; estimate "moy.aj. trt1,cov=330" intercept 1 trt 1 0 cov 330/e; estimate "moy.aj. trt1-trt3" intercept 0 trt 1 -1 /e; output out=out r=resid p=pred; run;quit; ********************************************* * Si on n'a pas la covariable: * * Anova à un facteur * *********************************************; proc glm data=bouffe ; class trt; model rep=trt; run;quit; ***********************************************************************************************; ***********************************************************************************************; *********************** * ANALYSE DES RÉSIDUS * ***********************; * À chaque fois qu'un modèle est ajusté, il est nécessaire de faire l'analyse des résidus pour s'assurer que les tests utilisés pour conclure et passer à un autre modèle sont valides. On peut ajouter la ligne output ci-dessous dans Proc glm (après l'énoncé model) et utiliser les procédures univariate et gplot pour analyser les résidus; * output out=out r=resid p=pred; proc univariate plot normal data=out; var resid; histogram resid; run; title "Vérification de l'homoscédasticité"; proc sgplot data=out; scatter y=resid x=pred; refline 0; yaxis label="Résidus du modèle 2"; xaxis label="Valeurs prédites du modèle 2"; run; title; ***********************************************************************************************; ***********************************************************************************************;