/******************************************************************* ******************************************************************** Lec 17 MARCH 19 ******************************************************************* ********************************************************************/ /******************************************************************* CHAPTER 10, EXAMPLE 1 Illustration of - fitting both a full random coefficient model as in Chapter 9 and a and modified random coefficient model with intercepts random and slopes fixed for the dental data using PROC MIXED. - obtaining BLUPs of random effects and random intercepts (and slopes where applicable) for both models. The model for each child is assumed to be a straight line. The intercepts and slopes may have different means depending on gender. However, for the modified model, slopes are taken to be the SAME for all children within each gender. This assumption is probably not true, but is made for illustrative purposes to show how such a model may be specified in PROC MIXED. For both models, we take D to be common to both genders and take Ri = sigma^2_G I for girls and Ri = sigma^2_B for boys using the REPEATED statement. We use the RANDOM statement to specify how random effects enter the model AND to ask for the BLUPs of the bi to be printed in each case. We also use an option in the MODEL statement to ask for the BLUPs of the individual means at each time point for each child. *******************************************************************/ options ls=80 ps=59 nodate; run; /****************************************************************** Read in the data set (See Example 1 of Chapter 4) *******************************************************************/ data dent1; infile "dental.txt"; input obsno child age distance gender; time=age; run; /******************************************************************* Use PROC MIXED to fit the two linear mixed effects models. For all of the fits, we use usual normal ML rather than REML (the default). We call PROC MIXED twice to fit each model, for reasons described below. In all cases, we use the usual parameterization for the mean model. Here, we use the syntax for versions 7 and higher of SAS for outputting calculations to data sets from PROC MIXED. In the first call to PROC MIXED: To fit the full random coefficient model, we must specify that both intercept and slope are random in the RANDOM statement. To fit the modified model where slopes are taken to be constant across all children within a gender, we specify only that intercept is random in the RANDOM statement. *******************************************************************/ * MODEL (i) -- full random coefficient model; * Call to PROC MIXED to get the printed results; title 'FULL RANDOM COEFFICIENT MODEL WITH BOTH'; title2 'INTERCEPTS AND SLOPES RANDOM FOR EACH GENDER'; proc mixed method=ml data=dent1; class gender child; * OUTPRED=dataset option in the MODEL statement requests that individual mean ; * at each time point in the data set be printed; * SOLUTION option in the MODEL statement tells SAS to print estimates of the fixed effects; model distance = gender gender*age / noint solution outpred=pdata; * SOLUTION option in the RANDOM statement tells SAS to print BLUP of the random effects; random intercept age / type=un subject=child solution; repeated / group=gender subject=child; run; proc print data=pdata; run; /******************************************************************* ******************************************************************** Cov Parm Subject Group Estimate UN(1,1) child 3.1978 UN(2,1) child -0.1103 UN(2,2) child 0.01976 Residual child gender 0 0.4449 Residual child gender 1 2.6294 Fit Statistics -2 Log Likelihood 406.0 AIC (smaller is better) 424.0 AICC (smaller is better) 425.9 BIC (smaller is better) 435.7 Null Model Likelihood Ratio Test DF Chi-Square Pr > ChiSq 4 72.20 <.0001 Solution for Fixed Effects Standard Effect gender Estimate Error DF t Value Pr > |t| gender 0 17.3727 0.7386 54 23.52 <.0001 gender 1 16.3406 1.1114 54 14.70 <.0001 age*gender 0 0.4795 0.06180 54 7.76 <.0001 age*gender 1 0.7844 0.09722 54 8.07 <.0001 Solution for Random Effects Std Err Effect child Estimate Pred DF t Value Pr > |t| Intercept 1 -0.4853 1.1744 54 -0.41 0.6811 age 1 -0.06820 0.1017 54 -0.67 0.5052 Intercept 2 -1.1922 1.1744 54 -1.02 0.3146 age 2 0.1420 0.1017 54 1.40 0.1683 Intercept 3 -0.8535 1.1744 54 -0.73 0.4705 age 3 0.1773 0.1017 54 1.74 0.0869 Intercept 4 1.7024 1.1744 54 1.45 0.1530 age 4 0.04017 0.1017 54 0.40 0.6943 Intercept 5 0.9136 1.1744 54 0.78 0.4400 ----------------------------------------------------------------------------- Type 3 Tests of Fixed Effects Num Den Effect DF DF F Value Pr > F gender 2 54 384.72 <.0001 age*gender 2 54 62.66 <.0001 # We use the OUTPRED=dataset option in the MODEL statement. # The data set contains values : Xi beta.hat + Zi b.hat.i # # This requests that the (approximate) Best Linear Unbiased Predictors # for the individual means at each time point in the data set for # each child be put in dataset (along with the original data for comparison). # These may be printed with a print statement, as shown. FULL RANDOM COEFFICIENT MODEL WITH BOTH 5 INTERCEPTS AND SLOPES RANDOM FOR EACH GENDER S t d d i E s g r o c t e r A L U R b h a n t P P l o p e O s i a n d i r r p w p s b n l g c e m e e D h e e i s o d e e r e d d F a r r d 1 1 1 8 21.0 0 8 20.1783 0.43711 54 0.05 19.3019 21.0546 0.82175 2 2 1 10 20.0 0 10 21.0009 0.33796 54 0.05 20.3234 21.6785 -1.00095 3 3 1 12 21.5 0 12 21.8236 0.34908 54 0.05 21.1238 22.5235 -0.32365 4 4 1 14 23.0 0 14 22.6463 0.46259 54 0.05 21.7189 23.5738 0.35366 5 5 2 8 21.0 0 8 21.1527 0.43711 54 0.05 20.2763 22.0290 -0.15266 6 6 2 10 21.5 0 10 22.3957 0.33796 54 0.05 21.7181 23.0733 -0.89570 7 7 2 12 24.0 0 12 23.6387 0.34908 54 0.05 22.9389 24.3386 0.36126 8 8 2 14 25.5 0 14 24.8818 0.46259 54 0.05 23.9543 25.8092 0.61822 ******************************************************************* ********************************************************************/ /******************************************************************* In this PROC MIXED, we use the ODS statement to produce data sets containing the fixed effects estimates and the BLUPs for the random effects. We use the Output Delivery System in SAS, or ODS. The first ODS call with "listing exclude" suppresses printing of the fixed and random effects. The output data sets FIXED1 and RANDOM1 we ask PROC MIXED to create in the ODS statements contain the estimated fixed effects (betahats) and random effects (the BLUPs of bis), respectively. We now combine these into a single data set in order to compute the BLUPs of the individual betais. This is accomplished by manipulating the output data sets and then merging them. *******************************************************************/ * Call to PROC MIXED to produce the output data sets; proc mixed method=ml data=dent1; class gender child; * SOLUTION option in MODEL statement = fixed effects, beta; model distance = gender gender*age / noint solution; * SOLUTION option in RANDOM statement = random effects, b.i; random intercept age / type=un subject=child solution ; repeated / group=gender subject=child; ods listing exclude SolutionF; * SolutionF is a key word recognized by PROC MIXED as identifying this data set; ods output SolutionF=fixed1; * SolutionR is the key word identifying this data set. * For each child, there is a separate row in the file for the intercept BLUP and the slope BLUP (b0i and b1i); ods listing exclude SolutionR; ods output SolutionR=rand1; run; /******************************************************************* ******************************************************************** Type 3 Tests of Fixed Effects Num Den Effect DF DF F Value Pr > F gender 2 54 384.72 <.0001 age*gender 2 54 62.66 <.0001 ******************************************************************* ********************************************************************/ data fixed1; set fixed1; keep gender effect estimate; run; title3 'FIXED EFFECTS OUTPUT DATA SET'; proc print data=fixed1; run; /******************************************************************* ******************************************************************** FULL RANDOM COEFFICIENT MODEL WITH BOTH 10 INTERCEPTS AND SLOPES RANDOM FOR EACH GENDER FIXED EFFECTS OUTPUT DATA SET Obs Effect gender Estimate 1 gender 0 17.3727 2 gender 1 16.3406 3 age*gender 0 0.4795 4 age*gender 1 0.7844 ******************************************************************* ********************************************************************/ proc sort data=fixed1; by gender; run; data fixed12; set fixed1; by gender; retain fixint fixslope; if effect='gender' then fixint=estimate; if effect='age*gender' then fixslope=estimate; if last.gender then do; output; fixint=.; fixslope=.; end; drop effect estimate; run; title3 'RECONFIGURED FIXED EFFECTS DATA SET'; proc print data=fixed12; run; /******************************************************************* ******************************************************************** FULL RANDOM COEFFICIENT MODEL WITH BOTH 11 INTERCEPTS AND SLOPES RANDOM FOR EACH GENDER RECONFIGURED FIXED EFFECTS DATA SET Obs gender fixint fixslope 1 0 17.3727 0.47955 2 1 16.3406 0.78438 ******************************************************************* ********************************************************************/ data rand1; set rand1; gender=1; if child<12 then gender=0; keep child gender effect estimate; run; title3 'RANDOM EFFECTS OUTPUT DATA SET'; proc print data=rand1; run; /******************************************************************* ******************************************************************** # Estimate bi0 and bi1 for all the subjects FULL RANDOM COEFFICIENT MODEL WITH BOTH 12 INTERCEPTS AND SLOPES RANDOM FOR EACH GENDER RANDOM EFFECTS OUTPUT DATA SET Obs Effect child Estimate gender 1 Intercept 1 -0.4853 0 2 age 1 -0.06820 0 3 Intercept 2 -1.1922 0 4 age 2 0.1420 0 5 Intercept 3 -0.8535 0 6 age 3 0.1773 0 7 Intercept 4 1.7024 0 8 age 4 0.04017 0 9 Intercept 5 0.9136 0 10 age 5 -0.08680 0 11 Intercept 6 -0.6740 0 ******************************************************************* ********************************************************************/ proc sort data=rand1; by child; run; data rand12; set rand1; by child; retain ranint ranslope; if effect='Intercept' then ranint=estimate; if effect='age' then ranslope=estimate; if last.child then do; output; ranint=.; ranslope=.; end; drop effect estimate; run; proc sort data=rand12; by gender child; run; title3 'RECONFIGURED RANDOM EFFECTS DATA SET'; proc print data=rand12; run; /******************************************************************* ******************************************************************** # re - organizing of the random intercept and random slope FULL RANDOM COEFFICIENT MODEL WITH BOTH 14 INTERCEPTS AND SLOPES RANDOM FOR EACH GENDER RECONFIGURED RANDOM EFFECTS DATA SET Obs child gender ranint ranslope 1 1 0 -0.48526 -0.06820 2 2 0 -1.19224 0.14198 3 3 0 -0.85346 0.17726 4 4 0 1.70243 0.04017 5 5 0 0.91363 -0.08680 6 6 0 -0.67403 -0.07292 7 7 0 -0.05461 0.03641 8 8 0 1.93498 -0.11486 9 9 0 -0.21898 -0.11515 10 10 0 -2.99738 -0.09085 11 11 0 1.92494 0.15297 12 12 1 1.34688 0.08788 13 13 1 -0.86755 -0.04068 14 14 1 -0.35750 -0.02176 15 15 1 1.59462 -0.02772 16 16 1 -1.15811 -0.04153 17 17 1 0.89718 0.02260 ******************************************************************* ********************************************************************/ data both1; merge fixed12 rand12; by gender; beta0i=fixint+ranint; beta1i=fixslope+ranslope; run; title3 'RANDOM INTERCEPTS AND SLOPES'; proc print data=both1; run; /******************************************************************* ******************************************************************** # Merging all the ESTIMATES; and recovering the subject specific predicted trajectory FULL RANDOM COEFFICIENT MODEL WITH BOTH 15 INTERCEPTS AND SLOPES RANDOM FOR EACH GENDER RANDOM INTERCEPTS AND SLOPES Obs gender fixint fixslope child ranint ranslope beta0i beta1i 1 0 17.3727 0.47955 1 -0.48526 -0.06820 16.8875 0.41135 2 0 17.3727 0.47955 2 -1.19224 0.14198 16.1805 0.62152 3 0 17.3727 0.47955 3 -0.85346 0.17726 16.5193 0.65681 4 0 17.3727 0.47955 4 1.70243 0.04017 19.0752 0.51971 5 0 17.3727 0.47955 5 0.91363 -0.08680 18.2864 0.39274 6 0 17.3727 0.47955 6 -0.67403 -0.07292 16.6987 0.40662 7 0 17.3727 0.47955 7 -0.05461 0.03641 17.3181 0.51595 8 0 17.3727 0.47955 8 1.93498 -0.11486 19.3077 0.36469 9 0 17.3727 0.47955 9 -0.21898 -0.11515 17.1537 0.36440 10 0 17.3727 0.47955 10 -2.99738 -0.09085 14.3753 0.38869 11 0 17.3727 0.47955 11 1.92494 0.15297 19.2977 0.63251 12 1 16.3406 0.78438 12 1.34688 0.08788 17.6875 0.87225 13 1 16.3406 0.78438 13 -0.86755 -0.04068 15.4731 0.74369 14 1 16.3406 0.78438 14 -0.35750 -0.02176 15.9831 0.76262 15 1 16.3406 0.78438 15 1.59462 -0.02772 17.9352 0.75665 16 1 16.3406 0.78438 16 -1.15811 -0.04153 15.1825 0.74285 17 1 16.3406 0.78438 17 0.89718 0.02260 17.2378 0.80697 18 1 16.3406 0.78438 18 -0.68894 -0.02853 15.6517 0.75584 19 1 16.3406 0.78438 19 -0.14433 -0.07348 16.1963 0.71090 20 1 16.3406 0.78438 20 -0.12730 0.02544 16.2133 0.80981 21 1 16.3406 0.78438 21 2.53489 0.10877 18.8755 0.89315 22 1 16.3406 0.78438 22 -0.22609 -0.08535 16.1145 0.69903 23 1 16.3406 0.78438 23 -0.63735 0.00651 15.7033 0.79088 24 1 16.3406 0.78438 24 -1.70079 0.11392 14.6398 0.89830 25 1 16.3406 0.78438 25 0.23870 -0.03166 16.5793 0.75272 26 1 16.3406 0.78438 26 0.11799 0.06104 16.4586 0.84542 27 1 16.3406 0.78438 27 -0.82229 -0.07545 15.5183 0.70893 ******************************************************************* ********************************************************************/ * MODEL (ii) -- RANDOM INTERCEPT FIXED SLOPE, common within each gender; * Call to PROC MIXED to get the printed results; * To save space, we do not print the predicted values; title 'MODIFIED RANDOM COEFFICIENT MODEL WITH'; title2 'INTERCEPTS RANDOM, SLOPES FIXED'; proc mixed method=ml data=dent1; class gender child; model distance = gender gender*age / noint solution ; random intercept / type=un subject=child solution; repeated / group=gender subject=child; run; /******************************************************************* ******************************************************************** Covariance Parameter Estimates Cov Parm Subject Group Estimate UN(1,1) child 3.1405 Residual child gender 0 0.5920 Residual child gender 1 2.7286 Fit Statistics -2 Log Likelihood 409.4 AIC (smaller is better) 423.4 AICC (smaller is better) 424.5 BIC (smaller is better) 432.4 Null Model Likelihood Ratio Test DF Chi-Square Pr > ChiSq 2 68.89 <.0001 Solution for Fixed Effects Standard Effect gender Estimate Error DF t Value Pr > |t| gender 0 17.3727 0.7903 79 21.98 <.0001 gender 1 16.3406 1.1272 79 14.50 <.0001 age*gender 0 0.4795 0.05187 79 9.24 <.0001 age*gender 1 0.7844 0.09234 79 8.49 <.0001 Solution for Random Effects Std Err Effect child Estimate Pred DF t Value Pr > |t| Intercept 1 -1.2154 0.6434 79 -1.89 0.0626 Intercept 2 0.3364 0.6434 79 0.52 0.6025 Intercept 3 1.0527 0.6434 79 1.64 0.1058 Intercept 4 2.1270 0.6434 79 3.31 0.0014 Intercept 5 -0.02170 0.6434 79 -0.03 0.9732 Intercept 6 -1.4542 0.6434 79 -2.26 0.0266 Intercept 7 0.3364 0.6434 79 0.52 0.6025 Intercept 8 0.6945 0.6434 79 1.08 0.2837 Intercept 9 -1.4542 0.6434 79 -2.26 0.0266 Intercept 10 -3.9611 0.6434 79 -6.16 <.0001 Intercept 11 3.5595 0.6434 79 5.53 <.0001 ---------------------------------------------------------------------------- Type 3 Tests of Fixed Effects Num Den Effect DF DF F Value Pr > F gender 2 79 346.69 <.0001 age*gender 2 79 78.81 <.0001 ******************************************************************* ********************************************************************/ * Call to PROC MIXED to get the output data sets; proc mixed method=ml data=dent1; class gender child; model distance = gender gender*age / noint solution; random intercept / type=un subject=child solution; repeated / group=gender subject=child; ods listing exclude SolutionF; ods output SolutionF=fixed2; ods listing exclude SolutionR; ods output SolutionR=rand2; run; /******************************************************************* ******************************************************************** Covariance Parameter Estimates Cov Parm Subject Group Estimate UN(1,1) child 3.1405 Residual child gender 0 0.5920 Residual child gender 1 2.7286 Fit Statistics -2 Log Likelihood 409.4 AIC (smaller is better) 423.4 AICC (smaller is better) 424.5 BIC (smaller is better) 432.4 Null Model Likelihood Ratio Test DF Chi-Square Pr > ChiSq 2 68.89 <.0001 Type 3 Tests of Fixed Effects Num Den Effect DF DF F Value Pr > F gender 2 79 346.69 <.0001 age*gender 2 79 78.81 <.0001 # the ODS statement requests tables to be produced for each of fixed and random effects ******************************************************************* ********************************************************************/ data fixed2; set fixed2; keep gender effect estimate; run; title3 'FIXED EFFECTS OUTPUT DATA SET'; proc print data=fixed2; run; proc sort data=fixed2; by gender; run; /******************************************************************* ******************************************************************** MODIFIED RANDOM COEFFICIENT MODEL WITH 21 INTERCEPTS RANDOM, SLOPES FIXED FIXED EFFECTS OUTPUT DATA SET Obs Effect gender Estimate 1 gender 0 17.3727 2 gender 1 16.3406 3 age*gender 0 0.4795 4 age*gender 1 0.7844 ******************************************************************* ********************************************************************/ data fixed22; set fixed2; by gender; retain fixint fixslope; if effect='gender' then fixint=estimate; if effect='age*gender' then fixslope=estimate; if last.gender then do; output; fixint=.; fixslope=.; end; drop effect estimate; run; title3 'RECONFIGURED FIXED EFFECTS DATA SET'; proc print data=fixed22; run; /******************************************************************* ******************************************************************** MODIFIED RANDOM COEFFICIENT MODEL WITH 22 INTERCEPTS RANDOM, SLOPES FIXED RECONFIGURED FIXED EFFECTS DATA SET Obs gender fixint fixslope 1 0 17.3727 0.47955 2 1 16.3406 0.78438 ******************************************************************* ********************************************************************/ data rand2; set rand2; gender=1; if child<12 then gender=0; keep child gender effect estimate; run; title3 'RANDOM EFFECTS OUTPUT DATA SET'; proc print data=rand2; run; /******************************************************************* ******************************************************************** MODIFIED RANDOM COEFFICIENT MODEL WITH 23 INTERCEPTS RANDOM, SLOPES FIXED RANDOM EFFECTS OUTPUT DATA SET Obs Effect child Estimate gender 1 Intercept 1 -1.2154 0 2 Intercept 2 0.3364 0 3 Intercept 3 1.0527 0 4 Intercept 4 2.1270 0 5 Intercept 5 -0.02170 0 6 Intercept 6 -1.4542 0 7 Intercept 7 0.3364 0 8 Intercept 8 0.6945 0 9 Intercept 9 -1.4542 0 10 Intercept 10 -3.9611 0 11 Intercept 11 3.5595 0 ******************************************************************* ********************************************************************/ proc sort data=rand2; by child; run; data rand22; set rand2; by child; retain ranint ranslope; if effect='Intercept' then ranint=estimate; if last.child then do; output; ranint=.; end; drop effect estimate; run; proc sort data=rand22; by gender child; run; title3 'RECONFIGURED RANDOM EFFECTS DATA SET'; proc print data=rand22; run; /******************************************************************* ******************************************************************** MODIFIED RANDOM COEFFICIENT MODEL WITH 24 INTERCEPTS RANDOM, SLOPES FIXED RECONFIGURED RANDOM EFFECTS DATA SET Obs child gender ranint 1 1 0 -1.21545 2 2 0 0.33642 3 3 0 1.05266 4 4 0 2.12703 5 5 0 -0.02170 6 6 0 -1.45420 7 7 0 0.33642 8 8 0 0.69454 9 9 0 -1.45420 10 10 0 -3.96105 11 11 0 3.55952 --------------------------------------- ******************************************************************* ********************************************************************/ data both2; merge fixed22 rand22; by gender; beta0i=fixint+ranint; beta1i=fixslope; run; title3 'RANDOM INTERCEPTS AND FIXED SLOPES'; proc print data=both2; run; /******************************************************************* ******************************************************************** MODIFIED RANDOM COEFFICIENT MODEL WITH 25 INTERCEPTS RANDOM, SLOPES FIXED RANDOM INTERCEPTS AND FIXED SLOPES Obs gender fixint fixslope child ranint beta0i beta1i 1 0 17.3727 0.47955 1 -1.21545 16.1573 0.47955 2 0 17.3727 0.47955 2 0.33642 17.7091 0.47955 3 0 17.3727 0.47955 3 1.05266 18.4254 0.47955 4 0 17.3727 0.47955 4 2.12703 19.4998 0.47955 5 0 17.3727 0.47955 5 -0.02170 17.3510 0.47955 6 0 17.3727 0.47955 6 -1.45420 15.9185 0.47955 7 0 17.3727 0.47955 7 0.33642 17.7091 0.47955 8 0 17.3727 0.47955 8 0.69454 18.0673 0.47955 9 0 17.3727 0.47955 9 -1.45420 15.9185 0.47955 10 0 17.3727 0.47955 10 -3.96105 13.4117 0.47955 11 0 17.3727 0.47955 11 3.55952 20.9322 0.47955 12 1 16.3406 0.78438 12 2.28494 18.6256 0.78438 13 1 16.3406 0.78438 13 -1.30935 15.0313 0.78438 14 1 16.3406 0.78438 14 -0.59049 15.7501 0.78438 15 1 16.3406 0.78438 15 1.36069 17.7013 0.78438 16 1 16.3406 0.78438 16 -1.61743 14.7232 0.78438 17 1 16.3406 0.78438 17 1.15531 17.4959 0.78438 18 1 16.3406 0.78438 18 -1.00127 15.3394 0.78438 19 1 16.3406 0.78438 19 -0.89857 15.4421 0.78438 20 1 16.3406 0.78438 20 0.12837 16.4690 0.78438 21 1 16.3406 0.78438 21 3.72265 20.0633 0.78438 22 1 16.3406 0.78438 22 -1.10396 15.2367 0.78438 23 1 16.3406 0.78438 23 -0.59049 15.7501 0.78438 24 1 16.3406 0.78438 24 -0.59049 15.7501 0.78438 25 1 16.3406 0.78438 25 -0.07702 16.2636 0.78438 26 1 16.3406 0.78438 26 0.74453 17.0852 0.78438 27 1 16.3406 0.78438 27 -1.61743 14.7232 0.78438 ******************************************************************* ********************************************************************/