libname tab 'C:\GRAN\HARG';
options ps=650 ls=256 pageno=1;

title 'Data analysis (adjusted data)';

*Get Diameters and Volumes as provenance-means;
data work.data_new;
set tab.Harg_totdata_ba;
if hojd11 > 30 then hojd11=.; *Remove obviously erroneous height (200 m...);
if vol11>10000 then vol11=.;
if voltot11>10000 then voltot11=.;
if diam94 >= (diam11*1) then diam11=. ;
if diam94 >= (diam11*1) then voltot11=.; * Remove diameters that have not grown 94-11;
if diam11 > (1.5*diam94) then diam11=.;
if diam11 > (1.5*diam94) then voltot=.; * Remove diameters that have grown > 50% in 15 yrs;
run;

proc sort data=work.data_new;
by uppr fled;
run;

proc means data=work.data_new noprint;
var DIAM11 VOLTOT11;
by uppr fled;
output out=work.HARG11_DIA_VOL mean=DIAM11 VOLTOT11 n=n_dv;
run;

data tab.HARG11_DIA_VOL;
set work.HARG11_DIA_VOL;
volha_11=VOLTOT11*n_dv/10;
run;

*Get sorted values of BA;
proc sort data=tab.Harg_ba;
by uppr fled;
run;

* Now add these provenance means to the old data;
data tab.Datatot_all_means;
*merge tab.datatot (keep=uppr fled h81 volha grkval diam hojd);
merge tab.datatot
      tab.Nscore_means tab.R75 tab.HARG11_DIA_VOL tab.HARG_ba;
by UPPR FLED;
h75=h75/10;
hdif7581=h81-h75;
hdif8194=hojd-h81;
*Also add the group;
If fled eq 11 then grupp='sv';
If fled eq 121 then grupp='sv';
If fled eq 124 then grupp='sv';
If fled eq 162 then grupp='sv';
If fled eq 4 then grupp='sv';
If fled eq 12 then grupp='sv';
If fled eq 329 then grupp='ve';
If fled eq 123 then grupp='ve';
If fled eq 175 then grupp='ve';
If fled eq 109 then grupp='ve';
If fled eq 1 then grupp='oe';
If fled eq 79 then grupp='oe';
If fled eq 46 then grupp='oe';
If fled eq 77 then grupp='oe';
If fled eq 78 then grupp='oe';
If fled eq 2 then grupp='oe';
*Create the total volume i.e. volha_11 + vol_gall;
vol_tot_w_thinn=volha_11+vol_gall;
run;

*This cannot be done with the new data!!
* Get the correlations (old and new);
*proc corr data=tab.Datatot_all_means;
*var h75 h81 hojd hdif7581 hdif8194 diam volha grkval RAKH94 SPROT94 busk fr75 FR81 DIAM11 VOLHA_11 VOLTOT11 BA;
*with h75 h81 hojd hdif7581 hdif8194 diam volha grkval RAKH94 SPROT94 busk fr75 FR81 DIAM11 VOLHA_11 VOLTOT11 BA;
*run;

*************************************************;
* Only the group comparison is reliable
*************************************************;

* Statistical analysis - group comparison;

proc glm data=tab.Datatot_all_means;
class uppr grupp;
*model diam11=uppr grupp; *Use this line to check group differences for diameter 2011;
*model ba=uppr grupp; *Use this line to check group differences for basal area 2011;
model vol_tot_w_thinn=uppr grupp; *Use this line to check group differences for total volume (incl thinned volume 1995) 2011;
lsmeans grupp /pdiff lines;
RUN;
quit;

*This cannot be done with the new data!!
*proc glm data=tab.Datatot_all_means;
*class uppr fled;
*model vol_tot_w_thinn=uppr fled;
*means fled/tukey;
*lsmeans fled;
*RUN;
*quit;
