# R and Stan code for a research collaboration between
# St George's, University of London,
# Kingston University and
# Princess Alice Hospice
#
# Results contained in the paper:
# Reed E, Todd J, Lawton S, Grant RL, Sadler C, Berg J,
# Lucas C, Watson M. "A multi-professional educational
# intervention to improve and sustain respondents’
# confidence to deliver palliative care: A mixed-methods
# study". Palliative Medicine (2018); 32(2): 571-580.



############################################################################
# This is free and unencumbered software released into the public domain.
#
# Anyone is free to copy, modify, publish, use, compile, sell, or
# distribute this software, either in source code form or as a compiled
# binary, for any purpose, commercial or non-commercial, and by any
# means.
#
# In jurisdictions that recognize copyright laws, the author or authors
# of this software dedicate any and all copyright interest in the
# software to the public domain. We make this dedication for the benefit
# of the public at large and to the detriment of our heirs and
# successors. We intend this dedication to be an overt act of
# relinquishment in perpetuity of all present and future rights to this
# software under copyright law.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
# IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.
#
# For more information, please refer to <http://unlicense.org/>
############################################################################



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

# you can obtain this from my code repo robertgrantstats.co.uk
source('~/git/html-reports/html-reports.R')



# get ECEPC dataset
ec<-read.csv('ecepc-6mo.csv')
ec$ntime=6
ec<-ec[,c(dim(ec)[2],(1:(dim(ec)[2]-1)))]
ec<-filter(ec,consent!="No")
ec<-filter(ec,!is.na(candidate_number))

# append baseline and 3 month data
ec3<-read.csv('ecepc-baseline-3mo.csv')
# add empty columns for compatibility
ec3$role_changed<-rep("",dim(ec3)[1])
ec3$role_changed_detail<-rep("",dim(ec3)[1])
ec3$ebp_barriers<-rep("",dim(ec3)[1])
ec3<-ec3[,c(1:16,59,60,17:33,61,34:58)] # bit of hard-coding here

ec$done_online<-rep("",dim(ec)[1])
ec$done_distance<-rep("",dim(ec)[1])
ec$done_reflective<-rep("",dim(ec)[1])
ec$done_remote_support<-rep("",dim(ec)[1])
ec$done_blended<-rep("",dim(ec)[1])
ec$exp_online<-rep("",dim(ec)[1])
ec$exp_distance<-rep("",dim(ec)[1])
ec$exp_reflective<-rep("",dim(ec)[1])
ec$exp_remote_support<-rep("",dim(ec)[1])
ec$exp_blended<-rep("",dim(ec)[1])
ec$edu_text<-rep("",dim(ec)[1])
ec$journals_applied<-rep("",dim(ec)[1])
ec<-ec[,c(1:19,49:59,20,21,60,22:48)]
ec3<-ec3[,-(dim(ec3)[2])]
names(ec3)[1]<-'ntime'

# compare variable names side-by-side
mm<-max(length(names(ec)),length(names(ec3)))
cbind(c(names(ec3),rep("",mm-length(names(ec3)))),
      c(names(ec),rep("",mm-length(names(ec)))))
rm(mm)

# append
ec<-rbind(ec3,ec)
rm(ec3)

# make binary predictors
ec$female<-as.numeric(ec$sex=='Female')
ec$female[ec$sex=='']<-NA
ec$time3<-as.numeric(ec$ntime==3)
ec$time6<-as.numeric(ec$ntime==6) # this compares 6 to 0, but can change...
ec$role_doctor<-as.numeric(ec$role=="Doctor" | ec$role=="Specialist Doctor")
ec$role_other<-as.numeric(ec$role=="Occupational therapist" |
                            ec$role=="Pharmacist" |
                            ec$role=="Other*" |
                            ec$role=="")
ec$place_hospital<-as.numeric(ec$workplace=="Hospital")
ec$place_carehome<-as.numeric(ec$workplace=="Care home")
ec$place_primary<-as.numeric(ec$workplace=="Primary care/community")
ec$place_other<-as.numeric(ec$workplace=="Other*" |
                           ec$workplace=="")

# recoding time in role
temp<-as.numeric(ec$time_in_role)
ec$time_in_role<-temp
ec$time_in_role[temp==6]<-2
ec$time_in_role[temp<6 & temp>1]<-temp[temp<6 & temp>1]+1
ec$time_in_role[temp==8]<-NA

