library(mice)
library(foreign)
library(MASS)
library(writexl)
library(car)

##########################################################################################################################
########################################### Model without Dropouts #########################################################

###model with missing data to check VIF -> lm to use vif()
fit1_missing <- lm(vas_today_a ~ vas_retro_bl + study_center + Sex + Age + CASMIN + migration + Personshousehold + proxy +
                     Pre_fracture_dependency + Fracture_type + Charlson_commorbity_index_cat + polypharmacy + Social_support +
                     subjective_need + Surgery  + Referral_to_rehabilitation_facility + Malnutrition + PHQ_4_cat + OHS, data=dat_withoutDropout)
summary(fit1_missing)
## VIF
vif(fit1_missing)
#####
fit2_missing <- lm(EQ5Dindex_FU ~ EQ5Dindex + study_center + Sex + Age + CASMIN + migration + Personshousehold + proxy +
                     Pre_fracture_dependency + Fracture_type + Charlson_commorbity_index_cat + polypharmacy + Social_support +
                     subjective_need + Surgery + Referral_to_rehabilitation_facility + Malnutrition + PHQ_4_cat + OHS, data=dat_withoutDropout)
summary(fit2_missing)
## VIF
vif(fit2_missing)


###################### Imputation ##############################

ini <- mice(dat_withoutDropout, seed=1, maxit=0, m=5)

meth <- ini$method
meth
#### passive imputation for diff-variables
meth["vas_diff"] <- "~I(vas_today_a - vas_retro_bl)"
meth["EQ5Dindex_diff"] <- "~I(EQ5Dindex_FU - EQ5Dindex)"


predmat <- ini$predictorMatrix
# Predictormatrix
predmat
#### change Predictormatrix
predmat[,20] <- c(rep(0,25))
predmat[,21] <- c(rep(0,25))

###### compute imputation #######
imp <- mice(dat_withoutDropout, seed=1,method=meth,predictorMatrix=predmat, maxit=20, m=10)

imp$imp
imp$loggedEvents

##### Diagnostics ############
# for numerical variables
stripplot(imp)
bwplot(imp)
densityplot(imp)
#convergence?
plot(imp)

#####################################################################################################
############################ Analysis ##############################
##################################################
# linear regression
fit1 <- with(imp, lm(vas_today_a ~ vas_retro_bl + study_center + Sex + CASMIN + migration + Personshousehold + proxy +
                       Pre_fracture_dependency + Fracture_type + Charlson_commorbity_index_cat + polypharmacy + Social_support +
                       subjective_need + Surgery + ICU_episode + Referral_to_rehabilitation_facility + Malnutrition + PHQ_4_cat + OHS))


summary.fit1 <- summary(pool(fit1))
summary.fit1$p.value <- round(summary.fit1$p.value,3)
summary.fit1$df <- round(summary.fit1$df,1)
summary.fit1$statistic <- round(summary.fit1$statistic,2)
summary.fit1$std.error <- round(summary.fit1$std.error,1)
summary.fit1$estimate <- round(summary.fit1$estimate,2)
summary.fit1

pool.r.squared(fit1)
pool.r.squared(fit1, adjusted = T)

fit2 <- with(imp, lm(EQ5Dindex_FU ~ EQ5Dindex + study_center + Sex + Age + CASMIN + migration + Personshousehold + proxy +
                       Pre_fracture_dependency + Fracture_type + Charlson_commorbity_index_cat + polypharmacy + Social_support +
                       subjective_need + Surgery + ICU_episode + Referral_to_rehabilitation_facility + Malnutrition + PHQ_4_cat + OHS))


summary.fit2 <- summary(pool(fit2))
summary.fit2$p.value <- round(summary.fit2$p.value,3)
summary.fit2$df <- round(summary.fit2$df,1)
summary.fit2$statistic <- round(summary.fit2$statistic,2)
summary.fit2$std.error <- round(summary.fit2$std.error,3)
summary.fit2$estimate <- round(summary.fit2$estimate,4)
summary.fit2

pool.r.squared(fit2)
pool.r.squared(fit2, adjusted = T)

#### export output
#update path!
write_xlsx(x = as.data.frame(summary.fit1) , path = 
             "C:/..../results_QoL_vas_FU_sensitivity1_withoutDropout.xlsx")
write_xlsx(x = as.data.frame(summary.fit2) , path = 
             "C:/..../results_QoL_EQ5D_FU_sensitivity1_withoutDropout.xlsx")