/*------------------------------------ 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. */
/*-----------------------------------------------------------------------------------------------------*/
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
 );