Date:         Thu, 5 Jun 2003 10:15:16 +0200
Reply-To:     "Dr. Olaf Kruse" <olaf.kruse@VST-GMBH.DE>
Sender:       "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From:         "Dr. Olaf Kruse" <olaf.kruse@VST-GMBH.DE>
Subject:      MLE in SAS -- Please help!
Content-Type: text/plain; charset="iso-8859-1"

>Guys,
>
>I am trying to do a MLE analysis in SAS. Assume I have 2000
>observations, each containing a response Y and a X vector. Assume
>Y(i) follows Normal distribution with mean 0 and standard deviation
>X'Beta. Is there any SAS procedure that can fit the Beta vector?
>
>Thanks,
>FG

Below you'll find a sniplet from the SAS-sample-library. It's also part of
the SAS-documentaion for PROC GENMOD. Maybe you can modify it for your
purpouses. I would also recommend to take  look into the literature (as a
starter: M. Aitkin: Modelling Variance-Heterogeneity in Normal Regression
using GLIM, Jou.Roy.Stat.Soc. C 36, 1987, p. 332-229).

HTH Olaf




 /* Example 7 */
   *---------------Modeling Variance Heterogeneity-----*
   *---------------------------------------------------*;

data semi;
   drop j;
   do p = 'a','b';
      do x = 10 to 15;
         if p = 'a' then do;
            do j = 1 to 2;
               y = x + 2*rannor(1302);
               output;
            end;
         end;
         else if p = 'b' then do;
            do j = 1 to 2;
            y = 3*x + 5*rannor(4567);
            output;
            end;
         end;
      end;
   end;
   run;

%macro semimod;

ods listing close;
ods output;

data work;
* input data set;
   set semi;
   wgt = 1;
   sum = 0;
run;

data tmp;
   sum = 0;
run;

%let conv = 0;
%let iter = 0;
* iterate until convergence;
%do %while( &conv = 0 );

data _old;
   set tmp;
   devold = sum;
   keep devold;
run;

* mean model;
* set NOSCALE so that scale is not estimated;
* select OBSTATS to get residuals;
* SCWGT selects dispersion parameter weights;
proc genmod data=work;
   class p;
   scwgt wgt;
   make 'obstats' out=A;
   make 'parameterestimates' out=meanests;
   model y = p x p*x / noscale
                       obstats ;
run;

* this data set contains squared residuals;
data work;
   drop sum pred xbeta std hesswgt upper lower
        resraw reschi resdev;
   set A ;
   set work;
   rsquare = resraw*resraw;
run;

* variance model;
* set scoring=100 to get Fisher scores;
* set scale = .5 for 1 dof in gamma distribution;
proc genmod data=work;
   class p;
   make 'obstats' out=C;
   make 'parameterestimates' out=varests;
   model rsquare = p x p*x / dist=gamma
                             link=log
                             obstats
                             scoring=100
                             noscale
                             scale=.5;
   run;

* get weights for 1st model;
* compute sum = overall deviance;
data work;
   drop xbeta std hesswgt upper lower resraw reschi resdev;
   set work ;
   set C;
   wgt = 1./pred;
   sum + (rsquare/pred + log(pred) + log(6.283185));
run;

data tmp;
   set work nobs = last;
   keep sum;
   if( _n_ = last );
run;

%let iter=%eval( &iter+1 );
%put &iter;

* check convergence;
data _NULL_;
   set tmp;
   set _old;
   put sum devold;
   if( abs(sum - devold) <= 1.e-3 ) then conv = 1;
   else conv = 0;
   call symput( 'dev', left(put(sum,12.4)) );
   call symput( 'conv', left(put(conv,3.)) );
run;

%end;

title 'Mean Model';
ods listing;
proc print data=meanests; run;
title 'Variance Model';
proc print data=varests; run;
title;
%put Number of Iterations: &iter;
%put Overall Deviance: &dev;
%mend semimod;


%semimod;