Our shop does this a lot.  If you don't have an indicator variable for each dataset, you will need to create one.  Then you can do something like the following (NOTE: This is a different PROC, so look at the various ODS table names for PHREG to get what you want..): This part generates the by variables and merges them back onto the original dataset: /*Process anlsdat byvar*/ proc sort data=anlsdatin out=byvars (keep = paramid parameter) nodupkey; by paramid parameter; run; data byvars; set byvars end=eof; by paramid parameter; byvar = _N_; if eof then call symputx('LastByVar', _N_); run; proc sort data=anlsdatin; by paramid parameter; run; data anlsdat; merge anlsdatin byvars; by paramid parameter; run; This extracts the fit information from the PROC (GLIMMIX here):   %macro runtimes(); %do b=2 %to &LastByVar; proc sort data=anlsdat (where=(byvar=&b.)) out=anlsdat_&b.; by byvar paramid parameter grp_no sex animal; run; %macro runtests(covtype=, it=); proc glimmix data=anlsdat_&b. method = laplace; nloptions maxiter = 5000; by byvar paramid parameter; class grp_no sex animal; model value = grp_no sex grp_no*sex / dist = gamma ddfm=bw; random grp_no / type = &covtype. subject=animal; ods output fitstatistics = _rm_fit_&b._&it._1 convergencestatus=_convdat1_&b._&it._1; run; data _rm_fit_&b._&it._1; length covtype $10 value 8 descr $100; %if %sysfunc(exist(_rm_fit_&b._&it._1)) %then %do; set _rm_fit_&b._&it._1; %end; covtype = "&covtype."; Distribution="Gamma"; run; %mend runtests;/* %runtests(covtype=CSH, it=1); %runtests(covtype=CHOL, it=2); data fit&b.; set _rm_fit_&b.:; run; proc datasets nolist lib=work; delete _rm_fit_&b.:; run; proc sort data=fit&b. out=checkcov&b. (where=(lowcase(descr) = "aicc (smaller is better)" or lowcase(descr) = "aicc (smaller is better)")); by value; run; %let cs = 0; %let dist = ; /* data checkcov&b.; %if %sysfunc(exist(checkcov&b.)) %then %do; set checkcov&b.; %end; if _N_ = 1 then do; call symputx('cs', covtype); call symputx('dist', Distribution); selected = 1; end; run; */ %if "&cs" = "0" %then %do; %let cs = UN(1); %let dist = Gamma; %end; %mend runtimes   If you use this, be sure that all the %do loops are closed, as I have edited out some things we would rather keep proprietary.   SteveDenham       ... View more

Suppose for family_id=1 you have 3 individuals.  Currently, you would have each individual in a separate record.  By "roll-up", I would get the sum of those three individuals' case counts.  You would then analyze cases per family, rather than per individual.  This may introduce some overdispersion, so you might want to also consider count distributions that handle overdispersion, such as the negative binomial or generalized Poisson.   SteveDenham ... View more

Family_id will be nested within f_region, so try this:   Proc glimmix data=temp method=quad; class toxin yeargrp SES f_region family_id; model case=toxin yeargrp SES/ dist=poisson link=log offset=lnpyrs covb cl solution; random intercept / subject=f_region; random intercept / subject=family_id(f_region); run; Key here is that if you just have a few observations within family_id (like >95% are just one) then there may be convergence problems.  In this case, you may want to "roll up" the number of cases so that each family_id has a single observation.   SteveDenham   ... View more

Yes, the RANDOM statement does mean that you are accounting for the effect of f_region.  If you want to see what is going on with that statement, add SOLUTON as an option. It should display the empirical Bayes estimates for each of the levels of f_region.   SteveDenham ... View more

Not necessarily a bug, but a standard feature of the CLASS statement.  It does exactly what I would expect it to do.   SteveDenham ... View more

Following up on @Reeza 's comment - If I understand correctly you have two dummy variables - event1 and event2..  Event1 is 1 if it is 2020Q1, and 0 for all other quarters, while event2 is 1 if it is 2020Q2 and 0 for all other quarters.  Is that how the dummy variables are coded?  This is probably more critical to working out the effect of a level shift due to something changing from Q1 to Q2 in 2020.  I would code these somewhat differently - event1 is 1 for all quarters up to and including 2020Q1 and 0 for all quarters after that, and event2 is 0 for all quarters up to and including 2020Q1 and 1 for the single quarter after that.  If it is coded any way besides these two, then there will be missing values for the other quarters and the design matrix will not be full rank.  PANEL interprets this as multicollinearity, but it actually is misspecification of the model.   Now if the coding is as the first, and you treat event1 and event2 as continuous covariates (your first code), all will work out as the design matrix is still full rank.  However, once you move these two variables to class variables, the two columns in the design matrix will be a linear combination of one another, and you are now not full rank.  If you code in the second manner, I suspect the same will happen.  Here is where @PaigeMiller 's suggestion of the OUTTRANS= option will come in handy to check what is going on.  Additionally, use the CORR or COV option in the MODEL statement to see if there a columns/rows that are all zero--these are the culprits in the non-full rank/multicollinearity error.   So what to do with non-full rank designs?  Use methods that do not require full rank matrices and actually use reduced rank matrices.  I offer again PROC GLIMMIX. Or reconsider how to code the event.  It is possible to use a single variable to code what you have, rather than two.  That would solve the non-full rank issue. I suggest coding all of the quarters up to and including 2020Q1 as 0, and 2020Q2 as 1.  Call this variable level_event, and use it as a CLASS variable.   SteveDenham ... View more