# WPPS BSc data # robert.grant@sgul.kingston.ac.uk library(foreign) library(MASS) library(coin) library(dplyr) library(reshape2) library(stringr) library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) # prevent dplyr-MASS conflict select<-dplyr::select # 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-bsc.csv') # clean errors # add no to levels of M2 levels(students$M1.A2.B)<-c(levels(students$M1.A2.B),'no') # the order is diff to other vars but will still melt as character vectors OK # M1, assessment 2 should be unblinded university students$M1.A2.B[students$Module.1.code=='M1']<-'no' students$M1.A2.T[students$Module.1.code=='M1']<-'' students$M1.A2.U[students$Module.1.code=='M1']<-'yes' # student 106 should have withdrawal code 1 students$WD.code[students$Student.code==106]<-1 # clean coding levels(students$Year.of.Entry)<-c("05","06","07","08","09","09","10","10","11","12") yo<-2004+as.numeric(students$Year.of.Entry) png('year-of-entry-hist.png') barplot(table(yo),xlab='Year of entry',ylab='Number of students') #hist(yo,breaks=2005:2012,xlab='Year of entry',col='#cccccc',ylab='Number of students',main='') dev.off() students$Age.code<-as.factor(students$Age.code.U21.21..groups) levels(students$Age.code)<-c('<21','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 - Bangladeshi', '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', 'Not known or not given', 'Information refused') # Arab (17) not used 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)<10]<-2 students$ethnic[as.numeric(students$Ethnicity.code)>9 & as.numeric(students$Ethnicity.code)<13]<-3 students$ethnic[as.numeric(students$Ethnicity.code)>12 & as.numeric(students$Ethnicity.code)<18]<-4 students$ethnic[as.numeric(students$Ethnicity.code)>17]<-5 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 information refused') students$Qual.code<-as.factor(students$Qual.code) levels(students$Qual.code)<-c('A levels', 'Access / BTEC', 'Degree level', 'Other 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', '4 Small employers & own account workers', '5 Lower supervisory & technical occupations', '6 Semi-routine occupations', '7 Routine occupations', '8 Never worked & long-term unemployed', '9 Not classified', '10 Not available/blank') # fix missing level students$Socio.Economic.Class.Code.1[students$Socio.Economic.Class.Code.1=='9 Not classified']<-'10 Not available/blank' students$Socio.Economic.Class.Code.1[students$Socio.Economic.Class.Code.1=='8 Never worked & long-term unemployed']<-'9 Not classified' students$Parents.in.HE.code<-as.factor(students$Parents.in.HE.code) levels(students$Parents.in.HE.code)<-c('yes', 'no', 'prefer not to say/undisclosed', "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', 'Exclusion/professionally unsuitable/death', '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('1st', '2:1', '2:2', 'Still enrolled', 'Non-qualifying degree', 'No award', 'Certificate/Diploma') # preliminary cleaning checks only #source('preliminary-cleaning-checks.R') # generate BSc analysis output in HTML format tableno<-figureno<-1 con<-file('wpps-bsc.html',open='w') html_start(projecttitle='Widening participation in physiotherapy students: BSc 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="Final award level",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,(49:255))],id.vars='Student.code') # this differs from MSc code temp1$variable<-as.character(temp1$variable) temp1$value[temp1$value=='N/A']<-NA # change 'Module.' names into 'M.' for(i in 1:14) { 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:23])] # this line differs from MSc 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)[c(4:6,8,9)]<-c('attempt1','attempt2','blinded','timed','university') # different to MSc # drop empty variables - not in MSc assessments<-select(assessments,-one_of(c('B..1','U..1'))) # drop empty assessments assessments<-filter(assessments,!(blinded==""&university=="")) # some missing timed - not in MSc assessments$timed[assessments$Module.code=='M4' & assessments$blinded=='yes' & assessments$timed=='']<-'yes' assessments$timed[assessments$Module.code=='M11' & assessments$blinded=='yes' & assessments$timed=='']<-'no' # 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))) # merge in level of module modules<-read.csv('bsc-modules.csv') modules<-filter(modules,!is.na(Module.level)) assessments<-left_join(assessments,modules) # 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('bsc_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='bsc_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('bsc_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='bsc_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)) qmark_all<-quantile(as.numeric(assessments$attempt1), 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. ', 'Taking all assessments together, these statistics were ', qmark_all[2],' (',qmark_all[1],' to ',qmark_all[3],'). ')) # describe difference by WP factors: html_h('Crude associations with widening participation factors') png('bsc_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,5.7)) axis(1,las=2,labels=c('White', 'Asian', 'Black', 'Other', 'Unknown'), at=c(1,2,3,4,5)) dev.off() figureno<-html_img(img_file='bsc_box_all_attempt1_ethnic.png', caption='Boxplot of marks by ethnic group', figureno=figureno) png('bsc_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='bsc_box_all_attempt1_polar.png', caption='Boxplot of marks by POLAR', figureno=figureno) png('bsc_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='bsc_box_all_attempt1_disability.png', caption='Boxplot of marks by disabiliity', figureno=figureno) png('bsc_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='bsc_box_all_attempt1_parents.png', caption='Boxplot of marks by whether a parent attended HE', figureno=figureno) png('bsc_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('A levels', 'Access / BTEC', 'Degree level', 'Other or not known'), at=1:4) dev.off() figureno<-html_img(img_file='bsc_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('Given the large number of assessments and that many of these must be included as ', 'latent variables in the model, we consider only the factors that show some ', 'sign of association with marks in the boxplots above: ethnicity and disability. ', 'The table below shows the maximal model formed with all predictors. The ', 'Stan code for fitting the model can be read by ', 'clicking here.')) html_p('For student i taking module j:




