/*
    Data series with SD that is 2 times larger than what Poisson model,
    using the relation sd = sqrt[mean], would predict.
*/
options ps=55 ls = 75 nodate;

DATA series;

INPUT

id
n_cases /* note mean = 16, sd = 8 [a little less if use divisor of n-1]    */
        /* method of Maximum Likelihood uses divisor of n]                 */ 
        /* very lumpy distribution of chocolate chips or raisins or olives */
        /* assume denominator is constant, == 1                            */
        /* so ln[denominator] = 0, no need for offset                      */
;
LINES;
 1    8
 2   24
 3    8
 4   24
 5    8
 6   24
 7    8
 8   24
 9    8
10   24
;

title1 'null [mean- or intercept-only ] Gaussian model ';
title2 'SD = Root MSE is slightly > 8, when estimated with divisor of n-1 = 9  ';
title3 'SE[mean] = SE[intercept] = SD / sqrt[10] by usual SE formula   ';
proc reg;
  model  n_cases  = ;
run;

title1 'null [mean- or intercept-only ] Gaussian model ';
title2 'note SD = scale =  8, when estimated by Max. Likelihood  ';
title3 'SE[mean] = SE[intercept] = SD / sqrt[10] by usual SE formula   ';
proc genmod;
  model  n_cases  = ;
run;

title1 'null [mean- or intercept-only ] POISSON model ';
title2 'DEFAULT scale =  1 =>> no allowance for extra-Poisson variation   ';
title3 'again, SE[mean] = SE[intercept] should be SD/sqrt[10] ';
title4 'model-based SE uses Poisson-based SD ie SD[Poisson] = sqrt[mean] = 4 ';
proc genmod;
  model  n_cases  = /
            link = identity dist = poisson ;
run;

title1 'null [mean- or intercept-only ] POISSON model, BUT with SE adjusted ';
title2 'dscale  =>> estimate the extra-Poisson variation and adjust SE accordingly';
title3 'so, real SD = SD[Poisson] /sqrt[10] ';
title4 'so SE[mean] = scale * SD[Poisson] /sqrt[10] ';
proc genmod;
  model  n_cases  = /
            link = identity dist = poisson dscale ;
run;


title1 'null [mean- or intercept-only ] POISSON model, empirical SE ';
title2 'can use GEE approach even if no clusters';
title3 'SE[mean] = SD[empirical] /sqrt[10]';
title4 'cf. Gaussian model via ML. SD is estimated SEPARATELY from MEAN';
proc genmod;
  class id;
  model  n_cases  = /
            link = identity dist = poisson  ;
  repeated subject = id ;
run;