# WPPS MSc data

# robert.grant@sgul.kingston.ac.uk

library(foreign)
library(dplyr)
library(reshape2)
library(stringr)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# Mac
setwd('~/Dropbox/Papers pending/Widening participation in physio')
source("~/git/html-reports/html-reports.R")
# Windows
#setwd('C:/Users/ku52502/Dropbox/Papers pending/Widening participation in physio')
#source("C:/Users/ku52502/git/html-reports/html-reports.R")

# get student data
students<-read.csv('wpps-msc.csv')[1:268,]

# get new award codes
new_awards<-read.csv('WPPS-MSc-v3-awards.csv')
#replace award codes
students$Award.class.code<-new_awards$Award.class.code
rm(new_awards)

# clean coding
students$Year.of.Entry<-as.factor(as.numeric(students$Year.of.Entry))
levels(students$Year.of.Entry)<-c("07","08","09","10","10","11","12","13")
students$Age.code<-as.factor(students$Age.code)
levels(students$Age.code)<-c('21-29','30+')
students$Gender.code<-as.factor(students$Gender.code)
levels(students$Gender.code)<-c('male','female')
students$Ethnicity.code<-as.factor(students$Ethnicity.code)
levels(students$Ethnicity.code)<-c('White','White - British','White - Irish',
                                  'White - Other white background',
                                  'Asian or Asian British - Indian',
                                  'Asian or Asian British - Pakistani',
                                  'Asian - Chinese',
                                  'Asian - Other',
                                  'Black - African',
                                  'Black - Caribbean',
                                  'Black - Other black background',
                                  'Mixed - White and Asian',
                                  'White - Black/African',
                                  'White - Black/Caribbean',
                                  'Other mixed background',
                                  'Other ethnic background',
                                  'Information refused')
students$old.Ethnicity.group.code<-students$Ethnicity.group.code # coding problems in this
students$ethnic<-NA
students$ethnic[as.numeric(students$Ethnicity.code)<5]<-1
students$ethnic[as.numeric(students$Ethnicity.code)>4 &
                                as.numeric(students$Ethnicity.code)<9]<-2
students$ethnic[as.numeric(students$Ethnicity.code)>8 & 
                                as.numeric(students$Ethnicity.code)<12]<-3
students$ethnic[as.numeric(students$Ethnicity.code)>11 & 
                                as.numeric(students$Ethnicity.code)<17]<-4
students$ethnic[as.numeric(students$Ethnicity.code)>16]<-6
students$Ethnicity.group.code<-as.factor(students$ethnic)
levels(students$Ethnicity.group.code)<-c('White as per HESA sub group',
                                         'Asian as per HESA sub group',
                                         'Black as per HESA sub group',
                                         'Other or mixed as per HESA sub group',
                                         'Not known or not given',
                                         'Information refused')
students$Qual.code<-as.factor(students$Qual.code)
levels(students$Qual.code)<-c('1st degree - L6',
                              'PGCert or PGDip or Masters or Doctorate - Level 7 or above',
                              'Below L6 or not known')
students$Disability.code<-as.factor(students$Disability.code)
levels(students$Disability.code)<-c('none',
                                    'physical',
                                    'learning',
                                    'mental',
                                    'other',
                                    'not known or not given')
students$DSA.code<-as.factor(students$DSA.code)
levels(students$DSA.code)<-c('Receiving DSA',
                             'Not receiving DSA (but has disability)',
                             'Not known (but has disability)',
                             'No disability')
students$TT.Accom.code<-as.factor(students$TT.Accom.code)
levels(students$TT.Accom.code)<-c('on campus',
                                  'off campus (inc shared house, own home, parents home)',
                                  'not known/other')
students$Socio.Economic.Class.Code.1<-as.factor(students$Socio.Economic.Class.Code.1)
levels(students$Socio.Economic.Class.Code.1)<-c('1 Higher managerial & professional occupations',
                                                '2 Lower managerial & professional occupations',
                                                '3 Intermediate occupations',
                                                '5 Lower supervisory & technical occupations',
                                                '6 Semi-routine occupations',
                                                '7 Routine occupations',
                                                '9 Not classified',
                                                '10 Not available/blank')
students$Parents.in.HE.code<-as.factor(students$Parents.in.HE.code)
levels(students$Parents.in.HE.code)<-c('yes',
                                  'no',
                                  "not known/don't know",
                                  'Not given/blank/N/A')
students$Interrupt.Abeyance.code<-as.factor(students$Interrupt.Abeyance.code)
levels(students$Interrupt.Abeyance.code)<-c('yes',
                                            'no')
students$WD.code<-as.factor(students$WD.code)
levels(students$WD.code)<-c('Withdrawn by University (academic failure)',
                            'Student has withdrawn',
                            'Student successfully completed course',
                            'Student still enrolled',
                            'Not known')
students$Award.class.code[students$Award.class.code==9]<-4
students$Award.class.code[students$Award.class.code==6]<-5
students$Award.class.code[students$Award.class.code==8]<-6
students$Award.class.code<-as.factor(students$Award.class.code)

levels(students$Award.class.code)<-c('MSc Distinction',
                                     'MSc Merit',
                                     'MSc Pass',
                                     'MSc (not classified)',
                                     'Other PG qualification',
                                     'Still enrolled',
                                     'Withdrawn / failed')
                            

