/* multilogit.sas */
title2 'Multinomial Logit Model';
%include '/folders/myfolders/441s18/Lecture/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;

/* Multinomial Logit model for missused by outcome is

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


proc logistic data=mathex outest=ParmNames; 
     title3 'Misssused by Outcome with proc logistic';
     model outcome (ref='Pass') = missused / link = glogit;
     contrast 'Hsgpa, hscalc or precalc missing' missused 1;
run;
 
/* Find out the parameter names, written with outest. */

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

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

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

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

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

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

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

   We want to know whether failing is different from disappearing in terms 
   of their relationship to the explanatory variables. We are getting 
   advanced here. What is H0?

   Recall the response categories are 1=Fail 2=Gone 3=Pass.

Model (using b instead of beta) is 

     ln(pi1/pi3) = b01 + b11 hsgpa + b21 hscalc + b31 precalc + b41 mtongue
     ln(pi2/pi3) = b02 + b12 hsgpa + b22 hscalc + b32 precalc + b42 mtongue

     The null hypothesis is b11=b12, b21=b22, b31=b32, b41=b42

Parameter names are easy to guess.  */

proc logistic data = mathex;
     title3 'Different coefficients for Gone and Fail?';
     model outcome (ref='Pass') = hsgpa hscalc precalc mtongue / link = glogit;
     DiffOverall: test hsgpa_Fail = hsgpa_Gone, hscalc_Fail = hscalc_Gone,  
                       precalc_Fail = precalc_Gone, mtongue_Fail = mtongue_Gone;
     Diff_hsgpa:   test hsgpa_Fail = hsgpa_Gone;
     Diff_hscalc:  test hscalc_Fail = hscalc_Gone;
     Diff_precalc: test precalc_Fail = precalc_Gone;
     Diff_mtongue: test mtongue_Fail = mtongue_Gone;
run;

/************************** Replication ***********************
For interpretation, want to replicate 8 findings:
Gone vs. Pass and Fail vs. Pass for each explanatory variable. 
***************************************************************/

%include '/folders/myfolders/441s18/Lecture/readreplic.sas';
if (0<=mark<=49) then outcome = 'Fail';
   else if (50<=mark<=100) then outcome = 'Pass';
   else outcome = 'Gone';

proc logistic data = mathrep; 
     title2 'Replicate hsgpa hscalc precalc calc mtongue 0.05/8 = .00625';
     model outcome (ref='Pass') = hsgpa hscalc precalc mtongue / link = glogit;
     Diff_mtongue: test mtongue_Fail = mtongue_Gone;