') # 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. # 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$over20<-as.numeric(assessments$Age.code=='21-29' | 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_access<-as.numeric(as.numeric(assessments$Qual.code)==2) assessments$qual_degree<-as.numeric(as.numeric(assessments$Qual.code)==3) assessments$qual_other<-as.numeric(as.numeric(assessments$Qual.code)==4) 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-test11bsc.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 over201[Nmark]; int over200[N0]; int over2040[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_access1[Nmark]; int qual_access0[N0]; int qual_access40[N40]; int qual_degree1[Nmark]; int qual_degree0[N0]; int qual_degree40[N40]; int qual_other1[Nmark]; int qual_other0[N0]; int qual_other40[N40]; real mark[Nmark]; int stu1[Nmark]; int stu0[N0]; int stu40[N40]; int mod1[Nmark]; int mod0[N0]; int mod40[N40]; } parameters { real alpha; real beta[18]; real sigma; real latentmark0[N0]; real latentmark40[N40]; real ustu[Nstu]; real sigmastu; real umod[Nmod]; real 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]*over201[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])+ (beta[12]*polar1[i,1])+ (beta[13]*polar1[i,2])+ (beta[14]*polar1[i,3])+ (beta[15]*polar1[i,4])+ (beta[16]*qual_access1[i])+ (beta[17]*qual_degree1[i])+ (beta[18]*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]*blinded_untimed0[i])+ (beta[2]*blinded_timed0[i])+ (beta[3]*clinical0[i])+ (beta[4]*female0[i])+ (beta[5]*over200[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])+ (beta[12]*polar0[i,1])+ (beta[13]*polar0[i,2])+ (beta[14]*polar0[i,3])+ (beta[15]*polar0[i,4])+ (beta[16]*qual_access0[i])+ (beta[17]*qual_degree0[i])+ (beta[18]*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]*blinded_untimed40[i])+ (beta[2]*blinded_timed40[i])+ (beta[3]*clinical40[i])+ (beta[4]*female40[i])+ (beta[5]*over2040[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])+ (beta[12]*polar40[i,1])+ (beta[13]*polar40[i,2])+ (beta[14]*polar40[i,3])+ (beta[15]*polar40[i,4])+ (beta[16]*qual_access40[i])+ (beta[17]*qual_degree40[i])+ (beta[18]*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_access0=a0$qual_access, qual_access1=a1$qual_access, qual_access40=a40$qual_access, qual_degree0=a0$qual_degree, qual_degree1=a1$qual_degree, qual_degree40=a40$qual_degree, qual_other0=a0$qual_other, qual_other1=a1$qual_other, qual_other40=a40$qual_other, over200=a0$over20, over201=a1$over20, over2040=a40$over20, 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-test11bsc.stan', data=wppsdata, iter=5000, chains = 4, seed=7144) stanoutput<-as.matrix(summary(wppsfit)[[1]])[c('alpha', paste0('beta[',1:11,']'), '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('Baseline mean mark (%)', 'Blinded, untimed', 'Blinded, timed', 'Clinical', 'Female', 'Over 20', 'Asian ethnicity', 'Black ethnicity', 'Other ethnicity', 'Physical disability', 'Learning disability', 'Other disability', 'Residual SD', 'Student ability SD', 'Module difficulty SD') # other row names in test11 but removed when we rolled back to test10 #'POLAR 1', #'POLAR 2', #'POLAR 4', #'POLAR 5', #'Access/BTEC qualification', #'Previous degree', #'Other pre-entry qualification', 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. Bayesian ', 'Data Analysis 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, falling a little short of statistical significance. Students aged ', '21 or over when commencing their studies got significantly higher marks than ', 'those under 21, while blinded asessments, both timed and untimed, got ', 'significantly lower marks than unblinded ones.')) save(wppsfit,file='wpps-bsc-maineffects.Rdata') # interaction model stancon<-file('wpps-inter2.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 over201[Nmark]; int over200[N0]; int over2040[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_access1[Nmark]; int qual_access0[N0]; int qual_access40[N40]; int qual_degree1[Nmark]; int qual_degree0[N0]; int qual_degree40[N40]; int qual_other1[Nmark]; int qual_other0[N0]; int qual_other40[N40]; real mark[Nmark]; int stu1[Nmark]; int stu0[N0]; int stu40[N40]; int mod1[Nmark]; int mod0[N0]; int mod40[N40]; } parameters { real alpha; real beta[18]; matrix[3,3] gamma; real sigma; real latentmark0[N0]; real latentmark40[N40]; real ustu[Nstu]; real sigmastu; real umod[Nmod]; real 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]*over201[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])+ (beta[12]*polar1[i,1])+ (beta[13]*polar1[i,2])+ (beta[14]*polar1[i,3])+ (beta[15]*polar1[i,4])+ (beta[16]*qual_access1[i])+ (beta[17]*qual_degree1[i])+ (beta[18]*qual_other1[i])+ (gamma[1,1]*ethn_asian1[i]*blinded_untimed1[i])+ (gamma[2,1]*ethn_black1[i]*blinded_untimed1[i])+ (gamma[3,1]*ethn_other1[i]*blinded_untimed1[i])+ (gamma[1,2]*ethn_asian1[i]*blinded_timed1[i])+ (gamma[2,2]*ethn_black1[i]*blinded_timed1[i])+ (gamma[3,2]*ethn_other1[i]*blinded_timed1[i])+ (gamma[1,3]*ethn_asian1[i]*clinical1[i])+ (gamma[2,3]*ethn_black1[i]*clinical1[i])+ (gamma[3,3]*ethn_other1[i]*clinical1[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]*over200[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])+ (beta[12]*polar0[i,1])+ (beta[13]*polar0[i,2])+ (beta[14]*polar0[i,3])+ (beta[15]*polar0[i,4])+ (beta[16]*qual_access0[i])+ (beta[17]*qual_degree0[i])+ (beta[18]*qual_other0[i])+ (gamma[1,1]*ethn_asian0[i]*blinded_untimed0[i])+ (gamma[2,1]*ethn_black0[i]*blinded_untimed0[i])+ (gamma[3,1]*ethn_other0[i]*blinded_untimed0[i])+ (gamma[1,2]*ethn_asian0[i]*blinded_timed0[i])+ (gamma[2,2]*ethn_black0[i]*blinded_timed0[i])+ (gamma[3,2]*ethn_other0[i]*blinded_timed0[i])+ (gamma[1,3]*ethn_asian0[i]*clinical0[i])+ (gamma[2,3]*ethn_black0[i]*clinical0[i])+ (gamma[3,3]*ethn_other0[i]*clinical0[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]*over2040[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])+ (beta[12]*polar40[i,1])+ (beta[13]*polar40[i,2])+ (beta[14]*polar40[i,3])+ (beta[15]*polar40[i,4])+ (beta[16]*qual_access40[i])+ (beta[17]*qual_degree40[i])+ (beta[18]*qual_other40[i])+ (gamma[1,1]*ethn_asian40[i]*blinded_untimed40[i])+ (gamma[2,1]*ethn_black40[i]*blinded_untimed40[i])+ (gamma[3,1]*ethn_other40[i]*blinded_untimed40[i])+ (gamma[1,2]*ethn_asian40[i]*blinded_timed40[i])+ (gamma[2,2]*ethn_black40[i]*blinded_timed40[i])+ (gamma[3,2]*ethn_other40[i]*blinded_timed40[i])+ (gamma[1,3]*ethn_asian40[i]*clinical40[i])+ (gamma[2,3]*ethn_black40[i]*clinical40[i])+ (gamma[3,3]*ethn_other40[i]*clinical40[i])+ ustu[stu40[i]]+umod[mod40[i]]; latentmark40[i] ~ normal(mu40[i],sigma); } } ',file=stancon) close(stancon) wppsfit<-stan(file='wpps-inter2.stan', data=wppsdata, iter=5000, chains = 4, seed=7316) stanoutput<-as.matrix(summary(wppsfit)[[1]])[c('alpha', paste0('beta[',1:11,']'), paste0('gamma[',rep(1:3,3),',',rep(1:3,each=3),']'), '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('Baseline mean mark (%)', 'Blinded, untimed', 'Blinded, timed', 'Clinical', 'Female', 'Over 20', 'Asian ethnicity', 'Black ethnicity', 'Other ethnicity', 'Physical disability', 'Learning disability', 'Other disability', 'Asian ethnicity + blinded, untimed', 'Black ethnicity + blinded, untimed', 'Other ethnicity + blinded, untimed', 'Asian ethnicity + blinded, timed', 'Black ethnicity + blinded, timed', 'Other ethnicity + blinded, timed', 'Asian ethnicity + clinical', 'Black ethnicity + clinical', 'Other ethnicity + clinical', 'Residual SD', 'Student ability SD', 'Module difficulty SD') # other row names in test11 but removed when we rolled back to test10 #'POLAR 1', #'POLAR 2', #'POLAR 4', #'POLAR 5', #'Access/BTEC qualification', #'Previous degree', #'Other pre-entry qualification', 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 ', 'clicking here. ', 'To arrive at this final model, we removed POLAR and qualifications from ', 'the model above and then explored the following interaction terms: ', '
  • ethnicity with assessment type: includes significant ', 'interactions, see table below.
    • ', '
  • ethnicity with sex: none were significant or borderline
    • ', '
  • ethnicity with age could not be estimated because of low ', 'student numbers in combinations of groups
    • ', '
  • disability with assessment type: could not reliably be ', 'estimated because of low student numbers in combinations of groups
    • ', '
  • disability with sex: none were significant or borderline
    • ', '
  • disability with age could not be estimated because of low ', 'student numbers in combinations of groups
