$ontext

   Example of estimation of the Brody growth curve

   Reference:
        Miroslav Kaps, William R. Lamberson
        Biostatistics for Animal Science
        2004 CABI Publishing

   SAS Input:

   DATA a;
   INPUT age weight @@;
   DATALINES;
      8  280    12  340    24  340    36  480    48  550
      60 580    72  590    84  600    96  590   108  600
   ;
   PROC NLIN;
      PARMS A=600 weight0=280 k=0.05
      MODEL weight=A-(A-weight0)*exp(-k*(age-8));
      RUN;


   SAS Output

                    Dependent Variable weight
                       Method: Gauss-Newton

                         Iterative Phase
                                                     Sum of
        Iter          A      weight0        k        Squares
           0      600.0       280.0      0.0500       2540.5
           1      610.2       285.8      0.0355       1388.7
           2      612.2       283.7      0.0381        966.9
           3      612.9       283.9      0.0379        965.9
           4      612.9       283.9      0.0380        965.9
           5      612.9       283.9      0.0380        965.9

              NOTE: Convergence criterion met.

                           Sum of        Mean                  Approx
 Source             DF     Squares      Square     F Value     Pr > F
 Regression         3     2663434       887811      446.69     <.0001
 Residual           7       965.9       138.0
 Uncorrected Total 10     2664400
 Corrected Total    9      124240

                         Approx
 Parameter   Estimate   Std Error   Approximate 95% Confidence Limit
 A              612.9     9.2683      590.9       634.8
 weight0        283.9     9.4866      261.5       306.3
 k              0.0380    0.00383     0.0289      0.0470

             Approximate Correlation Matrix
                A            weight0             k

 A          1.0000000       0.2607907     -0.8276063
 weight0    0.2607907       1.0000000     -0.4940824
 k         -0.8276063      -0.4940824      1.0000000


$offtext


*-----------------------------------------------------------------------------
* data represents weight of an Angus cow at different ages
*-----------------------------------------------------------------------------

set i 'data points' /i1*i10/;

table data(i,*)
      weight    age
*       kg    months
  i1    280      8
  i2    340     12
  i3    430     24
  i4    480     36
  i5    550     48
  i6    580     60
  i7    590     72
  i8    600     84
  i9    590     96
  i10   600    108
;

parameter Weight(i), Age(i);
Weight(i) = data(i,'weight');
Age(i) = data(i,'age');


*-----------------------------------------------------------------------------
* statistical model
*-----------------------------------------------------------------------------

scalar Age0 'initial age (months)' /8/;

variables
    sse       'sum of squared errors'
    A         'Estimated asymptotic (mature) weight'
    Weight0   'Estimated initial weight at age0'
    k         'Estimated maturing weight index'
;

equations
    fit(i)     'the non-linear model'
    obj        'objective'
;

obj..     sse =n= 0;
fit(i)..  Weight(i) =e= A - (A-Weight0) * exp[-k*(Age(i)-Age0)];

option nlp=nls;
model nlfit /obj,fit/;

*-----------------------------------------------------------------------------
* initial values
*-----------------------------------------------------------------------------

A.l = 600;
weight0.l = 280;
k.l = 0.05;

solve nlfit using nlp minimizing sse;
display sse.l,A.l,weight0.l,k.l;



*-----------------------------------------------------------------------------
* load and display the confidence intervals
*-----------------------------------------------------------------------------


sets
   alpha /'95%'/
   v /A,Weight0,k/
   interval /lo,up/
;

parameter fitresult(v,*);
fitresult('A',      'Estimate') = A.l;
fitresult('Weight0','Estimate') = Weight0.l;
fitresult('k',      'Estimate') = k.l;
fitresult('A',      'Std.Err.') = A.m;
fitresult('Weight0','Std.Err.') = Weight0.m;
fitresult('k',      'Std.Err.') = k.m;
display fitresult;


parameter
   confint(alpha,v,interval) "confidence intervals"
   confint95(v,interval)     "95% confidence interval"
   covar(v,v)                "Covariance matrix"
   corr(v,v)                 "Correlation matrix"
;
execute_load 'nls.gdx',confint,covar;
confint95(v,interval) = confint('95%',v,interval);
display confint95;

alias(v,vv);
corr(v,vv) = covar(v,vv)/(sqrt(covar(v,v))*sqrt(covar(vv,vv)));
option decimals=7;
display corr;