# recoding confidence scale
for(i in 37:59) {
  temp<-ec[,i]
  temp<-as.character(temp)
  temp<-str_sub(temp,start=1,end=2)
  ec[,i]<-as.numeric(temp)+1
}
rm(temp)


# define complete dataset
comp<-filter(ec,!is.na(Age) &
               !is.na(female) &
               !is.na(d_effects_pt) &
               !is.na(time_in_role) &
               !is.na(d_effects_fam) &
               !is.na(d_dying) &
               !is.na(d_death_pt) &
               !is.na(d_death_fam) &
               !is.na(d_bereavement) &
               !is.na(d_prognosis) &
               !is.na(d_pain) &
               !is.na(assess) &
               !is.na(aetiology) &
               !is.na(symptoms) &
               !is.na(prescribe) &
               !is.na(effects) &
               !is.na(psychological) &
               !is.na(social) &
               !is.na(spiritual) &
               !is.na(mdt) &
               !is.na(ref_physio) &
               !is.na(ref_ot) &
               !is.na(ref_complementary) &
               !is.na(ref_lymphoedema) &
               !is.na(ref_psychiatric) &
               !is.na(ref_spiritual))

# introduce ID number
comp$nid<-as.numeric(as.factor(comp$candidate_number))

# sort
comp<-arrange(comp,nid,ntime)

#tag ID
comp$tag<-0
comp$tag[1]<-1
comp$tag[c(FALSE,comp$nid[-1]!=comp$nid[-dim(comp)[1]])]<-1

# read in data
ecepcdata<-list(K=11,
                D=3,
                Dp=9,
                N=dim(comp)[1],
                M=max(comp$nid),
                part=comp$nid,
                d_effects_pt=comp$d_effects_pt,
                d_effects_fam=comp$d_effects_fam,
                d_dying=comp$d_dying,
                d_death_pt=comp$d_death_pt,
                d_death_fam=comp$d_death_fam,
                d_bereavement=comp$d_bereavement,
                d_prognosis=comp$d_prognosis,
                d_pain=comp$d_pain,
                assess=comp$assess,
                aetiology=comp$aetiology,
                symptoms=comp$symptoms,
                prescribe=comp$prescribe,
                effects=comp$effects,
                psychological=comp$psychological,
                social=comp$social,
                spiritual=comp$spiritual,
                mdt=comp$mdt,
                ref_physio=comp$ref_physio,
                ref_ot=comp$ref_ot,
                ref_complementary=comp$ref_complementary,
                ref_lymphoedema=comp$ref_lymphoedema,
                ref_psychiatric=comp$ref_psychiatric,
                ref_spiritual=comp$ref_spiritual,
                x=cbind(rep(1,dim(comp)[1]),
                        comp$time3,
                        comp$time6),
                xp=cbind(comp$Age[comp$tag==1],
                        comp$female[comp$tag==1],
                        comp$time_in_role[comp$tag==1],
                        comp$role_doctor[comp$tag==1],
                        comp$role_other[comp$tag==1],
                        comp$place_hospital[comp$tag==1],
                        comp$place_carehome[comp$tag==1],
                        comp$place_primary[comp$tag==1],
                        comp$place_other[comp$tag==1]))
