/*------------------------------------ Problem 2. Orange Tree Data ------------------------------------*/ /* PROBLEM: Based on the following orange tree data, the problem is to run a nonlinear mixed-effects */ /* model so as to obtain estimates of the population growth rate parameters for a NLME model */ /* having a logistic growth curve structure similar to the NLMIXED or NLNMIX code below. The */ /* problem is to start with the SAS code below and with the aid of the %COVPARMS macro (you */ /* need to load and run this), determinine what parameters should or can be treated as random */ /* and which should be treated as strictly fixed-effects parameters (i.e. have no random */ /* effects associated with them. Specifically, we are looking to justify why only the */ /* limiting growth parameter (asymptote) should be selected as the only random effect. */ /* Then run both NLMIXED with METHOD=FIRO and with METHOD=GAUSS with QPOINTS=1 and compare */ /* those results to NLINMIX using EXPAND=ZERO and EXPAND=EBLUP. */ /* HINT: Use the COVPARMS statement from the CEFAMANDOLE data example and adapt it to this problem. */ /*-----------------------------------------------------------------------------------------------------*/ %include 'c:\edward\macro\nlinmix\covparms.sas' / nosource2;run; %include 'c:\edward\macro\nlinmix\nlinmix.sas' / NOSOURCE2; data otree; input TREE $ DAYS Y; INTERCEP=1; x=days; GROUP=1; cards; 1 118 30 1 484 58 1 664 87 1 1004 115 1 1231 120 1 1372 142 1 1582 145 2 118 33 2 484 69 2 664 111 2 1004 156 2 1231 172 2 1372 203 2 1582 203 3 118 30 3 484 51 3 664 75 3 1004 108 3 1231 115 3 1372 139 3 1582 140 4 118 32 4 484 62 4 664 112 4 1004 167 4 1231 179 4 1372 209 4 1582 214 5 118 30 5 484 49 5 664 81 5 1004 125 5 1231 142 5 1372 174 5 1582 177 ; proc sort DATA=otree; by tree; run; proc print data=otree; run; proc nlmixed data=otree qpoints=1; parms b1=150 b2=700 b3=350 sigma=60 psi21=0 psi31=0 psi32=0; num = b1+u1; e = exp(-(x-(b2+u2))/(b3+u3)); den = 1 + e; predmean = (num/den); predvar = sigma; resid = y - predmean; model y ~ normal(predmean, predvar); random u1 u2 u3 ~ normal([0,0,0], [psi11, psi21, psi22, psi31, psi32, psi33]) subject=tree; predict predmean out=predout der; id resid; run; %nlinmix(data=otree, procopt=method=ml empirical, model=%str( num = b1+u1; den = 1+exp(-(x-(b2+u2))/(b3+u3)); predv=num/den; ), parms=%str(b1=150 b2=700 b3=350), stmts=%str(class tree; model pseudo_y = d_b1 d_b2 d_b3 / noint notest s cl; random d_u1 d_u2 d_u3 / subject=tree type=un s cl; ), expand=zero ); /* START CODE TO GET REDUCED RANDOM-EFFECTS STRUCTURE */ %macro get_starting_values; ods output ParameterEstimates=OLSparms; proc nlmixed data=otree qpoints=1; parms b1=150 b2=700 b3=350 sigma=60; num = b1+u1; e = exp(-(x-(b2+u2))/(b3+u3)); den = 1 + e; predmean = (num/den); predvar = sigma; resid = y - predmean; model y ~ normal(predmean, predvar); random u1 u2 u3 ~ normal([0,0,0], [0, 0, 0, 0, 0, 0]) subject=tree; predict predmean out=predout der; id resid; run; proc print data=predout; run; %covparms(parms=OLSparms, predout=predout, resid=resid, method=ml, random=der_u1 der_u2 der_u3, subject=tree, type=un, covname=psi, output=MLEparms); proc nlmixed data=otree qpoints=1; parms /data=MLEparms; num = b1+u1; e = exp(-(x-(b2+u2))/(b3+u3)); den = 1 + e; predmean = (num/den); predvar = sigma; resid = y - predmean; model y ~ normal(predmean, predvar); random u1 u2 u3 ~ normal([0,0,0], [psi11, psi21, psi22, psi31, psi32, psi33]) subject=tree; run; data MLEparms1; set MLEparms; if Parameter in ('psi31' 'psi32' 'psi33') then delete; run; %covparms(parms=OLSparms, predout=predout, resid=resid, method=ml, random=der_u1 der_u2 , subject=tree, type=un, covname=psi, output=MLEparms1); proc nlmixed data=otree qpoints=1; parms /data=MLEparms1; num = b1+u1; e = exp(-(x-(b2+u2))/b3); den = 1 + e; predmean = (num/den); predvar = sigma; resid = y - predmean; model y ~ normal(predmean, predvar); random u1 u2 ~ normal([0,0], [psi11, psi21, psi22]) subject=tree; run; %covparms(parms=OLSparms, predout=predout, resid=resid, method=ml, random=der_u1 der_u2 der_u3, subject=tree, type=vc, covname=psi, output=MLEparms2); proc nlmixed data=otree qpoints=1; parms /data=MLEparms2; num = b1+u1; e = exp(-(x-(b2+u2))/(b3+u3)); den = 1 + e; predmean = (num/den); predvar = sigma; resid = y - predmean; model y ~ normal(predmean, predvar); random u1 u2 u3 ~ normal([0,0,0], [psiu1, 0 , psiu2, 0 , 0 , psiu3]) subject=tree; run; %covparms(parms=OLSparms, predout=predout, resid=resid, method=ml, random=der_u1, subject=tree, type=vc, covname=psi, output=MLEparms3); proc nlmixed data=otree qpoints=1; parms /data=MLEparms3; num = b1+u1; e = exp(-(x-b2)/b3); den = 1 + e; predmean = (num/den); predvar = sigma; resid = y - predmean; model y ~ normal(predmean, predvar); random u1 ~ normal(0, psiu1) subject=tree; run; %mend get_starting_values; %macro final_models; proc nlmixed data=otree method=firo; parms /data=MLEparms3; num = b1+u1; e = exp(-(x-b2)/b3); den = 1 + e; predmean = (num/den); predvar = sigma; resid = y - predmean; model y ~ normal(predmean, predvar); random u1 ~ normal(0, psiu1) subject=tree; run; proc nlmixed data=otree method=gauss; parms /data=MLEparms3; num = b1+u1; e = exp(-(x-b2)/b3); den = 1 + e; predmean = (num/den); predvar = sigma; resid = y - predmean; model y ~ normal(predmean, predvar); random u1 ~ normal(0, psiu1) subject=tree; run; %nlinmix(data=otree, procopt=method=ml empirical, model=%str( num = b1+u1; den = 1+exp(-(x-b2)/b3); predv=num/den; ), parms=%str(b1=150 b2=700 b3=350), stmts=%str(class tree; model pseudo_y = d_b1 d_b2 d_b3 / noint notest s cl; random d_u1 / subject=tree type=un s cl; ), expand=zero ); %nlinmix(data=otree, procopt=method=ml empirical, model=%str( num = b1+u1; den = 1+exp(-(x-b2)/b3); predv=num/den; ), parms=%str(b1=150 b2=700 b3=350), stmts=%str(class tree; model pseudo_y = d_b1 d_b2 d_b3 / noint notest s cl; random d_u1 / subject=tree type=un s cl; ), expand=eblup ); %mend final_models;