job 


" CHICKWEED GROWTH EXPT 1996
ANNA CHRYSTAL, SAC, CROP SYSTEMS, EDINBURGH
PROGRAMMED BY KERRY BROWN, BIOSS, EDINBURGH

  Example of Genstat 5 program to fit an individual Gompertz curve
  to growth data from each experimental unit and to save the estimated
  Gompertz parameters as 4 variates to which an analysis of variance is
  applied. 


SCALARS SET IN PROGRAM AS FOLLOWS :

     ND - No. of times on which growth was recorded
     NU - No. of experimental units i.e. 54
     NC - No. of coefficients in model i.e. 4 for gompertz

   Treatment factors :

     PN - No. of populations
     CN - No. of collections
      T - No. of temperatures
     RN - No. of replicates


DATA INPUT READ FROM FILE chick.dt IN FOLLOWING FORMAT : 

  Times of recording growth (in this case days from start of expt.)
     e.g. 1 2 3 4 5 ...., :

  Read (in parallel) :

    population number 
    collection number 
    temperature       
    replicate number  
    no. of seeds germinating at each recording [1...ND] 

 "      
set [outprint=dots ; diagnostic=faults ]

scalar [value=24] ND
scalar [value=54] NU
scalar [value=5] NC
scalar [value=3] PN
scalar [value=3] CN
scalar [value=2] T
scalar [value=3] RN

calc NC1=NC*NU
calc NU1=1-NU
calc AA=CN*T*RN
calc BB=T*RN
calc CC=PN*CN

variate [nvalues=NU] day[1...ND]
variate [nvalues=ND] newday[1...NU]
variate [nvalues=NC1] coeffs
variate [nvalues=ND] x

open 'chick.dt' ;chan=2
read [chan=2] x
read [chan=2] popn, collcn, temp, rep, day[1...ND]


equate [oldformat=!(((1,#NU1)#ND,-1)#NU)] day; newstructures=newday
calc newday2[1...NU]=cum(newday[1...NU])

for i=1...NU
  model newday2[i]; fittedvalues=fitted[i]
  fitcurve [curve=gompertz;p=*] x
  rkeep newday2[i]; estimate=coeff[i]
  graph[nrows=11; ncolumns=31] y=newday2[i], fitted[i]; x=x;method=curve,p
endfor

equate !P(coeff[1...NU]); newstructures=coeffs
variate [values=(1...#NC)#NU] ind ; decimals=0

for j=1...NC
  subset [condition=ind.eq.j] oldv=coeffs; newv=coef[j]
endfor

factor [levels=3;values=#AA(1...#PN)] popn2
factor [levels=3;values=#BB(1...#CN)#PN] collcn2
factor [levels=2;values=#RN(1...#T)#CC] temp2
factor [levels=3;values=(1...#RN)#AA] rep2

print rep2, temp2, popn2, collcn2, coef[1...NC] ; field=4(4),5(8) ; \
        deci=4(0),5(2)

for i=1...NC
  tabulate [class=temp2,popn2 ; margins=yes ; percent=50] coef[i] ; \
   quantiles=tab[i]
  print tab[i]

endfor

endjob
stop