Thank you Steve - my original approach works when I drop time effects, and the model also works in PROC BGLIMM. This follow-up question may be a better fit for a new thread, so please let me so if that is the case and I will mark your post as the accepted solution. After speaking with a colleague, it sounds like we are going to shift to a multinomial logit design as it better theoretically fits the question - i.e. with outcomes of 0-3 for no participation in any program/participation in related programs, rather than assessing each independently. When I do this (code below), the model fails to run because it cannot be fit in a reasonable amount of time, even when I aggregate schools up to the district level for random effects. The version of SAS I am using does not have multinomial modeling built into proc bglimm. Any suggestions? proc glimmix data=data NOCLPRINT NOITPRINT ;
class school;
model program = treat year time_since_treat x1 x2 x3
/ solution dist=multinomial link=glogit;
random intercept / subject=school;
run;
... View more