NOTE: Copyright (c) 2002-2008 by SAS Institute Inc., Cary, NC, USA. NOTE: SAS (r) Proprietary Software 9.2 (TS2M3) Licensed to UNIVERSIDAD NACIONAL DE ROSARIO, Site 70092939. NOTE: This session is executing on the W32_VSHOME platform. NOTE: SAS initialization used: real time 3.10 seconds cpu time 2.16 seconds 1 data base; 2 input aglo_id grupo tiempo vname $ val; 3 cards; NOTE: The data set WORK.BASE has 798 observations and 5 variables. NOTE: DATA statement used (Total process time): real time 0.04 seconds cpu time 0.03 seconds 802 ; 803 run; 804 805 data merlab; 806 set base; 807 run; NOTE: There were 798 observations read from the data set WORK.BASE. NOTE: The data set WORK.MERLAB has 798 observations and 5 variables. NOTE: DATA statement used (Total process time): real time 0.01 seconds cpu time 0.01 seconds 808 data merlab4; 809 set merlab; 810 if vname="Act" then vname1=1; 811 if vname="Des" then vname1=2; 812 if vname="Emp" then vname1=3; 813 *Efectos Aleatorios: solo pendiente; 814 if vname1=1 then slopeact=tiempo; else slopeact=0; 815 if vname1=2 then slopedes=tiempo; else slopedes=0; 816 if vname1=3 then slopeemp=tiempo; else slopeemp=0; 817 *Ordenadas al origen por variable respuesta y grupo; 818 if vname1=1 and grupo=1 then boact_1=1; else boact_1=0; 819 if vname1=1 and grupo=2 then boact_2=1; else boact_2=0; 820 if vname1=2 and grupo=1 then bodes_1=1; else bodes_1=0; 821 if vname1=2 and grupo=2 then bodes_2=1; else bodes_2=0; 822 if vname1=3 and grupo=1 then boemp_1=1; else boemp_1=0; 823 if vname1=3 and grupo=2 then boemp_2=1; else boemp_2=0; 824 *Pendiente según grupo y variable respuesta; 825 if vname1=1 and grupo=1 then slopeact_1=tiempo; else slopeact_1=0; 826 if vname1=1 and grupo=2 then slopeact_2=tiempo; else slopeact_2=0; 827 if vname1=2 and grupo=1 then slopedes_1=tiempo; else slopedes_1=0; 828 if vname1=2 and grupo=2 then slopedes_2=tiempo; else slopedes_2=0; 829 if vname1=3 and grupo=1 then slopeemp_1=tiempo; else slopeemp_1=0; 830 if vname1=3 and grupo=2 then slopeemp_2=tiempo; else slopeemp_2=0; 831 run; NOTE: There were 798 observations read from the data set WORK.MERLAB. NOTE: The data set WORK.MERLAB4 has 798 observations and 21 variables. NOTE: DATA statement used (Total process time): real time 0.06 seconds cpu time 0.06 seconds 832 833 proc sort data=merlab4; 834 by vname1 grupo aglo_id tiempo; 835 run; NOTE: There were 798 observations read from the data set WORK.MERLAB4. NOTE: The data set WORK.MERLAB4 has 798 observations and 21 variables. NOTE: PROCEDURE SORT used (Total process time): real time 0.03 seconds cpu time 0.03 seconds 836 %macro jointpair ( 837 data=merlab4, 838 respons=val, 839 fixed=%str(boact_1 boact_2 slopeact_1 slopeact_2 bodes_1 bodes_2 slopedes_1 slopedes_2 boemp_1 839! boemp_2 slopeemp_1 slopeemp_2 840 ), 841 random=%str(slopeact slopedes slopeemp 842 ), 843 outcome_ind=vname1, 844 id=aglo_id, 845 errorstructure=un, 846 timecl=tiempo, 847 mixedoptions=%str(method=ML) 848 ); 849 850 /* 851 0 = PREPARATION 852 853 I = ANALYSIS 854 855 II = CONSTRUCTION of COVARIANCE MATRICES of JOINT MODEL 856 A = COVARIANCE MATRIX of RANDOM EFFECTS 857 B = COVARIANCE MATRIX of ERROR COMPONENTS 858 859 III = MEAN EFFECTS of JOINT MODEL 860 861 IV = CALCULATION OF STANDARD ERRORS FOR FIXED EFFECTS 862 863 */ 864 865 866 /* 867 DATA=input dataset (REQUIRED). The structure of the dataset is of the PROC MIXED 868 datastructure type (thus a row for each observation) 869 OUTCOME_IND=(REQUIRED) variable in dataset with levels indicating 870 the various outcomes 871 RESPONS= (REQUIRED) respons variable 872 ID=(REQUIRED) variable indicating subjects. 873 874 FIXED=(REQUIRED, NUMERIC) string variable indicating ALL fixed effects of the joint model 875 (do not use comma's to separate entries). Since all models are fitted 876 in PROC MIXED using the noint option, the fixed effects should explicitly 877 contain intercepts. 878 879 RANDOM=(REQUIRED, NUMERIC) String variable indicating ALL random effects of the joint model. 880 881 TIMECL=(NUMERIC) variable indicating time. Used in the repeated statement to model 882 the variance of the error components (default=TIMECL) 883 884 ERRORSTRUCTURE=structure for 2x2covariance matrix of error components at specific 885 timepoint. Note that the correlation of the measurements within the 886 same outcome is modelled entirely by the random effects 887 (default=UN(1):uncorrelated error components with outcome specific variance). 888 correlated errors are obtained by specifying ERRORSTRUCTURE=UN 889 890 MIXEDOPTIONS=String variable with options used in the PROC MIXED statement 891 e.g. %str(method=ML) 892 893 /* created global macro variables */ 894 /*-----------------------------------------------------*/ 895 /* nobs = number of observations in dataset 896 /* noutcomes = number of outcomes in dataset 897 /* nfixed = number of fixed effects 898 /* nrandom = number of random effects 899 /* nfixedpair = number of fixed effects for current pair of outcomes 900 /* npair = number of pairs of outcomes 901 902 /*usefull output datasets*/ 903 /*-----------------------*/ 904 /* est_beta_robust= fixed effects (averaged over all pairs) of joint model, 905 including the robust standard errors (based on 905! pseudolikelihood theory) 906 /* est_d = covariance matrix (averaged over all pairs) of the 907 /* random effects of the joint model 908 /* est_dcor = correlation matrix (averaged over all pairs) of the 909 /* random effects of the joint model 910 /* est_r = covariance matrix (averaged over all pairs) of the 911 /* error components of the joint model 912 /* est_rcor = correlation matrix (averaged over all pairs) of the 913 /* error components of the joint model 914 /* 915 /* 916 /* 917 /* 918 /*intermediate datasets*/ 919 /***********************/ 920 /* _allfixed = dataset with a row for each fixed effect 921 /* _allfixedest = dataset with a row for each fixed effect 922 /* _allrandom = dataset with a row for each random effect 923 /* _analysis = dataset containing the observations of the input dataset 924 /* referring to the two outcomes of the current pair 925 /* _covariantie = dataset containing all estimates for the covariance parameters 926 /* for all pairs 927 /* _covariantieD = subset of _COVARIANTIE (covariance parameters of random efffects) 928 /* _covariantieError= subset of _COVARIANTIE (covariance parameters of error components) 929 /* 930 /* _fixedeffects = dataset with one row per outcome and a string variable indicating 931 /* all fixed effects for that outcome 932 /* _fixedest 933 /* 934 /* _input = copy of input dataset, adding an indicator for each outcome 934! starting with 1, 935 /* and increasing with one unit 936 /* _nfixedpair = 937 /* _noutcomespersubject= number of outcomes per subject in the original dataset 938 /* _nprevrandom = 939 /* _outcomes = dataset with a row for each outcome 940 /* _randomeffects = dataset with one row per outcome, a string variable indicating 941 /* all random effects for that outcome 942 /* _randomeffectsN = dataset with one row per outcome, a string variable indicating 943 /* all random effects for that outcome and NRANDOM, a 943! variable indicating 944 /* the number of random effects for that outcome 945 /* _resid = dataset with residuals (observed minus predicted mean) stacked per pair*/ 946 /* _tablepairs = dataset with one row for each possible pair of outcomes 947 /* _tmp = 948 /* _uit = 949 /* _uit1 = 950 /* _uitcov = dataset with estimates of covariance parameters of current pair 951 /* _uitF = dataset with estimates of fixed effects of current pair 952 /* 953 954 /* 955 /* Used sub macros */ 956 /*-----------------------------------------------------*/ 957 /***put elements from a string into a dataset***/ 958 %macro put_in_dataset(string=,n=,out=); 959 data &out;set _null_;run; 960 961 %local word count k; 962 %let count=1; 963 %let word=%qscan(&string,&count,%str( )); 964 %do k=1 %to %sysevalf(&n) %by 1; 965 %let WORD=%QSCAN(&string,&count,%str( )); 966 data _hulp; 967 length out $ 80; 968 out="&word";run; 969 data &out;set &out _hulp;run; 970 %let count=%EVAL(&count+1); 971 972 %end; 973 %mend; 974 975 976 /*****sort string variable****/ 977 978 %macro sorted(list=, maxItemLen=100); 979 980 %local i item mvars; 981 982 %let list=%lowcase(&list); 983 984 %do %while (1); 985 %let i = %eval(&i.+1); 986 %let item = %qscan(%superq(list),&i.); 987 %let len = %length(%superq(item)); 988 %if &len.=0 %then %goto out; 989 %local SORTED_&i.; 990 %let SORTED_&i. = &item.%qsysfunc(repeat(%str( ), &maxItemLen.-&len.- 991 1)); 992 %if &i.=1 %then 993 %let mvars = SORTED_1; 994 %else 995 %let mvars = &mvars., SORTED_&i.; 996 %end; 997 %out: 998 999 %syscall sortc(&mvars.); 1000 1001 %local r; 1002 %do i = 1 %to %eval(&i.-1); 1003 %let r = &r. &&SORTED_&i..; 1004 %end; 1005 1006 %*;&r. 1007 1008 %mend; 1009 /*****count number of elements in a string****/; 1010 %macro words(string); 1011 %local count word; 1012 %let count=1; 1013 %let word=%qscan(&string,&count,%str( )); 1014 %do %while(&word ne); 1015 %let count=%EVAL(&count+1); 1016 %let WORD=%QSCAN(&string,&count,%str( )); 1017 %end; 1018 %eval(&count-1) 1019 %mend; 1020 /***this macro selects out of a string of variable names a substring with 1021 names corresponding to nonzero columns in the dataset***/ 1022 %macro selectdesignmatrix(inputdata=,set=); 1023 /* inputdata= dataset */ 1024 /* set = string of variable names*/ 1025 /* substring with variable names is given in the macro variable EFFECTSSUB*/ 1026 1027 %global effectssub; 1028 ods select none; 1029 proc univariate data=&inputdata; 1030 var &set; 1031 ods output basicmeasures=_uit;run; 1032 ods select all; 1033 data _uit1(keep=varname number); 1034 set _uit; 1035 where locmeasure="Mean" and ((varvalue>0) or (varvalue=0 and locvalue^=0)); 1036 number=_N_; 1037 run; 1038 /*Find current number of effects*/ 1039 %let dsid=%sysfunc(open(_uit1)); 1040 %let neffectssub=%sysfunc(attrn(&dsid,nobs)); 1041 %let rc = %sysfunc (close (&dsid)); 1042 1043 %let effectssub=; 1044 %do j=1 %to &neffectssub; 1045 data hulpsub;set _uit1;where number=&j;call symput('add',varname);run; 1046 1047 %let effectssub=&effectssub%str( )&add; 1048 %end; 1049 %mend selectdesignmatrix; 1050 /*-----------------------------------------------------*/ 1051 /* 0. PREPARATION */ 1052 /*-----------------------------------------------------*/ 1053 %let startt = %sysfunc(time(),tod9.2); 1054 %put Macro Jointpair started at &sysday, &sysdate9, &startt; 1055 1056 %let timecl_Upper=%upcase(&timecl); 1057 1058 /*----------------------------------------------------------*/ 1059 /*0.1. check if all necessary input parameters are specified*/ 1060 /*----------------------------------------------------------*/ 1061 1062 %if %length(&data)=0 %then %do; 1063 %put YOU FORGOT TO SPECIFY AN INPUT DATASET; 1064 %goto endmac; %end; 1065 1066 %if %length(&outcome_ind)=0 %then %do; 1067 %put YOU FORGOT TO SPECIFY THE VARIABLE IN THE INPUT DATASET INDICATING; 1068 %put THE DIFFERENT OUTCOMES (OUTCOME_IND); 1069 %goto endmac; %end; 1070 1071 %if %length(&respons)=0 %then %do; 1072 %put YOU FORGOT TO SPECIFY THE RESPONS VARIABLE; 1073 %goto endmac; %end; 1074 1075 %if %length(&id)=0 %then %do; 1076 %put YOU FORGOT TO SPECIFY THE ID VARIABLE; 1077 %goto endmac; %end; 1078 1079 proc printto log='c:\abcdefg.tmp';run; 1080 proc sql; 1081 create table _nsubjects as 1082 select count(&id) as nsubjects from 1083 (select distinct &id from &data); 1084 quit; 1085 data _nsubjects;set _nsubjects; 1086 call symput('nsubjects', nsubjects); 1087 proc printto log=log;run; 1088 %put THERE ARE %trim(&nsubjects) SUBJECTS IN THE ORIGINAL DATASET; 1089 proc printto log='c:\abcdefg.tmp';run; 1090 1091 1092 1093 /*------------------------------------------------------------------*/ 1094 /* 0.2 REMOVE */ 1095 /* - all rows with missing info */ 1096 /* - all irrelevant columns 1096! */ 1097 /* - all rows corresponding to outcomes not involved */ 1098 /*------------------------------------------------------------------*/ 1099 1100 data _inputtmp (keep=&id &respons &outcome_ind &timecl &fixed &random); 1101 set &data; 1102 %let varlist=&id &outcome_ind &timecl &fixed &random; 1103 %put &varlist; 1104 %let vars=&varlist; 1105 %let var=; 1106 %let nextword=1; 1107 %let var=%scan(&vars, &nextword, %str( )); 1108 1109 %do %while(%length(&var) > 0); 1110 1111 if &var=. then delete; 1112 1113 %let nextword=%eval(&nextword+1); 1114 %let var=%scan(&vars,&nextword,%str( )); 1115 %end; 1116 1117 /*removes all rows corresponding to outcomes not involved*/ 1118 hulp=sum(of &fixed); 1119 if hulp=0 then delete; 1120 run; 1121 1122 1123 1124 1125 /*------------------------------------------------------------------*/ 1126 /* 0.3 create an indicator variable (_OUTCOMENUM) for the outcomes */ 1127 /*------------------------------------------------------------------*/ 1128 /*create an indicator for the outcomes; an indicator starting with 1, 1129 and increasing with one unit*/ 1130 proc sort data=_inputtmp;by &outcome_ind;run; 1131 data _inputtmp; 1132 set _inputtmp; 1133 retain _outcomenum(0); 1134 by &outcome_ind;if first.&outcome_ind then _outcomenum=_outcomenum+1; 1135 run; 1136 1137 /*Find number of outcomes*/ 1138 proc sql; 1139 create table _outcomes as 1140 select distinct _outcomenum from _inputtmp; 1141 quit; 1142 %let dsid=%sysfunc(open(_outcomes)); 1143 %let noutcomes=%sysfunc(attrn(&dsid,nobs)); 1144 %let rc = %sysfunc (close (&dsid)); 1145 1146 /*------------------------------------------------------------------*/ 1147 /* 0.4 Remove all subjects with more than 1 missing outcome */ 1148 /*------------------------------------------------------------------*/ 1149 proc sql; 1150 create table _noutcomespersubject as 1151 select &id, count(_outcomenum) as _noutcomespersubject from 1152 (select distinct &id, _outcomenum from _inputtmp 1153 where &respons^=.) 1154 group by &id; 1155 1156 proc sort data=_inputtmp;by &id;run; 1157 1158 data _inputtmp (drop=_noutcomespersubject); 1159 merge _inputtmp _noutcomespersubject; 1160 by &id; 1161 /*delete subject when more than 1 outcome is missing*/ 1162 if _noutcomespersubject<(&noutcomes-1) then delete; 1163 run; 1164 1165 /*-----------------------------------------------------*/ 1166 /* 0.5 create an indicator (_IDNUM) for the subjects */ 1167 /*-----------------------------------------------------*/ 1168 /*create an indicator for the patients; an indicator starting with 1, 1169 and increasing with one unit*/ 1170 proc sort data=_inputtmp;by &id;run; 1171 1172 data _input; 1173 set _inputtmp; 1174 retain _idnum (0); 1175 by &id;if first.&id then _idnum=_idnum+1; 1176 run; 1177 1178 proc sql; 1179 create table _nsubjectsR as 1180 select count(_idnum) as nsubjectsR from 1181 (select distinct _idnum from _input); 1182 quit; 1183 data _nsubjectsR;set _nsubjectsR; 1184 call symput('nsubjectsR', nsubjectsR); 1185 proc printto log=log;run; 1186 %put THERE ARE %trim(&nsubjectsR) SUBJECTS IN THE ANALYSED DATASET; 1187 proc printto log='c:\abcdefg.tmp';run; 1188 1189 1190 1191 1192 /*-----------------------------------------------------*/ 1193 /* 0.6 Create additional necessary information */ 1194 /*-----------------------------------------------------*/ 1195 1196 1197 /*Find number of observations in dataset*/ 1198 %let dsid=%sysfunc(open(_input)); 1199 %let nobs=%sysfunc(attrn(&dsid,nobs));/*number*/ 1200 %let rc = %sysfunc (close (&dsid)); 1201 1202 /*count total number of fixed effects and total number of random effects*/ 1203 %let nfixed=%words(&fixed); 1204 %let nrandom=%words(&random); 1205 %let randoms=%sorted(list=&random); 1206 1207 %put_in_dataset(string=&randoms,n=&nrandom,out=idrandom); 1208 data idrandom; 1209 set idrandom; 1210 idrandom=_N_; 1211 run; 1212 1213 1214 /*create a dataset TABLEPAIRS with all possible pairs of outcomes, 1215 columnames are OUTCOME1 and OUTCOME2*/ 1216 proc iml symsize=10000 worksize=10000; 1217 use _outcomes; 1218 read all var{_outcomenum} into _outcomenum; 1219 close _outcomes; 1220 N=nrow(_outcomenum); 1221 Npair=N*(N-1)/2; 1222 pairs=J(Npair,2,0); 1223 hulp=1; 1224 do i=1 to (N-1); 1225 do j=i+1 to N; 1226 pairs[hulp,1]=_outcomenum[i]; 1227 pairs[hulp,2]=_outcomenum[j]; 1228 hulp=hulp+1; 1229 end; 1230 end; 1231 c_name = { "outcome1", "outcome2" } ; 1232 create _tablepairs from pairs[colname=c_name]; 1233 append from pairs; 1234 close _tablepairs; 1235 quit; 1236 data _tablepairs;set _tablepairs;pair=_N_;run; 1237 1238 /*count number of pairs*/ 1239 %let dsid=%sysfunc(open(_tablepairs)); 1240 %let npair=%sysfunc(attrn(&dsid,nobs)); 1241 %let rc = %sysfunc (close (&dsid)); 1242 1243 /*create table FIXEDEFFECTS and RANDOMEFFECTS from input string fixed and random*/ 1244 1245 data _fixedeffects; 1246 set _null_;run; 1247 1248 data _randomeffects; 1249 set _null_;run; 1250 1251 %do i=1 %to &noutcomes; 1252 1253 data hulpcreate; 1254 set _input; 1255 where _outcomenum=&i; 1256 run; 1257 1258 %selectdesignmatrix(inputdata=hulpcreate, set=&fixed); 1259 data _tmp; 1260 format fixedeffects $250.; 1261 _outcomenum=&i; 1262 fixedeffects="&effectssub"; 1263 run; 1264 data _fixedeffects;set _fixedeffects _tmp;run; 1265 1266 %selectdesignmatrix(inputdata=hulpcreate, set=&random); 1267 data _tmp; 1268 format randomeffects $250.; 1269 _outcomenum=&i; 1270 randomeffects="&effectssub"; 1271 run; 1272 data _randomeffects;set _randomeffects _tmp;run; 1273 1274 %end; 1275 1276 /*count number of random effects for each outcome=NRANDOM)*/ 1277 data _randomeffectsN (drop=hulp hulp2 n); 1278 set _randomeffects; 1279 by _outcomenum; 1280 n=1;hulp2=1; 1281 do while (hulp2>0); 1282 hulp=scan(randomeffects,n); 1283 hulp2=(length(hulp)>1); 1284 n=n+hulp2; 1285 end; 1286 nrandom=n-1; 1287 run; 1288 1289 /*create dataset with a new observation for each random effect;ALLRANDOM*/ 1290 /*e.g. if two random effect for each outcome and 10 outcomes: 20 observations*/ 1291 data _allrandom (drop=hulp hulp2 n); 1292 set _randomeffects; 1293 by _outcomenum; 1294 n=1;hulp2=1; 1295 do while (hulp2>0); 1296 hulp=scan(randomeffects,n); 1297 hulp2=(length(hulp)>1); 1298 n=n+hulp2; 1299 allrandom=compress(hulp); 1300 output; 1301 end; 1302 run; 1303 data _allrandom;set _allrandom;by _outcomenum;if last._outcomenum then delete;run; 1304 data _allrandom;set _allrandom;numberrandom=_N_;run; 1305 1306 /*create dataset with a new observation for each fixed effect;_ALLFIXED*/ 1307 data _allfixed (drop=hulp hulp2 n); 1308 set _fixedeffects; 1309 by _outcomenum; 1310 n=1;hulp2=1; 1311 do while (hulp2>0); 1312 hulp=scan(fixedeffects,n); 1313 hulp2=(length(hulp)>1); 1314 n=n+hulp2; 1315 allfixed=compress(hulp); 1316 output; 1317 end; 1318 run; 1319 data _allfixed;set _allfixed;by _outcomenum;if last._outcomenum then delete;run; 1320 data _allfixed;set _allfixed;numberfixed=_N_;run; 1321 1322 data _nfixedpair; 1323 set _null_;run; 1324 1325 1326 /*-----------------------------------------------------*/ 1327 /* I. ANALYSIS */ 1328 /*-----------------------------------------------------*/ 1329 1330 /*loop over analyses of all possible pairs*/ 1331 1332 data _covariantie;set _null_;run; /*dataset with covariance parameters*/ 1333 1334 data _fixedest;set _null_;run; /*dataset with fixed estimates*/ 1335 1336 data _resid;set _null_;run; /*dataset with residuals*/ 1337 1338 data _randomlist;set _null_;run; 1339 1340 data _status;set _null_;run;/*information on convergence status*/ 1341 1342 %let residstring=; 1343 1344 1345 %do i=1 %to &npair; 1346 1347 data hulp;set _tablepairs;where pair=&i; 1348 call symput('outcome1',outcome1); 1349 call symput('outcome2',outcome2);run; 1350 %put outcome1=&outcome1; 1351 %put outcome2=&outcome2; 1352 1353 data _analysis; 1354 set _input; 1355 where _outcomenum in (&outcome1,&outcome2); 1356 run; 1357 1358 %selectdesignmatrix(inputdata=_analysis, set=&fixed); 1359 %let fixedpair=&effectssub; 1360 %selectdesignmatrix(inputdata=_analysis, set=&random); 1361 %let randompair=%sorted(list=&effectssub);/*sort random effects within a pair*/ 1362 %let nfixedpair=%words(&fixedpair); 1363 %let nrandompair=%words(&randompair); 1364 1365 proc printto log=log;run; 1366 %put pair &i on &npair; 1367 %put FIXEDPAIR &fixedpair; 1368 %put Number of FIXED effects in current PAIR: &nfixedpair; 1369 %put Number of RANDOM effects in current PAIR: &nrandompair; 1370 %put RANDOMPAIR &randompair; 1371 proc printto log='c:\abcdefg.tmp';run; 1372 1373 /*create dataset with random effect and pair as variables*/ 1374 /*put elements of string in dataset*/ 1375 %put_in_dataset(string=&randompair,n=&nrandompair,out=hulpje); 1376 data hulpje;set hulpje;pair=&i;run; 1377 data _randomlist;set _randomlist hulpje;run; 1378 1379 1380 ods noresults; 1381 ods listing close; 1382 proc mixed data=_analysis covtest noclprint &mixedoptions; 1383 class _idnum &TIMECL _outcomenum ; 1384 model &respons=&fixedpair/noint solution covb outpm=_uitresid; 1385 random &randompair/subject=_idnum type=un; 1386 repeated _outcomenum/subject=_idnum*&TIMECL type=&ERRORSTRUCTURE; 1387 ods output covparms=_uitCov solutionF=_uitF convergencestatus=uitstatus/*covB=covB&i */; 1388 run; 1389 1390 1391 data _uitCov;set _uitCov;outcome1=&outcome1;outcome2=&outcome2;pair=&i;run; 1392 data _covariantie;set _covariantie _uitcov;run; 1393 1394 data _uitF;set _uitF;outcome1=&outcome1;outcome2=&outcome2;pair=&i;run; 1395 data _fixedest;set _fixedest _uitF;run; 1396 1397 data _uitresid&i;set _uitresid;outcome1=&outcome1;outcome2=&outcome2;pair=&i;run; 1398 /*%let residstring=%str(&residstring%str( )%str(_uitresid)&i);*/ 1399 1400 data uitstatus;set uitstatus;outcome1=&outcome1;outcome2=&outcome2;pair=&i;run; 1401 data _status;set _status uitstatus;run; 1402 1403 data hulp (keep=effect numberfixed); 1404 set _allfixed; 1405 effect=allfixed;run; 1406 proc sort data=hulp;by effect;run; 1407 1408 data _hulpje;pair=&i;nfixedpair=&nfixedpair;run; 1409 data _nfixedpair; 1410 set _nfixedpair _hulpje;run; 1411 1412 /*voeg een rij toe voor ontbrekende eerste outcome of voor 1413 ontbrekende tweede outcome 1414 het zou efficienter zijn om dit te implementeren in de 1415 grotere input dataset*/ 1416 proc sort data=_uitresid&i;by _idnum &timecl _outcomenum;run; 1417 1418 data _addresid; 1419 set _uitresid&i; 1420 by _idnum &timecl _outcomenum; 1421 /*create empty record for first outcome if only second is present*/ 1422 if first.&timecl and _outcomenum=&outcome2 then do; 1423 _outcomenum=&outcome1;output;end; 1424 /*create empty record for second outcome if only first is present*/ 1425 else if last.&timecl and _outcomenum=&outcome1 then do; 1426 _outcomenum=&outcome2;output;end; 1427 run; 1428 1429 data _uitresid&i;set _uitresid&i _addresid;run; 1430 1431 proc sort data=_uitresid&i;by _idnum &timecl _outcomenum;run; 1432 1433 1434 %end; /*end loop over all pairs*/ 1435 1436 proc printto log=log;run; 1437 %put ALL PAIRWISE MODELS HAVE BEEN FITTED; 1438 proc printto log='c:\abcdefg.tmp';run; 1439 1440 1441 1442 /*------------------------------------------------------------*/ 1443 /* II.CONSTRUCTING COVARIANCE MATRICES OF JOINT MODEL */ 1444 /* 1444! */ 1445 /* (Matrices of covariance parameters of random effects and */ 1446 /* error components */ 1447 /* as averages of matrices of different pairs) */ 1448 /*------------------------------------------------------------*/ 1449 1450 proc sort data=_randomlist;by out;run; 1451 proc sort data=idrandom;by out;run; 1452 data _randomlist;merge _randomlist idrandom; 1453 by out;rename out=effect;run; 1454 proc sort;by pair effect;run; 1455 1456 1457 data _covariantie;set _covariantie;covnummer=_N_;run; 1458 1459 data hulp1 (keep=outcome1 nrandom1); 1460 set _randomeffectsN; 1461 outcome1=_outcomenum; 1462 nrandom1=nrandom;run; 1463 1464 data hulp2 (keep=outcome2 nrandom2); 1465 set _randomeffectsN; 1466 outcome2=_outcomenum; 1467 nrandom2=nrandom;run; 1468 1469 proc sort data=_covariantie;by outcome1;run; 1470 data _covariantie;merge _covariantie (in=x) hulp1;by outcome1;if x;run; 1471 proc sort data=_covariantie;by outcome2;run; 1472 data _covariantie;merge _covariantie (in=x) hulp2;by outcome2;if x;run; 1473 proc sort data=_covariantie;by outcome1 outcome2;run; 1474 proc sort data=_tablepairs;by outcome1 outcome2;run; 1475 data _covariantie;merge _covariantie (in=x) _tablepairs;by outcome1 outcome2;if x;run; 1476 1477 1478 /*determination of the number of random effects used for the previous 1479 outcome. This needs to be done separately for the first and second 1480 outcome of each pair (NPREVRANDOM1 NPREVRANDOM2*/ 1481 /*data _Nprevrandom (keep=_outcomenum Nprevrandom); 1482 set _randomeffectsN; 1483 retain nprevrandom (0); 1484 hulp=lag1(nrandom); 1485 if _N_=1 then nprevrandom=0;else Nprevrandom=Nprevrandom+hulp; 1486 run; 1487 data hulp1 (keep=outcome1 Nprevrandom1); 1488 set _Nprevrandom; 1489 outcome1=_outcomenum; 1490 Nprevrandom1=Nprevrandom; 1491 run; 1492 data hulp2 (keep=outcome2 Nprevrandom2); 1493 set _Nprevrandom; 1494 outcome2=_outcomenum; 1495 Nprevrandom2=Nprevrandom; 1496 run; 1497 proc sort data=_covariantie;by outcome1;run; 1498 data _covariantie;merge _covariantie (in=x) hulp1;by outcome1;if x;run; 1499 */ 1500 proc sort data=_covariantie;by outcome2;run; 1501 1502 /*extracting the row and column number(within ONE proc mixed!) 1503 of the covariance estimates, e.g. the 2 from UN(2,1)*/ 1504 data _covariantie(drop=hulp1 hulp2);merge _covariantie (in=x) hulp2;by outcome2;if x; 1505 r=0; 1506 c=0; 1507 hulp1=substr(covparm,4,1)+0; 1508 hulp2=substr(covparm,6,1)+0; 1509 r=r+hulp1; 1510 c=c+hulp2; 1511 run; 1512 1513 /*splitting up the covariance estimates in two parts: 1514 - COVARIANTIED=covariance matrix of random effects 1515 - COVARIANTIEERROR=covariance matrix of error components*/ 1516 1517 proc sort data=_covariantie;by covnummer;run; 1518 1519 data _covariantieD; 1520 set _covariantie; 1521 subject1=upcase(subject); 1522 if subject1^="_IDNUM" then delete;run; 1523 1524 1525 data _randomlist;set _randomlist;by pair; 1526 retain r;/*r=number within pair*/ 1527 if first.pair then r=1; 1528 else r=r+1; 1529 run; 1530 1531 proc sort data=_covariantieD;by pair r;run; 1532 proc sort data=_randomlist;by pair r;run; 1533 1534 data _covariantieD; 1535 merge _covariantieD _randomlist; 1536 by pair r; 1537 rename idrandom=idrandomr; 1538 run; 1539 1540 data hulpje (keep=pair c idrandomc); 1541 set _randomlist; 1542 idrandomc=idrandom;c=r;run; 1543 1544 proc sort data=_covariantieD;by pair c;run; 1545 proc sort data=hulpje;by pair c;run; 1546 data _covariantieD; 1547 merge _covariantieD hulpje; 1548 by pair c; 1549 run; 1550 1551 1552 data _covariantieError; 1553 set _covariantie; 1554 subject1=upcase(subject); 1555 if subject1^="_IDNUM*&timecl_upper" then delete;run; 1556 1557 1558 1559 /*II.A. Constructing covariance matrix for random effects*/ 1560 /*-------------------------------------------------------*/ 1561 proc iml symsize=10000 worksize=10000; 1562 use idrandom;read all var{out} into labels;close idrandom; 1563 use _covariantieD; 1564 read all var{estimate pair covnummer idrandomc idrandomr}; 1565 close _covariantieD; 1566 nbi=nrow(labels); 1567 covmatrix=J(nbi,nbi,0); 1568 max=nrow(estimate); 1569 1570 /*determining for each covariance estimate its place in the covmatrix*/ 1571 /*result is: ROW,COL*/ 1572 row=idrandomr; 1573 col=idrandomc; 1574 1575 1576 /*putting the estimates in the covariance matrix*/ 1577 count=J(nbi,nbi,0); /*matrix with counts of estimates*/ 1578 1579 do i=1 to max; 1580 covmatrix[row[i],col[i]]=covmatrix[row[i],col[i]]+estimate[i]; 1581 count[row[i],col[i]]=count[row[i],col[i]]+1; 1582 end; 1583 1584 /*filling in upper diagonal part*/ 1585 covmatrix=covmatrix+t(covmatrix)-diag(covmatrix); 1586 count=count+t(count)-diag(count); 1587 1588 /*obtaining mean estimates*/ 1589 est_D=covmatrix/count; 1590 1591 /*calculating correlation matrix of random effects*/ 1592 est_Dcor=J(nbi,nbi,0); 1593 do i=1 to nbi; 1594 do j=1 to nbi; 1595 est_Dcor[i,j]=est_D[i,j]/(sqrt(est_D[i,i]#est_D[j,j])); 1596 end; 1597 end; 1598 1599 tlabels=t(labels); 1600 1601 create est_Dcor from est_Dcor[colname=labels rowname=tlabels]; 1602 append from est_Dcor[rowname=tlabels]; 1603 create est_D from est_D[colname=labels rowname=tlabels]; 1604 append from est_D[rowname=tlabels]; 1605 quit; 1606 1607 /*II. B. constructing covariance matrix for error components*/ 1608 /*----------------------------------------------------------*/ 1609 proc iml symsize=10000 worksize=10000; 1610 use _covariantieError; 1611 read all var{estimate stderr outcome1 outcome2 pair covnummer r c}; 1612 close _covariantieError; 1613 use _outcomes; 1614 read all var{_outcomenum}; 1615 close _outcomes; 1616 1617 nerror=max(outcome2); 1618 covmatrix=J(nerror,nerror,0); 1619 1620 max=nrow(estimate); 1621 1622 /*determining for each covariance estimate its place in the covmatrix*/ 1623 /*result is: ROW,COL*/ 1624 row=J(max,1,0); 1625 col=J(max,1,0); 1626 do i=1 to max; 1627 if r[i]=1 then do;row[i]=outcome1[i];col[i]=outcome1[i];end; 1628 else if c[i]=2 then do;row[i]=outcome2[i];col[i]=outcome2[i];end; 1629 else do;row[i]=outcome2[i];col[i]=outcome1[i];end; 1630 end; 1631 1632 /*putting the estimates in the covariance matrix of the error 1633 components */ 1634 count=J(nerror,nerror,0); /*matrix with counts of estimates*/ 1635 do i=1 to max; 1636 covmatrix[row[i],col[i]]=covmatrix[row[i],col[i]]+estimate[i]; 1637 count[row[i],col[i]]=count[row[i],col[i]]+1; 1638 end; 1639 1640 /*filling in upper diagonal part*/ 1641 covmatrix=covmatrix+t(covmatrix)-diag(covmatrix); 1642 count=count+t(count)-diag(count); 1643 1644 1645 /*obtaining mean estimates */ 1646 est_R=covmatrix/count; 1647 1648 est_Rcor=J(nerror,nerror,0); 1649 do i=1 to nerror; 1650 do j=1 to nerror; 1651 est_Rcor[i,j]=est_R[i,j]/(sqrt(est_R[i,i]#est_R[j,j])); 1652 end; 1653 end; 1654 1655 labels=char(_outcomenum); 1656 outcome=t(labels); 1657 create est_Rcor from est_Rcor[colname=labels rowname=outcome]; 1658 append from est_Rcor[rowname=outcome]; 1659 create est_R from est_R[colname=labels rowname=outcome]; 1660 append from est_R[rowname=outcome]; 1661 quit; 1662 1663 /*-------------------------------------------------------------*/ 1664 /* III. CONSTRUCTING MEAN FIXED EFFECTS OF JOINT MODEL */ 1665 /*-------------------------------------------------------------*/ 1666 data _fixedest_sorted;set _fixedest;run; 1667 proc sort data=_fixedest_sorted;by effect;run; 1668 data _allfixedest; 1669 set _fixedest_sorted; 1670 by effect;if first.effect then output;run; 1671 1672 /*construction of mean effects as averages over all pairs */ 1673 /*-------------------------------------------------------------*/ 1674 1675 proc iml symsize=10000 worksize=10000; 1676 use _fixedest_sorted; 1677 read all var{effect estimate outcome1 outcome2}; 1678 close _fixedest_sorted; 1679 use _allfixedest; 1680 read all var{effect} into label; 1681 close _allfixedest; 1682 n_effects=nrow(label); 1683 n_estimates=nrow(estimate); 1684 1685 hulp1=J(n_effects,1,0); /*matrix for sum of fixed estimates*/ 1686 count1=J(n_effects,1,0); 1687 do i=1 to n_effects; 1688 do j=1 to n_estimates; 1689 if label[i]=effect[j] then do; 1690 hulp1[i]=hulp1[i]+estimate[j]; 1691 count1[i]=count1[i]+1; 1692 1693 end; 1694 end; 1695 end; 1696 est_beta=hulp1/count1; 1697 1698 1699 1700 create est_beta from est_beta[colname={'est'} rowname=label]; 1701 append from est_beta[rowname=label]; 1702 quit; 1703 1704 1705 /*----------------------------------------------------------*/ 1706 /* IV. Calculation of hessian and gradient */ 1707 /*----------------------------------------------------------*/ 1708 1709 1710 1711 /*give a unique number (seqnr) and a common number (nr_parameter)*/ 1712 data _fixedest1 ; 1713 set _fixedest ; 1714 seqnr=_N_;run; 1715 proc sort data=_fixedest1; 1716 by effect;run; 1717 data nr_parameter; 1718 set _fixedest1 (keep=effect seqnr); 1719 by effect; 1720 if first.effect then output;run; 1721 proc sort data=nr_parameter;by seqnr;run; 1722 data nr_parameter;set nr_parameter;nr_parameter=_N_;run; 1723 proc sort data=_fixedest1;by effect;run; 1724 proc sort data=nr_parameter;by effect;run; 1725 data _fixedest; 1726 merge _fixedest1 nr_parameter; 1727 by effect; 1728 run; 1729 proc sort data=_fixedest;by seqnr;run; 1730 proc sort data=nr_parameter;by nr_parameter;run; 1731 1732 1733 proc printto log=log;run; 1734 proc iml symsize=10000 worksize=10000; 1735 reset print; 1736 free gradient; 1737 1738 /*begin loop over pairs*/ 1739 1740 %do k=1 %to &npair;/*number of pairs*/ 1741 1742 %put pair &k; 1743 1744 free gradient&k; 1745 %let index=%sysevalf(&k); 1746 1747 1748 /*---------------------------------------*/ 1749 /* Construction of pair specific D */ 1750 /*---------------------------------------*/ 1751 /*Dpair(k)*/ 1752 1753 use _covariantieD where (pair=&k); 1754 read all var{estimate pair r c}; 1755 close _covariantieD; 1756 1757 1758 nbi=(-1+sqrt(1+8*nrow(estimate)))/2; /*number of random effects (quadratic equation)*/ 1759 covmatrix=J(nbi,nbi,0); 1760 max=nrow(estimate); 1761 1762 do i=1 to max; 1763 covmatrix[r[i],c[i]]=covmatrix[r[i],c[i]]+estimate[i]; 1764 end; 1765 1766 /*filling in upper diagonal part*/ 1767 hulp=diag(covmatrix); 1768 covmatrix=covmatrix+t(covmatrix)-diag(covmatrix); 1769 Dpair&index=covmatrix; 1770 1771 /*-------------------------------------------------------------*/ 1772 /* Construction of pair specific R (variance error components */ 1773 /*-------------------------------------------------------------*/ 1774 /*errorpair(k)*/ 1775 1776 use _covariantieerror where (pair=&k); 1777 read all var{estimate pair r c}; 1778 close _covariantieerror; 1779 1780 1781 nbi=(-1+sqrt(1+8*nrow(estimate)))/2; /*number of random effects (quadratic equation)*/ 1782 covmatrix=J(nbi,nbi,0); 1783 max=nrow(estimate); 1784 1785 do i=1 to max; 1786 covmatrix[r[i],c[i]]=covmatrix[r[i],c[i]]+estimate[i]; 1787 end; 1788 1789 /*filling in upper diagonal part*/ 1790 hulp=diag(covmatrix); 1791 covmatrix=covmatrix+t(covmatrix)-diag(covmatrix); 1792 Errorpair&k=covmatrix; 1793 1794 /*-------------------------------------------------------------*/ 1795 /* Construction of Subject-specific contributions */ 1796 /* to Gradient and Hessian */ 1797 /*-------------------------------------------------------------*/ 1798 /* for gradient: cross products BETWEEN pairs*/ 1799 /* for hessian: can be constructed within a pair*/ 1800 1801 /* Selecting subjectspecific Xi and Zi matrices and residi vector*/ 1802 /*-----------------------------------------------------------------*/ 1803 /*select levels of two outcomes: &one, &two*/ 1804 use _tablepairs where (pair=&k); 1805 read var{outcome1 outcome2}; 1806 close _tablepairs; 1807 %let one=outcome1; 1808 %let two=outcome2; 1809 1810 1811 /*select from dataset with residuals the information 1812 corresponding to current pair*/ 1813 1814 use _uitresid&k; 1815 read all var{&respons _outcomenum &timecl resid pair}; 1816 read all var{_idnum} into id; 1817 read all var{&fixed} into X; 1818 read all var{&random} into Z; 1819 close _uitresid&k; 1820 1821 1822 /*select fixed effects design matrix of current pair*/ 1823 test=X[##,]; /*calculation of cross-products*/ 1824 X=X[,loc(test>0)];/*selection of columns in X matrix, not containing all zero's*/ 1825 1826 /*select random effects design matrix of current pair*/ 1827 test=Z[##,]; /*calculation of cross-products*/ 1828 Z=Z[,loc(test>0)];/*selection of columns in Z matrix, not containing all zero's*/ 1829 1830 numobs=nrow(X); 1831 1832 /*Count number of observations per subject and store in VECLEN.*/ 1833 count=1; 1834 free veclen; /*empty veclen*/ 1835 do i=1 to (numobs-1); 1836 IF id[i]=id[i+1] then 1837 do; 1838 count=count+1; 1839 END; 1840 ELSE IF id[i]^=id[i+1] then 1841 DO; 1842 veclen=veclen//count; 1843 count=1; 1844 END; 1845 IF i=numobs-1 THEN veclen=veclen//count; 1846 END; 1847 1848 1849 y=&respons; 1850 nsubjects=nrow(veclen); 1851 1852 /*number of fixed effects*/ 1853 1854 p=ncol(X); 1855 1856 Hessian&k=J(p,p,0); 1857 1858 /*begin loop over subjects*/ 1859 do i=1 to nsubjects; 1860 1861 if i=1 then pos=1; 1862 1863 Zi=Z[pos:pos+veclen[i]-1,]; 1864 Xi=X[pos:pos+veclen[i]-1,]; 1865 yi=y[pos:pos+veclen[i]-1]; 1866 residi=resid[pos:pos+veclen[i]-1]; 1867 1868 /*create per subject its Ri matrix (dimension=ni x ni)*/ 1869 ni=nrow(yi); /*number of measurements for specific subject*/ 1870 hulp=diag(J(ni/2,1,1)); 1871 Ri=hulp@Errorpair&k;/*covariance matrix of error components*/ 1872 1873 /*row vector indicating rows corresponding to nonmissing y values*/ 1874 notmissy=loc(residi^=.); 1875 1876 /*row vector indicating rows corresponding to nonmissing X values*/ 1877 test=ncol(Xi); 1878 check=(Xi^=.); 1879 check2=check[,+]; 1880 notmissX=loc(check2=test); 1881 1882 1883 /*when for each observation (row) of a subject its y or x value is 1884 missing, go to next subject*/ 1885 if (ncol(notmissy)>0 & ncol(notmissx)>0) 1886 then do; 1887 notmiss=xsect(notmissy, notmissX); /*intersection of both row vectors*/ 1888 Zi=Zi[notmiss,];/*deleting rows*/ 1889 Xi=Xi[notmiss,]; 1890 yi=yi[notmiss,]; 1891 Residi=Residi[notmiss,]; 1892 Ri=Ri[notmiss,notmiss];/*delete rows and columns*/ 1893 1894 Wi=ginv(Zi*Dpair&k*t(Zi)+Ri); 1895 1896 XiWi=t(Xi)*Wi; 1897 XiWiXi=t(Xi)*Wi*Xi; 1898 1899 contributionHessian=t(Xi)*Wi*Xi; 1900 Hessian&k=Hessian&k+contributionHessian; /*p*p matrix*/ 1901 1902 /*create subject-specific vector with first derivatives*/ 1903 grad_ik=t(Xi)*Wi*Residi; /*pair and subjectspecific vector with first derivatives*/ 1904 gradient&k=gradient&k||grad_ik;/*horizontal concatenation: put subjects side by side: 1904! p_k x N matrix*/ 1905 end; 1906 pos=pos+veclen[i]; 1907 1908 end;/*end loop over subjects*/ 1909 1910 /*collect pairspecific p_k x ni matrices into one matrix: vertical concatenation*/ 1911 gradient=gradient//gradient&k; 1912 1913 /*collect pairspecific Hessians into one blockdiagonal matrix*/ 1914 if &k=1 then HESSIAN=hessian&k; 1915 else HESSIAN=block(HESSIAN,hessian&k); 1916 1917 1918 1919 /*end loop over pairs*/ 1920 %end; 1921 %put end loop over pairs; 1922 1923 1924 /*calculate crossproducts for the gradients*/ 1925 ptotal=nrow(gradient);/*total number of parameters*/ 1926 CPgradient=J(ptotal,ptotal,0); /*matrix with crossproducts of subjectspecific gradient vectors*/ 1927 1928 nsubjects=ncol(gradient); 1929 do i=1 to nsubjects; 1930 contribution=gradient[,i]*t(gradient[,i]); 1931 CPgradient=CPgradient+contribution; 1932 end; 1933 1934 CPgradient=CPgradient#1/nsubjects; /*expectation*/ 1935 Hessian=Hessian#1/nsubjects; /*expectation*/ 1936 1937 JKJ=ginv(hessian)*CPgradient*ginv(Hessian);/*covariance of sqrt(n)*(theta-theta_0)*/ 1938 cov=JKJ#1/nsubjects; /*covariance of (theta-theta_0)*/ 1939 1940 1941 /*combine the fixed effects parameters obtained from the different pairs 1942 and compute its covariance matrix*/ 1943 use _fixedest; 1944 read all var{effect estimate seqnr nr_parameter}; 1945 close _fixedest; 1946 1947 p=max(nr_parameter); /*totaal aantal definitieve parameters*/ 1948 pseq=max(seqnr);/*totaal aantal parameters komende uit de paren*/ 1949 1950 free A; /*A: p rows, pseq columns*/ 1951 do i=1 to p; 1952 currentnumber=0; 1953 hulp=J(pseq,1,0); 1954 do j=1 to pseq; 1955 if nr_parameter[j,]=i then do; 1956 currentnumber=currentnumber+1; /*aantal keren dat parameter is geschat over alle 1956! paren*/ 1957 hulp[j,]=1;end; 1958 end; 1959 hulp=hulp#1/currentnumber; 1960 A=A//t(hulp); 1961 end; 1962 mean=A*estimate; /*mean over fixed effects obtained in different pairs*/ 1963 print "A" A; 1964 print "estimate" estimate; 1965 1966 use nr_parameter; 1967 read all var{effect} into label; 1968 close nr_parameter; 1969 1970 covrobust=A*cov*t(A); 1971 serobust=vecdiag(covrobust)##(1/2); 1972 1973 B=mean||serobust; 1974 tlabel=t(label); 1975 create covariance_robust from covrobust[rowname=tlabel]; 1976 append from covrobust[rowname=tlabel]; 1977 create est_beta_robust from B[colname={beta se_robust} rowname=tlabel]; 1978 append from B[rowname=tlabel]; 1979 quit; 1980 1981 %endmac: 1982 1983 ods listing; 1984 ods results; 1985 1986 1987 proc printto log=log;run; 1988 1989 %let endt = %sysfunc(time(),tod9.2); 1990 %put Macro Jointpair ended at &sysday, &sysdate9, &endt; 1991 1992 %mend jointpair; 1993 1994 %jointpair ( 1995 data=merlab4, 1996 respons=val, 1997 fixed=%str(boact_1 boact_2 slopeact_1 slopeact_2 bodes_1 bodes_2 slopedes_1 slopedes_2 boemp_1 1997! boemp_2 slopeemp_1 slopeemp_2 1998 ), 1999 random=%str(slopeact slopedes slopeemp 2000 ), 2001 outcome_ind=vname1, 2002 id=aglo_id, 2003 errorstructure=un, 2004 timecl=tiempo, 2005 mixedoptions=%str(method=ML) 2006 ); Macro Jointpair started at Sunday, 01MAY2011, 0:00:41.1 NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.01 seconds cpu time 0.01 seconds THERE ARE 19 SUBJECTS IN THE ORIGINAL DATASET NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.00 seconds cpu time 0.00 seconds THERE ARE 19 SUBJECTS IN THE ANALYSED DATASET NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.00 seconds cpu time 0.00 seconds pair 1 on 3 FIXEDPAIR boact_1 boact_2 slopeact_1 slopeact_2 bodes_1 bodes_2 slopedes_1 slopedes_2 Number of FIXED effects in current PAIR: 8 Number of RANDOM effects in current PAIR: 2 RANDOMPAIR slopeact slopedes NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.00 seconds cpu time 0.00 seconds pair 2 on 3 FIXEDPAIR boact_1 boact_2 slopeact_1 slopeact_2 boemp_1 boemp_2 slopeemp_1 slopeemp_2 Number of FIXED effects in current PAIR: 8 Number of RANDOM effects in current PAIR: 2 RANDOMPAIR slopeact slopeemp NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.00 seconds cpu time 0.00 seconds pair 3 on 3 FIXEDPAIR bodes_1 bodes_2 slopedes_1 slopedes_2 boemp_1 boemp_2 slopeemp_1 slopeemp_2 Number of FIXED effects in current PAIR: 8 Number of RANDOM effects in current PAIR: 2 RANDOMPAIR slopedes slopeemp NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.00 seconds cpu time 0.00 seconds ALL PAIRWISE MODELS HAVE BEEN FITTED NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.00 seconds cpu time 0.00 seconds NOTE: Worksize = 10240000 NOTE: Symbol size = 10240000 NOTE: IML Ready pair 1 ERROR: ESTIMATE is not in the scope of variables for the data set. statement : READ at line 1994 column 1 NOTE: Writing HTML Body file: sashtml3.htm ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : nbi, nbi, *LIT1007 nbi 1 row 1 col (numeric) 0 nbi 1 row 1 col (numeric) 0 *LIT1007 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : DIAG at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : T at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 WARNING: Data set WORK._COVARIANTIEERROR is empty. statement : USE at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : nbi, nbi, *LIT1015 nbi 1 row 1 col (numeric) 0 nbi 1 row 1 col (numeric) 0 *LIT1015 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : DIAG at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : T at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 WARNING: Data set WORK._UITRESID1 is empty. statement : USE at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : [ at line 1994 column 1 operands : X, $SUB0008, X 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : > at line 1994 column 1 operands : test, *LIT1023 test 0 row 0 col (type ?, size 0) *LIT1023 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : [ at line 1994 column 1 operands : Z, $SUB0008, Z 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : > at line 1994 column 1 operands : test, *LIT1024 test 0 row 0 col (type ?, size 0) *LIT1024 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : val val 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : p, p, *LIT1033 p 1 row 1 col (numeric) 0 p 1 row 1 col (numeric) 0 *LIT1033 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : // at line 1994 column 1 operands : gradient, gradient1 gradient 0 row 0 col (type ?, size 0) gradient1 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : hessian1 Hessian1 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 pair 2 ERROR: ESTIMATE is not in the scope of variables for the data set. statement : READ at line 1994 column 1 ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : nbi, nbi, *LIT1056 nbi 1 row 1 col (numeric) 0 nbi 1 row 1 col (numeric) 0 *LIT1056 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : DIAG at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : T at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 WARNING: Data set WORK._COVARIANTIEERROR is empty. statement : USE at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : nbi, nbi, *LIT1064 nbi 1 row 1 col (numeric) 0 nbi 1 row 1 col (numeric) 0 *LIT1064 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : DIAG at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : T at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 WARNING: Data set WORK._UITRESID2 is empty. statement : USE at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : [ at line 1994 column 1 operands : X, $SUB0008, X 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : > at line 1994 column 1 operands : test, *LIT1072 test 0 row 0 col (type ?, size 0) *LIT1072 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : [ at line 1994 column 1 operands : Z, $SUB0008, Z 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : > at line 1994 column 1 operands : test, *LIT1073 test 0 row 0 col (type ?, size 0) *LIT1073 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : val val 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : p, p, *LIT1082 p 1 row 1 col (numeric) 0 p 1 row 1 col (numeric) 0 *LIT1082 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : // at line 1994 column 1 operands : gradient, gradient2 gradient 0 row 0 col (type ?, size 0) gradient2 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : BLOCK at line 1994 column 1 operands : HESSIAN, hessian2 HESSIAN 0 row 0 col (type ?, size 0) Hessian2 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 pair 3 ERROR: ESTIMATE is not in the scope of variables for the data set. statement : READ at line 1994 column 1 ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : nbi, nbi, *LIT1105 nbi 1 row 1 col (numeric) 0 nbi 1 row 1 col (numeric) 0 *LIT1105 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : DIAG at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : T at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 WARNING: Data set WORK._COVARIANTIEERROR is empty. statement : USE at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : nbi, nbi, *LIT1113 nbi 1 row 1 col (numeric) 0 nbi 1 row 1 col (numeric) 0 *LIT1113 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : DIAG at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : T at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : covmatrix covmatrix 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 WARNING: Data set WORK._UITRESID3 is empty. statement : USE at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : [ at line 1994 column 1 operands : X, $SUB0008, X 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : > at line 1994 column 1 operands : test, *LIT1121 test 0 row 0 col (type ?, size 0) *LIT1121 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : [ at line 1994 column 1 operands : Z, $SUB0008, Z 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : > at line 1994 column 1 operands : test, *LIT1122 test 0 row 0 col (type ?, size 0) *LIT1122 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MOVE at line 1994 column 1 operands : val val 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : p, p, *LIT1131 p 1 row 1 col (numeric) 0 p 1 row 1 col (numeric) 0 *LIT1131 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : // at line 1994 column 1 operands : gradient, gradient3 gradient 0 row 0 col (type ?, size 0) gradient3 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : BLOCK at line 1994 column 1 operands : HESSIAN, hessian3 HESSIAN 0 row 0 col (type ?, size 0) Hessian3 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 end loop over pairs ERROR: (execution) Invalid operand to operation. operation : J at line 1994 column 1 operands : ptotal, ptotal, *LIT1148 ptotal 1 row 1 col (numeric) 0 ptotal 1 row 1 col (numeric) 0 *LIT1148 1 row 1 col (numeric) 0 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : # at line 1994 column 1 operands : CPgradient, *LIT1150 CPgradient 0 row 0 col (type ?, size 0) *LIT1150 1 row 1 col (numeric) 1 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : # at line 1994 column 1 operands : Hessian, *LIT1151 HESSIAN 0 row 0 col (type ?, size 0) *LIT1151 1 row 1 col (numeric) 1 statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : GINV at line 1994 column 1 operands : hessian HESSIAN 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : # at line 1994 column 1 operands : JKJ, *LIT1152 JKJ 0 row 0 col (type ?, size 0) *LIT1152 1 row 1 col (numeric) 1 statement : ASSIGN at line 1994 column 1 WARNING: Data set WORK._FIXEDEST is empty. statement : USE at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MAX at line 1994 column 1 operands : nr_parameter nr_parameter 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : MAX at line 1994 column 1 operands : seqnr seqnr 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : * at line 1994 column 1 operands : A, estimate A 0 row 0 col (type ?, size 0) ESTIMATE 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: Matrix A has not been set to a value. statement : PRINT at line 1994 column 1 ERROR: Matrix ESTIMATE has not been set to a value. statement : PRINT at line 1994 column 1 WARNING: Data set WORK.NR_PARAMETER is empty. statement : USE at line 1994 column 1 WARNING: End of File reached. statement : READ at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : * at line 1994 column 1 operands : A, cov A 0 row 0 col (type ?, size 0) cov 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : VECDIAG at line 1994 column 1 operands : covrobust covrobust 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : || at line 1994 column 1 operands : mean, serobust mean 0 row 0 col (type ?, size 0) serobust 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: (execution) Matrix has not been set to a value. operation : T at line 1994 column 1 operands : label label 0 row 0 col (type ?, size 0) statement : ASSIGN at line 1994 column 1 ERROR: Matrix covrobust has not been set to a value. statement : CREATE at line 1994 column 1 ERROR: No data set is currently open for output. statement : APPEND at line 1994 column 1 ERROR: Matrix B has not been set to a value. statement : CREATE at line 1994 column 1 ERROR: No data set is currently open for output. statement : APPEND at line 1994 column 1 NOTE: Exiting IML. NOTE: The SAS System stopped processing this step because of errors. NOTE: PROCEDURE IML used (Total process time): real time 1.05 seconds cpu time 0.90 seconds NOTE: PROCEDURE PRINTTO used (Total process time): real time 0.00 seconds cpu time 0.00 seconds Macro Jointpair ended at Sunday, 01MAY2011, 0:00:51.2 2007 2008 title 'Results for Fixed effects'; 2009 proc print data=est_beta_robust noobs; ERROR: File WORK.EST_BETA_ROBUST.DATA does not exist. 2010 run; NOTE: The SAS System stopped processing this step because of errors. NOTE: PROCEDURE PRINT used (Total process time): real time 0.01 seconds cpu time 0.01 seconds 2011 2012 title 'Results for random effects'; 2013 title2 'Covariance Matrix'; 2014 proc print data=est_D noobs; 2015 run; NOTE: Writing HTML Body file: sashtml4.htm NOTE: There were 3 observations read from the data set WORK.EST_D. NOTE: PROCEDURE PRINT used (Total process time): real time 0.07 seconds cpu time 0.07 seconds 2016 title2 'Correlation Matrix'; 2017 proc print data=est_Dcor noobs; 2018 run; NOTE: Writing HTML Body file: sashtml5.htm NOTE: There were 3 observations read from the data set WORK.EST_DCOR. NOTE: PROCEDURE PRINT used (Total process time): real time 0.06 seconds cpu time 0.04 seconds 2019 2020 title 'Results for error components'; 2021 title2 'Covariance Matrix'; 2022 proc print data=est_R noobs; ERROR: File WORK.EST_R.DATA does not exist. 2023 run; NOTE: The SAS System stopped processing this step because of errors. NOTE: PROCEDURE PRINT used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 2024 title2 'Correlation Matrix'; 2025 proc print data=est_Rcor noobs; ERROR: File WORK.EST_RCOR.DATA does not exist. 2026 run; NOTE: The SAS System stopped processing this step because of errors. NOTE: PROCEDURE PRINT used (Total process time): real time 0.00 seconds cpu time 0.00 seconds 2027 2028 /*---------------------------------------------------------------------------*/ 2029 /* */ 2030 /* EXAMPLE 2031 /* */ 2032 /*---------------------------------------------------------------------------*/ 2033 /* 2034 2035 data example_hearing; 2036 infile 'c:\example_hearing.txt'; 2037 input id TRESH INTR500 R500v1 SLOPER500 INTR750 R750v1 SLOPER750 INTR1000 R1000v1 SLOPER1000 2037! INTR1500 R1500v1 SLOPER1500 INTR2000 R2000v1 SLOPER2000 time nr 2038 ; 2039 run; 2040 2041 %jointpair( 2042 data=basenueva4, 2043 respons=val, 2044 fixed=%str( 2045 2046 ), 2047 random=%str( 2048 2049 ), 2050 outcome_ind=vname, 2051 timecl=time, 2052 id=id, 2053 errorstructure=un, 2054 mixedoptions=%str(method=ML) 2055 ); 2056 2057 title 'Results for Fixed effects'; 2058 proc print data=est_beta_robust noobs; 2059 run; 2060 2061 title 'Results for random effects'; 2062 title2 'Covariance Matrix'; 2063 proc print data=est_D noobs; 2064 run; 2065 title2 'Correlation Matrix'; 2066 proc print data=est_Dcor noobs; 2067 run; 2068 2069 title 'Results for error components'; 2070 title2 'Covariance Matrix'; 2071 proc print data=est_R noobs; 2072 run; 2073 title2 'Correlation Matrix'; 2074 proc print data=est_Rcor noobs; 2075 run; 2076 2077 2078 /* 2079 In the example, the dataset has the following structure for one subject. 2080 2081 id nr time intR500 R500v1 slopeR500 intR1000 R1000v1 slopeR1000 ... 2082 ------------------------------------------------------------------- 2083 1 500 7 1 1 7 0 0 0 ... 2084 1 500 12 1 0 12 0 0 0 ... 2085 1 500 16 1 0 16 0 0 0 ... 2086 .... 2087 1 1000 7 0 0 0 1 1 7 ... 2088 1 1000 12 0 0 0 1 0 12 ... 2089 1 1000 16 0 0 0 1 0 16 ... 2090 .... 2091 2092 */