$ontext

  Non-Linear Least Squares Regression example

  Erwin Kalvelagen,  may 2008

  Reference:
   http://www.statsci.org/data/general/rabbit.html


Age and Eye Lens Weight for Rabbits in Australia
--------------------------------------------------------------------------------

Description
The European rabbit Oryctolagus cuniculus is a major pest in Australia. A reliable method of
age determination for rabbits caught in the wild would be of importance in ecological studies.
In this study, the dry weight of the eye lens was measured for 71 free-living wild rabbits of
known age. Eye lens weight tends to vary much less with environmental conditions than does
total body weight, and therefore may be a much better indicator of age

The rabbits were born and lived free in an experimental 1.7 acre enclosure at Gungahlin, ACT.
The birth data and history of each individual were accurately known. Rabbits in the enclosure
depended on the natural food supply. In this experiment, 18 of the eye lenses were collected
from rabbits that died in the course of the study from various causes such as coccidiosis,
bird predation or starvation. The remaining 53 rabbits were deliberately killed, immediately
after being caught in the enclosure or after they had been kept for some time in cages.
The lenses were preserved and their dry weight determined.



--------------------------------------------------------------------------------
 Variable        Description
--------------------------------------------------------------------------------
 Age             Age of rabbit in days
 Lens            Dry weight of eye lens in milligrams
--------------------------------------------------------------------------------


Source
Dudzinski, M. L., and Mykytowycz, R. (1961). The eye lens as an indicator of age in the wild rabbit in Australia. CSIRO Wildlife Research, 6, 156-159.

Ratkowsky, D. A. (1983). Nonlinear Regression Modelling. Marcel Dekker, New York.

Wei, B.-C. (1998). Exponential Family Nonlinear Models. Springer, Singapore. Examples 2.4 and 6.


Analysis
Dudzinski and Mykytowycz (1961) concluded the eye lens weight is a reliable indicator of the rabbit's
age up to about 150 days. The deterministic component of their model is

    Lens = a exp{ -b / (Age + g) }

Since the variance of Lens increases with Age, it may be appropriate to fit the nonlinear model

    log(Lens) = a - b / (Age + g)


$offtext

set i /i1*i71/;

table data(i,*)

         Age        Lens
 i1      15        21.66
 i2      15        22.75
 i3      15        22.30
 i4      18        31.25
 i5      28        44.79
 i6      29        40.55
 i7      37        50.25
 i8      37        46.88
 i9      44        52.03
 i10     50        63.47
 i11     50        61.13
 i12     60        81.00
 i13     61        73.09
 i14     64        79.09
 i15     65        79.51
 i16     65        65.31
 i17     72        71.90
 i18     75        86.10
 i19     75        94.60
 i20     82        92.50
 i21     85        105.00
 i22     91        101.70
 i23     91        102.90
 i24     97        110.00
 i25     98        104.30
 i26     125       134.90
 i27     142       130.68
 i28     142       140.58
 i29     147       155.30
 i30     147       152.20
 i31     150       144.50
 i32     159       142.15
 i33     165       139.81
 i34     183       153.22
 i35     192       145.72
 i36     195       161.10
 i37     218       174.18
 i38     218       173.03
 i39     219       173.54
 i40     224       178.86
 i41     225       177.68
 i42     227       173.73
 i43     232       159.98
 i44     232       161.29
 i45     237       187.07
 i46     246       176.13
 i47     258       183.40
 i48     276       186.26
 i49     285       189.66
 i50     300       186.09
 i51     301       186.70
 i52     305       186.80
 i53     312       195.10
 i54     317       216.41
 i55     338       203.23
 i56     347       188.38
 i57     354       189.70
 i58     357       195.31
 i59     375       202.63
 i60     394       224.82
 i61     513       203.30
 i62     535       209.70
 i63     554       233.90
 i64     591       234.70
 i65     648       244.30
 i66     660       231.00
 i67     705       242.40
 i68     723       230.77
 i69     756       242.57
 i70     768       232.12
 i71     860       246.70
;

parameters
   Age(i)         'Age of rabbit in days '
   Lens(i)        'Dry weight of eye lens in milligrams'
;

Age(i) = data(i,'Age');
Lens(i) = data(i,'Lens');


variables
   a1  'parameter to be estimated'
   a2  'parameter to be estimated'
   a3  'parameter to be estimated'
   sse 'sum of squared errors'
;

equation
   fit(i)     'equation to fit'
   sumsq      'dummy objective'
;

sumsq..    sse =n= 0;
fit(i)..   log(Lens(i)) =e= a1 - a2/(Age(i)+a3);

option nlp=nls;

models m /sumsq,fit/;
solve m minimizing sse using nlp;



*
* plot results
*

file pltdat /rabbit.dat/;
loop(i,
    put pltdat Age(i):17:5,Lens(i):17:5/;
);
putclose;


file plt /rabbit.plt/;
put plt;
put "a1=",a1.l:0:16/;
put "a2=",a2.l:0:16/;
put "a3=",a3.l:0:16/;
put "fit1(x)=exp(a1-a2/(x+a3))"/;
put "fit2(x)=a1-a2/(x+a3)"/;
putclose 'set term png'/
         'set key off'/
         'set output "rabbit.png"'/
         'set multiplot layout 1,2 scale 1,0.5'/
         'set xlabel "Lens"'/
         'set ylabel "Age"'/
         'plot "rabbit.dat" using 1:2,fit1(x)'/
         'set xlabel "Lens"'/
         'set ylabel "log(Age)"'/
         'plot "rabbit.dat" using 1:(log($2)),fit2(x)'/
         'unset multiplot'/
;

execute 'type rabbit.plt';

execute '=wgnuplot.exe rabbit.plt';
execute '=shellexecute rabbit.png';