/* GPL-2 License: Copyright (C) 2015 Kalema George This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. */ %macro CorrPoisson( CovData=, /* dataset containing covariates for the mean*/ ID=, /* identification variable in 'CovData' dataset */ OrderVar=, /* variable with ordering of observations within subject */ Xcov=, /*covariates for the desired or target marginal mean */ Xcov2=, /*covariates for the hierarchical mean */ Alpha=, /* marginal mean parameters (without the reference categories) */ Class=, /* classification covariates among those specified in the Xcov macro argument*/ Class2=, /* classification covariates among those specified in the Xcov2 macro argument*/ outData=out, /* final dataset with generated Poisson counts */ param=glm, /* Parameterization for the classification variables in marginal means model. Options are EFFECT, GLM, ORDINAL, THERMOMETER, POLYNOMIAL, POLY, REFERENCE, REF, ORTHEFFECT, ORTHORDINAL, ORTHOTHERM, ORTHPOLY, ORTHREF. See documentation of param (logistic or genmod procedure) at http://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/ viewer.htm#statug_logistic_sect006.htm */ random=, /* covariates for the random effects e.g., for random intercept and slope model, specify 'random=time'. Random intercept model fit by default*/ meanNormalRE=, /* mean of the random effects. 0 by default */ seed=123, /* set seed */ desiredVarCov=, /* desired marginal variance-covariance matrix e.g., d11 d12 d22 or d11 d12 d13 d22 d23 d33 */ GammaRandEff=2, /* Either '0', '1' or '2' for the Gamma random effects {0=No Gamma RE, 1=Independent Gamma RE, 2=Correlated Gamma RE} */ NormalRandEff=2, /* Either '0', '1' or '2' for the normal random effects {0=No normal RE, 1=Independent normal RE, 2=Correlated normal RE} */ Estimates=1, /* Either '0'=NO or '1'=YES, to print the estimated unknowns */ INCLUDE= /* Path to the RinSASIML.sas file, must be specified */ ); %if %length(&include) = 0 or %sysfunc(fileexist("&include\RinSASIML.sas")) = 0 %then %do; %put ERROR: The path to RinSASIML.sas must to be correctly specified.; %put ERROR: Macro will not execute without the correct path specification.; %put ERROR: Please specify the correct path to the file RinSASIML.sas using the INCLUDE macro argument.; %goto endmac; %end; options minoperator mindelimiter=','; %if %length(&CovData)=0 or not(%sysfunc(exist(&CovData))) %then %do; %put ERROR: Macro argument 'CovData' needs to be correctly specified.; %put ERROR: Macro will stop executing NOW.....; %put ERROR: Please specify the dataset containing subject or cluster data; %goto endmac; %end; %put CovData=&CovData NormalRandEff=&NormalRandEff and GammaRandEff=&GammaRandEff; %let NormRE = Intercept &random; %let NormalRE=; %if &NormalRandEff in 0,1,2 and &GammaRandEff in 0,1,2 %then %do; %if &GammaRandEff=0 %then %let GammaRE=No; %if &GammaRandEff=1 %then %let GammaRE=Independent; %if &GammaRandEff=2 %then %let GammaRE=Correlated; %if &NormalRandEff=0 %then %let NormalRE=No; %if %length(&random) > 0 %then %do; %if &NormalRandEff=1 %then %let NormalRE=Independent; %if &NormalRandEff=2 %then %let NormalRE=Correlated; %end; %let text=&NormalRE Normal and &GammaRE Gamma random effects; %end; %else %do; %put ERROR: Permitted values for the macro arguments 'NormalRandEff' and 'GammaRandEff' are strictly {0, 1 or 2}; %put ERROR: Macro will stop executing NOW.....; %put ERROR: Please specify the correct values for these macro arguments.; %goto endmac; %end; ods exclude all; proc sort data=&CovData; by &id &OrderVar; run; /* number of observations per subject */ proc freq data=&CovData; tables &id / out=freq(drop=percent rename=(count=n)) noprint; run; /* Create X-tilde design matrix */ %desmat(data=&covdata, Xdes=&Xcov, Clss=&Class, prm=¶m, seedn=&seed, idn=&id, outdes=Xtilde); %if %length(&Xcov2)>0 %then %do; %desmat(data=&covdata, Xdes=&Xcov2, Clss=&Class2, prm=¶m, seedn=&seed, idn=&id, outdes=X); proc contents data=X out=checkx noprint varnum; run; proc sort data=checkx(keep=name varnum); by varnum;run; %end; proc contents data=Xtilde out=check noprint varnum; run; proc sort data=check(keep=name varnum); by varnum;run; proc sql noprint; select '"'||compress(name)||'"' into: vardesmat separated by ' ' from WORK.check; %if %length(&Xcov2)>0 %then %do; select '"'||compress(name)||'"' into: vardesmat2 separated by ' ' from WORK.checkx;;%end; select min(n), max(n) into: minni , : maxni from WORK.freq; select distinct "'"||"Y"||strip(LEFT(PUT(&OrderVar,3.)))||"'" into :varnames separated by ' ' from &CovData; select distinct "'"||"Gammacov_&OrderVar"||strip(LEFT(PUT(&OrderVar,3.)))||"'" into :Gammacov separated by ' ' from &CovData; select distinct "'"||"Gammacorr_&OrderVar"||strip(LEFT(PUT(&OrderVar,3.)))||"'" into :Gammacorr separated by ' ' from &CovData; quit; proc iml; covvec={&desiredVarCov}`; alpha={&alpha}`; seed=&seed; call randseed(seed); varNames = {&varnames}; printcov={&vardesmat}`; %if %length(&Xcov2)>0 %then %do; printcov2={&vardesmat2}`;;%end; use &CovData NOBS nobs; read all var{&OrderVar} into OrderVar; read all var{&id} into id; %if %length(&random) > 0 %then %do; read all var{&random} into random; %end; close &CovData; use Xtilde; read all into Xtilde; read all var{intercept} into intercept; close Xtilde; %if %length(&Xcov2)>0 %then %do; use X; read all into X; close X; %end; %else %do; X=Xtilde; %end; use freq NOBS K; read all var{n} into niVec; close freq; call symput('K',left(char(K))); %let K=&K; * Create desired marginal variance covariance matrix from input vector ; start Vmat(CovVec)global(Vmat); nRows=nrow(CovVec); *nRows=dimension(Vmat){dimension(Vmat) + 1}/2; * Dim(dimVmat)=solving 'nRows=dimension(Vmat){dimension(Vmat) + 1}/2'; dimVmat=0.5 # (-1 + sqrt(1 - 4 # (-2)# nRows)); Vmat=j(dimVmat,dimVmat); l = 1; do s = 1 to dimVmat; do t = s to dimVmat; Vmat[s,t] = CovVec[l]; if s 0 %then %do; || random %end;; q=ncol(Z); p=ncol(Xtilde); lambdaStar=Xtilde * alpha;; *dim(lambdaStar)=nobs X 1; M=diag(exp(lambdaStar)); invMi=solve(M, i(nrow(M)));* inv(M); Mx1=invMi*(Vblock-M)*invMi; invZtZ=ginv(Z`*Z); %if &NormalRandEff in 1,2 %then D=invZtZ*Z`*log(Mx1 + j(nobs,nobs,1))*Z*invZtZ;; %if &NormalRandEff=1 %then D=diag(D);; beta=ginv(X`*X)*X`*(lambdaStar %if &NormalRandEff in 1,2 %then - 0.5#((Z*D)#Z)[,+];); %if %length(&Xcov2) = 0 %then diff=alpha-beta;;; %if &NormalRandEff in 1,2 and &GammaRandEff in 0 %then %do; %if %length(&meanNormalRE) > 0 %then %do; meanNormalRE={&meanNormalRE}`; %end; %else %do; meanNormalRE=j(ncol(D),1,0); * mean of the random effects is set to 0 ; %end; bi = RANDNORMAL( K, meanNormalRE, D); %if &syserr %then %do; quit; %goto endmac; %end; b = shape(bi, q*K, 1); first=1; random=j(nobs,q); do i = 1 to K; last=cumm[i]; z_i = z[first:last,]; bigz=i(K)[i,]@z_i; Zblock=Zblock//bigz; random[first:last,]=bi[i,] @ j(niVec[i],1,1); first=last + 1; end; mu=exp(X*beta + Zblock*b); Y=ranpoi(&seed,mu); GammaCov=j(nobs,1,.); GammaCorr=j(nobs,1,.); shape=j(nobs,1,.); scale=j(nobs,1,.); gamV=j(nobs,1,.); %let GammaCov=GammaCov; %let GammaCorr=GammaCorr; %end; %else %if &NormalRandEff in 0,1,2 and &GammaRandEff in 1,2 %then %do; %if &NormalRandEff in 1,2 %then %do; first=1; last=cumm[1]; expZDZTblock=exp(-1#(Z[first:last, ]*D*Z[first:last, ]`)); do i = 2 to K; first=last + 1; last=cumm[i]; expZDZTblock = block(expZDZTblock, exp(-1#(Z[first:last, ]*D*Z[first:last, ]`)) ); end ; covGam=expZDZTblock#(Mx1 + JBlock) - JBlock; * covariance Gamma random effects; %end; %else covGam=(Mx1 + JBlock) - JBlock;; * covariance Gamma random effects; %if &GammaRandEff=1 %then %do; covGam=diag(covGam); scale=vecdiag(covGam); shape=1/scale; gamV=rangam(seed,shape); gamV=scale#gamV; CorGam=i(nobs); %end; %if &NormalRandEff in 1,2 %then %do; %if %length(&meanNormalRE) > 0 %then %do; meanNormalRE={&meanNormalRE}; %end; %else %do; meanNormalRE=j(ncol(D),1,0); * mean of the random effects is set to 0 ; %end; bi = RANDNORMAL( K, meanNormalRE, D); %if &syserr %then %do; quit; %goto endmac; %end; b = shape(bi, q*K, 1); %end; %if &GammaRandEff=2 %then %do; * calculate correlation from covariance matrix ; StdCovGam = sqrt(diag(CovGam));* standard deviations of each time point; StdCovGamInv = inv(StdCovGam); CorGam = StdCovGamInv * CovGam * StdCovGaminv; * correlation matrix ; %end; first=1; %if &NormalRandEff in 1,2 %then random=j(nobs,q);; GammaCov=j(nobs,&maxni,.); GammaCorr=j(nobs,&maxni,.); do i = 1 to K; last=cumm[i]; z_i = z[first:last,]; bigz=i(K)[i,]@z_i; Zblock=Zblock//bigz; %if &NormalRandEff in 1,2 %then random[first:last,]=bi[i,] @ j(niVec[i],1,1);; GammaCov[first:last,1:niVec[i]]=CovGam[first:last, first:last]; GammaCorr[first:last,1:niVec[i]]=CorGam[first:last, first:last]; first=last + 1; end; %if &GammaRandEff=2 %then %do; * create vector containing all the correlation parameters [output=Rvec]; niChoose2=niVec#(niVec-1)/2; Rvec=j(niChoose2[+],1); count=0; first=1; do y=1 to K; last=cumm[y]; do i=1 to (niVec[y]-1); do j=i+1 to niVec[y]; count=count+1; Rvec[count]=CorGam[first:last,first:last][i,j]; end; end; first=last + 1; end; * generate the gamma part from MGamma(shape=1/CovGam, scale=covGam) ; * Sampling from a Multivariate Gamma distribution done in R (copulas); scale=vecdiag(CovGam); shape=1/scale; *From restriction that the mean of the Gammas is 1.; call ExportMatrixToR(seed, "seed"); call ExportMatrixToR(Rvec, "Rvec"); call ExportMatrixToR(scale, "scale"); call ExportMatrixToR(shape, "shape"); call ExportMatrixToR(niVec, "niVec"); %inc "&include.\RinSASIML.sas";/* RUN R IN SAS IML */ call ImportMatrixFromR(gamV, "dat"); %end; mu=gamV # exp(X*beta %if &NormalRandEff in 1,2 %then + Zblock*b;); %end; %if &NormalRandEff in 0 and &GammaRandEff in 0 %then mu=lambdaStar;; Y=ranpoi(&seed,mu); data = id || OrderVar || %if &GammaRandEff in 1,2 %then GammaCov || GammaCorr || shape || scale || gamV ||; mu || Y %if &NormalRandEff in 1,2 %then || random;; create &outData from data[colname={"&id" "&OrderVar" %if &GammaRandEff in 1,2 %then &GammaCov &GammaCorr "Gamma_shape" "Gamma_scale" "theta_ij"; "mu" "Y" %if &NormalRandEff in 1,2 %then "randomIntercept";}]; append from data; ods exclude none; %if &estimates %then %do; *StdVmat = sqrt(diag(Vmat)); *StdVmatInv = inv(StdVmat); *CorrVmat = StdVmatInv * Vmat * StdVmatinv; * desired correlation matrix ; print , "Generate Correlated Poisson data from CM with", "&text", "****************************************************"; print ,"Sample size (K) = &K ", "minimum number of measurements per subject = %trim(&minni)", "maximum number of measurements per subject = %trim(&maxni)", "Normal random effects covariates = &normRE"; print "Given Mean parameters are:" printcov[label=''] "=" alpha[label=''],; print "Given variance-covariance matrix = " vmat[label=''][c=varNames r=varNames],,; %if %length(&Xcov2) = 0 %then %do; print printcov[label='Parameter'] alpha beta Diff %if &NormalRandEff in 1,2 %then %do; D %end;;; %end; %else %do; print printcov[label='X-tilde'] alpha printcov2[label='X'] beta %if &NormalRandEff in 1,2 %then D ;; %end;; %end; quit; data &outData; set &outData; format %if &GammaRandEff in 1,2 %then Gamma_shape Gamma_scale theta_ij; mu 8.2; run; data &outData; merge &CovData &outData; by &id &OrderVar; run; /* proc datasets nolist; delete check Xdes Xtilde %if %length(&Xcov2) > 0 %then X; freq names solutionuni test2; run; */ quit; %endmac: *quiting macro; %mend; /* Macro to create the X-tilde (and X) design matrices */ %macro desmat(data=, Xdes=, Clss=, prm=, seedn=, idn=, outdes=Xtilde); data &data; set &data; resp=ranpoi((&seedn + &idn),10); /* add pseudo response variable (resp) to aid creation of the design matrix*/ output; run; proc logistic data=&data outdesign=Xdes(drop=resp) outdesignonly; %if %length(&Clss)> 0 %then %do; Class &Clss / param=&prm;;%end; model resp=&Xdes; run; proc datasets lib=work memtype=data; modify Xdes; attrib _all_ label=' '; quit; proc genmod data=&data; /*obtain initial values*/ %if %length(&Clss)> 0 %then class &Clss / param=&prm;; model resp=&Xdes/ dist=poisson link=log; ods output ParameterEstimates=SolutionUni(drop=df LowerWaldCL UpperWaldCL ChiSq ProbChiSq where=(parameter ^= "Scale")); run;quit; data &data(drop=resp); set &data; run; %if %length(&Xdes) > 0 and %length(&Clss)>0 %then %do; data solutionuni; set solutionuni; varnum = _n_; run; proc contents data=Xdes order=varnum out=names(keep = varnum name) noprint; run; proc sort data=names; by varnum; run; data test2; merge names solutionuni; by varnum; run; proc sql noprint; select name into :droplist separated by ' ' from test2 where estimate=0; quit; data &outdes; * design matrix without columns for reference categories ; set Xdes; drop &droplist; run; %end; %else %do; data &outdes; * design matrix without columns for reference categories ; set Xdes; run; %end; %mend; /**************************************************************************************************************************** TEST CALL procedure: 1. Save RinSASIML.sas in a directory, e.g., C:=\TEMP. Then, specify this path in the macro using the 'include' argument 2. Ensure that SAS is launched with the -RLANG system option and that the connection between R and SAS is okay 3. Generate data for the covariates, for example, using our little macro below 4. Call macro corrpoisson to generate the Poisson counts %macro gentestdata(data=test, n=100, ni=4, equalni=1, seed=123); data &data; do id=1 to &n; if id<&n/2 then trt=0; else trt=1; do time=1 to %if &equalni=1 %then %do; ∋ %end; %else %do; (2+ &ni*ranuni(&seed + id)); %end; ; output; end; end; run; %mend; %gentestdata(data=temp, n=100, ni=4, equalni=1, seed=123); %CorrPoisson(CovData=temp, ID=id, OrderVar=time, Xcov=trt time trt*time, Alpha=1.521 0.437 -0.254 0.145, Class=trt, outData=out1, param=glm, random=, meanNormalRE=, seed=123, desiredVarCov=256 128 144 224 208 228 172 299 296 567, GammaRandEff=2, NormalRandEff=1, Estimates=1, include=C:\TEMP); %gentestdata(data=temp, n=3, ni=3, equalni=1, seed=123); %CorrPoisson(CovData=temp, ID=id, OrderVar=time, Xcov=trt time trt*time, Alpha=1.521 0.437 -0.254 0.145, Class=trt, outData=out1, param=glm, random=, meanNormalRE=, seed=123, desiredVarCov=256 128 144 224 208 228, GammaRandEff=2, NormalRandEff=1, Estimates=1, include=C:\TEMP); %gentestdata(data=temp, n=3, ni=3, equalni=1, seed=123); %CorrPoisson(CovData=temp, ID=id, OrderVar=time, Xcov=trt time trt*time, Alpha=1.521 0.437 -0.254 0.145, Class=trt, outData=out1, param=glm, random=time, meanNormalRE=, seed=123, desiredVarCov=256 128 144 224 208 228, GammaRandEff=2, NormalRandEff=2, Estimates=1, include=C:\TEMP); *****************************************************************************************************************************/