data insomnia; infile "/courses/ddf5e9e5ba27fe300/STAT770/Insomnia.txt"; input case treat time outcome; y1=0; y2=0; y3=0; y4=0; if outcome=1 then y1=1; if outcome=2 then y2=1; if outcome=3 then y3=1; if outcome=4 then y4=1; run; proc nlmixed qpoints=40; bounds i2 > 0; bounds i3 > 0; eta1 = i1 + treat*beta1 + time*beta2 + treat*time*beta3 + u; eta2 = i1 + i2 + treat*beta1 + time*beta2 + treat*time*beta3 + u; eta3 = i1 + i2 + i3 + treat*beta1 + time*beta2 + treat*time*beta3 + u; p1 = exp(eta1)/(1 + exp(eta1)); p2 = exp(eta2)/(1 + exp(eta2)) - exp(eta1)/(1 + exp(eta1)); p3 = exp(eta3)/(1 + exp(eta3)) - exp(eta2)/(1 + exp(eta2)); p4 = 1 - exp(eta3)/(1 + exp(eta3)); ll = y1*log(p1) + y2*log(p2) + y3*log(p3) + y4*log(p4); model y1 ~ general(ll); estimate 'interc2' i1+i2; * this is alpha_2 in model, and i1 is alpha_1; estimate 'interc3' i1+i2+i3; * this is alpha_3 in model; estimate 'd vs p at 2 weeks ' exp(beta1+beta3); estimate 'd vs p at baseline' exp(beta1); random u ~ normal(0, sigma*sigma) subject=case; run; *Treating initial outcome as covariate; data insomnia; input treat initial outcome count @@; datalines; 1 1 1 7 1 1 2 4 1 1 3 1 1 1 4 0 1 2 1 11 1 2 2 5 1 2 3 2 1 2 4 2 1 3 1 13 1 3 2 23 1 3 3 3 1 3 4 1 1 4 1 9 1 4 2 17 1 4 3 13 1 4 4 8 0 1 1 7 0 1 2 4 0 1 3 2 0 1 4 1 0 2 1 14 0 2 2 5 0 2 3 1 0 2 4 0 0 3 1 6 0 3 2 9 0 3 3 18 0 3 4 2 0 4 1 4 0 4 2 11 0 4 3 14 0 4 4 22 ; run; proc logistic; class initial / param=ref; model outcome = initial treat initial*treat; freq count; contrast 'sleep=1' treat 1 initial*treat 1 0 0 / estimate=exp; contrast 'sleep=2' treat 1 initial*treat 0 1 0 / estimate=exp; contrast 'sleep=3' treat 1 initial*treat 0 0 1 / estimate=exp; contrast 'sleep=4' treat 1 initial*treat 0 0 0 / estimate=exp; run;