stancon<-file('ecepc-sem7a.stan',open='w')
write('
data {
      int<lower=2> K; // number of categories
      int<lower=0> N; // number of responses
      int<lower=0> M; // number of participants
      int<lower=0> D; // number of predictors at response level
      int<lower=0> Dp; // number of predictors at person level
      int<lower=1,upper=M> part[N];
      int<lower=1,upper=K> d_effects_pt[N];
      int<lower=1,upper=K> d_effects_fam[N];
      int<lower=1,upper=K> d_dying[N];
      int<lower=1,upper=K> d_death_pt[N];
      int<lower=1,upper=K> d_death_fam[N];
      int<lower=1,upper=K> d_bereavement[N];
      int<lower=1,upper=K> d_prognosis[N];
      int<lower=1,upper=K> d_pain[N];
      int<lower=1,upper=K> assess[N];
      int<lower=1,upper=K> aetiology[N];
      int<lower=1,upper=K> symptoms[N];
      int<lower=1,upper=K> prescribe[N];
      int<lower=1,upper=K> effects[N];
      int<lower=1,upper=K> psychological[N];
      int<lower=1,upper=K> social[N];
      int<lower=1,upper=K> spiritual[N];
      int<lower=1,upper=K> mdt[N];
      int<lower=1,upper=K> ref_physio[N];
      int<lower=1,upper=K> ref_ot[N];
      int<lower=1,upper=K> ref_complementary[N];
      int<lower=1,upper=K> ref_lymphoedema[N];
      int<lower=1,upper=K> ref_psychiatric[N];
      int<lower=1,upper=K> ref_spiritual[N];
      row_vector[D] x[N]; // predictors at response level, including a column of ones
      row_vector[Dp] xp[M]; // predictors at person level
  }
  parameters {
      vector[M] u_part; // person random effect for all confidence
      real<lower=0> tau_part; // SD of person random effect
      vector[18] u_topic; // topic random effect
      real<lower=0> tau_topic; // SD
      vector[D] beta; // regression coefficients at response level
      vector[Dp] beta_p; // regression coefficients at person level
      ordered[K-1] c;
      real pt; // additional effect of talking to the patient rather than family
  }
  model {
      real mu[M];
      beta ~ normal(0,20); // priors
      beta_p ~ normal(0,20); // priors
      tau_part ~ normal(0,20); // prior
      u_part ~ normal(0,tau_part);
      tau_topic ~ normal(0,20); // prior
      u_topic ~ normal(0,tau_topic);
      c ~ normal(0,20); // priors
      for (i in 1:M) {
        mu[i] <- u_part[i] + xp[i]*beta_p;
      }
      for (i in 1:N) {
        d_effects_pt[i] ~ ordered_logistic(u_topic[1]+pt+mu[part[i]]+(x[i]*beta), c);
        d_effects_fam[i] ~ ordered_logistic(u_topic[1]+mu[part[i]]+(x[i]*beta), c);
        d_dying[i] ~ ordered_logistic(u_topic[2]+mu[part[i]]+(x[i]*beta), c);
        d_death_pt[i] ~ ordered_logistic(u_topic[3]+pt+mu[part[i]]+(x[i]*beta), c);
        d_death_fam[i] ~ ordered_logistic(u_topic[3]+mu[part[i]]+(x[i]*beta), c);
        d_bereavement[i] ~ ordered_logistic(u_topic[4]+mu[part[i]]+(x[i]*beta), c);
        d_prognosis[i] ~ ordered_logistic(u_topic[5]+mu[part[i]]+(x[i]*beta), c);
        d_pain[i] ~ ordered_logistic(u_topic[6]+mu[part[i]]+(x[i]*beta), c);
        assess[i] ~ ordered_logistic(u_topic[7]+mu[part[i]]+(x[i]*beta),c);
        aetiology[i] ~ ordered_logistic(u_topic[8]+mu[part[i]]+(x[i]*beta),c);
        symptoms[i] ~ ordered_logistic(u_topic[9]+mu[part[i]]+(x[i]*beta),c);
        prescribe[i] ~ ordered_logistic(u_topic[10]+mu[part[i]]+(x[i]*beta),c);
        effects[i] ~ ordered_logistic(u_topic[1]+mu[part[i]]+(x[i]*beta),c);
        psychological[i] ~ ordered_logistic(u_topic[11]+mu[part[i]]+(x[i]*beta),c);
        social[i] ~ ordered_logistic(u_topic[12]+mu[part[i]]+(x[i]*beta),c);
        spiritual[i] ~ ordered_logistic(u_topic[13]+mu[part[i]]+(x[i]*beta),c);
        mdt[i] ~ ordered_logistic(u_topic[14]+mu[part[i]]+(x[i]*beta),c);
        ref_physio[i] ~ ordered_logistic(u_topic[15]+mu[part[i]]+(x[i]*beta),c);
        ref_ot[i] ~ ordered_logistic(u_topic[16]+mu[part[i]]+(x[i]*beta),c);
        ref_complementary[i] ~ ordered_logistic(u_topic[17]+mu[part[i]]+(x[i]*beta),c);
        ref_lymphoedema[i] ~ ordered_logistic(u_topic[18]+mu[part[i]]+(x[i]*beta),c);
        ref_psychiatric[i] ~ ordered_logistic(u_topic[11]+mu[part[i]]+(x[i]*beta),c);
        ref_spiritual[i] ~ ordered_logistic(u_topic[13]+mu[part[i]]+(x[i]*beta),c);
      }
  }
      ',file=stancon)
close(stancon)
ecepcfit<-stan(file='ecepc-sem7a.stan',
              data=ecepcdata,
              iter=4000,
              warmup=2000,
              chains = 4,
              seed=3871)
print(ecepcfit)
save(ecepcfit,file='ecepc-sem7a.Rdata')