# preliminary cleaning checks only
#source('preliminary-cleaning-checks.R')

# generate MSc analysis output in HTML format
tableno<-figureno<-1
con<-file('wpps-msc.html',open='w')
html_start(projecttitle='Widening participation in physiotherapy students: MSc results',handle=con)

html_h('Demographics')
tableno<-html_unitab(students$Age.code,handle=con,caption="Age",tableno=tableno)
tableno<-html_unitab(students$Gender.code,handle=con,caption="Sex",tableno=tableno)
tableno<-html_unitab(students$Ethnicity.group.code,handle=con,caption="Ethnicity",tableno=tableno)
tableno<-html_unitab(students$Disability.code,handle=con,caption="Disability",tableno=tableno)
tableno<-html_unitab(students$DSA.code,handle=con,caption="Receipt of DSA",tableno=tableno)
tableno<-html_unitab(students$Socio.Economic.Class.Code.1,handle=con,
                     caption="Socio-economic class",tableno=tableno)
tableno<-html_unitab(students$Parents.in.HE.code,handle=con,
                     caption="Parents took part in Higher Education",tableno=tableno)
tableno<-html_unitab(students$POLAR.3,handle=con,
                     caption="POLAR social participation score",tableno=tableno)


html_h('Details of study')
tableno<-html_unitab(students$Year.of.Entry,handle=con,caption="Year of entry",tableno=tableno)
tableno<-html_unitab(students$Qual.code,handle=con,caption="Qualification at entry",tableno=tableno)
tableno<-html_unitab(students$TT.Accom.code,handle=con,caption="Accommodation",tableno=tableno)
tableno<-html_unitab(students$Interrupt.Abeyance.code,handle=con,
                     caption="Interruption to studies",tableno=tableno)

html_h('Outcomes')
tableno<-html_unitab(students$WD.code,handle=con,
                     caption="Completion or withdrawal from study",tableno=tableno)
tableno<-html_unitab(students$Award.class.code,handle=con,
                     caption="Award outcome",tableno=tableno)

html_p(paste("Note that institutions classify degrees in different ways, in particular",
             "one participating HEI issues MSc degrees but does not classify them at all."))



##################################################################################################
##############   manipulate data into one row per module & one row per assessment   ##############
##################################################################################################

# make long data frame of results as first stepping stone...
temp1<-melt(students[c(1,(46:249))],id.vars='Student.code')
temp1$variable<-as.character(temp1$variable)
temp1$value[temp1$value=='N/A']<-NA
# change 'Module.' names into 'M.'
for(i in 1:17) {
  temp1$variable[temp1$variable==paste0('Module.',i,'.code')]<-paste0('M',i,'.code')
}
# move module order number into its own variable
single_digit<-(str_sub(temp1$variable,3,3)==".")
temp1$module_order<-str_sub(temp1$variable,
                                  start=2,
                                  end=3-as.numeric(single_digit))
temp1$variable<-str_sub(temp1$variable,
                              start=5-as.numeric(single_digit))
rm(single_digit)
# move module code into its own variable
temp1$module_code<-''
temp1<-arrange(temp1,Student.code,module_order)
temp1$module_code[temp1$variable=='code']<-temp1$value[temp1$variable=='code']

###############   obsolete but included as a C++ idea   ##################
# LOCF for module codes: this is unpleasantly slow but zoo::na.locf() doesn't do character vectors
# could usefully convert to Rcpp
#pb<-txtProgressBar(min=1,max=dim(temp1)[1],style=3)
#for(i in 1:dim(temp1)[1]) {
#  setTxtProgressBar(pb,i)
#  if(temp1$module_code[i]!='') {
#    thecode<-temp1$module_code[i]
#  }
#  else {
#    temp1$module_code[i]<-thecode
#  }
#}
#close(pb)
#rm(thecode)
#########################################################################
mcodeindex<-which(temp1$module_code!='')
mcodetag<-as.numeric(temp1$module_code!='')
mcodenum<-cumsum(mcodetag)
mcoderow<-mcodeindex[mcodenum]
mcode<-temp1$module_code[mcoderow]
temp1$module_code<-mcode
rm(mcodeindex)
rm(mcodetag)
rm(mcodenum)
rm(mcoderow)
rm(mcode)

# move back to wide format as a stepping stone...
temp2<-dcast(temp1,Student.code+module_order ~ variable)
  # drop unnecessary variables (they are still in 'students')
temp2<-select(temp2,-one_of(c('module_order','.credit.value','credit.value','mark','P.F.plment.result')))
  # move module code to the 2nd column
temp2<-temp2[,c('Student.code','code',names(temp2)[2:28])]
names(temp2)[2]<-'Module.code'
# now go back to long assessment-level data as the third stepping stone...
temp3<-melt(temp2,id.vars=c('Student.code','Module.code'))
  # drop empty modules
temp3<-filter(temp3,Module.code!='')
  # sort
temp3<-arrange(temp3,Student.code,Module.code)
  # move assessment number into its own variable
single_digit<-(str_sub(temp3$variable,3,3)==".")
temp3$a_number<-str_sub(temp3$variable,
                        start=2,
                        end=3-as.numeric(single_digit))
temp3$variable<-str_sub(temp3$variable,
                        start=5-as.numeric(single_digit))
