* use the analyze.smoking.sas program for linear mixed modeling 
  of the smoking data via the LCVmixed macro after running the
  load.smoking.sas code to load in the smoking data set

  last changed 7/20/12

;

options pagesize=54 linesize=76 pageno=1 nodate;
title "Analyses of Smoking Data";
* switch to the second or third set of options for
  macro testing purposes
;
options nomlogic nomprint nosymbolgen nosource2;
* options mlogic mprint symbolgen nosource2;
* options mlogic mprint symbolgen source2;

%let macrolib=c:\macros;
%include "&macrolib\LCVmixed.20120720.sas";

%let kfold=15;
/*
* setting the number of folds;
* 15 folds generates the first local max in LCV scores over
  multiples of 5 folds for model based on either CS or AR(1)
  together with the full factorial fixed effects model
;
* subject 12 was first to have all 7 measurements;
* id of matchvar and time of withnvrs automatically added to
  class statment so not needed in classvrs
;
%let xvars=smoker time smoker*time;
%let ecovtype=ar(1);
%*let ecovtype=cs;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,
          ecovtype=&ecovtype,ecvindxs=12,rintrcpt=n,rvars=,
          procmod=y);
*/
/*
* alternate tests for zero fixed effects;
%let xvars=smoker time smoker*time;
%let ddfm=betwithin;
%let ecovtype=ar(1);
%*let ecovtype=cs;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=&ecovtype,
          ecvindxs=12,ddfm=&ddfm,rintrcpt=n,rvars=,procmod=y);
%let ddfm=contain;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=&ecovtype,
          ecvindxs=12,ddfm=&ddfm,rintrcpt=n,rvars=,procmod=y);
%let ddfm=kenwardroger;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=&ecovtype,
          ecvindxs=12,ddfm=&ddfm,rintrcpt=n,rvars=,procmod=y);
%let ddfm=residual;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=&ecovtype,
          ecvindxs=12,ddfm=&ddfm,rintrcpt=n,rvars=,procmod=y);
%let ddfm=satterth;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=&ecovtype,
          ecvindxs=12,ddfm=&ddfm,rintrcpt=n,rvars=,procmod=y);
*/
/*
* use to compute likelihood ratio tests with/without interaction;
%let xvars=smoker time smoker*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          ecvindxs=12,procmod=y);
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          ecvindxs=12,procmod=y);
*/
/*
* generating Table 1;
* basic covariance models with full factorial fixed effects model;
%let xvars=smoker time smoker*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker,withnvrs=time,matchvar=id,ecovtype=csh,
          ecvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          ecvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=arh(1),
          ecvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,
          ecovtype=sp(pow)(time),ecvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=un,
          ecvindxs=12,procmod=y);
*/
/*
* generating Table 2;
* choosing fixed effects using CS;
%let xvars=smoker time smoker*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%let xvars=smoker;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%let xvars=time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
*/
/*
* using regression fixed coefficient;
%let xvars=smoker time smoker*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%let xvars=time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%let xvars=smoker time time*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
%let xvars=time time*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=cs,
          ecvindxs=12,procmod=y);
*/
/*
* generating Table 3;
* using the additive fixed effects model with the best LCV score
  so far and choosing a random effects model considering only
  random coefficient models since random factor models
  are degenerate
;
* random effects with independent errors;
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=indep,
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,withnvrs=time,matchvar=id,
          rintrcpt=y,rvars=time2,ecovtype=indep,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=indep,
          rintrcpt=y,rvars=time2,rcovtype=un,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,withnvrs=time,matchvar=id,
          rintrcpt=y,rvars=time2 time2*time2,ecovtype=indep,rcvindxs=12,
          procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,withnvrs=time,matchvar=id,
          rintrcpt=y,rvars=time2 time2*time2,ecovtype=indep,rcvindxs=12,
          rcovtype=un,procmod=y);
*/
/*
* generating Table 4;
* random effects with dependent errors;
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=csh,
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=arh(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=sp(pow)(time),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=un,
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);

%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=csh,
          rintrcpt=y,rvars=time2,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=time2,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=arh(1),
          rintrcpt=y,rvars=time2,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=sp(pow)(time),
          rintrcpt=y,rvars=time2,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=un,
          rintrcpt=y,rvars=time2,rcvindxs=12,procmod=y);

%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=csh,
          rintrcpt=y,rvars=time2,rcovtype=un,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=time2,rcovtype=un,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=arh(1),
          rintrcpt=y,rvars=time2,rcovtype=un,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=sp(pow)(time),
          rintrcpt=y,rvars=time2,rcovtype=un,rcvindxs=12,procmod=y);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=un,
          rintrcpt=y,rvars=time2,rcovtype=un,rcvindxs=12,procmod=y);
*/
/*
* not reported in paper;
* assessing somemiss effects on random component;
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y,rcgrpvar=somemiss);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y,ecgrpvar=somemiss);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y,ecgrpvar=somemiss,
          rcgrpvar=somemiss);
*/
/*
* not reported in paper;
* assessing smoker effects on random component;
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y,rcgrpvar=smoker);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y,ecgrpvar=smoker);
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y,ecgrpvar=smoker,
          rcgrpvar=somemiss);
*/
/*
* generating Table 5;
* choosing fixed effects using the best covariance model so far;
%let xvars=smoker time smoker*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%let xvars=smoker;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%let xvars=time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
* using regression fixed coefficient;
%let xvars=smoker time smoker*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%let xvars=time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%let xvars=smoker time time*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
%let xvars=time time*time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y);
*/
/*
* generating estimates for the final selected model;
%let xvars=smoker time;
%LCVmixed(datain=smoking,yvar=FEV1,xvars=&xvars,foldcnt=&kfold,
          classvrs=smoker time,matchvar=id,withnvrs=time,ecovtype=ar(1),
          rintrcpt=y,rvars=,rcvindxs=12,procmod=y,predlist=y,
          plistby=smoker time);
*/
quit;