Reducing Parameter Variance in a Tree Biomass Model

/*--------------------------------------------------------------

                    SAS Sample Library

        Name: modex22.sas
 Description: Example program from SAS/ETS User's Guide,
              The MODEL Procedure
       Title: Reducing Parameter Variance in a Tree Biomass Model
     Product: SAS/ETS Software
        Keys: SUR, systems of equations, heteroscedasticity
        PROC: MODEL
       Notes:

--------------------------------------------------------------*/

data trees;
   input tree Dbh Height LCL Age Wood Bark Crown total;
datalines;
    1   73.2      29.0    16.2      93    4463.4    572.9     186.4   5222.7
    2   30.5      18.3    10.4      40     550.7     83.9      44.9    679.5
    3   48.3      22.9    11.9      69    1689.2    225.0      93.4   2007.6
    4   69.6      27.4    18.3      74    3441.5    435.5     178.3   4055.3
    5   28.7      19.8    11.6      38     482.2     75.3      45.4    602.9
    6   53.1      32.0    14.9      78    2281.6    307.5      57.2   2646.3
    7   45.7      32.0    17.7      79    1771.3    230.4      24.9   2026.6
    8   46.5      30.5    17.4      83    1611.6    228.2      17.2   1857.0
    9   33.8      25.9    10.1      68     861.8    122.5      20.0   1004.3
   10   55.4      30.5    13.1      70    2952.9    367.4      65.3   3385.6
   11   30.5      24.4     9.1      75     679.9    125.2      17.7    822.8
   12   70.1      30.5    13.4      81    3867.8    546.1      79.4   4493.3
   13   41.9      24.4    12.2      79    1289.1    185.5      36.7   1511.3
   14   42.9      25.9    13.7      76    1495.5    201.4      71.7   1768.6
   15   66.3      35.1    17.7      81    4091.0    413.2      99.3   4603.5
   16   40.6      25.9    11.6      64    1264.2    175.5      12.7   1452.4
   17   28.7      21.3     9.8      75     485.8     67.6      19.1    572.5
   18   80.5      32.0    17.4      93    5782.0    657.7     186.9   6626.6
   19   66.5      33.8    16.8     106    4085.1    524.4     112.0   4721.5
   20   61.2      29.0    17.1      71    2621.4    225.4      82.6   2929.4
   21   21.8      25.9    11.3      35     292.6     51.7       5.0    349.3
   22   29.0      24.4     8.8      35     616.4     94.3      21.8    732.5
   23   49.8      27.4    16.2      35    1757.2    220.0      63.5   2040.7
   24   34.3      27.4    10.4      37     902.2    144.7      21.8   1068.7
   25   43.2      27.4    15.5      41    1251.5    212.3      20.9   1484.7
   26   25.1      24.4     8.8      38     437.7     63.5       6.4    507.6
   27   29.7      25.9    13.7      40     704.4     89.4      10.0    803.8
   28   34.8      25.9    13.1      40     906.7    117.5      23.1   1047.3
   29   38.6      27.4    13.4      41    1309.5    148.8      26.8   1485.1
   30   49.8      22.9    14.3      40    1497.3    160.1      49.0   1706.4
   31   36.1      25.9    16.5      34     794.3    116.6      20.9    931.8
   32   37.1      21.3    11.0      67     846.9    120.7      26.8    994.4
   33   33.0      21.3     9.8      67     635.5     89.8      27.7    753.0
   34   57.7      25.9    18.6      84    2545.1    371.0      70.3   2986.4
   35   53.8      25.9    12.2      87    2275.7    359.3      32.7   2667.7
   36   57.9      27.4    10.1      89    2822.3    379.2      61.2   3262.7
   37   75.4      27.4    13.7      91    3782.1    579.2      61.7   4423.0
   38   57.2      25.9    12.2      86    2055.7    362.0      45.4   2463.1
   39   69.1      27.4    14.6      87    3618.4    498.1      87.1   4203.6
;

proc model data=trees;
   endo wood bark crown total;

   vol = dbh**2*height;
   cvol = vol*lcl/1000;

   wood  = b10 + b11*vol;
   bark  = b20 + b21*vol;
   crown = b30 + b31*cvol + b32*height;
   total = b40 + b41*vol + b42*cvol + b43*height;

   fit / ols outs=stree out=otree outresid;
quit;

proc model data=trees;
ods output ParameterEstimates=ineff;
   endo wood bark crown total;

   vol = dbh**2*height;
   cvol = vol*lcl/1000;

   wood  = b10 + b11*vol;
   bark  = b20 + b21*vol;
   crown = b30 + b31*cvol + b32*height;
   total = b40 + b41*vol + b42*cvol + b43*height;

   fit / nools sur sdata=stree;

   restrict b10 + b20 + b30 = b40,
            b11 + b21 = b41,
            b42 = b31,
            b43 = b32;
quit;

proc model data=otree;
ods output parameterestimates=vparms;
   vol = dbh**2*height;
   cvol = vol*lcl/1000;

   eq.varwood  = log(wood*wood) - (w1 + w2*log(vol));
   eq.varbark  = log(bark*bark) - (k1 + k2*log(vol));
   eq.varcrown = log(crown*crown) - (c1 + c2*log(cvol) - c3*height**2);
   eq.vartotal = log(total*total) - (t1 + t2*log(vol));

   fit varwood varbark varcrown vartotal;
quit;

data _null_;
   set vparms;
   call symput(Parameter,Estimate);
run;

proc model data=trees;
ods output ParameterEstimates=eff;
   endo wood bark crown total;

   vol = dbh**2*height;
   cvol = vol*lcl/1000;

   wood  = b10 + b11*vol;
   bark  = b20 + b21*vol;
   crown = b30 + b31*cvol + b32*height;
   total = b40 + b41*vol + b42*cvol + b43*height;

   h.wood = exp(&w1)*vol**&w2;
   h.bark = exp(&k1)*vol**&k2;
   h.crown = exp(&c1)*cvol**&c2 * exp (-&c3*height*height);
   h.total = exp(&t1)*vol**&t2;

   fit / sur;

   restrict b10 + b20 + b30 = b40,
            b11 + b21 = b41,
            b42 = b31,
            b43 = b32;
quit;

data varred;
   keep parameter stderrchange;
   merge ineff eff(rename=(estimate=effestimate stderr=efferr));
   format stderrchange percent.;
   label stderrchange = "StdErr Rel. Change";
   stderrchange = (efferr - stderr)/stderr;
run;

proc print data=varred(obs=11) noobs label;
run;