3/25/2010 36-402/608 ADA-II H. SeltmanBreakout #17 CommentsThis problem comes from http://www.statsci.org/data/oz/wallaby.html.The data give growth measurements on Tammar wallabies (Macropus eugenii). We willfocus on the pattern of change in weight in grams (original variable is tenths of grams)between ages 100 and 300 days. Other potential explanatory variables are gender andlocation.Here is some EDA using separate splines for each animal: ! "# $ % & ' ( ) * + , -%. + / 0 1 %$ # -%$ 2 3 (45 1 6%47 %8 /+ 9 ) 2 $ # : 5 ( %47 %; + , )2 $ # !"##"$%&!'()*+&,-./+*+ 7 49< < < = < > = # = = = ? = > ? = ? !> 2 > # > = > ? > ! > @ ! 2 ! $ ! >! ! @ $ @ # @ < @ ? @ ! 2 2 2 2 $2 < 2 = 2 ! 2 2 $ 2 2 # 2 2 < 2 2 = 2 2 ? 2 2 @2 $ 2 $ $ 2 $ # 2 $ < 2 $ = 2 $ ? 2 $ > 2 # 2 2 # <2 # = 2 # ? 2 # ! 2 # @ 2 < 2 2 < $ 2< ! 2 < @ 2 = 2 = $ 2 = < 2 = > 2 ? = 2 ? ?Question 1: What fixed effect model would you fit? What random effects would you con-sider including? Using your knowledge of biology and statistics, why would this analysisbe harder for birth through death?Clearly there is at least one bend in the weight vs. age graph, and perhaps two, so (unre-lated to mixed modeling) at least a square and possibly a cubic polynomial term is need tofit this EDA plot. (Transformations other than polynomial might also work.) Due to thecommon prevalence of sexual dimorphism across mammalian species, allowing separatemale and female polynomial shapes is a good idea, so a gender term plus interaction ofgender with the polynomial terms is worth while.The individual curves do NOT look parallel, so in addition to a random intercept, weshould try random slopes and curvatures. And since this is truly a repeated measuresanalysis, the benefit of a within-subject serial correlation model should be checked.The curves we see are for a period of rapid growth (presumable adolescence). A full lifeto death growth curve would need a much higher order polynomial (or maybe splines) toaccommodate high newborn growth rates, then a slower growth rate, then an adolescentspurt, then leveling off or slow adult growth.2Here are SAS results for a rich fixed model with just a random intercept:title "Handout 17 Wallaby Data";data wallaby;infile "wallaby.dat" firstobs=2;input anim sex loca$ leng head ear arm leg pes tail weight age;grams = weight/10;male = 1; <<< create male indicator variableif sex=2 then male=0;drop leng head ear arm leg pes tail sex weight; <<< drop unneeded columnsif age<100 OR age>300 then delete; <<< drop unneeded rowsdaysC = age-100; <<< put the intercept into the datadaysC2 = daysC*daysC; <<< compute polynomial termsdaysC3 = daysC*daysC2;run;proc print data=here.wallaby(obs=5);run;title2 "EDA";proc freq;tables loca male;run;proc univariate;var age grams;run;title2 "Rich fixed effects + random intercept";proc mixed covtest;class loca male;model grams = daysC|male daysC2|male daysC3|male loca;/* expands to male + days + male:days + days2 + male:days2 + days3 + male:days3 + loca */random int / subject=anim;run;### The log:NOTE: Convergence criteria met.NOTE: The PROCEDURE MIXED printed pages 3-4.3### Key results:Rich fixed effects + random interceptThe Mixed ProcedureModel InformationData Set HERE.WALLABYDependent Variable gramsCovariance Structure Variance Components >> meaningless for RI onlySubject Effect anim << one random intercept per animalEstimation Method REML << Unbiased for random effectsResidual Variance Method Profile >>Fixed Effects SE Method Model-Based >> Highly technical infoDegrees of Freedom Method Containment >>Class Level InformationClass Levels Valuesloca 12 "G" "H1" "H12" "H2" "H3" "H7""H8" "H9" "Ha" "Hb" "K" "W"male 2 0 1DimensionsCovariance Parameters 2 << intercept variance and residual varianceColumns in X 24 << useful betas plus some overparameterizationColumns in Z Per Subject 1 << just a random intercept (group indicator var.)Subjects 59Max Obs Per Subject 16Number of ObservationsNumber of Observations Read 600Number of Observations Used 600Number of Observations Not Used 04Iteration HistoryIteration Evaluations -2 Res Log Like Criterion0 1 7598.187086721 2 7357.02282769 0.00000000Convergence criteria met. << We really need this!Covariance Parameter EstimatesStandard ZCov Parm Subject Estimate Error Value Pr > ZIntercept anim 11806 2939.27 4.02 <.0001Residual 11624 710.23 16.37 <.0001>>p-values are no very reliable. Estimates are variances. Corresponding>>square roots reflect the size of subject-to-subject variability (Intercept)>>and within-subject variability (residual).Fit Statistics-2 Res Log Likelihood 7357.0AIC (smaller is better) 7361.0AICC (smaller is better) 7361.0BIC (smaller is better) 7365.2 << Smaller is better. Compare REML to REML.<< No meaning except comparing models.Type 3 Tests of Fixed EffectsNum DenEffect DF DF F Value Pr > FdaysC 1 534 24.07 <.0001male 1 534 0.04 0.8344daysC*male 1 534 1.90 0.1681daysC2 1 534 89.82 <.0001daysC2*male 1 534 5.47 0.0197daysC3 1 534 2.25 0.1345daysC3*male 1 534 8.91 0.0030 << Useless fixed effects will beloca 11 534 0.74 0.6991 << dropped later.Question 2: Explain everything except “highly technical” and AICC.I note that the BIC is much smaller than for the fixed effects only model (not shown;generated by dropping the RANDOM statement).5/* With more than one random effect (here, random int. and slope) useTYPE=UN(STRUCTURED) to allow correlated random effects. */title2 "Rich fixed effects + random intercept + random time";proc mixed covtest;class loca male;model grams = daysC|male daysC2|male daysC3|male loca;random int daysC/ subject=anim type=UN;run;### The log:NOTE: Convergence criteria met.### Selected results:Model InformationCovariance Structure UnstructuredDimensionsCovariance Parameters 4Columns in X 24Columns in Z Per Subject 2Iteration HistoryConvergence criteria met.Covariance Parameter EstimatesStandard ZCov Parm Subject Estimate Error Value Pr ZUN(1,1) anim 1293.45 553.72 2.34 0.0097 << rand. int.UN(2,1) anim -53.1466 14.6732 -3.62 0.0003 << covUN(2,2) anim 2.2014 0.4606 4.78 <.0001 << rand. slopeResidual 3950.47 250.82 15.75 <.0001Fit StatisticsBIC (smaller is better) 6832.4Question 3: Compare the models. Calculate the estimated correlation of theintercept and slope: UN(2,1)/sqrt(UN(1,1))/sqrt(UN(2,2)).This model has a 2-by-2 random effect covariance matrix (for each animal) that hasintercept variance 1293, slope variance 2.20, and correlation of the random intercept andslope