/*--------------------------------------------------------------*/ /* */ /* %mvarlst */ /* Make a variable list from the class list, model */ /* specification, and random specification. */ /* */ /*--------------------------------------------------------------*/ %macro mvarlst; %let varlst =; %let mdllst = &mdlspec; /*---get response variable---*/ %if %index(&response,/) %then %let varlst = %scan(&response,1,/) %scan(&response,2,/) &varlst; %else %let varlst = &response &varlst; /*---get fixed effects---*/ %if %index(&mdllst,@) %then %do; %let j = 1; %let mdl = &mdllst; %let mdllst=; %do %while(%length(%scan(&mdl,&j,' '))); %let var=%scan(&mdl,&j,' '); %if %index(&var,@) %then %do; %let b = %eval(%index(&var,@)-1); %let mdllst = &mdllst %substr(%quote(&var),1,&b); %end; %else %let mdllst = &mdllst &var; %let j = %eval(&j+1); %end; %end; %let iv = 1; %do %while (%length(%scan(&mdllst,&iv))); %let varlst = &varlst %scan(&mdllst,&iv); %let iv = %eval(&iv + 1); %end; /*---get random effects---*/ %let iv = 1; %do %while (%length(%scan(&rndlst,&iv))); %let temp = %scan(&rndlst,&iv); %if &temp ne INT and &temp ne INTERCEPT %then %let varlst = &varlst &temp; %let iv = %eval(&iv + 1); %end; /*---get repeated effects---*/ %let iv = 1; %do %while (%length(%scan(&replst,&iv))); %let temp = %scan(&replst,&iv); %if &temp ne DIAG %then %let varlst = &varlst &temp; %let iv = %eval(&iv + 1); %end; %let varlst = &varlst &class &id &freq &weight; %mend mvarlst; /*--------------------------------------------------------------*/ /* */ /* %trimlst */ /* Get rid of repetitions in a list */ /* */ /*--------------------------------------------------------------*/ %macro trimlst(name,lst); %let i1 = 1; %let tname =; %do %while (%length(%scan(&lst,&i1,%str( )))); %let first = %scan(&lst,&i1,%str( )); %let i2 = %eval(&i1 + 1); %do %while (%length(%scan(&lst,&i2,%str( )))); %let next = %scan(&lst,&i2,%str( )); %if %quote(&first) = %quote(&next) %then %let i2=10000; %else %let i2 = %eval(&i2 + 1); %end; %if (&i2<10000) %then %let tname = &tname &first; %let i1 = %eval(&i1 + 1); %end; %let &name = &tname; %mend trimlst; /*--------------------------------------------------------------*/ /* */ /* %remove */ /* Remove a word from a string */ /* */ /*--------------------------------------------------------------*/ %macro remove(sntnc,wrd); %let sentence=%str( )%nrbquote(&sntnc); %if &sentence^=%str( ) %then %do; %let word=%str( )%nrbquote(&wrd); %let answer=; %let i=%index(&sentence,&word); %if &i and &word^=%str( ) %then %do; %if &i>1 %then %let answer=%qsubstr(&sentence,1,&i-1); %let j=%eval(&i+%index(%qsubstr(&sentence,&i+1),%str( ))); %if &j>&i %then %let answer=&answer%qsubstr(&sentence,&j); %end; %else %let answer=&sentence; %unquote(&answer) %end; %mend remove; /*--------------------------------------------------------------*/ /* */ /* %errlink */ /* Create macro variables that contain the link, invlink, */ /* derivative, and the variance funtion for the following */ /* error distributions and link functions: */ /* */ /* error distn: normal, binomial, poisson, gamma, and */ /* inverse gaussian */ /* link func: logit, probit, complementary log-log, */ /* log-log, identity, power(), log, exp, and */ /* reciprocal. */ /* */ /* The user-defined specification is given by leaving the */ /* error distribution field blank and then giving the link, */ /* the derivative of the link, the inverse link, and the */ /* variance function. The parameters for each are: */ /* */ /* mu: variance function, link, and the derivative */ /* of the link */ /* eta: inverse link; */ /* */ /*--------------------------------------------------------------*/ %macro errlink; /*---error distributions: set variance and deviance functions---*/ %let exiterr = 0; %if %length(&error)=0 %then %do; %if %length(&errvar) and %length(&errdev) %then %let error=USER; %else %let error=BINOMIAL; %end; %if %length(&linkn) and %length(&link)=0 %then %let link=NLIN; %if %length(&linku) and %length(&link)=0 %then %let link=USER; %if &error=BINOMIAL or &error=B %then %do; %let errorfn=BINOMIAL; %let varform=mu*(1-mu); %let devform=_wght*2*(_y*log(_y/mu) + (1-_y)*log((1-_y)/(1-mu))); %if %length(&link)=0 %then %let link=LOGIT; %end; %else %if &error=POISSON or &error=P %then %do; %let errorfn=POISSON; %let varform=mu; %let devform=_wght*2*_y*log(_y/mu); %if %length(&link)=0 %then %let link=LOG; %end; %else %if &error=NORMAL or &error=N %then %do; %let errorfn=NORMAL; %let varform=1; %let devform=_wght*(_y-mu)**2; %if %length(&link)=0 %then %let link=IDENTITY; %end; %else %if &error=GAMMA or &error=G %then %do; %let errorfn=GAMMA; %let varform=mu**2; %let devform=_wght*-2*log(_y/mu); %if %length(&link)=0 %then %let link=RECIPROCAL; %end; %else %if &error=INVGAUSSIAN or &error=IG %then %do; %let errorfn=INVGAUSSIAN; %let varform=mu**3; %let devform=_wght*(((_y-mu)**2)/(_y*mu*mu)); %if %length(&link)=0 %then %let link=INVGAUSSIAN; %end; %else %if &error=USER or &error=U %then %do; %let errorfn=USER; %if %length(&errvar) %then %let varform=&errvar; %else %let exiterr = 1; %if %length(&errdev) %then %let devform=&errdev; %else %let exiterr = 1; %if %length(&link)=0 %then %let link=LOGIT; %end; /*---truncate link function, so we can match if a power link---*/ %if %length(&link)>5 %then %let trlink=%substr(&link,1,5); %else %let trlink=&link; /*---link functions; set eta, mu, and derivative formulas---*/ %if &trlink=LOGIT %then %do; %let linkfn=LOGIT; %let etaform=log(mu/(1-mu)); %let detaform=1/(mu*(1-mu)); %let muform=exp(eta)/(1+exp(eta)); %let dmuform=exp(eta)/(1+exp(eta))**2; %end; %else %if &trlink=PROBI %then %do; %let linkfn=PROBIT; %let etaform=probit(mu); %let detaform=(probit(mu+&numder)-probit(mu-&numder))/ (2*&numder); /* %let detaform=(probit(min(1,mu+&numder*(1+abs(mu)))) - probit(max(1,mu-&numder*(1+abs(mu))))) / (2*&numder*(1+abs(mu))); */ %let muform=probnorm(eta); %let dmuform=(probnorm(eta+&numder)-probnorm(eta-&numder))/ (2*&numder); %end; %else %if &trlink=CLOGL %then %do; %let linkfn=COMPLEMENTARY LOG LOG; %let etaform=log(-log(1-mu)); %let detaform=-1/((1-mu)*log(1-mu)); %let muform=1-exp(-exp(eta)); %let dmuform=exp(-exp(eta))*exp(eta); %end; %else %if &trlink=LOGLO %then %do; %let linkfn=LOG LOG; %let etaform=-log(-log(mu)); %let detaform=-1/(mu*log(mu)); %let muform=exp(-exp(-eta)); %let dmuform=exp(-exp(-eta))*exp(-eta); %end; %else %if &trlink=IDENT %then %do; %let linkfn=IDENTITY; %let etaform=mu; %let detaform=1; %let muform=eta; %let dmuform=1; %end; %else %if &trlink=LOG %then %do; %let linkfn=LOG; %let etaform=log(mu); %let detaform=1/mu; %let muform=exp(eta); %let dmuform=exp(eta); %end; %else %if &trlink=EXP %then %do; %let linkfn=EXPONENTIAL; %let etaform=exp(mu); %let detaform=exp(mu); %let muform=log(eta); %let dmuform=1/eta; %end; %else %if &trlink=RECIP %then %do; %let linkfn=INVERSE; %let etaform=1/mu; %let detaform=-1/mu**2; %let muform=1/eta; %let dmuform=-1/eta**2; %end; %else %if &trlink=POWER %then %do; %let linklen = %eval(%length(&link)-7); %let expon=%substr(&link,7,&linklen); %let linkfn=POWER(&expon); %let etaform=mu**(&expon); %let detaform=(&expon)*mu**(&expon-1); /* %let detaform=((mu+&numder)**(&expon)-(mu-&numder)**(&expon))/ (2*&numder); %let detaform=((mu+&numder*(1+abs(mu)))**(&expon) - (mu-&numder*(1+abs(mu)))**(&expon))/ (2*&numder*(1+abs(mu))); */ %let muform=eta**(1/(&expon)); %let dmuform=(1/(&expon))*eta**(1/(&expon)-1); %end; %else %if &trlink=INVGA %then %do; %let linkfn=POWER(-2); %let etaform=mu**(-2); %let detaform=-2*mu**(-3); %let muform=eta**(-1/2); %let dmuform=(-1/2)*eta**(-3/2); %end; %else %if &trlink=BOXCO %then %do; %let linkfn=BOX-COX; %let linklen = %eval(%length(&link)-8); %let expon=%substr(&link,8,&linklen); %let etaform=(mu**(&expon)-1)/(&expon); %let detaform=mu**((&expon)-1); %let muform=((&expon)*eta + 1)**(1/(&expon)); %let dmuform=((&expon)*eta + 1)**(1/(&expon)-1); %end; %else %if &trlink=USER %then %do; %let linkfn=USER; %if %length(&linku) %then %let etaform=&linku; %else %let exiterr = 1; %if %length(&linkud) %then %let detaform=&linkud; %else %let exiterr = 1; %if %length(&linkui) %then %let muform=&linkui; %else %let exiterr = 1; %if %length(&linkuid) %then %let dmuform=&linkuid; %else %let exiterr = 1; %end; %else %if &trlink=NLIN %then %do; %let linkfn=NONLINEAR; %if %length(&linkn) %then %let nlinform=&linkn; %else %let exiterr = 1; %if %length(&linknd) %then %let nlinder=&linknd; %else %let exiterr = 1; %end; %if %index(&options,DEBUG) %then %do; %put options = &options; %put intopt = &intopt; %put varlst = &varlst; %put error = &errorfn; %put variance = &varform; %put deviance = &devform; %put link: eta = &etaform; %put dlink: deta = &detaform; %put invlink: mu = &muform; %end; %mend errlink; /*--------------------------------------------------------------*/ /* */ /* %init */ /* Sets the initial values for the iterations. */ /* */ /*--------------------------------------------------------------*/ %macro init; %if %index(&options,DEBUG) %then %put Initializing.; %else %if not %index(&options,NOTES) %then %do; options nonotes nodate nonumber; %end; /*---turn off printing to log and output---*/ %global _print_; %if not %index(&options,PRINTALL) %then %let _print_ = off; %let off = &offset; %if %index(&intopt,NLIN) %then %do; /*---determine number of parameters---*/ %let nb = 0; %let i = 1; %do %while(%index(&linkni,B&i)); %let nb = %eval(&nb + 1); %let i = %eval(&i + 1); %end; %let nu = 0; %let ns = 0; %let nus = 0; %let varlst = &varlst one mu; %end; data _ds; set %unquote(&data); %if not %index(&options,INITIAL) %then %do; /*---move away from parameter space boundary for the binomial error situation---*/ %if %index(&response,/) %then %do; mu = (%scan(&response,1,/) + &cf)/(%scan(&response,2,/) + 2*&cf); _wght = %scan(&response,2,/) ; %end; %else %if &errorfn=BINOMIAL %then %do; mu = (&response + &cf)/(1 + 2*&cf); _wght = 1; %end; %else %do; mu = &response + &cf; _wght = 1; %end; %if %length(&weight) %then %do; _wght = &weight * _wght; %end; _y = &response; var = &varform; _offset = &off; %if %index(&intopt,NLIN) %then %do; array b{&nb} b1-b&nb; array db{&nb} db1-db&nb; one = 1; %do i = 1 %to &nb; %let idx = %index(&linkni,B&i); b&i = %scan(%substr(&linkni,&idx),2,'=' ' '); %end; &nlinform _z = _y - mu; &nlinder do i = 1 to &nb; _z = _z + db{i}*b{i}; end; _w = _wght / var; %end; %else %do; eta = &etaform ; deta = &detaform ; _w = _wght / ((deta**2)*(var)); _z = (_y-mu)*deta + eta - _offset; %end; %if %length(&freq) %then %do; do i = 1 to &freq; if i=1 then _orig='y'; else _orig='n'; output; end; %end; %else %do; _orig='y'; %end; if (_w = .) then _w = 1; /* keep _y _z _w _offset _wght _orig &varlst; */ %end; run; %if %index(&options,PRINTDATA) %then %do; proc print; run; %end; %let iter = 0; %mend init; /*--------------------------------------------------------------*/ /* */ /* %newdata */ /* Create the new data set with the updated values */ /* */ /*--------------------------------------------------------------*/ %macro newdata; /*---save previous parameter estimates---*/ data _oldsoln; set _soln; run; /*---save previous estimates of covariance matrix---*/ data _oldcov; set _cov; %let covsaved = 1; run; %if %index(&options,DEBUG) %then %put Creating new pseudo data.; /*---make new data set---*/ data _ds; %if %index(&intopt,NLIN) %then %do; set _ds; array b{&nb} b1-b&nb; array db{&nb} db1-db&nb; /* array u{&nu,&ns} u1-u&nus; array du{&nu,&ns} du1-du&nus; array _r{&nu} _r1-_rν */ &nlinform _z = _y - mu; &nlinder do i = 1 to &nb; _z = _z + db{i}*b{i}; end; /* do i = 1 to ν _r{i} = du{i,subject}; _z = _z + _r{i}*u{i,subject}; end; end; */ var = &varform; _w = _wght / var; /* keep _y _z _w _offset _wght _orig &varlst db1-db&nb one x; */ %end; %else %do; set _pred; eta = _pred_ + _offset; mu = &muform; deta = &detaform; var = &varform; %if %index(&options,HYBRID) %then %do; eta = _predm_ + _offset; mu = &muform; /* deta = &detaform; var = &varform; eta = _pred_ + _offset; mu = &muform; */ %end; _w = _wght / ((deta**2)*(var)); _z = (_y - mu)*deta + eta - _offset; /* keep _y _z _w _offset _wght _orig &varlst; */ %end; run; %if %index(&options,PRINTDATA) %then %do; proc print; run; %end; %mend newdata; /*--------------------------------------------------------------*/ /* */ /* %mixed */ /* Calculate parameter estimates using PROC MIXED. */ /* */ /*--------------------------------------------------------------*/ %macro mixed; %if %index(&options,DEBUG) %then %put Calling PROC MIXED.; %let mivque0 = 0; %again: /*---get rid of predicted data set---*/ proc datasets lib=work nolist; delete _pred; run; /*---use mivque0 if did not converge the first time---*/ %if (&mivque0 = 1) %then %do; /*---save the original method---*/ %let procopt0 = &procopt; %let procopt = &procopt METHOD=MIVQUE0; %end; proc mixed data=_ds &procopt; %if %length(&class) %then %do; class &class ; %if not %index(&procopt,NOCLPRINT) %then %do; make 'classlevels' out=_class; %end; %end; model _z = %unquote(&mdlspec) %unquote(&mdlopt); %if %length(%scan(&mdlspec,1)) and not %index(&mdlopt,NOTEST) %then %do; make 'tests' out=_tests; %end; make 'solutionf' out=_soln; %if %index(&options,HYBRID) %then %do; make 'predmeans' out=_predm noprint; make 'predicted' out=_pred noprint; %end; %else %if %index(&options,MQL) %then %do; make 'predmeans' out=_pred noprint; %end; %else %do; make 'predicted' out=_pred noprint; %end; make 'covparms' out=_cov; %if %index(&options,FITTING) %then %do; make 'fitting' out=_fitted; %end; weight _w; %if %length(&spec) %then %do; %unquote(&spec) %if %index(&intopt,SOLNR) %then %do; make 'solutionr' out=_solnr noprint; %end; %if %index(&spec,ESTIMATE) %then %do; make 'estimate' out=_est; %end; %if %index(&spec,CONTRAST) %then %do; make 'contrast' out=_con; %end; %if %index(&spec,LSMEANS) %then %do; make 'lsmeans' out=_lsm; %if %index(&intopt,LSMDIFF) %then %do; make 'diffs' out=_diff; %end; %if %index(&intopt,LSMSLICE) %then %do; make 'slices' out=_slice; %end; %end; %end; %if &covsaved=1 and not %index(&options,NOPREV) and not %index(&procopt,MIVQUE0) and (&mivque0 = 0) %then %str(parms / pdata=_oldcov &parmopt2;); %else %if (%length(&parmspec) or %length(&parmopt)) and (&mivque0 = 0) %then %do; parms &parmspec &parmopt ; %end; %if %index(&procopt,INFO) %then %do; make 'model' out=_model; make 'dimensions' out=_dim; %end; id _y _offset _wght _orig &varlst; run; %if %index(&options,PRINTDATA) %then %do; proc print data=_pred; run; %end; /*---check for convergence, if not, then run again with method=mivque0---*/ %if (&mivque0 = 1) %then %do; %let procopt = &procopt0; %let mivque0 = 0; %end; %else %do; %let there = no; data _null_; set _pred; call symput('there','yes'); run; %if ("&there" = "no") %then %do; %put Computing MIVQUE0 estimates in iteration &iter because; %put %str( )PROC MIXED did not converge.; %let mivque0 = 1; %goto again; %end; %end; /*---set up for hybrid Taylor series---*/ %if %index(&options,HYBRID) %then %do; data _predm; set _predm; _predm_ = _pred_; keep _predm_; run; data _pred; merge _predm _pred; run; %end; /*---merge in new estimates of b and u for nonlinear link---*/ %if %index(&intopt,NLIN) %then %do; %if (&nb) %then %do; proc transpose data=_soln out=_beta; var _est_; run; data _beta; set _beta; array b{&nb} b1-b&nb; array col{&nb} col1-col&nb; do i = 1 to &nb; b{i} = col{i}; end; one = 1; keep one b1-b&nb; run; data _ds; set _ds; drop b1-b&nb; run; data _ds; merge _ds _beta; by one; run; %end; /*---this isn't finished---*/ /* %if (&nu) %then %do; proc transpose data=_solnr out=_blup; var _est_; run; data _blup; set _blup; array u{&nus} u1-u&nus; array col{&nus} col1-col&nus; do i = 1 to &nus; u{i} = col{i}; end; one = 1; keep one u1-u&nus; run; %end; */ %end; %mend mixed; /*--------------------------------------------------------------*/ /* */ /* %compare */ /* Compare the last two parameter estimates to check for */ /* convergence if no random components, else compare */ /* estimates of covariance matrix. */ /* */ /*--------------------------------------------------------------*/ %macro compare; /*---Use relative difference of parameter estimates and covariance matrix as a measure of convergence---*/ %let crit = 0; %compit(soln,_est_); %compit(cov,est); data _null_; crit = &crit; if (crit<&converge) then conv = 1; else conv = 0; call symput('conv',conv); run; %mend compare; %macro compit(type,est); /*---save convergence information in conv and crit---*/ data _compare; merge _old&type(rename=(&est=oldest)) _&type end=last; retain crit &crit; denom = (abs(oldest) + abs(&est))/2; if (denom > &converge) then do; reldiff = abs(oldest - &est) / denom; crit = max(crit,reldiff); end; output; if last then do; call symput('crit',left(crit)); end; run; %if %index(&options,DEBUG) %then %do; proc print data=_compare; run; %end; %mend compit; /*--------------------------------------------------------------*/ /* */ /* %iterate */ /* Iteration process */ /* */ /*--------------------------------------------------------------*/ %macro iterate; %let conv = 0; %let iter = 1; %do %while(&iter <= &maxit); %newdata %mixed %compare %if not %index(&options,NOITPRINT) %then %do; %if (&iter=1) %then %do; %put %str( ) GLIMMIX Iteration History; %put; %put Iteration Convergence criterion; %end; %if (&iter<10) %then %put %str( ) &iter &crit; %else %if (&iter<100) %then %put %str( ) &iter &crit; %else %put %str( ) &iter &crit; %end; %let iter = %eval(&iter+1); %if (&conv=1) %then %let iter=%eval(&maxit+1); %end; %mend iterate; /*--------------------------------------------------------------*/ /* */ /* %compile */ /* Compile the macro results. */ /* */ /*--------------------------------------------------------------*/ %macro compile; /*---get variance estimate---*/ %let scale = 1; data _null_; set _cov; if covparm='Residual' then call symput('scale',est); run; /*---calculate deviance and Pearson Chi-Squared---*/ %if %index(&intopt,NLIN) %then %do; data _pred; set _pred; one = 1; run; data _pred; merge _pred _beta; by one; run; %end; data _stats; set _pred end=last; retain deviance 0 pearson 0 nn 0; if ((_y ne .) and (_pred_ ne .)) then do; _y = _y + 1e-10*(_y=0) - 1e-10*(_y=1); eta = _pred_ + _offset; %if %index(&intopt,NLIN) %then %do; &nlinform; %end; %else %do; mu = &muform; %end; deviance + &devform; pearson + _wght * ((_y-mu)**2/(&varform)); nn + 1; end; if last; if last then do; call symput('deviance',deviance); call symput('n',nn); end; keep deviance pearson; run; data _stats; length descript $35; set _stats; descript = 'Deviance'; value = deviance; output; descript = 'Scaled Deviance'; value = deviance / &scale; output; descript = 'Pearson Chi-Square'; value = pearson; output; descript = 'Scaled Pearson Chi-Square'; value = pearson / &scale; output; descript = 'Extra-Dispersion Scale'; value = &scale; output; keep descript value; label descript = 'Description' value = 'Value'; run; /*---ESTIMATE statement results---*/ %if %index(&spec,ESTIMATE) %then %do; data _est; set _est; %if not %index(&intopt,NLIN) %then %do; eta = est; mu = &muform; label mu = 'Mu'; drop eta; %if %index(&intopt,ESTCL) %then %do; eta = lower; lowermu = &muform; eta = upper; uppermu = &muform; label lowermu = 'LowerMu' uppermu = 'UpperMu'; %end; %end; run; %end; /*---least squares means---*/ %if %index(&spec,LSMEANS) %then %do; data _lsm; set _lsm; %if not %index(&intopt,NLIN) %then %do; eta = _lsmean_; mu = &muform; label mu = 'Mu'; %if %index(&intopt,LSMCL) %then %do; eta = _lower_; lowermu = &muform; eta = _upper_; uppermu = &muform; label lowermu = 'LowerMu' uppermu = 'UpperMu'; %end; drop eta; %end; run; %end; %mend compile; /*--------------------------------------------------------------*/ /* */ /* %printout */ /* Print out the macro results. */ /* */ /*--------------------------------------------------------------*/ %macro printout; %if %index(&procopt,INFO) %then %do; &titlen 'GLIMMIX Model Information'; proc print data=_model noobs label; run; %end; %if %length(&class) and not %index(&procopt,NOCLPRINT) %then %do; &titlen 'Class Level Information'; proc print data=_class noobs label; run; %end; %if %index(&procopt,INFO) %then %do; &titlen 'Dimensions'; proc print data=_dim noobs label; run; %end; %if %index(&intopt,COVMOD) %then %do; &titlen 'Covariance Parameter Estimates'; proc print data=_cov noobs label; where covparm ne 'Residual'; run; %end; &titlen 'GLIMMIX Model Statistics'; proc print data=_stats noobs label; format value 10.4; run; %if %index(&intopt,SOLNF) %then %do; &titlen 'Parameter Estimates'; proc print data=_soln noobs label; format _est_ 10.4 _se_ 10.4; run; %end; %if %index(&intopt,SOLNR) %then %do; &titlen 'Random Effects Estimates'; proc print data=_solnr noobs label; format _est_ 10.4 _sepred_ 10.4; run; %end; %if %index(&options,FITTING) %then %do; &titlen 'Model Fitting Information for _Z Weighted by _W'; proc print data=_fitted noobs label; run; %end; %if %length(%scan(&mdlspec,1)) and not %index(&mdlopt,NOTEST) %then %do; &titlen 'Tests of Fixed Effects'; proc print data=_tests noobs label; run; %end; %if %index(&spec,CONTRAST) %then %do; &titlen 'CONTRAST Statement Results'; proc print data=_con noobs label; run; %end; %if %index(&spec,ESTIMATE) %then %do; &titlen 'ESTIMATE Statement Results'; proc print data=_est noobs label; format est 10.4 se 10.4; %if not %index(&intopt,NLIN) %then %do; format mu 10.4; %if %index(&intopt,ESTCL) %then %do; format lower 10.4 upper 10.4 lowermu 10.4 uppermu 10.4; %end; %end; run; %end; %if %index(&spec,LSMEANS) %then %do; &titlen 'Least Squares Means'; proc print data=_lsm noobs label; format _lsmean_ 10.4 _se_ 10.4; %if not %index(&intopt,NLIN) %then %do; format mu 10.4; %if %index(&intopt,LSMCL) %then %do; format _lower_ 10.4 _upper_ 10.4 lowermu 10.4 uppermu 10.4; %end; %end; run; %if %index(&spec,DIFF) or %index(&spec,ADJUST) %then %do; &titlen 'Differences of Least Squares Means'; proc print data=_diff noobs label; format _diff_ 10.4 _se_ 10.4; run; %end; %if %index(&spec,SLICE) %then %do; &titlen 'Tests of Effect Slices'; proc print data=_slice noobs label; run; %end; %end; %mend printout; /*--------------------------------------------------------------*/ /* */ /* %outinfo */ /* Make an output data set. */ /* */ /*--------------------------------------------------------------*/ %macro outinfo(output); data &outfile; set _pred; eta = _pred_ + _offset; stdereta = _sepred_; %if (&outalpha = 0.05) %then %do; lowereta = _l95_ + _offset; uppereta = _u95_ + _offset; %end; %else %do; lowereta = _lower_ + _offset; uppereta = _upper_ + _offset; %end; mu = &muform; dmu = &dmuform; stdermu = abs(dmu)*stdereta; eta = lowereta; lowermu = &muform; eta = uppereta; uppermu = &muform; eta = _pred_ + _offset; var = &varform; resraw = _y - mu; reschi = (_y-mu)/sqrt(&scale * var); deta = &detaform; _w = _wght /((deta**2)*(var)); _z = (_y - mu)*deta + eta - _offset; if _orig='y'; run; %mend outinfo; /*--------------------------------------------------------------*/ /* */ /* %glimmix */ /* Put it all together */ /* */ /*--------------------------------------------------------------*/ %macro glimmix(data=,procopt=,stmts=,weight=,freq=, error=BINOMIAL,errvar=,errdev=,link=,linku=,linkud=, linkui=,linkuid=,linkn=,linknd=,linkni=,numder=1e-5,cf=0.5, converge=1e-8,maxit=20,offset=0,out=,outalpha=0.05,options=); /*---default data set---*/ %if %bquote(&data)= %then %let data=&syslast; %let exiterr = 0; %let covsaved = 0; %let there = no; /*---check that it is there---*/ data _null_; set &data; call symput('there','yes'); run; %if ("&there" = "no") %then %do; %let exiterr = 1; %goto exit; %end; /*---change to uppercase---*/ %let data = %qupcase(&data); %let procopt = %qupcase(&procopt); %let weight = %qupcase(&weight); %let freq = %qupcase(&freq); %let stmts = %qupcase(&stmts); %let error = %qupcase(&error); %let errvar = %qupcase(&errvar); %let errdev = %qupcase(&errdev); %let link = %qupcase(&link); %let linkn = %qupcase(&linkn); %let linknd = %qupcase(&linknd); %let linkni = %qupcase(&linkni); %let linku = %qupcase(&linku); %let linkud = %qupcase(&linkud); %let linkui = %qupcase(&linkui); %let linkuid = %qupcase(&linkuid); %let offset = %qupcase(&offset); %let outfile = %qupcase(&out); %let outalpha = %qupcase(&outalpha); %let options = %qupcase(&options); /*---title---*/ %let ntitle=0; %let titlesas=; data _null_; set sashelp.vtitle; if (number=1) and (text="The SAS System") then call symput("titlesas",right(text)); else call symput("ntitle",left(put(number,2.))); if number=10 then call symput("title10",text); run; %if &ntitle=10 %then %let titlen = title10; %else %let titlen = title%eval(&ntitle + 1); /*---loop through statements and extract information---*/ %let spec = ; %let class = ; %let parms = ; %let id = ; %let rndlst = ; %let replst = ; %let intopt = ; %let iv = 1; %do %while (%length(%scan(&stmts,&iv,;))); %let stmt = %qscan(&stmts,&iv,%str(;)); %let first = %qscan(&stmt,1); %let fn = %eval(%index(&stmt,&first) + %length(&first)); /*---check RANDOM options and extract random effects list---*/ %if %index(&first,RANDOM) or %index(&first,REPEATED) %then %do; %let intopt = &intopt COVMOD; %end; %if %index(&first,RANDOM) %then %do; %let i = %index(&stmt,/); %if &i = 0 %then %let i = %length(&stmt); %else %do; %let rndopt = %substr(&stmt,&i); %if %index(&rndopt,%str( S )) or %index(&rndopt,SOLUTION) or %index(&rndopt,CL) or %index(&rndopt,ALPHA) %then %do; %let intopt = &intopt SOLNR; %end; %end; %let rndlst = %substr(&stmt,&fn,%eval(&i-&fn+1)); %end; %if %index(&first,REPEATED) %then %do; %let i = %index(&stmt,/); %if &i = 0 %then %let i = %length(&stmt); %else %do; %let repopt = %substr(&stmt,&i); %let j = %index(&repopt,EXP); %if &j ne 0 %then %do; %let k = %index(&repopt,EXP); %let repexp = %bquote(%substr(&repopt,&k)); %let k1 = %index(&repexp,%str(%()); %let k1 = %eval(&k1 + 1); %let repexp1 = %bquote(%substr(&repexp,&k1)); %let k2 = %index(&repexp1,%str(%))); %let replst = &replst %substr(&repexp,&k1, %eval(&k2-1)); %end; %let j = %index(&repopt,TYPE=SP); %if &j ne 0 %then %do; %let k = %index(&repopt,TYPE=SP); %let repexp = %bquote(%substr(&repopt,&k)); %let k1 = %eval(%index(&repexp,%str(%))) + 1); %let repexp = %bquote(%substr(&repexp,&k1)); %let k2 = %eval(%index(&repexp, %str(%()) + 1); %let k3 = %eval(%index(&repexp, %str(%))) - 1); %let replst = &replst %substr(&repexp,&k2, %eval(&k3-&k2+1)); %end; %end; %let j = %eval(&i-&fn); %if &j > 0 %then %let replst = &replst %substr(&stmt,&fn,&j); %end; /*---check ESTIMATE and LSMEANS options---*/ %if %index(&first,ESTIMATE) and (%index(&stmt,CL) or %index(&stmt,ALPHA)) %then %do; %let intopt = &intopt ESTCL; %end; %if %index(&first,LSMEANS) %then %do; %if %index(&stmt,CL) or %index(&stmt,ALPHA) %then %do; %let intopt = &intopt LSMCL; %end; %if %index(&stmt,DIFF) or %index(&stmt,ADJ) %then %do; %let intopt = &intopt LSMDIFF; %end; %if %index(&stmt,SLICE) %then %do; %let intopt = &intopt LSMSLICE; %end; %end; /*---save statements---*/ %if %index(&first,CLASS) %then %do; %let class = %qsubstr(&stmt,&fn); %end; %else %if %index(&first,MODEL) %then %do; %let model = %qsubstr(&stmt,&fn); %end; %else %if %index(&first,ID) %then %do; %let id = %qsubstr(&stmt,&fn); %end; %else %if %index(&first,PARMS) %then %do; %let parms = %qsubstr(&stmt,&fn); %end; %else %let spec = &spec &stmt %str(;) ; %let iv = %eval(&iv + 1); %end; /*---get response, model specification, and model options---*/ %let response = %scan(&model,1,=); %let eqidx = %eval(%index(&model,=)+1); %if (&eqidx > %length(&model)) %then %let mdl = %str(); %else %let mdl = %str( ) %qsubstr(&model,&eqidx); %if %index(&mdl,/) %then %do; %let mdlspec = %qscan(&mdl,1,/); %let mdlopt = / %qscan(&mdl,2,/); %if %index(&mdlopt,%str( S )) or %index(&mdlopt,SOLUTION) or %index(&mdlopt,CL) or %index(&mdlopt,ALPHA) %then %do; %let intopt = &intopt SOLNF; %end; %if %index(&options,HYBRID) %then %do; %let mdlopt = &mdlopt S PM P; %end; %else %if %index(&options,MQL) %then %do; %let mdlopt = &mdlopt S PM; %end; %else %do; %let mdlopt = &mdlopt S P; %end; %end; %else %do; %let mdlspec = &mdl; %if %index(&options,HYBRID) %then %do; %let mdlopt = / S PM P; %end; %else %if %index(&options,MQL) %then %do; %let mdlopt = / S PM; %end; %else %do; %let mdlopt = / S P; %end; %end; %let mdlopt = &mdlopt alphap = &outalpha; /*---add an @ sign if it is missing---*/ %if %index(&mdlspec,|) %then %do; %let mdl=&mdlspec; %let mdlspec=; %let i=1; %do %while(%length(%scan(&mdl,&i,' '))); %let mdlterm = %scan(&mdl,&i,' '); %if %index(&mdlterm,|) and %index(&mdlterm,@)=0 %then %do; %let j=1; %do %while(%length(%scan(&mdlterm,&j,|))); %let j = %eval(&j+1); %end; %let atvalue = %eval(&j-1); %let mdlterm = &mdlterm.@&atvalue; %end; %let mdlspec = &mdlspec &mdlterm; %let i = %eval(&i+1); %end; %end; /*---get parms specification, and parms options---*/ %let parmopt2 = ; %if %index(&parms,/) %then %do; %let parmspec = %scan(&parms,1,/); %let parmopt = / %scan(&parms,2,/); %if not %index(&options,NOPREV) %then %do; %let parmopt2 = %scan(&parms,2,/); %let parmopt2 = %remove(%quote(&parmopt2),PARMSDATA); %let parmopt2 = %remove(%quote(&parmopt2),PDATA); %let parmopt2 = %remove(%quote(&parmopt2),OLS); %let parmopt2 = %remove(%quote(&parmopt2),RATIOS); %end; %end; %else %do; %let parmspec = %qupcase(&parms); %let parmopt = ; %end; %if %length(&linkn) %then %let intopt = &intopt NLIN; /*---create local variables---*/ %local varlst errorfn linkfn varform devform etaform detaform muform dmuform nlinform nlinder deviance scale n nb nu ns nus crit conv cf intopt iter; /*---get variable list and trim it---*/ %mvarlst %trimlst(varlst,&varlst); /*---set error and link function macro variables---*/ %errlink /*---print header---*/ %if not %index(&options,NOPRINT) %then %do; %if %index(&data,.)=0 %then %let data=WORK.&data; %put; %put %str( ) The GLIMMIX Macro; %put; %put Data Set : &data; %put Error Distribution : &errorfn; %put Link Function : &linkfn; %put Response Variable : &response; %if %length(&weight) %then %put Weight : &weight; %if %length(&freq) %then %put Frequency : &freq; %put; %put; %end; /*---initialize iteration starting values---*/ %init /*---run first estimates for convergence tests---*/ %mixed /*---iterate until convergence---*/ %iterate /*---turn the printing back on---*/ %let _print_ = on; /*---final PROC MIXED run---*/ %if %index(&options,PRINTLAST) %then %do; %put; %put Output from last PROC MIXED run:; %put; %mixed %end; /*---stop if did not converge---*/ %if &conv ne 1 %then %do; %put GLIMMIX did not converge.; %end; /*---otherwise compile and print results---*/ %else %do; %compile %if not %index(&options,NOPRINT) %then %do; %printout %end; /*---create output data set---*/ %if %length(&outfile) %then %do; %outinfo(%quote(&outfile)); %end; %end; %exit:; %if not %index(&options,NOTES) %then %do; options notes date number; %end; /*---title cleanup---*/ &titlen; %if %length(&titlesas) %then %do; title "The SAS System"; %end; %else %if &ntitle=10 %then title10 &title10; %mend glimmix;