$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';