rm(single_digit)
# now reshape to wide assessment-level data
assessments<-dcast(temp3,Student.code+Module.code+a_number ~ variable)
names(assessments)[4:8]<-c('attempt1','attempt2','blinded','timed','university')
  # drop empty assessments
assessments<-filter(assessments,!(blinded==""&university==""))

# merge in student details
studentvars<-c("Student.code",
               "Year.of.Entry",
               "Age.code",
               "Gender.code",
               "Ethnicity.group.code",
               "Qual.code",
               "FT.PT.code",
               "Disability.code",
               "DSA.code",
               "TT.Accom.code",
               "Socio.Economic.Class.Code.1",
               "Parents.in.HE.code",
               "POLAR.3",
               "Interrupt.Abeyance.code",
               "WD.code",
               "Award.class.code")
assessments<-left_join(assessments,select(students,one_of(studentvars)))

# assessments dplyr-grouped by module could be useful
g_modules<-group_by(assessments,Module.code)




##################################################################################################
###################################   summarise assessments   ####################################
##################################################################################################

n_assess<-dim(assessments)[1]
n_pf<-sum(is.na(as.numeric(assessments$attempt1)))
html_p(paste0('We collected data on ',n_assess,' assessments, ',n_pf,' of which were marked only ',
              'as pass or fail (mostly clinical practice assessments. For the remainder, ',
              'the distribution of the percentage marks at first attempt are shown below.'))

  # distribution of % marks
png('msc_hist_all_attempt1.png',width=560)
hist(as.numeric(assessments$attempt1),
     breaks=seq(from=0,to=100,by=5),
     col='#f2f2f2',
     xlab='Marks at first attempt',
     ylab='Number of assessments',
     main='')
dev.off()
figureno<-html_img(img_file='msc_hist_all_attempt1.png',
                   caption='Percentage marks at first attempt over all assessments',
                   figureno=figureno)

  # variation across modules
mymedian<-function(x) {
  y<-as.numeric(x)
  return(median(y,na.rm=TRUE))
}
myn<-function(x) {
  y<-as.numeric(x)
  return(sum(!is.na(y)))
}
g_modules_marks<-summarise(g_modules,medmarks=mymedian(attempt1),nmarks=myn(attempt1))
g_modules_marks$module_difficulty<-NA
g_modules_marks$module_difficulty[g_modules_marks$nmarks>9]<-
      rank(g_modules_marks$medmarks[g_modules_marks$nmarks>9],ties.method='random')
g_modules_marks<-arrange(g_modules_marks,module_difficulty)
assessments<-left_join(assessments,g_modules_marks)
png('msc_box_all_attempt1_module.png',width=560)
boxplot(as.numeric(assessments$attempt1)~assessments$module_difficulty,
        xlab='Modules ordered by median mark',
        ylab='Mark at first attempt')
dev.off()
figureno<-html_img(img_file='msc_box_all_attempt1_module.png',
                   caption='Boxplot of marks per module',
                   figureno=figureno)

  # describe difference by blinded or not, etc
qmark_b_no<-quantile(as.numeric(assessments$attempt1[assessments$blinded=='no']),
                     na.rm=TRUE,
                     probs=c(0.25,0.5,0.75))
qmark_b_yes<-quantile(as.numeric(assessments$attempt1[assessments$blinded=='yes']),
                      na.rm=TRUE,
                      probs=c(0.25,0.5,0.75))
qmark_t_no<-quantile(as.numeric(assessments$attempt1[assessments$timed=='no']),
                     na.rm=TRUE,
                     probs=c(0.25,0.5,0.75))
qmark_t_yes<-quantile(as.numeric(assessments$attempt1[assessments$timed=='yes']),
                      na.rm=TRUE,
                      probs=c(0.25,0.5,0.75))
qmark_u_no<-quantile(as.numeric(assessments$attempt1[assessments$university=='no']),
                     na.rm=TRUE,
                     probs=c(0.25,0.5,0.75))
qmark_u_yes<-quantile(as.numeric(assessments$attempt1[assessments$university=='yes']),
                      na.rm=TRUE,
                      probs=c(0.25,0.5,0.75))
html_p(paste0('The median (and inter-quartile range) for first attempt marks among blinded assessments was ',
              qmark_b_yes[2],' (',qmark_b_yes[1],' to ',qmark_b_yes[3],') ',
              'and among unblinded assessments was ',
              qmark_b_no[2],' (',qmark_b_no[1],' to ',qmark_b_no[3],'). ',
              'The same statistics for timed and untimed blinded assessments were ',
              qmark_t_yes[2],' (',qmark_t_yes[1],' to ',qmark_t_yes[3],') ',
              'and ',
              qmark_t_no[2],' (',qmark_t_no[1],' to ',qmark_t_no[3],') respectively. ',
              'The same statistics for university and clinical assessments were ',
              qmark_u_yes[2],' (',qmark_u_yes[1],' to ',qmark_u_yes[3],') ',
              'and ',
              qmark_u_no[2],' (',qmark_u_no[1],' to ',qmark_u_no[3],') respectively.'))

# describe difference by WP factors:
html_h('Crude associations with widening participation factors')

png('msc_box_all_attempt1_ethnic.png')
boxplot(as.numeric(assessments$attempt1)~assessments$Ethnicity.group.code,
        xaxt='n',ylab='% marks at first attempt',xlim=c(0.3,4.7))
