library(lcmm)
data2<-read.csv(file.choose())
data2$day<-as.factor(data2$day)
View(data2)

md1<- hlme(fixed = GDF15 ~ -1 + day,
           random=~- 1 + day,
           ng=1,
           idiag=FALSE,
           data=data2,
           subject='id')
BIC<-c(md1$ng,md1$BIC)
models<-list()

for(i in 2:6) {
  mdi <- hlme(fixed = GDF15 ~ - 1 + day,
              mixture = ~ - 1 + day,
              random = ~ -1 + day,
              ng=i,
              nwg = TRUE,
              idiag = FALSE,
              data = data2,
              subject = "id",
              B = md1
  )
  print(paste('number of group:',i,'BIC',mdi$BIC))
  BIC<-rbind(BIC,c(i, mdi$BIC))
  models[[i]]<-mdi
}
best_model1 = models[[1]]
best_model2  = models[[2]]
best_model3  = models[[3]]
best_model4 = models[[4]]
best_model5 = models[[5]]
best_model6 = models[[6]]

summarytable(md1,md2,md3,md4,md5,md6, which = c("G", "loglik", "conv", "npm", "AIC", "BIC", "SABIC", "entropy", "%class"))
summarytable(best_model,best_model2,best_model3,best_model4,best_model5,best_model6,
             which = c("G", "loglik", "conv", "npm", "AIC", "BIC", "SABIC", "entropy", "%class"))

best_model<-models[[4]]
summary(best_model)
best_model$BIC
best_model$pprob
summarytable(best_model)
library(dplyr)
par(lwd = 2 , cex.lab = 1.8 , cex.axis = 1.8 , mai = c(1,1,1,1)) 
newdata<-data.frame(day=c(1,3,7))
plotpred<-predictY(best_model,newdata,
                   var.time='day',draws = TRUE)
plot(plotpred,lwd=3,lty=1,xlab='day',ylab='GDF15',
     legend.loc='topleft',cex=0.75 , shades = true)

plot(plotpred, lty=1, xlab='day', ylab='aim',
     legend.loc='topleft', cex=0.75,  xlim=c(1, 7))

plotpred <- predictY(best_model, newdata, var.time = 'Day', draws = FALSE)
plot(plotpred, col=c("#97CBB7","#DFB58F","#A3836D","#F3A7AD") ,lwd=3, lty = 1, xlab = 'Day', ylab = 'GDF15',
     legend.loc = FALSE, cex = 0.75, shade = TRUE)



summary(predictY(best_model, newdata, var.time = 'day', draws = TRUE))

summary(best_model)

x<-best_model$pprob[,1:2]
newdata<-merge(data,x,by = "id" )
View(newdata)
write.csv(newdata,file = "nwedata.csv")