/************************************************************************** * * * GLS analysis of the subject 5 indomethacin data * * for variance = 2 times a power theta of the mean. * * * * Here, define a SAS macro to perform the 3-step GLS algorithm * * * **************************************************************************/ options ps=55 ls=80 nodate; /************************************************************************** * * * Define a sas macro called "glsalg" to control the GLS iteration. * * The macro takes as input the name of the SAS working data set, the * * names of the x and y variables (here, x is a scalar, but the program * * may be modified easily to pass more than one explanatory variable), * * and starting values for each WLS calculation. * * * **************************************************************************/ %macro glsalg(dset,xvar,yvar,b1,b2,b3,b4); /************************************************************************** * * * Set up the data set for the first pass through the algorithm, which * * will compute the OLS estimate. The variable "pred1" containing the * * predicted values from the previous iteration is set equal to 1 here, * * so that the first pass with weights based on pred1^theta will actually * * be OLS (weights all = 1). * * * **************************************************************************/ data &dset; set &dset; pred1=1; if &xvar=. or &yvar=. then delete; b1_new=&b1; b2_new=&b2; b3_new=&b3; b4_new=&b4; /************************************************************************** * * * Set the initial value for theta (OLS) * * * **************************************************************************/ theta1=0.0; /************************************************************************** * * * Call the macro "step3," defined below, that actually implements the * * call to PROC NLIN to do WLS with the current set of fixed weights. * * The argument sets options for supressing the printing of the results * * of each call to PROC NLIN. * * * **************************************************************************/ %step3(1) %step3(1) %step3(1) %step3(1) %step3(1) %step3(1) %step3(1) %step3(1) %step3(1) %step3(1) %step3(0) /************************************************************************** * * * After the final iteration, compute the final estimate of sigma * * using the newest values of the weights * * * **************************************************************************/ data sigma; set outnlin(keep = resid pred); data sigma2; set new(keep = theta1); data sigma; merge sigma sigma2; data sigma; set sigma; wresid=resid/(pred**theta1); proc means data=sigma noprint; var wresid; output out=sigout uss=rss; data sigout; set sigout; sigma=sqrt(rss/7); data sigout; set sigout(keep = sigma); if _n_>1 then delete; proc print data=sigout; title2 "Final estimate of sigma"; %mend glsalg; /************************************************************************** * * * Define the macro "step3." The argument "foo" controls options in * * calling PROC NLIN (see below). * * * **************************************************************************/ %macro step3(foo); /************************************************************************** * * * PROC NLIN prints out a summary by default. We print out this full * * output only for the final GLS fit; the "noprint" option is invoked * * for all other fits. * * * **************************************************************************/ %if &foo>0 %then %do; %let opt=noprint; %end; %else %do; %let opt= ; %end; /************************************************************************** * * * The call to PROC NLIN to implement the current WLS fit with theta * * equal to theta1, the previous estimate * * * **************************************************************************/ proc nlin data=&dset method=gauss &opt; parms b1=&b1 b2=&b2 b3=&b3 b4=&b4; if _iter_=-1 then do; b1=b1_new; b2=b2_new; b3=b3_new; b4=b4_new; end; /************************************************************************** * * * Programming statements to define the mean function and its derivatives * * * **************************************************************************/ eb1=exp(b1); eb2=exp(b2); eb3=exp(b3); eb4=exp(b4); f=eb1*exp(-eb2*&xvar) + eb3*exp(-eb4*&xvar); db1=eb1*exp(-eb2*&xvar); db2=-eb1*eb2*&xvar*exp(-eb2*&xvar); db3=eb3*exp(-eb4*&xvar); db4=-eb3*eb4*&xvar*exp(-eb4*&xvar); /************************************************************************** * * * The model statement and derivative statements to tell PROC NLIN the * * form of the model and derivatives * * * **************************************************************************/ model &yvar = f; der.b1=db1; der.b2=db2; der.b3=db3; der.b4=db4; /************************************************************************** * * * The "_weight_" statement defines the values of the weights to use for * * WLS. Here, we use weights computed using the fixed value of theta1 * * and the values of the predicted values from the previous iteration. * * Thus, the weights are FIXED (do not depend on current predicted * * values). PROC NLIN thus knows to do WLS with FIXED weights rather * * than IRWLS. * * * **************************************************************************/ _weight_ = (1/abs(pred1))**(2*theta1); /************************************************************************** * * * Form an output data set containing everything in the input data set * * plus the predicted values (to set up weights for next time), the * * residuals (so we can calculate sigma on the final iteration), and the * * parameter estimates from this iteration (so we can print them). * * * **************************************************************************/ output out=outnlin p=pred r=resid sse=sigma parms=b1 b2 b3 b4; run; /************************************************************************** * * * Now estimate theta by PL with the "Trick", defining a "dummy" * * nonlinear regression problem. Here, fdot = geometric mean of the * * current estimated means, dummy = dummy responses all = 0 * * * **************************************************************************/ data outnlin; set outnlin; if pred<=0 then delete; sigma=sqrt(sigma/7); lpred=log(pred); proc means data=outnlin noprint; var lpred; output out=outmeans mean=mlog; data outnlin; merge outnlin outmeans(drop=_type_ _freq_); retain xx 0; if _n_=1 the xx=mlog; if mlog=. then mlog=xx; data outnlin; set outnlin; dummy=0; fdot=exp(mlog); /* relax the convergence criterion for variance parameters */ proc nlin data=outnlin method=gauss converge=0.00001 &opt; parms theta = 0 to 1.5 by 0.05; trk = resid*((fdot/pred)**theta); model dummy = trk; der.theta = resid*((fdot/pred)**theta)*log(fdot/pred); output out=new parms=theta; run; /************************************************************************** * * * Set up the data set containing the new predicted values, etc, for use * * on the next iteration. Also, create the data set "parms" containing * * only the parameter estimates from the iteration just completed for * * printing * * * **************************************************************************/ data new; set new(keep = &xvar &yvar b1 b2 b3 b4 pred theta sigma rename=(b1=b1_new b2=b2_new b3=b3_new b4=b4_new pred=pred1 theta=theta1 sigma=sigma1)); data parms; set new(drop = &xvar &yvar); if _n_>1 then delete; proc print data=parms; data &dset; set new; %mend step3; /************************************************************************** * * * End of the macro definitions. Now the program to read the data and * * call the macros begins. * * * **************************************************************************/ data indo; input time conc; cards; 0.25 2.05 0.50 1.04 0.75 0.81 1.00 0.39 1.25 0.30 2.00 0.23 3.00 0.13 4.00 0.11 5.00 0.08 6.00 0.10 8.00 0.06 ; /************************************************************************** * * * Call the macro "glsalg" with the name of the indomethacin data set, * * the names of x (time) and y (conc), and the starting values. * * * **************************************************************************/ title1 "3-STEP GLS ALGORITHM APPLIED TO THE INDOMETHACIN DATA"; %glsalg(indo,time,conc,0.69,0.69,-1.6,-1.6);