axis(1,las=2,labels=c('White',
                      'Asian',
                      'Black',
                      'Other'),
     at=c(1,2,3,4))
dev.off()
figureno<-html_img(img_file='msc_box_all_attempt1_ethnic.png',
                   caption='Boxplot of marks by ethnic group',
                   figureno=figureno)

png('msc_box_all_attempt1_polar.png')
boxplot(as.numeric(assessments$attempt1)~assessments$POLAR.3,
        xaxt='n',ylab='% marks at first attempt')
axis(1,las=2,labels=c('Missing',
                      '1',
                      '2',
                      '3',
                      '4',
                      '5',
                      'Non UK',
                      'N/A'),
     at=1:8)
dev.off()
figureno<-html_img(img_file='msc_box_all_attempt1_polar.png',
                   caption='Boxplot of marks by POLAR',
                   figureno=figureno)

png('msc_box_all_attempt1_disability.png')
boxplot(as.numeric(assessments$attempt1)~assessments$Disability.code,
        xaxt='n',ylab='% marks at first attempt')
axis(1,las=2,labels=c('None',
                      'Physical',
                      'Learning',
                      'Mental',
                      'Other',
                      'Missing'),
     at=1:6)
dev.off()
figureno<-html_img(img_file='msc_box_all_attempt1_disability.png',
                   caption='Boxplot of marks by disabiliity',
                   figureno=figureno)

png('msc_box_all_attempt1_parents.png')
boxplot(as.numeric(assessments$attempt1)~assessments$Parents.in.HE.code,
        xaxt='n',ylab='% marks at first attempt',xlim=c(0.3,4.7))
axis(1,las=2,labels=c('Yes',
                      'No',
                      'Undisclosed',
                      'Missing'),
     at=1:4)
dev.off()
figureno<-html_img(img_file='msc_box_all_attempt1_parents.png',
                   caption='Boxplot of marks by whether a parent attended HE',
                   figureno=figureno)

png('msc_box_all_attempt1_qual.png')
boxplot(as.numeric(assessments$attempt1)~assessments$Qual.code,
        xaxt='n',ylab='% marks at first attempt')
axis(1,las=2,labels=c('L6',
                      'L7+',
                      'Other / NK'),
     at=1:3)
dev.off()
figureno<-html_img(img_file='msc_box_all_attempt1_qual.png',
                   caption='Boxplot of marks by entry qualification',
                   figureno=figureno)


##################################################################################################
####################################   fit Bayesian model   ######################################
##################################################################################################

html_h("Bayesian regression models")
html_p(paste0('In these models, we examine the effect of widening participation factors ',
              'on attainment - assessment marks as well as award grade. ',
              'The other demographic factors and characteristics of the assessment are of ',
              'interest only insofar as they interact with WP factors, but we must consider adjusting ',
              'for them anyway. Coarsened marks (fails recorded as zero, as well as practice ',
              'assessments recorded only as pass or fail) were imputed from the distribution of ',
              'all marks, conditional on the covariates, and truncated to the pass/fail region. ',
              'Random effects were included for students and modules, reflecting the fact that ',
              'each assessment belongs to a student and a module and that they are ',
              'not independent events as classical statistical tests assume. These complexities ',
              'are one reason for using a Bayesian model, the other being that we can readily ',
              'obtain intuitive statements in terms of probabilities, such as the probability ',
              'that covariate X has an effect bigger than Y - which cannot be obtained by other ',
              'methods.'))
html_p(paste0('The table below shows the maximal model formed with all predictors. The ',
              'Stan code for fitting the model can be read by ',
              '<a href="wpps-test10.stan">clicking here</a>.'))