.', 'On the basis of these analyses, we retained ethnicity vs assessment type.')) tableno<-html_anytab(stanoutput,tableno=tableno,caption='Model with significant WP interactions') save(wppsfit,file='wpps-bsc-interactions.Rdata') ############################################################################### ############## analyse student-level associations with award ############## ############################################################################### html_h('Student-level predcitors of award level') html_p(paste0('We analysed associations with award level using ordinal logistic regression. In the ', 'statistics that follow, the effects are shown as odds ratios for a one step ', 'worse award in the following classification: 1st, 2:1, 2:2, non-qualifying degree, ', 'certificate or diploma, no award (means withdrawal for most students). ', 'So, a value greater than 1 represents worse results for the group concerned, and ', 'values less that 1 represent better results.')) # make numeric covariates students$n_polar<-as.numeric(students$POLAR.3) # needs not applicable and maybe non-UK to be treated as missing students$polar1<-as.numeric(students$n_polar==2) # POLAR.3 1st quintile has numeric code 2... students$polar2<-as.numeric(students$n_polar==3) # use 3rd quintile as baseline students$polar4<-as.numeric(students$n_polar==5) students$polar5<-as.numeric(students$n_polar==6) students$polar1[!(students$n_polar>1 & students$n_polar<7)]<-NA students$polar2[!(students$n_polar>1 & students$n_polar<7)]<-NA students$polar4[!(students$n_polar>1 & students$n_polar<7)]<-NA students$polar5[!(students$n_polar>1 & students$n_polar<7)]<-NA students$polar_all<-students$n_polar-1 students$polar_all[students$polar_all>5]<-NA students$polar_all[students$polar_all==0]<-NA students$female<-as.numeric(students$Gender.code)-1 students$over20<-as.numeric(students$Age.code=='21-29' | students$Age.code=='30+') students$ethn_asian<-as.numeric(as.numeric(students$Ethnicity.group.code)==2) students$ethn_black<-as.numeric(as.numeric(students$Ethnicity.group.code)==3) students$ethn_other<-as.numeric(as.numeric(students$Ethnicity.group.code)==4) students$ethn_all<-as.numeric(students$Ethnicity.group.code) students$ethn_all[students$ethn_all>4]<-NA students$disab_phys<-as.numeric(as.numeric(students$Disability.code)==2) students$disab_learn<-as.numeric(as.numeric(students$Disability.code)==3) students$disab_other<-as.numeric(as.numeric(students$Disability.code)==5) students$disab_all<-as.numeric(students$Disability.code) students$disab_all[students$disab_all==4]<-NA students$disab_all[students$disab_all==5]<-4 students$disab_all[students$disab_all>4]<-NA students$qual_access<-as.numeric(as.numeric(students$Qual.code)==2) students$qual_degree<-as.numeric(as.numeric(students$Qual.code)==3) students$qual_other<-as.numeric(as.numeric(students$Qual.code)==4) students$qual_all<-as.numeric(students$Qual.code) students$qual_all[students$qual_all>4]<-NA students$award<-as.numeric(students$Award.class.code) students$award[students$award==4]<-NA # remove 'still enrolled' students$award[students$award==5]<-4 # non-qualifying degree is the 4th category students$award[students$award==7]<-5 # certificate/diploma is the 5th category # all the 'no award' students have legitimate reasons under WD.code so no exclusions needed cor.test(as.numeric(students$award),students$polar_all,method='kendall',use='complete.obs') kruskal.test(students$award,as.factor(students$disab_all)) kruskal.test(students$award,as.factor(students$ethn_all)) wilcox.test(students$award,students$female) wilcox.test(students$award,students$over30) kruskal.test(students$award,as.factor(students$qual_all)) students$award<-as.factor(students$award) # the age effect is NS adjusting for sex. # only learning disability shows a sig/borderline effect # degree entry qualification NS ologmod<-polr(award~female+over20+ethn_asian+ethn_black+ethn_other+ disab_learn+disab_phys+disab_other+qual_access+qual_degree+qual_other, data=students,Hess=TRUE) ologresults <- coef(summary(ologmod)) ologresults <- ologresults[-((dim(ologresults)[1]-4):(dim(ologresults)[1])),] colnames(ologresults)[1]<-'OR' # convert to OR ologlogors<-ologresults[,1] ologresults[,1]<-exp(ologlogors) ologresults <- cbind(ologresults, "p value"=2*pnorm(abs(ologresults[,"t value"]),lower.tail=FALSE), "95% CI"=paste(round(exp(ologlogors-1.96*ologresults[,'Std. Error']), digits=4), 'to', round(exp(ologlogors+1.96*ologresults[,'Std. Error']), digits=4))) ologresults[,1:4]<-as.character(round(as.numeric(ologresults[,1:4]),digits=4)) tableno<-html_anytab(ologresults,tableno=tableno,caption='Significant predictors of award level') html_p(paste0('Female students and mature students had markedly better awards.. ', 'Learning disability was associated with worse awards, while physical and ', 'other disabilities were not significant. All non-White ethnic groups had ', 'significantly worse awards than White students, adjusting for the other factors. ', 'Access-level, and other entry qualifications were associated with significantly ', 'worse awards, while a prior degree was not significantly different to A-levels ', '(the baseline for comparison).')) # same but logistic regression for abeyance students$abeyance<-2-as.numeric(students$Interrupt.Abeyance.code) logmod<-glm(abeyance~female+over20+ethn_asian+ethn_black+ethn_other+ disab_learn+disab_phys+disab_other+qual_access+qual_degree+qual_other, data=students,family=binomial) logresults <- coef(summary(logmod)) # convert to OR logresults[,1]<-exp(logresults[,1]) colnames(logresults)[1]<-'OR' colnames(logresults)[4]<-'p value' logresults <- cbind(logresults, "95% CI"=paste(round(exp(logresults[,'OR']-1.96*logresults[,'Std. Error']), digits=4), 'to', round(exp(logresults[,'OR']+1.96*logresults[,'Std. Error']), digits=4))) logresults[,1:4]<-as.character(round(as.numeric(logresults[,1:4]),digits=4)) tableno<-html_anytab(logresults,tableno=tableno,caption='Significant predictors of abeyance') html_p(paste0('Female students had slightly better awards, while age was not significant. ', 'Learning disability was associated with worse awards, while physical and ', 'other disabilities were not significant. All non-White ethnic groups had ', 'significantly worse wards than White students, adjusting for the other factors. ', 'Access-level, and other entry qualifications were associated with significantly ', 'worse awards, while a prior degree was not significantly different to A-levels ', '(the baseline for comparison).')) close(con)