Transcript PROC MIXED
Weekend Workshop I PROC MIXED Random or Fixed ? RANDOM FIXED Levels: Selected at random from infinite population Finite number of possibilities Another Experiment Different selections from same population Same Levels Estimate variance components Compare means All levels in population Only levels used in the experiment. Goal Inference Twins: One gets SAS training method 1, the other gets method 2 Response Y = programming times PROC MIXED Model Model 38.25 33.00 28.75 30.00 46.50 55.25 Y 1 1 1 1 1 1 = ; Random 0 1 1 0 0 1 M1 1 0 M 2 1 0 0 1 Xb 1 0 0 Variance of g is G = 2 0 1 0 F 0 0 1 1 1 0 0 0 0 ; Repeated ; 0 0 e1 e 0 0 2 F1 e3 1 0 F 2 1 0 e4 F3 e5 0 1 0 1 e6 Zg ,Variance of e is R = + e 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 2 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 PROC MIXED DATA=TWINS; CLASS FAMILY METHOD; MODEL TIME = METHOD; * fixed; RANDOM FAMILY; *<- family ~ N(0, 2F) ; Covariance Parameter Estimates Cov Parm family Residual Estimate 21.2184 40.8338 Type 3 Tests of Fixed Effects Effect method Num DF 1 Den DF F Value 19 9.60 Pr > F 0.0059 Intraclass correlation (related to heritability) 2F /(2F + 2) Estimated as 21.2/62 or about 1/3. Q: Why not usual (Pearson) correlation? Demo Get_Twins.sas Twins_MIXED.sas BLUP Yij = + Fi + eij Di = Family mean – Fi + ei. best estimate of Fi = ? Variance of (Fi – b Di) is (1-b)22F + b2 2/2 Use b = 2F /(2F + 2/2) Estimate: b = 21.2/(21.2 + 40.8/2) = 0.510 Overall mean + 0.510(Family i mean – Overall mean) PROC GLM DATA=TWINS; CLASS FAMILY METHOD; MODEL TIME = FAMILY METHOD; LSMEANS FAMILY; PROC MIXED DATA=TWINS; CLASS FAMILY METHOD; MODEL TIME = METHOD; RANDOM FAMILY; ESTIMATE "1 " intercept 1 | family 1; ESTIMATE "2 " intercept 1 | family 0 1; MEANS (GLM) and BLUPs (MIXED) Demo Twins_BLUP.sas Twins_TEST.sas ML Estimation Search over all (fixed and random) parameters Estimates of variances biased low! REML Estimation (1) Regress out fixed effects (2) Maximze likelihood of residuals (mean known: 0) (3) Variance estimates less biased (unbiased in some simple cases) Unbalanced Data SUBJ Ear plug A I 25 (L) 19 (L) C 22 (R) D E 29 (R) 8 (R) 7 (L) II III B 7 (R) F G 16 (R) 25 (L) 23 (L) 16 (R) 14 (L) 12 (L) I vs. III free of subject effects for red data. Misses info in other data. 24 (R) proc glm; class plug worker; model loss = worker plug; Random Worker; Estimate "I vs III - GLM" Plug -1 0 1; run; proc mixed; class plug worker; model Loss=Plug; Random Worker; Estimate "I vs III - Mixed" Plug -1 0 1; run; Covariance Parameter Estimates Cov Parm Estimate worker Residual 37.578 6.1674 GLM Source worker plug DF 6 2 Parameter I vs III - GLM Type III SS 451.9062500 62.6562500 Estimate -4.8125 F Value Pr > F 12.21 0.0074 5.08 0.0625 Standard Error t Value Pr > |t| 1.9635 -2.45 0.0579 Type 3 Tests of Fixed Effects Num Den Effect DF DF F Value Pr > F plug 2 5 5.79 0.0499 Estimates Standard Label Estimate Error DF t Value Pr > |t| I vs III - Mixed -5.2448 1.9347 5 -2.71 0.0422 Demo Earplugs.sas SPLIT PLOT six dishes / aquarium one plant / dish soil x variety combinations 4 Aquariums, 2 aerated 2 not ANOVA Soil 1 1 2 2 3 3 Variety 1 2 1 2 1 2 Source Air Error A V S VA VS AS AVS Error B Compare Air to No Air within soil 1 Variance of this contrast is hard to figure out: (1/3)[MS(A)+2 MS(B)] Need Satterthwaite df AUTOMATIC IN MIXED!!! PROC MIXED; CLASS VAR AQUARIUM SOIL AIR; MODEL YIELD = AIR SOIL VAR SOIL*VAR AIR*SOIL AIR*VAR AIR*SOIL*VAR / DDFM=SATTERTHWAITE; RANDOM AQUARIUM(AIR); ESTIMATE "SOIL 1: AIR EFFECT" AIR -1 1 AIR*SOIL -1 1 0 0 0 0; RUN; Covariance Parameter Estimates Cov Parm Estimate AQUARIUM(AIR) 2.1833 Residual 7.7333 Type 3 Tests of Fixed Effects Num Den DF DF F Value Pr > F AIR 1 2 16.20 0.0565 SOIL 2 10 7.87 0.0088 VAR 1 10 24.91 0.0005 VAR*SOIL 2 10 0.04 0.9631 SOIL*AIR 2 10 1.08 0.3752 VAR*AIR 1 10 4.22 0.0669 VAR*SOIL*AIR 2 10 0.23 0.7973 Effect Standard Label SOIL 1: AIR EFFECT Estimate Error DF t Value Pr > |t| 5.2500 2.4597 5.47 2.13 0.0812 Demo Aquarium.sas Random Coefficient Models the basic idea mistakes Dave Average programmer a0 + b0 t Program writing time Line for individual j: (a0 + aj) + ( b0 + bj )t 0 A2 AB aj ~ N , 2 b j 0 AB B Hierarchial Models (1) Same as split plot - almost (2) Whole and split level continuous predictor variables (typically) Yij = [a0 + a1 pHi + b0 Nij + b1 pHi Nij] + [ai* +bi* Nij+eij] random Aquariumfixed level (level i): pH (1) i (2) Dish level: Soil nitrogen test (Nij) YPROC + bi NijDATA +eij = UNDERWATER; ij = ai MIXED MODEL = Ni +Pai*N*P; (3) Idea: ai =GROWTH a0 + a1 pH * RANDOM / SUBJECT = TANK TYPE=UN; bi = INTERCEPT b0 + b1 pHi N+ b i Yij = a ai0++bai 1NpH + ai* + (b bi 0N+ij+e b1ijpHi + bi* ) Nij+eij ij+e i ij aquarium N 1 1 1 1 1 1 2 2 2 2 2 2 3 3 3 3 3 3 4 4 4 4 4 4 2.21 1.25 4.36 7.14 8.61 6.53 6.58 3.12 5.28 1.09 4.83 9.61 7.99 7.79 8.32 2.53 6.85 4.73 0.95 2.00 9.99 0.23 0.13 1.17 pH growth 5.5 5.5 5.5 5.5 5.5 5.5 4.7 4.7 4.7 4.7 4.7 4.7 4.2 4.2 4.2 4.2 4.2 4.2 5.1 5.1 5.1 5.1 5.1 5.1 27.05 25.92 30.09 33.66 36.13 33.00 35.72 31.17 34.35 28.34 33.56 40.25 47.04 46.56 48.27 34.20 44.59 39.29 24.94 27.33 43.84 23.54 23.56 25.68 Num Effect DF N 1 pH 1 N*pH 1 Den DF 2 2.05 2 F Value 3.50 6.76 1.31 Pr > F 0.2018 0.1186 0.3702 Num Effect DF N 1 pH 1 Den DF 3 2.03 F Value 50.19 14.68 Pr > F 0.0058 0.0603 p Cov Parm UN(1,1) UN(2,1) UN(2,2) Residual pH N Estimate 1.8976 -0.5563 0.2596 0.0286 Demo Hierarchial.sas Next: Repeated Measures Notes in pdf from NCSU experimental design class (ST 711) Demo SURGERY.sas