html_p('<img src="mathtex1.gif"><br>
      <img src="mathtex2.gif"><br>
      <img src="mathtex3.gif"><br>
      <img src="mathtex4.gif">')
# note that the Student.code numbers work because those missing assessment/module data
# are at the end of the data frame - so the numbers are consecutive from 1.

#fix a couple of assessment types
assessments$blinded[assessments$Module.code=='M1' & assessments$a_number==1]<-'no'
assessments$blinded[assessments$Module.code=='M3' & assessments$a_number==2]<-'no'

# make numeric covariates
assessments$n_polar<-as.numeric(assessments$POLAR.3) # needs not applicable and maybe non-UK to be treated as missing
assessments$polar1<-as.numeric(assessments$n_polar==2) # POLAR.3 1st quintile has numeric code 2...
assessments$polar2<-as.numeric(assessments$n_polar==3) 
  # use 3rd quintile as baseline
assessments$polar4<-as.numeric(assessments$n_polar==5)
assessments$polar5<-as.numeric(assessments$n_polar==6)
assessments$polar1[!(assessments$n_polar>1 & assessments$n_polar<7)]<-NA
assessments$polar2[!(assessments$n_polar>1 & assessments$n_polar<7)]<-NA
assessments$polar4[!(assessments$n_polar>1 & assessments$n_polar<7)]<-NA
assessments$polar5[!(assessments$n_polar>1 & assessments$n_polar<7)]<-NA
assessments$female<-as.numeric(assessments$Gender.code)-1
assessments$over30<-as.numeric(assessments$Age.code=='30+')
assessments$ethn_asian<-as.numeric(as.numeric(assessments$Ethnicity.group.code)==2)
assessments$ethn_black<-as.numeric(as.numeric(assessments$Ethnicity.group.code)==3)
assessments$ethn_other<-as.numeric(as.numeric(assessments$Ethnicity.group.code)==4)
assessments$disab_phys<-as.numeric(as.numeric(assessments$Disability.code)==2)
assessments$disab_learn<-as.numeric(as.numeric(assessments$Disability.code)==3)
assessments$disab_other<-as.numeric(as.numeric(assessments$Disability.code)==5)
assessments$qual_level7<-as.numeric(as.numeric(assessments$Qual.code)==2)
assessments$qual_other<-as.numeric(as.numeric(assessments$Qual.code)==3)
assessments$blinded_untimed<-as.numeric(assessments$blinded=='yes' & (!assessments$timed=='yes'))
assessments$blinded_timed<-as.numeric(assessments$blinded=='yes' & assessments$timed=='yes')
assessments$clinical<-as.numeric(assessments$university=='no')
assessments$mark<-as.numeric(assessments$attempt1)
assessments$mod<-as.numeric(as.factor(assessments$Module.code))



# make data frame that contains known % marks
a1<-filter(assessments,!is.na(as.numeric(attempt1))&(as.numeric(attempt1)!=0))
# and one that contains zero marks and 'F'
a0<-filter(assessments,(attempt1=='F')|(as.numeric(attempt1)==0))
# and one that contains 'P' results
a40<-filter(assessments,attempt1=='P')

# reduce these data frames to the non-missing observations for the model
a0<-filter(a0,!is.na(polar1))
a1<-filter(a1,!is.na(polar1))
a40<-filter(a40,!is.na(polar1))

stancon<-file('wpps-test10.stan',open='w')
write('
data {
	int Nmark;
  int N0;
  int N40;
  int Nstu;
  int Nmod;
  real polar1[Nmark,4];
  real polar0[N0,4];
  real polar40[N40,4];
  int blinded_untimed1[Nmark];
  int blinded_untimed0[N0];
  int blinded_untimed40[N40];
  int blinded_timed1[Nmark];
  int blinded_timed0[N0];
  int blinded_timed40[N40];
  int clinical1[Nmark];
  int clinical0[N0];
  int clinical40[N40];
  int female1[Nmark];
  int female0[N0];
  int female40[N40];
  int over301[Nmark];
  int over300[N0];
  int over3040[N40];
  int ethn_asian1[Nmark];
  int ethn_asian0[N0];
  int ethn_asian40[N40];
  int ethn_black1[Nmark];
  int ethn_black0[N0];
  int ethn_black40[N40];
  int ethn_other1[Nmark];
  int ethn_other0[N0];
  int ethn_other40[N40];
  int disab_phys1[Nmark];
  int disab_phys0[N0];
  int disab_phys40[N40];
  int disab_learn1[Nmark];
  int disab_learn0[N0];
  int disab_learn40[N40];
  int disab_other1[Nmark];
  int disab_other0[N0];
  int disab_other40[N40];
  int qual_level71[Nmark];
  int qual_level70[N0];
  int qual_level740[N40];
  int qual_other1[Nmark];
  int qual_other0[N0];
  int qual_other40[N40];
  real mark[Nmark];
  int<lower=0> stu1[Nmark];
  int<lower=0> stu0[N0];
  int<lower=0> stu40[N40];
  int<lower=0> mod1[Nmark];
  int<lower=0> mod0[N0];
  int<lower=0> mod40[N40];
}
parameters {
  real alpha;
  real beta[17];
  real<lower=0> sigma;
  real<lower=0, upper=39.5> latentmark0[N0];
  real<lower=39.5, upper=100> latentmark40[N40];
  real ustu[Nstu];
  real<lower=0> sigmastu;
  real umod[Nmod];
  real<lower=0> sigmamod;
}
model {
  real mu[Nmark];
  real mu0[N0];
  real mu40[N40];
  sigma ~ normal(0,30);
  alpha ~ normal(50,30);
  beta ~ normal(0,20);
  sigmastu ~ normal(0,30);
  ustu ~ normal(0,sigmastu);
  sigmamod ~ normal(0,30);
  umod ~ normal(0,sigmamod);
  for (i in 1:Nmark) {
    mu[i] <- alpha+(beta[1]*polar1[i,1])+
            (beta[2]*polar1[i,2])+
            (beta[3]*polar1[i,3])+
            (beta[4]*polar1[i,4])+
            (beta[5]*blinded_untimed1[i])+
            (beta[6]*blinded_timed1[i])+
            (beta[7]*clinical1[i])+
            (beta[8]*female1[i])+
            (beta[9]*over301[i])+
            (beta[10]*ethn_asian1[i])+
            (beta[11]*ethn_black1[i])+
            (beta[12]*ethn_other1[i])+
            (beta[13]*disab_phys1[i])+
            (beta[14]*disab_learn1[i])+
            (beta[15]*disab_other1[i])+
            (beta[16]*qual_level71[i])+
            (beta[17]*qual_other1[i])+
            ustu[stu1[i]]+umod[mod1[i]];
    mark[i] ~ normal(mu[i],sigma);
  }
  for (i in 1:N0) {
    mu0[i] <- alpha+(beta[1]*polar0[i,1])+
              (beta[2]*polar0[i,2])+
              (beta[3]*polar0[i,3])+
              (beta[4]*polar0[i,4])+
              (beta[5]*blinded_untimed0[i])+
              (beta[6]*blinded_timed0[i])+
              (beta[7]*clinical0[i])+
              (beta[8]*female0[i])+
              (beta[9]*over300[i])+
              (beta[10]*ethn_asian0[i])+
              (beta[11]*ethn_black0[i])+
              (beta[12]*ethn_other0[i])+
              (beta[13]*disab_phys0[i])+
              (beta[14]*disab_learn0[i])+
              (beta[15]*disab_other0[i])+
              (beta[16]*qual_level70[i])+
              (beta[17]*qual_other0[i])+
              ustu[stu0[i]]+umod[mod0[i]];
    latentmark0[i] ~ normal(mu0[i],sigma);
  }
  for (i in 1:N40) {
    mu40[i] <- alpha+(beta[1]*polar40[i,1])+
                (beta[2]*polar40[i,2])+
                (beta[3]*polar40[i,3])+
                (beta[4]*polar40[i,4])+
                (beta[5]*blinded_untimed40[i])+
                (beta[6]*blinded_timed40[i])+
                (beta[7]*clinical40[i])+
                (beta[8]*female40[i])+
                (beta[9]*over3040[i])+
                (beta[10]*ethn_asian40[i])+
                (beta[11]*ethn_black40[i])+
                (beta[12]*ethn_other40[i])+
                (beta[13]*disab_phys40[i])+
                (beta[14]*disab_learn40[i])+
                (beta[15]*disab_other40[i])+
                (beta[16]*qual_level740[i])+
                (beta[17]*qual_other40[i])+
                ustu[stu40[i]]+umod[mod40[i]];
    latentmark40[i] ~ normal(mu40[i],sigma);
  }
}
',file=stancon)
close(stancon)
wppsdata<-list(Nmark=dim(a1)[1],
               N0=dim(a0)[1],
               N40=dim(a40)[1],
               Nstu=length(table(assessments$Student.code)),
               Nmod=length(table(assessments$Module.code)),
               polar1=as.matrix(a1[,c('polar1','polar2','polar4','polar5')]),
               polar0=as.matrix(a0[,c('polar1','polar2','polar4','polar5')]),
               polar40=as.matrix(a40[,c('polar1','polar2','polar4','polar5')]),
               blinded_untimed0=a0$blinded_untimed,
               blinded_untimed1=a1$blinded_untimed,
               blinded_untimed40=a40$blinded_untimed,
               blinded_timed0=a0$blinded_timed,
               blinded_timed1=a1$blinded_timed,
               blinded_timed40=a40$blinded_timed,
               clinical0=a0$clinical,
               clinical1=a1$clinical,
               clinical40=a40$clinical,
               female0=a0$female,
               female1=a1$female,
               female40=a40$female,
               ethn_asian0=a0$ethn_asian,
               ethn_asian1=a1$ethn_asian,
               ethn_asian40=a40$ethn_asian,
               ethn_black0=a0$ethn_black,
               ethn_black1=a1$ethn_black,
               ethn_black40=a40$ethn_black,
               ethn_other0=a0$ethn_other,
               ethn_other1=a1$ethn_other,
               ethn_other40=a40$ethn_other,
               disab_phys0=a0$disab_phys,
               disab_phys1=a1$disab_phys,
               disab_phys40=a40$disab_phys,
               disab_learn0=a0$disab_learn,
               disab_learn1=a1$disab_learn,
               disab_learn40=a40$disab_learn,
               disab_other0=a0$disab_other,
               disab_other1=a1$disab_other,
               disab_other40=a40$disab_other,
               qual_level70=a0$qual_level7,
               qual_level71=a1$qual_level7,
               qual_level740=a40$qual_level7,
               qual_other0=a0$qual_other,
               qual_other1=a1$qual_other,
               qual_other40=a40$qual_other,
               over300=a0$over30,
               over301=a1$over30,
               over3040=a40$over30,
               stu0=a0$Student.code,
               stu1=a1$Student.code,
               stu40=a40$Student.code,
               mod0=a0$mod,
               mod1=a1$mod,
               mod40=a40$mod,
               mark=a1$mark)
wppsfit<-stan(file='wpps-test10.stan',
              data=wppsdata,
              iter=5000,
              chains = 4,
              seed=7144)
stanoutput<-as.matrix(summary(wppsfit)[[1]])[c('alpha',
                                               paste0('beta[',1:17,']'),
                                               'sigma',
                                               'sigmastu',
                                               'sigmamod'),c('mean',
                                                             '2.5%',
                                                             '97.5%',
                                                             'n_eff',
                                                             'Rhat')]
stanoutput<-round(stanoutput,digits=3)
stanoutput<-cbind(stanoutput,paste(stanoutput[,2],"to",stanoutput[,3]))
stanoutput<-stanoutput[,c(1,6,4,5)]
rownames(stanoutput)<-c('Grand mean',
                        'POLAR quintile 1',
                        'POLAR quintile 2',
                        'POLAR quintile 4',
                        'POLAR quintile 5',
                        'Blinded, untimed',
                        'Blinded, timed',
                        'Clinical',
                        'Female',
                        'Over 30',
                        'Asian ethnicity',
                        'Black ethnicity',
                        'Other ethnicity',
                        'Physical disability',
                        'Learning disability',
                        'Other disability',
                        'Entry qualification L7',
                        'Entry qualification <L6, other',
                        'Residual SD',
                        'Student ability SD',
                        'Module difficulty SD')
colnames(stanoutput)<-c('Coefficient',
                        '95% CI',
                        'Effective sample size',
                        'Convergence')
tableno<-html_anytab(stanoutput,tableno=tableno,caption='Model with all predictors but no interactions')
html_p(paste0('We formed this model by adding covariates one at a time over ten steps. ',
              'For the record, we used Stan software through the rstan package in R ',
              '(mc-stan.org ; r-project.org). ',
              'Weakly informative priors were used throughout and four chains run for 5000 ',
              'iterations each, discarding the first half as warm-up (reference: ',
              'Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. <em>Bayesian ',
              'Data Analysis</em> 3rd edition, 2014, CRC Press, Boca Raton, Florida USA). ',
              'The "Effective sample size" column shows how the 10,000 iterations compares to ',
              'independent random numbers from the posterior distributions of the parameters we ',
              'want to learn about. A higher number is better. The r-hat statistic is shown ',
              ' in the "Convergence" column. This should be in the range 0.9 to 1.1 to indicate ',
              'that the simulation has converged to the desired posterior distribution. ',
              'Whether parent attended HE was omitted because the great majority of students ',
              'did not have this recorded. The ',
              'coefficients reflect the change in the mean mark that results from that covariate ',
              'being present. Each is measured relative to its own baseline level, for example ',
              'POLAR 3rd quintile is the baseline because it is not listed in the regression ',
              'coefficients. As the model was constructed, none of the coefficients changed much, ',
              'suggesting that there is relatively little confounding. Next, we consider whether ',
              'they interact with (also known as "effect modification" in the epidemiological',
              ' literature) the widening participation factors.' ))
html_p(paste0('The individual widening participation factors for further consideration can be ',
              'selected on the basis of the model above. Ethnicity and disability are ',
              'interesting, while POLAR and entry qualifications are not, based on the ',
              'coefficients and their confidence intervals.'))
html_p(paste0('Among the demographic and assessment covariates, female students had slightly ',
              'higher marks, with borderline statistical significance. Age bands split at 30 ',
              'showed no significant difference, while blinded asessments, both timed and ',
              'untimed, got significantly lower marks than unblinded ones'))


# interaction model


stancon<-file('wpps-inter7.stan',open='w')
write('
    data {
      int Nmark;
      int N0;
      int N40;
      int Nstu;
      int Nmod;
      int blinded_untimed1[Nmark];
      int blinded_untimed0[N0];
      int blinded_untimed40[N40];
      int blinded_timed1[Nmark];
      int blinded_timed0[N0];
      int blinded_timed40[N40];
      int clinical1[Nmark];
      int clinical0[N0];
      int clinical40[N40];
      int female1[Nmark];
      int female0[N0];
      int female40[N40];
      int over301[Nmark];
      int over300[N0];
      int over3040[N40];
      int ethn_asian1[Nmark];
      int ethn_asian0[N0];
      int ethn_asian40[N40];
      int ethn_black1[Nmark];
      int ethn_black0[N0];
      int ethn_black40[N40];
      int ethn_other1[Nmark];
      int ethn_other0[N0];
      int ethn_other40[N40];
      int disab_phys1[Nmark];
      int disab_phys0[N0];
      int disab_phys40[N40];
      int disab_learn1[Nmark];
      int disab_learn0[N0];
      int disab_learn40[N40];
      int disab_other1[Nmark];
      int disab_other0[N0];
      int disab_other40[N40];
      real mark[Nmark];
      int<lower=0> stu1[Nmark];
      int<lower=0> stu0[N0];
      int<lower=0> stu40[N40];
      int<lower=0> mod1[Nmark];
      int<lower=0> mod0[N0];
      int<lower=0> mod40[N40];
    }
    parameters {
      real alpha;
      real beta[11];
      matrix[2,1] gamma;
      real<lower=0> sigma;
      real<lower=0, upper=39.5> latentmark0[N0];
      real<lower=39.5, upper=100> latentmark40[N40];
      real ustu[Nstu];
      real<lower=0> sigmastu;
      real umod[Nmod];
      real<lower=0> sigmamod;
    }
    model {
      real mu[Nmark];
      real mu0[N0];
      real mu40[N40];
      sigma ~ normal(0,30);
      alpha ~ normal(50,30);
      beta ~ normal(0,20);
      sigmastu ~ normal(0,30);
      ustu ~ normal(0,sigmastu);
      sigmamod ~ normal(0,30);
      umod ~ normal(0,sigmamod);
      for (i in 1:Nmark) {
        mu[i] <- alpha+
                (beta[1]*blinded_untimed1[i])+
                (beta[2]*blinded_timed1[i])+
                (beta[3]*clinical1[i])+
                (beta[4]*female1[i])+
                (beta[5]*over301[i])+
                (beta[6]*ethn_asian1[i])+
                (beta[7]*ethn_black1[i])+
                (beta[8]*ethn_other1[i])+
                (beta[9]*disab_phys1[i])+
                (beta[10]*disab_learn1[i])+
                (beta[11]*disab_other1[i])+
                (gamma[1,1]*ethn_other1[i]*blinded_untimed1[i])+
                (gamma[2,1]*ethn_other1[i]*blinded_timed1[i])+
                ustu[stu1[i]]+umod[mod1[i]];
        mark[i] ~ normal(mu[i],sigma);
      }
      for (i in 1:N0) {
        mu0[i] <- alpha+
                (beta[1]*blinded_untimed0[i])+
                (beta[2]*blinded_timed0[i])+
                (beta[3]*clinical0[i])+
                (beta[4]*female0[i])+
                (beta[5]*over300[i])+
                (beta[6]*ethn_asian0[i])+
                (beta[7]*ethn_black0[i])+
                (beta[8]*ethn_other0[i])+
                (beta[9]*disab_phys0[i])+
                (beta[10]*disab_learn0[i])+
                (beta[11]*disab_other0[i])+
                (gamma[1,1]*ethn_other0[i]*blinded_untimed0[i])+
                (gamma[2,1]*ethn_other0[i]*blinded_timed0[i])+
                ustu[stu0[i]]+umod[mod0[i]];
        latentmark0[i] ~ normal(mu0[i],sigma);
      }
      for (i in 1:N40) {
        mu40[i] <- alpha+
                    (beta[1]*blinded_untimed40[i])+
                    (beta[2]*blinded_timed40[i])+
                    (beta[3]*clinical40[i])+
                    (beta[4]*female40[i])+
                    (beta[5]*over3040[i])+
                    (beta[6]*ethn_asian40[i])+
                    (beta[7]*ethn_black40[i])+
                    (beta[8]*ethn_other40[i])+
                    (beta[9]*disab_phys40[i])+
                    (beta[10]*disab_learn40[i])+
                    (beta[11]*disab_other40[i])+
                    (gamma[1,1]*ethn_other40[i]*blinded_untimed40[i])+
                    (gamma[2,1]*ethn_other40[i]*blinded_timed40[i])+
                    ustu[stu40[i]]+umod[mod40[i]];
        latentmark40[i] ~ normal(mu40[i],sigma);
      }
    }
',file=stancon)
close(stancon)
wppsfit<-stan(file='wpps-inter7.stan',
              data=wppsdata,
              iter=5000,
              chains = 4,
              seed=7316)
stanoutput<-as.matrix(summary(wppsfit)[[1]])[c('alpha',
                                               paste0('beta[',1:11,']'),
                                               paste0('gamma[',1:2,',1]'),
                                               'sigma',
                                               'sigmastu',
                                               'sigmamod'),c('mean',
                                                             '2.5%',
                                                             '97.5%',
                                                             'n_eff',
                                                             'Rhat')]
stanoutput<-round(stanoutput,digits=3)
stanoutput<-cbind(stanoutput,paste(stanoutput[,2],"to",stanoutput[,3]))
stanoutput<-stanoutput[,c(1,6,4,5)]
rownames(stanoutput)<-c('Grand mean',
                        'Blinded, untimed',
                        'Blinded, timed',
                        'Clinical',
                        'Female',
                        'Over 30',
                        'Asian ethnicity',
                        'Black ethnicity',
                        'Other ethnicity',
                        'Physical disability',
                        'Learning disability',
                        'Other disability',
                        'Other ethnicity + blinded, untimed',
                        'Other ethnicity + blinded, timed',
                        'Residual SD',
                        'Student ability SD',
                        'Module difficulty SD')
colnames(stanoutput)<-c('Coefficient',
                        '95% CI',
                        'Effective sample size',
                        'Convergence')
# note that interaction results are hard-typed below...
html_p(paste0('The table below shows the model with significant interactions among WP ',
              'factors and other covariates. The ',
              'Stan code for fitting the model can be read by ',
              '<a href="wpps-inter7.stan">clicking here</a>. ',
              'To arrive at this final model, we removed POLAR and qualifications from ',
              'the model above and then explored the following interaction terms: ',
              '<ul><li><strong>ethnicity with assessment type</strong>: the only terms that were ',
              'significant or borderline were ',
              'for "other" ethnicity combined with blinded assessments - either timed or untimed - ',
              'which had lower marks than would have been predicted from the previous model: ',
              'for timed assessments the estimated difference was -10.508, 95% CI -21.949 to 0.882, and ',
              'for untimed assessments the estimated difference was -4.66, 95% CI -10.698 to 1.349.</li> ',
              '<li><strong>ethnicity with sex</strong>: none were significant or borderline</li> ',
              '<li><strong>ethnicity with age</strong> could not be estimated because of no over ',
              '30 students in the Asian or Other ethnic groups</li>',
              '<li><strong>disability with assessment type</strong>: the only significant term was ',
              'for other disability with timed, blinded assessment, but this had a very wide ',
              'confidence interval, being based on one student.</li>',
              '<li><strong>disability with sex</strong>: none were significant or borderline</li>',
              '<li><strong>disability with age</strong> could not be estimated because of no over ',
              '30 students in the physical or other disability groups</li></ul>.',
              'On the basis of these analyses, we retained:',
              '<ul><li>other ethnicity + blinded, untimed</li>',
              '<li>other ethnicity + blinded, timed</li></ul>'))

tableno<-html_anytab(stanoutput,tableno=tableno,caption='Model with significant WP interactions')



close(con)