################################
#########  IMPUTATION  #########
################################
library(openxlsx)
rawdata  <- read.xlsx("HPF dataset_JEEHP.xlsx", sheet = 1)

# Reverse Coding #
rawdata$DREEM4_R <- 4 - rawdata$DREEM4
rawdata$DREEM8_R <- 4 - rawdata$DREEM8
rawdata$DREEM9_R <- 4 - rawdata$DREEM9
rawdata$DREEM17_R <- 4 - rawdata$DREEM17
rawdata$DREEM25_R <- 4 - rawdata$DREEM25
rawdata$DREEM35_R <- 4 - rawdata$DREEM35
rawdata$DREEM39_R <- 4 - rawdata$DREEM39
rawdata$DREEM48_R <- 4 - rawdata$DREEM48
rawdata$DREEM50_R <- 4 - rawdata$DREEM50

library(mice)
library(miceadds)

# two-way imputation for Perception items (DREEM) #
rawdata1 <- rawdata[, c(6:55, 68:76)] 
rawdata2 <- rawdata1[, -c(4,8, 9,17, 25, 35, 39, 48, 50)]
rawdata2.imp <- tw.imputation(rawdata2, integer=FALSE)

# two-way imputation for Intrinsic Motivation items #
rawdata3 <- rawdata[, c(56:67)]
rawdata3.imp <- tw.imputation(rawdata3, integer=FALSE)

# Analysis Dataset #
dataP1 <- rawdata[,1:5]
dataP2 <- cbind(rawdata2.imp, rawdata3.imp)
preprocessed.data <- cbind(dataP1, dataP2)
preprocessed.data$total.DREEM <- rowSums(preprocessed.data[, c(6:55)])
preprocessed.data$total.IM <- rowSums(preprocessed.data[, c(56:67)])


################################
########### ANALYSIS ########### 
################################
HPF.data  <- read.xlsx("HPF dataset_JEEHP.xlsx", sheet = 2)
HPF.data$grading.system <- as.integer(HPF.data$grading.system)

# Participants #
sexdistr <- table(HPF.data$grading.system, HPF.data$sex)
pathdistr <- table(HPF.data$grading.system, HPF.data$admission.type) 
chisq.test(sexdistr);chisq.test(pathdistr)

describeBy(HPF.data$total.DREEM, group = HPF.data$grading.system)
describeBy(HPF.data$total.IM, group = HPF.data$grading.system)
describeBy(HPF.data$performance, group = HPF.data$grading.system)

var.test(total.DREEM ~ grading.system, data=HPF.data)
var.test(total.IM ~ grading.system, data=HPF.data)
var.test(performance~ grading.system, data=HPF.data)

t.test(total.DREEM ~ grading.system, data=HPF.data, var.equal=TRUE)
t.test(total.IM ~ grading.system, data=HPF.data, var.equal=TRUE)
t.test(performance ~ grading.system, data=HPF.data, var.equal=TRUE)

# Main results - 1. Correlation #
cor.test(HPF.data$grading.system, HPF.data$performance)
cor.test(HPF.data$grading.system, HPF.data$total.IM)
cor.test(HPF.data$grading.system, HPF.data$total.DREEM)
cor.test(HPF.data$total.DREEM, HPF.data$performance)
cor.test(HPF.data$total.DREEM,HPF.data$total.IM)
cor.test(HPF.data$total.IM, HPF.data$performance)  

# Main results - 2. Serial mediation analysis #
process(data=HPF.data,
        y="performance",
        x="grading.system",
        m=c("total.DREEM","total.IM"),
        cov=c( "sex", "admission.type"),
        model=6,
        effsize=1,
        total=1,
        stand=1,
        boot=10000,
        modelbt=1,
        seed=123456)