/* multilogit.sas */
title2 'Multinomial Logit Model';
%include '/home/brunner0/441s20/mathread2.sas';
/* mathread2.sas creates missused (any of hsgpa hscalc precalc missing) 
   and outcome (Pass-Fail-Gone)*/

proc freq data=mathex;
     title3 'Outcome by passed just to check';
     tables outcome*passed 
            / norow nocol nopercent missing;
run;

proc freq data=mathex;
     title3 'Missused by Outcome with proc freq';
     tables missused * outcome / nocol nopercent chisq;
run;

/* Probabilities
     pi1 = P(Gone)
     pi2 = P(Pass)
     pi3 = P(Fail)

   Multinomial Logit model for missused by outcome is

     ln(pi1/pi3) = beta01 + beta11 x             Gone vs. Fail
     ln(pi2/pi3) = beta02 + beta12 x             Pass vs. Fail
*/     


proc logistic data=mathex outest=ParmNames; 
     title3 'Misssused by Outcome with proc logistic';
     model outcome (ref='Fail') = missused / link = glogit;
     contrast 'Hsgpa, hscalc or precalc missing' missused 1;
run;


proc transpose data=ParmNames;
run;
proc print; 
run;

proc logistic data = mathex;
     title3 'Contrast versus test';
     model outcome (ref='Fail') = missused / link = glogit;
     contrast 'Missused method 1' missused 1;
     MissusedMethod2: test missused_Pass = missused_Gone = 0;
run;

proc iml;
     title3 'Estimate Probabilities using output from proc logistic';
      b01 = 0.3878; b11 = 0.7209;   /* Gone */
      b02 = 1.4123; b12 = -0.6461;  /* Pass */
      missused = 0;
      L1 = b01 + b11*missused;
      L2 = b02 + b12*missused;
      denom = 1 + exp(L1) + exp(L2);
      Gone = exp(L1)/denom; Pass = exp(L2)/denom; Fail = 1/denom;
      print "Not missing:" Fail Gone Pass;   
      missused = 1;
      L1 = b01 + b11*missused;
      L2 = b02 + b12*missused;
      denom = 1 + exp(L1) + exp(L2);
      Gone = exp(L1)/denom; Pass = exp(L2)/denom; Fail = 1/denom;
      print "Yes Missing:" Fail Gone Pass;   
run;

proc freq  data = mathex;
     title3 'Missused by outcome again for comparison';
     tables missused * outcome / nocol nopercent;
run;

/* Lots of exploration is not shown here. */

proc logistic data = mathex;
     title3 'hsgpa hscalc precalc mtongue';
     model outcome (ref='Fail') = 
                   hsgpa hscalc precalc mtongue / link = glogit;
run;

/* Allowing for academic background, students whose first language is English
are less likely to disappear as opposed to failing the course, and less likely to
Pass as opposed to failing. If this is replicated, it will be very
interesting. Explore further.

/************************** Replication ***********************
For interpretation, want to replicate 8 findings:
Gone vs. Fail and Pass vs. Fail for each explanatory variable. 
***************************************************************/
%include '/home/brunner0/441s20/readreplic.sas';
if (0<=mark<=49) then outcome = 'Fail';
   else if (50<=mark<=100) then outcome = 'Pass';
   else outcome = 'Gone';

proc logistic data = replic; 
     title2 'Replicate hsgpa hscalc precalc mtongue 0.05/8 = .00625';
     model outcome (ref='Fail') = hsgpa hscalc precalc mtongue / link = glogit;